Re: questions about cz2ccm

From: <yihuiw_at_nyahnyahspammersnyahnyah>
Date: Mon, 8 Sep 2008 22:35:23 -0700 (PDT)

Hi,
Thank you for your help. I use the function cz2ccm to calculate
geopotential height successfully.
Can I have the original codes of function cz2ccm? I wrote a script to
calculate geopotential height by myself. The pattern of geopotential
height is almost the same while its value is 5% larger than the result by
using cz2ccm. So I'd like to know how does cz2ccm calculate geopotential
height. Thank you again.
All the best,
Yi-Hui

> 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

_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Sep 08 2008 - 23:35:23 MDT

This archive was generated by hypermail 2.2.0 : Wed Sep 10 2008 - 07:38:41 MDT