; vap_gts_6190.pro
; FINAL SYNTHETIC VAPOUR PRESSURE PROGRAM
; written by Tim Mitchell, 06.01.03 following Mark New's vap_gts.pro
; program to generate a set of synthetic anomaly grids for vapour pressure
; fixed ingredients:
;	1961-90 normals on 2.5deg grid for tmp,dtr,vap
; variable ingredients:
;	any period anomalies (v.1961-90) on 2.5deg grid for tmp,dtr 

pro vap_gts_6190,year1,year2,vap_out,dtr_ano,tmp_ano, $
	dumpbin=dumpbin,dumpglo=dumpglo,dumpclim=dumpclim,dumpraw=dumpraw

if n_params() lt 1 then begin
  print,'pro vap_gts_6190,year1,year2,vap_out,dtr_ano,tmp_ano,'
  print,'                 dumpbin=dumpbin,dumpglo=dumpglo'
  return
endif

close,/all
nor1=1961 & nor2=1990
climmonth=fltarr(144,72)
monthname=strarr(12)
monthname[*]=['01','02','03','04','05','06','07','08','09','10','11','12']

print, ''
print, '***** LOAD 2.5 NORMALS *****'
rdbin,dtrnor,'/cru/mark1/f080/gts/dtr/glo25/glo25.dtr.6190',gridsize=2.5,/quiet
lan25=where(dtrnor gt -9999, countland) & sea25=where(dtrnor le -9999, countsea)
dtrnor(lan25)=dtrnor(lan25)/10.0

rdbin,tmpnor,'/cru/mark1/f080/gts/tmp/glo25/glo25.tmp.6190',gridsize=2.5,/quiet
tmpnor(lan25)=tmpnor(lan25)/10.0

rdbin,vapnor,'/cru/tyn1/f014/_keep/glo25.vap.6190',gridsize=2.5,/quiet
vapnor(lan25)=vapnor(lan25)/10.0

print, '2.5deg grid: land,sea: ', countland/12,countsea/12
dtrmom=moment(dtrnor(lan25))
dtrmax=max(dtrnor(lan25),min=dtrmin)
print, 'dtr  (x,s2,<<,>>):',dtrmom[0],dtrmom[1],dtrmin,dtrmax
tmpmom=moment(tmpnor(lan25))
tmpmax=max(tmpnor(lan25),min=tmpmin)
print, 'tmp  (x,s2,<<,>>):',tmpmom[0],tmpmom[1],tmpmin,tmpmax
vapmom=moment(vapnor(lan25))
vapmax=max(vapnor(lan25),min=vapmin)
print, 'vap  (x,s2,<<,>>):',vapmom[0],vapmom[1],vapmin,vapmax

print, ''
print, '***** DEDUCE 2.5 NORMALS *****'

tmnnor=float(vapnor)*0.0 & tmnnor(sea25)=-999.0
tdwnor=float(vapnor)*0.0 & tdwnor(sea25)=-999.0
synnor=float(vapnor)*0.0 & synnor(sea25)=-999.0
tadj=float(vapnor)*0.0
v=float(vapnor)*0.0
nonzero=float(vapnor)*0.0

nonzero(lan25)=vapnor(lan25)>0.1
tdwnor(lan25)=tvap(nonzero(lan25))
tdwmom=moment(tdwnor(lan25))
tdwmax=max(tdwnor(lan25),min=tdwmin)
print, 'tdw  (x,s2,<<,>>):',tdwmom[0],tdwmom[1],tdwmin,tdwmax

tmnnor(lan25)=tmpnor(lan25)-(dtrnor(lan25)/2.0)
tadj(lan25)=tdwnor(lan25)-tmnnor(lan25)
tmnnor(lan25)=tmnnor(lan25)+tadj(lan25)
tmnmom=moment(tmnnor(lan25))
tmnmax=max(tmnnor(lan25),min=tmnmin)
print, 'tmn  (x,s2,<<,>>):',tmnmom[0],tmnmom[1],tmnmin,tmnmax

