Hello,
[1] It is likely a memory issue. You script ran just fine on our systems.
Try the attached script. If it fails you may have to run on a few
months, create a netCDF file [time =>unlimites); then use
ncrcat to concatenate the files.
[2] You could improve your memory management by
reusing variable space.
[3] FYI: you do not need to read netCDF, you can read the grib
directly.
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
begin
f = addfile("U78614.grb","r")
qf = addfile("U78615.grb","r")
P0 = 100000.
Rd = 287.
Rv = 461.
T = f->T_GDS4_HYBL_123
Q = qf->Q_GDS4_HYBL_123
T = T*(1+(Rv/Rd-1)*Q) ; compute virtual temperature
T_at_long_name = "Temperature (virtual)"
delete(Q) ; delete to save memory
Zs = f->Z_GDS4_HYBL_123
ps = f->LNSP_GDS4_HYBL_123
ps = exp(ps) ; replace
hyai = f->lv_HYBL_i4_a ;interface hybrid coefficient A
hybi = f->lv_HYBL_i4_b ;interface hybrid coefficient B
hyam = f->lv_HYBL3_a ;mid-layer hybrid coefficient A
hybm = f->lv_HYBL3_b ;mid-layer hybrid coefficient B
z2 = cz2ccm(ps,Zs,T,P0,hyam(::-1),hybm(::-1),hyai(::-1),hybi(::-1))
copy_VarCoords(T,z2)
z2_at_long_name = "geopotential height"
z2_at_units = "gpm"
;---- Create 'nominal' pressure level coordinate and
; replace the integer hybrid level coordinates.
; This is what the CAM model does.
delete(z2&lv_HYBL3) ; delete the integer values
lev = (hyam+hybm)*P0*0.01 ; calculate nominal pressure levels
lev!0= "lev"
lev_at_units = "hPa"
z2!1 = "lev" ; assign new name and p values
z2&lev = lev
printVarSummary( z2 )
; print one vertical profile over the ocean
; Use nominal pressure levels for p
print (z2&lev+" "+z2(0,:,{45},{180}))
end
yihuiw_at_uci.edu wrote:
> Hi,
> I try to calculate geopotential height by using function cz2ccm but I have
> some questions about it.
> I used ERA40 data and there are four hybrid coefficients, hybrid A and B
> coefficient at layer midpoints, and hybrid A and B coefficient at layer
> interface. In the ncl documentation, it requires hyam and hybm, which are
> one-dimensional array of hybrid coordinate coefficients for base pressures
> (layer midpoints) ordered bottom-to-top, and hyai and hybi, which are
> one-dimensional array of hybrid A and B coefficients (layer midpoints)
> ordered bottom-to-top. All hybrid coefficients in ncl documentation need
> to be at layer midpoints while two coefficients in ERA40 are at layer
> interface. Do conefficients in ERA40 satisfy the requirements in ncl
> documentation?
> I used hybrid A and B coefficient at layer midpoints in ERA40 as hyam and
> hybm, and hybrid A and B coefficient at layer interface as hyai and hybi.
> But I am getting a “segmentation fault error. The file is monthly mean
> data in one year with spatial resolution 128*256. Is it the problem of
> memory? How can I get past it? I attached the script below.
> Can I have the script of function cz2ccm because I am interested in how to
> calculate geopotential height in ncl? Thank you so much.
> All the best,
> Yi-Hui
>
> **************************************************************
>
> load "/usr/lib/ncarg/ncarg/nclscripts/csm/contributed.ncl"
> begin
> f=addfile("U78614.nc","r")
> qf=addfile("U78615.nc","r")
> P0=100000
> Rd=287
> Rv=461
>
> T=f->T_GDS4_HYBL_123(:,:,:,:)
> Zs=f->Z_GDS4_HYBL_123(:,:,:)
> Q=qf->Q_GDS4_HYBL_123(:,:,:,:)
> lnps=f->LNSP_GDS4_HYBL_123(:,:,:)
> ps=exp(lnps)
> hyai=f->lv_HYBL_i4_a ;interface hybrid coefficient A
> hybi=f->lv_HYBL_i4_b ;interface hybrid coefficient B
> hyam=f->lv_HYBL3_a ;mid-layer hybrid coefficient A
> hybm=f->lv_HYBL3_b ;mid-layer hybrid coefficient B
>
> Tv=new((/12,60,128,256/),float)
> Tv=T*(1+(Rv/Rd-1)*Q)
> copy_VarCoords(T,Tv)
> z2=cz2ccm(ps,Zs,Tv,P0,hyam(::-1),hybm(::-1),hyai(::-1),hybi(::-1))
>
> end
>
>
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk_at_ucar.edu
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
-- ====================================================== Dennis J. Shea tel: 303-497-1361 | P.O. Box 3000 fax: 303-497-1333 | Climate Analysis Section | Climate & Global Dynamics Div. | National Center for Atmospheric Research | Boulder, CO 80307 | USA email: shea 'at' ucar.edu | ====================================================== _______________________________________________ ncl-talk mailing list ncl-talk_at_ucar.edu http://mailman.ucar.edu/mailman/listinfo/ncl-talkReceived on Mon Sep 08 2008 - 14:43:14 MDT
This archive was generated by hypermail 2.2.0 : Mon Sep 08 2008 - 14:49:02 MDT