unable to run a ncl script

Dear Sir,

I am trying to calculate :Equivalent Pot. Temp variation". when I am running
script like: $ncl test-eq.ncl It is hanging and hanging. Nothing is
coming out.

Kindly help me to sort it.

Following is the script (file name: test-eq.ncl):
- ---- ---- ---- ------ ------ ------- ----------- ---------
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin

datadir = "./wrfout/"
FILES = systemfunc ("ls -1 " + datadir + "wrfout_d01* " )
numFILES = dimsizes(FILES)

type = "ps"

res = True
res@MainTitle = "Eq. Pot. Temp. Variation"
pltres = True
mpres = True

opts = res
opts@cnFillOn = True

do ifil = 0,numFILES-1
   a = addfile(FILES(ifil)+".nc","r")

; What times and how many time steps are in the data set?
  times = wrf_user_list_times(a) ; get times in the file
  ntimes = dimsizes(times) ; number of times in the file

  do it = 0,ntimes-1,2 ; TIME LOOP

    print("Working on time: " + times(it) )
    res@TimeLabel = times(it) ; Set Valid time to use on plots

day = systemfunc("echo "+times(it) +"|cut -c6-7")
yr = systemfunc("echo "+times(it) +"|cut -c1-4")
mon = systemfunc("date '+%b'|tr [a-z] [A-Z]")
hr = systemfunc("echo "+times(it) +"|cut -c12-13")
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need
    T = wrf_user_getvar(a,"T",it)
    P = wrf_user_getvar(a,"P",it)
    PB = wrf_user_getvar(a,"PB",it)
    qv = wrf_user_getvar(a,"QVAPOR",it)

P = P + PB

eth = wrf_eth(qv,T,P)
plevs = 850.
p = wrf_user_getvar(a,"pressure",it) ; pressure is our vertical coordinate
eth_plane = wrf_user_intrp3d(eth,p,"h",plevs,0.,False)
wks = gsn_open_wks(type,"eth_"+"hPa"+"_"+day+mon+yr+"_"+hr+"Z")

; Plotting options for Equivalent Pot Temp.
        opts = res
        opts@cnFillOn = True
  res@cnLevelSelectionMode = "ManualLevels" ; manually set the contour
levels with the following 3 resources
  res@cnMinLevelValF = -60. ; set the minimum contour level
  res@cnMaxLevelValF = 60. ; set the maximum contour level
  res@cnLevelSpacingF = 5. ; set the interval between
contours

        opts@mpOutlineBoundarySets = "AllBoundaries"
        mpres@mpOutlineBoundarySets = "National"
        mpres@mpOutlineBoundarySets = "AllBoundaries"
        mpres@mpNationalLineThicknessF = 1.5
        mpres@mpNationalLineColor = "black"
        mpres@mpGeophysicalLineColor = "black"
        mpres@mpGeophysicalLineThicknessF = 1.5

contour_eth = wrf_contour(a,wks,eth_plane,opts)
 if ( pressure .eq. 850 ) then
          plot = wrf_map_overlays(a,wks,contour_eth,pltres,mpres)
        end if

end do
end do
end
 - - - - - -- - -- - -- - -- - - -- - - - - - - - - - - - - - - - - - - - -
- - - - - - - - -

 My WRF out files are look like this:

$ncdump -h wrfout/wrfout_d01_2010-02-16_00\:00\:00 |more

netcdf wrfout_d01_2010-02-16_00\:00\:00 {
dimensions:
        Time = UNLIMITED ; // (1 currently)
        DateStrLen = 19 ;
        west_east = 369 ;
        south_north = 307 ;
        west_east_stag = 370 ;
        bottom_top = 26 ;
        south_north_stag = 308 ;
        bottom_top_stag = 27 ;
        soil_layers_stag = 5 ;
variables:
        char Times(Time, DateStrLen) ;
float P(Time, bottom_top, south_north, west_east) ;
                P:FieldType = 104 ;
                P:MemoryOrder = "XYZ" ;
                P:description = "perturbation pressure" ;
                P:units = "Pa" ;
                P:stagger = "" ;
                P:coordinates = "XLONG XLAT" ;
        float PB(Time, bottom_top, south_north, west_east) ;
                PB:FieldType = 104 ;
                PB:MemoryOrder = "XYZ" ;
                PB:description = "BASE STATE PRESSURE" ;
                PB:units = "Pa" ;
                PB:stagger = "" ;
                PB:coordinates = "XLONG XLAT" ;

etc. etc

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Feb 25 00:21:24 2010