pro mean_gts_tdm,prefix,year1,year2,ave=ave,sd=sd,mmperday=mmperday,$
 outprefix=outprefix,uneven=uneven,seas=seas
;calculates the bias in 1961-1990 anomalies, 
; then adjusts individual years by bias to produce
; anomalies that have a 1961-1990 mean of zero
;
if n_params() lt 1 then begin
 print,'mean_gts,prefix,year1,year2,ave=ave,sd=sd,mmperday=mmperday,'+$
  'outprefix=outprefix,uneven=uneven,seas=seas'
 return
endif
;
if keyword_set(seas) then seas=1 else seas=0
count=0.0
for i=year1,year2 do begin
 anomfl=strip(string(prefix,round(i)))
 rdbin,anom,anomfl,seas=seas
 if i ne year1 and keyword_set(uneven) eq 0 then goto,skip
 n_s=n_elements(anom(0,0,*))
 nsea=where(anom eq -9999,nnsea)
 nland=where(anom ne -9999,nnland)
 if i ne year1 then goto,skip
 sum=float(anom*0)
 sum2=float(anom*0)
 count=sum2
 skip:
 anom=float(anom)
 if keyword_set(mmperday) then $
  for im=0,n_s-1 do anom(*,*,im)=10.0*anom(*,*,im)/float(days(im))
 sum(nland)=sum(nland)+anom(nland)
 sum2(nland)=sum2(nland)+anom(nland)^2
 count(nland)=count(nland)+1.0
endfor
anom=0
;
;calculate ave and sd
ave=sum*0.0
sd=ave
nland=where(count gt round(0.75*(year2-year1+1)))
nsea=where(count le round(0.75*(year2-year1+1)),nnsea)
sum=sum(nland)
sum2=sum2(nland)
count=count(nland)
ave(nland)=sum/count
sd(nland)=sum2/(count-1) - $
 (count/(count-1)) * $
 (sum/count)^2
sum=0
sum2=0
count=0
sd(nland)=sqrt(sd(nland))
if nnsea gt 0 then begin
 ave(nsea)=-9999
 sd(nsea)=-9999
endif
;
;write binary grids
y1=year1-1900 & y2=year2-1900
if n_elements(outprefix) eq 0 then outprefix=prefix
avefl=strip(string(outprefix)+string(y1,y2,format='(2i2.2)'))
sdfl=strip(string(outprefix)+string(y1,y2,format='(2i2.2)'))+'.sd'
print,"***** mean_gts_tdm : "+strip(avefl)+" *****"
wrbin,ave,avefl,/compress
print,"***** mean_gts_tdm : "+strip(sdfl)+" *****"
wrbin,sd,sdfl,/compress
end