pro vap_gts_tdm,year1,year2, $
	dtr_ano=dtr_ano,tmp_ano=tmp_ano,dtr_abs=dtr_abs,tmp_abs=tmp_abs,$
    	nor1=nor1,nor2=nor2,outprefix=outprefix

; 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 used to be in percent, whole numbers.
; they are now in hPa ??

if n_params() lt 1 then begin
 print,' vap_gts_tdm,year1,year2,dtr_ano=dtr_ano,tmp_ano=tmp_ano,'
 print,'         dtr_abs=dtr_abs,tmp_abs=tmp_abs,'
 print,'         nor1=nor1,nor2=nor2,outprefix=outprefix'
 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_ano) ne 0 then begin
  rdbin,dtrnor,'/cru/mark1/f080/gts/dtr/means/glo.dtr.6190'
  dtrnor=dtrnor/10.0
endif

if keyword_set(tmp_ano) ne 0 then begin
  rdbin,tmpnor,'/cru/mark1/f080/gts/tmp/means/glo.tmp.6190'
  tmpnor=tmpnor/10.0
endif

rdbin,vapnor,'/cru/mark1/f080/gts/vap/means/glo.vap.6190'

vapnor=vapnor/10.0
;
;
; calculate 1961-1990 tmn normal from monthly tmp and dtr values
;
print,'@@@@@ calc tmn 6190 normal from tmp and dtr @@@@@'
for iy=nor1,nor2 do begin
 print,iy
 tmpfl=strip(string(tmp_abs,iy))
 dtrfl=strip(string(dtr_abs,iy))
 rdbin,tmpgrd,tmpfl,gridsize=2.5
 rdbin,dtrgrd,dtrfl,gridsize=2.5
 
 if iy eq 1961 then begin
  nsea=where(tmpgrd le -9999)
  nland=where(tmpgrd gt -9999)
  tmnnor=float(tmpgrd)*0.0
 endif
 tmnnor(nland)=tmnnor(nland)+(tmpgrd(nland)-(0.5*dtrgrd(nland)))/10.0
endfor
tmnnor(nland)=tmnnor(nland)/(nor2-nor1+1)
vapnor=vapnor(nland)>0.1
tdwnor=tvap(vapnor)
tadj=tdwnor-tmnnor(nland)

; calculate 1961-1990 synthetic normal from adjusted tmn
;
print,'@@@@@ calc tmn 6190 synthetic normal @@@@@'
for iy=nor1,nor2 do begin
 print,iy
 tmpfl=strip(string(tmp_abs,iy))
 dtrfl=strip(string(dtr_abs,iy))
 rdbin,tmpgrd,tmpfl,gridsize=2.5
 rdbin,dtrgrd,dtrfl,gridsize=2.5
 if iy eq 1961 then begin
  vapsyn=float(tmpgrd)*0.0
 endif
 tmn=(tmpgrd(nland)-(0.5*dtrgrd(nland)))/10.0 + tadj
 dtrgrd=1
 v=6.108*exp( (17.27*tmn) /(237.3+tmn) )>0.05
 v=v<esat(tmpgrd(nland)/10.0)
 vapsyn(nland)=vapsyn(nland)+v
 tmpgrd=1 
endfor
vapsyn(nland)=vapsyn(nland)/(nor2-nor1+1)>0.01
;
;
;  Calculate synthetic vap from tmin, convert to anomalies 
;  relative to synthetic mean vap, and apply to normal vap
third:
print,'@@@@@ calc synthetic anomalies @@@@@'
for iy=year1,year2 do begin
 print,iy
 tmpfl=strip(string(tmp_abs,iy))
 dtrfl=strip(string(dtr_abs,iy))
 rdbin,tmpgrd,tmpfl,gridsize=2.5
 rdbin,dtrgrd,dtrfl,gridsize=2.5
 if iy eq year1 then begin
  vapgrd=float(tmpgrd)*0.0
  nsea=where(tmpgrd le -9999)
  nland=where(tmpgrd gt -9999)
 endif
 tmn=tmpgrd(nland)/10.0-0.5*dtrgrd(nland)/10.0 + tadj
 dtrgrd=1
 v=6.108*exp( (17.27*tmn) /(237.3+tmn) )>0.05
 v=v<esat(tmpgrd(nland)/10.0)
 tmpgrd=1
 vapgrd(nland)=v
;
; this is the output percentage anomalies
;   vapgrd(nland)=100*(vapgrd(nland)-vapsyn(nland))/vapsyn(nland)
;
; this is the output hPa anomalies
 vapgrd(nland)=vapgrd(nland)-vapsyn(nland)
; 
 vapgrd(nland)=vapgrd(nland)>(-9998)
 vapgrd(nland)=vapgrd(nland)<9998
 if (nsea ne -1) then vapgrd(nsea)=-9999

 wrbin,fix(round(vapgrd)),strip(string(outprefix,iy)) ; write anomalies

 vapgrd(nland)=(vapgrd(nland)/100.0)*(vapnor*10.0)+(vapnor*10.0)
 wrbin,fix(round(vapgrd)),strip(string(outprefix2,iy))  ; write absolutes
endfor

end