Re: Help on computing equvalent potential temperature (fwd)

From: Dennis Shea <shea_at_nyahnyahspammersnyahnyah>
Date: Tue, 24 Apr 2007 11:54:51 -0600 (MDT)

>
> I want to calculate equivalent potential temperature with NCEP-NCAR
; 4 dimension dataset (longitude, latitude, level, time).
>
> I want to know if there is a NCL function which can be directly
> use to perform this computation ?
> ---------------------------------

There is no user accessible built-in function in NCL.

As I recall, NCEP-NCAR Reanalysis contain specific humidity.
I think the following should work. this is not tested.

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

function theta_eqv(T:numeric, X:numeric, p0[1], p[*]:numeric)
;
; calculate equivalent potential temperature
;
; T : Temperature [K] (time,lev,lat,lon)
; X : mixing ratio of water vapor in air [kg/kg] (time,lev,lat,lon)
; p0 : reference pressure [ usually 1000 mb or 100000 Pa ]
; p ; pressure at each level [must be same units as p0 ]
local T_e, L_v, c_pd, R, Rcpd, dimT, dimX, P
begin
   L_v = 2400. ; latent heat of evaporation [approx 2400 kJ/kg at 25C]
   c_pd = 1004. ; specific heat at constant pressure for air [approx 1004 J/(kg-K)]
   R = 287. ; specific gas constant for air [J/(kg-K)]
   Rcpd = R/c_pd

   dimT = dimsizes(T)
   dimX = dimsizes(X)
   if (dimsizes(dimT) .ne. dimsizes(dimX)) then
       print("theta_eqv: rank of T .ne. rank of X")
       exit
   end if

   P = conform(T,p,1) ; make p same shape/rank/size as T
   T_e = T + (L_v/c_pd)*X ; common approximation
   theta_e = T_e*(p0/P)^Rcpd

   theta_e_at_long_name = "equivalent potential temperature"
   theta_e_at_units = "K"
   copy_VarCoords(T, theta_e) ; assign coordinates

   return(theta_e)
end

begin
   f = addfile (...)
   p = f->level ; [*] .... not sure of name
   T = f->T
   X = f->Q ; spc humidity
   X = X/(1-X) ; convert spc humidity to mixing ratio
   X_at_long_name = "mixing ratio"

   p0 = 1000
   p0_at_units = "hPa" ; must match units of p

   THETAE = theta_eqv(T, X, p0, p)

   printVarSummary( THETAE )
   printMinMax ( THETAE, True )

end

Good luck
D
_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Tue Apr 24 2007 - 11:54:51 MDT

This archive was generated by hypermail 2.2.0 : Thu Apr 26 2007 - 08:55:37 MDT