; vap_gts_anom.pro is based on Mark New's vap_gts.pro
; this version corrects some deficiencies in the Mark New version

; required inputs are:
; ** vapour pressure and temperature normals on 2.5deg grid
;	(these come ready-supplied for a 1961-90 normal period)
; ** temp and dtr monthly anomalies on 2.5deg grid, including normal period 

; outputs to binary file or .glo files
; binary is for subsequent input to quick_interp_tdm2.pro, using synth_prefix
; .glo is for checking purposes, using plotspatialauto.pro

; modified Mark New headers:
; calculates vapour pressure anomaly time series from dtr and tmp, using tmn 
;  as proxy for tdew
; this version also calculates the tdew normal from the vapour pressure normal
;  and then adjusts the tmin to have an average equal to tdew normal
; synthetic values are contrained to be between 0.1 and saturated vap
;  pressure at mean temperature
; Note that anomalies are in hPa*10 (bin) or hPa (glo)

pro vap_gts_anom,dtr_prefix=dtr_prefix,tmp_prefix=tmp_prefix,year1,year2,$
    nor1=nor1,nor2=nor2,outprefix=outprefix, $
    dumpbin=dumpbin,dumpglo=dumpglo,dumpclim=dumpclim

if n_params() lt 1 then begin
 print,' vap_gts_anom,dtr_prefix=dtr_prefix,tmp_prefix=tmp_prefix,year1,year2,'
 print,'         nor1=nor1,nor2=nor2,outprefix=outprefix,'
 print,'         dumpbin=dumpbin,dumpglo=dumpglo,dumpclim=dumpclim'
 return
endif
close,/all
if keyword_set(nor1) eq 0 then nor1=1961
if keyword_set(nor2) eq 0 then nor2=1990
if keyword_set(outprefix) eq 0 then outprefix='glo.vap.ano.'
outprefix2='glo.vap.'
if keyword_set(dtr_prefix) eq 0 then dtr_prefix='../../dtr/glo/glo.dtr.'
if keyword_set(tmp_prefix) eq 0 then tmp_prefix='../../tmp/glo/glo.tmp.'
;
monthname=strarr(12)
monthname[*]=['01','02','03','04','05','06','07','08','09','10','11','12']
vapmonth=fltarr(144,72)
tmpdtrfac=10.0
;
rdbin,vapnor,'/cru/tyn1/f014/_keep/glo25.vap.6190',gridsize=2.5,/quiet
nsea =where(vapnor le -9999, countsea)
nland=where(vapnor gt -9999, countland)
print,'Land,sea:',countland,countsea
vapnor=vapnor/10.0
;
; load tmp normal, used for checking only
rdbin,tmpnor,'/cru/tyn1/f014/_keep/glo25.tmp.6190',gridsize=2.5,/quiet
tmpnor(nland)=tmpnor(nland)/10.0
;
tadj  =float(vapnor)*0.0 & tadj  (nsea)=-9999
tdwnor=float(vapnor)*0.0
tmn   =float(vapnor)*0.0 & tmn   (nsea)=-9999
vapsyn=float(vapnor)*0.0
vapgrd=float(vapnor)*0.0
v     =float(vapnor)*0.0
tmnnor=float(vapnor)*0.0
;
; calculate 1961-1990 tmn normal from monthly tmp and dtr values
print,'Calculating tmn normal'
for iy=nor1,nor2 do begin
 tmpfl=strip(string(tmp_prefix,iy))
 dtrfl=strip(string(dtr_prefix,iy))
 rdbin,tmpgrd,tmpfl,gridsize=2.5,/quiet
 rdbin,dtrgrd,dtrfl,gridsize=2.5,/quiet
 tmnnor(nland)=tmnnor(nland)+tmpgrd(nland)-(0.5*dtrgrd(nland))
endfor
tmnnor(nland)=tmnnor(nland)/((nor2-nor1+1)*tmpdtrfac)

vapnor(nland)=vapnor(nland)>0.01
tdwnor(nland)=tvap(vapnor(nland))
tadj(nland)=tdwnor(nland)-tmnnor(nland) ; adj 1961-90 tmn norm anom --> abs