v(lan25)=6.108*exp( (17.27*tmnnor(lan25) /(237.3+tmnnor(lan25)) ))>0.05
synnor(lan25)=v(lan25)<esat(tmpnor(lan25))
synnor(lan25)=synnor(lan25)>0.01
synmom=moment(synnor(lan25))
synmax=max(synnor(lan25),min=synmin)
print, 'syn  (x,s2,<<,>>):',synmom[0],synmom[1],synmin,synmax

if keyword_set(dumpclim) then begin
   for imon=0,11 do begin
     climmonth=synnor(*,*,imon)
     
     SaveFile=vap_out+monthname(imon)+'.6190.syn.glo'
     SaveTitle=strip('syn 1961-90 normal for month '+monthname(imon))
     SaveGlo,23,climmonth,CallFile=Savefile,CallTitle=SaveTitle
   endfor
endif

print, ''
print, '***** CALC SYNTHETIC VAP *****'

vapyear=float(vapnor)*0.0
tmpyear=float(vapnor)*0.0
dtryear=float(vapnor)*0.0
tmnyear=float(vapnor)*0.0

vapmonth=fltarr(144,72)

for iy=year1,year2 do begin
 yeartext=string(iy,form='(i4)')

 tmpfl=strip(string(tmp_ano,iy))
 rdbin,tmpyear,tmpfl,gridsize=2.5,/quiet
 
 dtrfl=strip(string(dtr_ano,iy))
 rdbin,dtryear,dtrfl,gridsize=2.5,/quiet
 
; anyzero = where (tmpyear(lan25) eq 0.0 or dtryear(lan25) eq 0.0, nanyzero)
; if (nanyzero gt 0) then begin
;   tmpyear(anyzero)=0 & dtryear(anyzero)=0
; endif
 
 tmpyear(lan25)=tmpnor(lan25)+(tmpyear(lan25)/10.0)
 dtryear(lan25)=dtrnor(lan25)+(dtryear(lan25)/10.0)
 tmnyear(lan25)=tmpyear(lan25)-(dtryear(lan25)/2.0)+tadj(lan25)
 
 tmpyear(sea25)=-999.0
 dtryear(sea25)=-999.0
 tmnyear(sea25)=-999.0
 tmnnor (sea25)=   0.0
 
 if keyword_set(dumpraw) then begin
   for imon=0,11 do begin
     vapmonth=tmnyear(*,*,imon)-tmnnor(*,*,imon)
     
     SaveFile=vap_out+'tmn.'+monthname(imon)+'.'+yeartext+'.glo'
     SaveTitle=strip('raw tmn anom for month '+monthname(imon)+' in year '+yeartext)
     SaveGlo,23,vapmonth,CallFile=Savefile,CallTitle=SaveTitle
   endfor
 endif

 vapyear(lan25)=6.108*exp( (17.27*tmnyear(lan25)) /(237.3+tmnyear(lan25)) )>0.05
 vapyear(lan25)=vapyear(lan25)<esat(tmpyear(lan25))	; best est of vap

 vapmom=moment(vapyear(lan25))				; best est moments
 vapmax=max(vapyear(lan25),min=vapmin)			; best est min/max
 print, yeartext,' (x,s2,<<,>>):',vapmom[0],vapmom[1],vapmin,vapmax
 
 vapyear(lan25)=10.0*(vapyear(lan25)-synnor(lan25))	; anomalise, *10
 
 vapyear(lan25)=vapyear(lan25)>(-9998)			; retain limits
 vapyear(lan25)=vapyear(lan25)<9998
 vapyear(sea25)=0.0
 
 if keyword_set(dumpbin) then $
 	wrbin,fix(round(vapyear)),strip(string(vap_out,iy)),/compre
 
 if keyword_set(dumpglo) then begin
   for imon=0,11 do begin
     vapmonth=vapyear(*,*,imon)/10.0
     
     SaveFile=vap_out+monthname(imon)+'.'+yeartext+'.glo'
     SaveTitle=strip('synthetic vap for month '+monthname(imon)+' in year '+yeartext)
     SaveGlo,23,vapmonth,CallFile=Savefile,CallTitle=SaveTitle
   endfor
 endif
endfor

print,''

end