call of function

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
Received on Fri Dec 19 2008 - 22:25:55 MST