pro rd0_gts,year1,year2,nor1,nor2,outprefix=outprefix,$
 prenorm=prenorm,rd0norm=rd0norm,synnorm,pre_prefix=pre_prefix
 
; calculates wet-day frequency anomaly time series as function of
;  monthly precip, and precip amd rd0 normal climatologies
if n_params() lt 4 then begin
 print,'rd0_gts,year1,year2,nor1,nor2,outprefix=outprefix'
 print,'prenorm=prenorm,rd0norm=rd0norm,synnorm,pre_prefix=pre_prefix'
 return
endif
close,/all
if keyword_set(outprefix) eq 0 then outprefix='glo.rd0.ano.'
if keyword_set(pre_prefix) eq 0 then pre_prefix='~/m1/gts/pre/glo2/glo.pre.'
;
;read in pre and rd0 normals
print,'Reading precip and rd0 normals'
if (keyword_set(prenorm) eq 0 or (keyword_set(prenorm) ne 0 and $
 n_elements(prenorm) ne 720*360.0*12)) then $    
 rdbin,prenorm,'~/m1/gts/pre/glo2/glo.pre.norm'
if (keyword_set(rd0norm) eq 0 or (keyword_set(rd0norm) ne 0 $
 and n_elements(rd0norm) ne 720*360.0*12)) then $    
 rdbin,rd0norm,'~/m1/gts/rd0/glo2/glo.rd0.norm'
nland=where(rd0norm gt -9999 and prenorm gt -9999)
nsea=where(rd0norm eq -9999 or prenorm eq -9999)
rd0norm(nland)=(rd0norm(nland)/10)>0.49
prenorm(nland)=prenorm(nland)>5.0
monthgrd=rd0norm
for im=0,11 do monthgrd(*,*,im)=days(im)
;
;
if n_elements(synnorm) gt 1 then begin
 rd0syn=synnorm
 goto,third
endif
; second, calculate 1961-1990 synthetic normal
print,'Calculating synthetic Rd0 normal'
for iy=nor1,nor2 do begin
 print,iy
 prefl=strip(string(pre_prefix,iy))
 rdbin,pregrd,prefl
 if iy eq nor1 then begin
  rd0syn=float(pregrd)*0.0
 endif
 npre=where(pregrd ge 0)
 rd=(rd0cal(pregrd(npre),rd0norm(npre),prenorm(npre)))>0.4
 rd0syn(npre)=rd0syn(npre)+(rd<monthgrd(npre))
endfor
rd0syn(nland)=rd0syn(nland)/(nor2-nor1+1)
rd0syn(nsea)=-999.9
;synnorm=rd0syn
;wrbin,synnorm*10,'synnorm',/compr
;return
;
;
; Third, calculate synthetic rd0, convert to anomalies 
;  relative to synthetic mean rd0, and apply to normal rd0
third:
print,'Calculating synthetic anomalies'
for iy=year1,year2 do begin
 print,iy
 prefl=strip(string(pre_prefix,iy))
 rdbin,pregrd,prefl
 rd0grd=float(pregrd)*0.0
 npre=where(pregrd ge 0)
 rd=(rd0cal(pregrd(npre),rd0norm(npre),prenorm(npre)))>0.4
 rd0grd(npre)=rd<monthgrd(npre)
 rd0grd(npre)=(rd0grd(npre)-rd0syn(npre))/rd0syn(npre)
; rd0grd=rd0grd<monthgrd
; rd0grd=rd0grd>0
  rd0grd(nsea)=-99.99
 wrbin,rd0grd*100,strip(string(outprefix,iy))
endfor
end