;
; Tests Ivanov'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
;
; Read in the land vapour pressure dataset
;
print,'Reading vapour pressures'
ncid=ncdf_open('/cru/u2/f055/data/obs/grid/surface/vap_new_19011995.mon.nc')
vmon=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,
; for temperature and vapour pressure
;
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(*,*,*)
;
vmon=reform(vmon,g.nx,g.ny,nmon,nyr)
vdtot=total(vmon,/nan,4)
vdnum=total(finite(vmon),4)
ml=where(vdnum lt 20.,nmiss)
vdtot=vdtot/float(vdnum)
if nmiss gt 0 then vdtot(ml)=!values.f_nan
;
; 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,*))
    onevapo=reform(vdtot(ix,iy,*))
    if finite(total(onetemp+onevapo)) then begin
      monpe=ivanov(onetemp,onevapo,/vapour)
      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() & map.origx=90. & map.limit=[30.,20.,85.,180.]
;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.
;
!p.multi(0)=1
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