; calculate 1961-1990 synthetic normal from adjusted tmn
print,'Calculating synthetic vap normal'
for iy=nor1,nor2 do begin
 tmpfl=strip(string(tmp_prefix,iy))
 dtrfl=strip(string(dtr_prefix,iy))
 rdbin,tmpgrd,tmpfl,gridsize=2.5,/quiet
 rdbin,dtrgrd,dtrfl,gridsize=2.5,/quiet
 tmpgrd(nland)=tmpgrd(nland)/tmpdtrfac
 dtrgrd(nland)=dtrgrd(nland)/tmpdtrfac
 
; anyzero = where (tmpgrd eq 0.0 or dtrgrd eq 0.0, nanyzero)
; if (nanyzero gt 0) then begin
;   tmpgrd(anyzero)=0.0 & dtrgrd(anyzero)=0.0
; endif
 
 tmn(nland)=tmpgrd(nland)-(dtrgrd(nland)/2.0) + tadj(nland)
 v(nland)=6.108*exp( (17.27*tmn(nland)) /(237.3+tmn(nland)) )>0.05
 tmpgrd(nland)=tmpgrd(nland)+tmpnor(nland)
 v(nland)=v(nland)<esat(tmpgrd(nland))
 vapsyn(nland)=vapsyn(nland)+v(nland)
endfor
vapsyn(nland)=vapsyn(nland)/(nor2-nor1+1)>0.01
vapsyn(nsea) =-9999

;
;
;  Calculate synthetic vap from tmin, convert to anomalies 
;  relative to synthetic mean vap, and apply to normal vap
third:
print,'Calculating synthetic anomalies'

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

 tmpfl=strip(string(tmp_prefix,iy))
 dtrfl=strip(string(dtr_prefix,iy))
 rdbin,tmpgrd,tmpfl,gridsize=2.5,/quiet
 rdbin,dtrgrd,dtrfl,gridsize=2.5,/quiet
 tmpgrd(nland)=tmpgrd(nland)/tmpdtrfac
 dtrgrd(nland)=dtrgrd(nland)/tmpdtrfac
 
; anyzero = where (tmpgrd eq 0.0 or dtrgrd eq 0.0, nanyzero)
; if (nanyzero gt 0) then begin
;   tmpgrd(anyzero)=0.0 & dtrgrd(anyzero)=0.0
; endif
 
 tmn(nland)=tmpgrd(nland)-(dtrgrd(nland)/2.0) + tadj(nland)
 v(nland)=6.108*exp( (17.27*tmn(nland)) / (237.3+tmn(nland)) )>0.05
 tmpgrd(nland)=tmpgrd(nland)+tmpnor(nland)
 v(nland)=v(nland)<esat(tmpgrd(nland))
 v(nland)=v(nland)-vapsyn(nland)
 
 vapmom=moment(v(nland))				; best est moments
 vapmax=max(v(nland),min=vapmin)			; best est min/max
 print, yeartext,' vap',' (x,s2,<<,>>):',vapmom[0],vapmom[1],vapmin,vapmax
 
 vapgrd(nland)=10.0*v(nland)
 vapgrd(nland)=vapgrd(nland)>(-9998)
 vapgrd(nland)=vapgrd(nland)<9998
 vapgrd(nsea)=-9999

 if keyword_set(dumpbin) then $
 	wrbin,fix(round(vapgrd)),strip(string(outprefix,iy))
 
 if keyword_set(dumpglo) then begin
   for imon=0,11 do begin
     vapmonth=vapgrd(*,*,imon)
     if (iy eq year1 and imon eq 0) then monsea=where(vapmonth eq -9999)
     vapmonth=vapmonth/10.0
     vapmonth(monsea)=-999.0
     
     SaveFile=outprefix+monthname(imon)+'.'+yeartext+'.glo'
     SaveTitle='synthetic vap for month '+monthname(imon)+' in year '+yeartext
     SaveGlo,23,vapmonth,CallFile=Savefile,CallTitle=SaveTitle
   endfor
 endif
endfor

if keyword_set(dumpclim) then begin
   for imon=0,11 do begin
     vapmonth=tadj(*,*,imon)
     if (imon eq 0) then monsea=where(vapmonth eq -9999)
     vapmonth(monsea)=-999.0
     
     SaveFile=outprefix+monthname(imon)+'.6190.adj.glo'
     SaveTitle='adj 1961-90 MNew normal for month '+monthname(imon)
     SaveGlo,23,vapmonth,CallFile=Savefile,CallTitle=SaveTitle
   endfor
endif

end