pro rd0_gts_anom,year1,year2,nor1=nor1,nor2=nor2,outprefix=outprefix,$
 	pre_prefix=pre_prefix,dumpbin=dumpbin,dumpglo=dumpglo
 
; calculates wet-day frequency anomaly time series as function of
;  monthly precip, and precip amd rd0 normal climatologies
if n_params() lt 2 then begin
 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.rd0.ano.'
if keyword_set(pre_prefix) eq 0 then pre_prefix='~/m1/gts/pre/glo2/glo.pre.'
;
monthname=strarr(12)
monthname[*]=['01','02','03','04','05','06','07','08','09','10','11','12']
rd0month=fltarr(144,72)
;
;read in pre and rd0 normals
print,'Reading precip and rd0 normals'
rdbin,prenorm,'/cru/tyn1/f014/_keep/glo25.pre.6190',gridsize=2.5,/quiet
rdbin,rd0norm,'/cru/tyn1/f014/_keep/glo25.rd0.6190',gridsize=2.5,/quiet

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)

rd0syn=float(prenorm)*0.0 & rd0syn(nsea)=-9999
rd    =float(prenorm)*0.0

; calculate 1961-1990 synthetic normal
print,'Calculating synthetic rd0 normal'

for iy=nor1,nor2 do begin
 prefl=strip(string(pre_prefix,iy))
 rdbin,pregrd,prefl,gridsize=2.5,/quiet
 pregrd(nland)=((pregrd(nland)/100.0)+1.0)*prenorm(nland) ; make pre anom into abs
 pregrd(nsea) =-9999
 npre=where(pregrd ge 0)
 
 rd(npre)=(rd0cal(pregrd(npre),rd0norm(npre),prenorm(npre)))>0.4
 rd0syn(npre)=rd0syn(npre)+(rd(npre)<monthgrd(npre))
endfor

landzeros=where(rd0syn eq 0,nlandzeros)
if (nlandzeros gt 0) then rd0syn(landzeros)=1

rd0syn(nland)=rd0syn(nland)/(nor2-nor1+1)
rd0syn(nsea)=-9999

; calculate synthetic rd0, convert to anomalies 
;  relative to synthetic mean rd0, and apply to normal rd0
print,'Calculating synthetic rd0 anomalies'

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

 prefl=strip(string(pre_prefix,iy))
 rdbin,pregrd,prefl,gridsize=2.5,/quiet
 pregrd(nland)=((pregrd(nland)/100.0)+1.0)*prenorm(nland) ; make pre anom into abs
 pregrd(nsea) =-9999
 npre=where(pregrd ge 0)
 
 rd0grd=float(pregrd)*0.0 & rd0grd(nsea)=-9999
 rd(npre)=(rd0cal(pregrd(npre),rd0norm(npre),prenorm(npre)))>0.4
 rd0grd(npre)=rd(npre)<monthgrd(npre)
 
 rd0grd(nland)=100.0*(rd0grd(nland)-rd0syn(nland))/rd0syn(nland)

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

end