Re: call of function

What is your error?

On Sat, Dec 20, 2008 at 12:25 AM, Jin Young Kim <Jin-Young.Kim_at_noaa.gov>wrote:

> Greeting..
>
> I've define a function and the used it later in attached scripts..
>
> But it's showing following fatal error ..
>
> Would you look at it?
> Do you anybody know its problem?
>
> Best regards,
> Jin-Young Kim
>
> ;*************************************************
> ; traj_1.ncl from www.ncl.ucar.edu
> ;
> ; 2008.12.8 (Thu) PM05:17
> ; 2008.12.19 (Fri) PM01:27 (modified adding tracker infor)
> ;*************************************************
> load "$NCARG_LIB/ncarg/nclscripts/csm/contributed.ncl"
> load "$NCARG_LIB/ncarg/nclscripts/csm/gsn_code.ncl"
> load "$NCARG_LIB/ncarg/nclscripts/csm/gsn_csm.ncl"
> load "$NCARG_LIB/ncarg/nclscripts/csm/shea_util.ncl"
>
> model = getenv ("MID")
> domain = getenv ("DOMAIN")
>
> ;*************************************************
> functionn theta_eqv(T:numeric,X:numeric,p0[1],p:numeric)
> ;
> ; calculate equivalent potential temperature
> ;
> ; T: temperature [K] (lev,lat,lon)
> ; X: mixing ratio of water vapor in air [kg/kg] (lev,lat,lon)
> ; p0: reference pressure [usually 1000mb or 100000Pa]
> ; 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 = 10004 ; specific heat at constant pressure for ait [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 t3
> 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
>
> ;************************************************
> ; read in grib file
> ;************************************************
>
> ifim = sprinti("wrfprs_d02.%0.2i_mod",ispan(0,72,6))
> ofim = sprinti("Wind_v_d02h_%0.2i",ispan(0,72,6))
> titl = sprinti("%0.2i",ispan(0,72,6))
>
> path = "./obs-tracker/"
> f_1 = asciiread(path + "obs.txt",(/12+1,4+1/),"float")
> f_2 = asciiread(path + "gfs.txt",(/12+1,4+1/),"float")
> f_3 = asciiread(path + "hrs.txt",(/12+1,15+1/),"float")
> path = "./obs-diapost/"
> f_4 = asciiread(path + "hrs.txt",(/12+1,7+1/),"float")
>
> olat = f_1(:,4)
> olon = f_1(:,3)
> oprs = f_1(:,2)
> ownd = f_1(:,1)
> glat = f_2(:,4)
> glon = f_2(:,3)
> gprs = f_2(:,2)
> gwnd = f_2(:,1)
> hour = f_3(:,3)
> hlat = f_3(:,4)
> hlon = f_3(:,5)
> hprs = f_3(:,7)
> hwnd = f_4(:,6)
> ; this is for atlancitc ocean hurricane
> do j=0,8
> if ( olon(j) .gt. 0 ) then
> olon = olon * (-1.0)
> ownd = ownd
> end if
> if ( glon(j) .gt. 0 ) then
> glon = glon * (-1.0)
> gwnd = gwnd
> end if
> if ( hlon(j) .gt. 0 ) then
> hlon = hlon * (-0.1)
> hlat = hlat * (0.1)
> hwnd = hwnd
> end if
> end do
> hour = hour*0.01
>
> print (hour)
> print ("lon="+ hlon + "lat="+hlat + "pressure="+hprs + "wind="+hwnd)
> print (hprs(0)-oprs(0))
> print (hwnd(0)-ownd(0))
> print (gprs(0)-oprs(0))
> print (gwnd(0)-ownd(0))
>
>
> do i= 0, dimsizes(ifim)-1
> grb_file = addfile(ifim(i)+".grb","r")
> wks =gsn_open_wks("ps",ofim(i))
>
> print (ifim)
>
> ;************************************************
> ; read in 1d lat/lon arrays
> ;************************************************
>
> lat1d = grb_file->g0_lat_0
> lon1d = grb_file->g0_lon_1
>
> lat_dim = dimsizes(lat1d)
> lon_dim = dimsizes(lon1d)
>
> minlat = min(lat1d)
> maxlat = max(lat1d)
> minlon = min(lon1d)
> maxlon = max(lon1d)
>
> print (min(lat1d))
> print (max(lat1d))
> print (min(lon1d))
> print (max(lon1d))
>
> ;************************************************
> ; Make sure data is ordered from N->S and W->E
> ;************************************************
>
> lat_flag = 0
> lon_flag = 0
>
> if (lat1d(0).gt.lat1d(1)) then
> ; data ordered from S->N
> ; flip lat to go from N->S
> lat1d = lat1d(::-1)
> lat_flag = 1
> end if
>
> if (lon1d(0).gt.lat1d(1)) then
> ; data ordered from E->W
> ; flip lon to go from W->E
> lon1d = lon1d(::-1)
> lon_flag = 1
> end if
> ;************************************************
> ; Create workstation and assign the colormap
> ;************************************************
>
> ; gsn_define_colormap(wks,"gui_default") ;blue - green - red
>
> ; color_maps = NhlPalGetDefined
> ; print(color_maps)
>
> gsn_define_colormap(wks,"BlWhRe") ; blue - whilete - red
> ; gsn_define_colormap(wks,"nrl_sirkes") ; blue - whilete - red
>
> ;************************************************
> ; Get Sfc Wind from grib file
> ;************************************************
>
> sh = grb_file->SPF_H_GDS0_ISBL
> rh = grb_file->R_H_GDS0_ISBL
> T = grb_file->TMP_GDS0_ISBL
> p = T
> do i=0,45
> p(i,:,:) = grb_file->lv_ISBL3(i)
> end do
> ;************************************************
> ; mixhum_ptd, mixhum_ptrh
> ;************************************************
>
>
> ; Reorder met data if done above for lat or lon
> if (lat_flag.eq.1) then
> ; flip lat dimension to go from N->S
> sh = sh(::-1,:)
> rh = rh(::-1,:)
> T = T(::-1,:)
> end if
>
> if (lon_flag.eq.1) then
> ; flip lon dimension to go from W->E
> sh = sh(:,::-1)
> rh = rh(:,::-1)
> T = T(:,::-1)
> end if
>
> X=sh
> X=X/(1-X); convert spc humidity to mixing ratio
> X_at_long_name="mixing ratio"
>
> ; X = mixhum_ptrh(p,T,rh,-1) ; mix ratio (g/kg)
> p0 = 1000
> p0_at_units="hPa"
>
> THETAE=theta_eqv(T,X,p0,p)
>
> ;************************************************
> ; Set map resources for plot
> ;************************************************
>
> res = True
>
> res_at_gsnFrame = False
> res_at_gsnAddCyclic = False
> res_at_trYReverse = True
> res_at_gsnMaximize = True
>
> res_at_gsnRightString = "Equiv. Potential Temperature (K)"
> res_at_tiMainString = titl(i) + "h Fcst.";title
> ;************************************************
> ; Set isotac contour resources.
> ;************************************************
> res_at_cnFillOn = True ; color plot desired
> res_at_gsnSpreadColors = True
> res_at_gsnScalarContour = True ; contours desired
> res_at_cnInfoLabelOn = False
>
> ; res_at_cnLevelSelectionMode = "ManualLevels"
> ; res_at_cnMinLevelValF = -70
> ; res_at_cnMaxLevelValF = 70
> ; res_at_cnLevelSpacingF = 10
>
> ;************************************************
> ; Draw everything for final image
> ;************************************************
>
> tc= hlat(i)
>
> map = gsn_csm_scalar(wks,THETAE(:,{tc},:),res)
>
> ;************************************************
> ; Draw everything for final image
> ;************************************************
> draw(map)
> end do
> end
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>
>

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Fri Dec 19 2008 - 22:33:22 MST