CCSM4 geopotential heights

Does the script below for computing 500hPa geopotential heights from CCSM4 data
look ok? The resulting values are slightly higher than I expected, by maybe 100 m,
but perhaps I am wrong in my concern.

Thanks, Michael

nlat=192
nlon=288

mns=(/"JAN","FEB","MAR","APR","MAY","JUN","JUL","AUG","SEP","OCT","NOV","DEC"/)

a=addfile("/data/Share/regcm_inputdata/CCSM4/historical/ccsm4.shumJAN.1971.nc","r")
lat=a->lat
lon=a->lon
lev=a->lev(::-1)
coefa=a->a(::-1)
coefb=a->b(::-1)
p0=a->p0

b=addfile("/data/Share/regcm_inputdata/CCSM4/historical/orog_fx_CCSM4_historical_r0i0p0.nc","r")
topo=b->orog ; height above sea level (m)

data=new((/32,12,nlat,nlon/),float) ; 1973-2004

do iy=0,31
 do im=0,11

 a=addfile("/data/Share/regcm_inputdata/CCSM4/historical/ccsm4.shum"+mns(im)+"."+(iy+1973)+".nc","r")
 q=a->hus(:,::-1,:,:) ; kg/kg spec hum
 w=q/(1.-q) ; convert spec hum to mixing ratio
 q=(/w/)

 time=a->time
 b=addfile("/data/Share/regcm_inputdata/CCSM4/historical/ccsm4.pres"+mns(im)+"."+(iy+1973)+".nc","r")
 ps=b->ps
 ps=ps/100.
 copy_VarCoords(b->ps,ps)
 c=addfile("/data/Share/regcm_inputdata/CCSM4/historical/ccsm4.air"+mns(im)+"."+(iy+1973)+".nc","r")
 t=c->ta(:,::-1,:,:)

 copy_VarCoords(t,q)

 zsfc=new((/dimsizes(time),nlat,nlon/),float)
 do it=0,dimsizes(time)-1
 zsfc(it,:,:)=topo
 end do
 copy_VarCoords(ps,zsfc)
 t=t*(1.+0.61*q) ; virtual T
 copy_VarCoords(q,t)

 p=new((/dimsizes(time),dimsizes(coefa),nlat,nlon/),float)
 do i=0,dimsizes(coefa)-1
 p(:,i,:,:)=doubletofloat(coefa(i))*p0+doubletofloat(coefb(i))*ps
 end do
 copy_VarCoords(t,p) ; suppose to be hPa

 z=p(time|:,lat|:,lon|:,lev|::-1)
 z=hydro(p(time|:,lat|:,lon|:,lev|::-1),t(time|:,lat|:,lon|:,lev|::-1),zsfc)
 z=z(:,:,:,::-1)

 p1=p(time|:,lat|:,lon|:,lev|:)
 h500=int2p(p1,z,500.,2)
 h500_less=h500(:,:,:,0)
 copy_VarCoords(ps,h500_less)
 v=h500_less(lat|:,lon|:,time|:)
 data(iy,im,:,:)=dim_avg_Wrap(v)

 delete(time)
 delete(a)
 delete(b)
 delete(c)
 delete(q)
 delete(ps)
 delete(t)
 delete(zsfc)
 delete(p)
 delete(z)
 delete(p1)
 delete(h500)
 delete(h500_less)
 delete(v)
 delete(w)

 end do
end do
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Mar 7 21:54:46 2013