;
; Tests Thornthwaite's potential evaporation formula across all of
; Russia using gridded monthly data.  For 1931-1960 means!
;
; Read in the land air temperature dataset (has to be Mark New's, since
; we need absolute values!)
;
print,'Reading temperature baseline'
restore,'/cru/u2/f055/data/obs/grid/surface/lat_new_19611990.idlsave'
;  g,ltm5,nmon
;
print,'Reading temperature anomalies'
ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/lat_new_19011995.mon.nc')
tmon=crunc_rddata(ncid,tst=[1931,0],tend=[1960,11],grid=g)
ncdf_close,ncid
;
; Compute long term mean anomalies for each month of the 1931-60 period
;
print,'Computing long term means'
nmon=12
nyr=g.nt/nmon
tmon=reform(tmon,g.nx,g.ny,nmon,nyr)
fdtot=total(tmon,/nan,4)
fdnum=total(finite(tmon),4)
ml=where(fdnum lt 20.,nmiss)
fdtot=fdtot/float(fdnum)
if nmiss gt 0 then fdtot(ml)=!values.f_nan
;
; Add back on the 1961-90 means (month by month)
;
fdtot(*,*,*)=fdtot(*,*,*)+ltm5(*,*,*)
;
; Put into Thornthwaite (if no missing data)
;
annpe=fltarr(g.nx,g.ny)*!values.f_nan
for ix = 0 , g.nx-1 do begin
  print,ix,format='($,I2)'
  for iy = 0 , g.ny-1 do begin
    onetemp=reform(fdtot(ix,iy,*))
    if finite(total(onetemp)) then begin
      monpe=thornthwaite(onetemp,g.y(iy))
      annpe(ix,iy)=total(monpe)
    endif
  endfor
endfor
;
; Plot russian region map
;
loadct,39
multi_plot,nrow=2
if !d.name eq 'X' then begin
  window,ysize=800
  !p.font=-1
endif else begin
  !p.font=0
  device,/helvetica,/bold,font_size=9
endelse
def_1color,80,color='deepblue'
def_1color,84,color='lgreen'
def_1color,88,color='red'
def_smearcolor,fromto=[80,84]
def_smearcolor,fromto=[84,88]
map=def_map(/npolar) & map.limit(0)=25.
coast=def_coast(/get_device) & coast.coasts=1 & coast.on=2
;
levs=[0.,125.,250.,375.,500.,625.,875.,1250.,1750.,2250.]
cols=findgen(9)+80.
;
scale_vert,levels=levs,c_colors=cols
inter_confd,annpe,g.x,g.y,map=map,coast=coast,$
  /cell_fill,levels=levs,miss_grey='white',$
  /hi_on,c_colors=cols
scale_vert,/off
;
end