Re: Script for cross section

Hi Jonathan,

I think the error probably has something to do with the WRFUserARW.ncl
you are using, because wrf_contour is defined in the WRFUserARW.ncl
that comes installed with NCL (as of V4.3.1).

Is there a reason why you need to load your own WRFUserARW.ncl file?

--Mary

On Sep 26, 2008, at 4:48 PM, Jonathan Smith wrote:

> To whom this may concern:
>
> When I run wrf_CrossSection4.ncl at line 121 I recieve an error.
> The error says that the function wrf_contour is not defined. I am
> not sure if it has something to do with the WRFUserARW.ncl that I am
> using. Can anyone take a look the scripts. I using a wrf output in
> this process.
>
> Thanks
>
> Jonathan
> ; Example script to produce plots for a WRF real-data run,
> ; with the ARW coordinate dynamics option.
> ; Plot data on a cross section
> ; This script will plot data at a set angle through a specified
> point
>
> load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
> ;load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
> load "/ptmp/jwsmith/NCLscripts/WRFUserARW.ncl"
>
> begin
> ;
> ; The WRF ARW input file.
> ; This needs to have a ".nc" appended, so just do it.
> dirname = "/ptmp/jwsmith/WRFChem3/run/"
> filename = "wrfout_d01_2006-05-01_00:00:00.nc"
> print(" reading from file "+dirname+filename)
> a = addfile (dirname + filename, "r")
>
>
> ; We generate plots, but what kind do we prefer?
> ; type = "x11"
> ; type = "pdf"
> type = "ps"
> ; type = "ncgm"
> wks = gsn_open_wks(type,"plt_CrossSection4")
>
>
> ; Set some basic resources
> res = True
> res_at_MainTitle = "REAL-TIME WRF"
> res_at_Footer = False
>
> pltres = True
>
>
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>
> ; What times and how many time steps are in the data set?
> FirstTime = True
> FirstTimeMap = True
> times = wrf_user_list_times(a) ; get times in the file
> ntimes = dimsizes(times) ; number of times in the file
>
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>
> xlat = wrf_user_getvar(a, "XLAT",0)
> xlon = wrf_user_getvar(a, "XLONG",0)
> ter = wrf_user_getvar(a, "HGT",0)
>
> do it = 0,ntimes-1,2 ; TIME LOOP
>
> print("Working on time: " + times(it) )
> res_at_TimeLabel = times(it) ; Set Valid time to use on plots
>
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
> ; First get the variables we will need
>
> tc = wrf_user_getvar(a,"tc",it) ; T in C
> rh = wrf_user_getvar(a,"rh",it) ; relative humidity
> z = wrf_user_getvar(a, "z",it) ; grid point height
>
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>
> do ip = 1, 3 ; we are doing 3 plots
> ; all with the pivot point (plane) in the center of the domain
> ; at angles 0, 45 and 90
> ; |
> ; |
> ; angle=0 is |, angle=90 is ------
> ; |
> ; |
>
> ; Build plane (pivot point) through which the cross section will go
> ; OR set to zero, if start and end points are specified
> ; IF using plane and angle, set opts in wrf_user_intrp3d to False
> dimsrh = dimsizes(rh)
> plane = new(2,float)
> plane = (/ dimsrh(2)/2, dimsrh(1)/2 /) ; pivot point is
> center of domain
> opts = False
>
> if(ip .eq. 1) then
> angle = 90.
> X_plane = wrf_user_intrp2d(xlon,plane,angle,opts)
> X_desc = "longitude"
> end if
> if(ip .eq. 2) then
> angle = 0.
> X_plane = wrf_user_intrp2d(xlat,plane,angle,opts)
> X_desc = "latitude"
> end if
> if(ip .eq. 3) then
> angle = 45.
> X_plane = wrf_user_intrp2d(xlon,plane,angle,opts)
> X_desc = "longitude"
> end if
>
> ; X-axis lables
> dimsX = dimsizes(X_plane)
> xmin = X_plane(0)
> xmax = X_plane(dimsX(0)-1)
> xspan = dimsX(0)-1
> nxlabs = floattoint( (xmax-xmin)/2 + 1)
>
> if (FirstTimeMap) then
> lat_plane = wrf_user_intrp2d(xlat,plane,angle,opts)
> lon_plane = wrf_user_intrp2d(xlon,plane,angle,opts)
> mpres = True
> mpres_at_mpGeophysicalLineColor = "Black"
> mpres_at_mpNationalLineColor = "Black"
> mpres_at_mpUSStateLineColor = "Black"
> mpres_at_mpGridLineColor = "Black"
> mpres_at_mpLimbLineColor = "Black"
> mpres_at_mpPerimLineColor = "Black"
> pltres = True
> pltres_at_FramePlot = False
> optsM = res
> optsM_at_NoHeaderFooter = True
> optsM_at_cnFillOn = True
> optsM_at_lbTitleOn = False
> contour = wrf_contour(a,wks,ter,optsM)
> plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)
> lnres = True
> lnres_at_gsLineThicknessF = 3.0
> lnres_at_gsLineColor = "Red"
> do ii = 0,dimsX(0)-2
> gsn_polyline(wks,plot,(/lon_plane(ii),lon_plane(ii+1)/),(/
> lat_plane(ii),lat_plane(ii+1)/),lnres)
> end do
> frame(wks)
> delete(lon_plane)
> delete(lat_plane)
> pltres_at_FramePlot = True
> end if
>
>
> if (FirstTime) then
> ; THIS IS NEEDED FOR LABLES - ALWAYS DO (Z axis only needed
> once. X everytime plane changes)
>
> ; Y-axis labels
> zmax = 6000. ; We only want to see the first 6 km
> zz = wrf_user_intrp3d(z,z,"v",plane,angle,opts)
> dims = dimsizes(zz)
> do imax = 0,dims(0)-1
> if ( .not.ismissing(zz(imax,0)) .and. zz(imax,0) .lt. zmax )
> then
> zmax_pos = imax
> end if
> end do
> zspan = zmax_pos
> zmin = z(0,0,0)
> zmax = zz(zmax_pos,0)
> zmin=zmin/1000.
> zmax=zmax/1000.
> nzlabs = floattoint(zmax + 1)
>
> FirstTime = False
> ; END OF ALWAYS DO
> end if
>
> ; Interpolate data vertically (in z)
>
> rh_plane = wrf_user_intrp3d(rh,z,"v",plane,angle,opts)
> tc_plane = wrf_user_intrp3d(tc,z,"v",plane,angle,opts)
>
> ; Options for XY Plots
> opts_xy = res
> opts_xy_at_tiXAxisString = X_desc
> opts_xy_at_tiYAxisString = "Height (km)"
> opts_xy_at_cnMissingValPerimOn = True
> opts_xy_at_cnMissingValFillColor = 0
> opts_xy_at_cnMissingValFillPattern = 11
> opts_xy_at_tmXTOn = False
> opts_xy_at_tmYROn = False
> opts_xy_at_tmXBMode = "Explicit"
> opts_xy_at_tmXBValues = fspan(0,xspan,nxlabs) ;
> Create the correct tick marks
> opts_xy_at_tmXBLabels = sprintf("%.
> 1f",fspan(xmin,xmax,nxlabs)); Create labels
> opts_xy_at_tmXBLabelFontHeightF = 0.015
> opts_xy_at_tmYLMode = "Explicit"
> opts_xy_at_tmYLValues = fspan(0,zspan,nzlabs) ;
> Create the correct tick marks
> opts_xy_at_tmYLLabels = sprintf("%.
> 1f",fspan(zmin,zmax,nzlabs)); Create labels
> opts_xy_at_tiXAxisFontHeightF = 0.020
> opts_xy_at_tiYAxisFontHeightF = 0.020
> opts_xy_at_tmXBMajorLengthF = 0.02
> opts_xy_at_tmYLMajorLengthF = 0.02
> opts_xy_at_tmYLLabelFontHeightF = 0.015
> opts_xy_at_PlotOrientation = tc_plane_at_Orientation
>
>
> ; Plotting options for RH
> opts_rh = opts_xy
> opts_rh_at_ContourParameters = (/ 10., 90., 10. /)
> opts_rh_at_pmLabelBarOrthogonalPosF = -0.1
> opts_rh_at_cnFillOn = True
> opts_rh_at_cnFillColors = (/"White","White","White", \
>
> "White","Chartreuse","Green", \
> "Green3","Green4", \
>
> "ForestGreen","PaleGreen4"/)
>
> ; Plotting options for Temperature
> opts_tc = opts_xy
> opts_tc_at_cnInfoLabelZone = 1
> opts_tc_at_cnInfoLabelSide = "Top"
> opts_tc_at_cnInfoLabelPerimOn = True
> opts_tc_at_cnInfoLabelOrthogonalPosF = -0.00005
> opts_tc_at_ContourParameters = (/ 5. /)
>
>
> ; Get the contour info for the rh and temp
> contour_tc = wrf_contour(a,wks,tc_plane(0:zmax_pos,:),opts_tc)
> contour_rh = wrf_contour(a,wks,rh_plane(0:zmax_pos,:),opts_rh)
>
>
> ; MAKE PLOTS
> plot = wrf_overlays(a,wks,(/contour_rh,contour_tc/),pltres)
>
> ; Delete options and fields, so we don't have carry over
> delete(opts_xy)
> delete(opts_tc)
> delete(opts_rh)
> delete(tc_plane)
> delete(rh_plane)
> delete(X_plane)
>
> end do ; make next cross section
>
> ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
>
> FirstTimeMap = False
> end do ; END OF TIME LOOP
>
> end
> ; function wrf_user_set_xy( var:numeric, xp:numeric, yp:numeric,
> x1:numeric, \
> ; y1:numeric,angle:numeric ,opts )
> ; function wrf_user_intrp3d( var3d:numeric, z:numeric,
> plot_type:string, \
> ; loc_param:numeric, angle:numeric,
> opts:logical )
> ; function wrf_user_intrp2d( var2d:numeric, \
> ; loc_param:numeric, angle:numeric,
> opts:logical )
> ; function wrf_user_getvar( nc_file:file, variable:string,
> time:integer )
> ; function wrf_user_list_times( nc_file:file )
> ; function wrf_user_latlon_to_ij( nc_file:file, latitude:numeric,
> ; longitude:numeric )
> ; function wrf_contour(nc_file:file,wks[1]: graphic, data[*]
> [*]:numeric, \
> ; opt_args[1]:logical)
> ; function wrf_vector(nc_file:file,wks[1]: graphic, data_u[*]
> [*]:numeric, \
> ; data_v[*][*]:numeric, opt_args[1]:logical)
> ; function wrf_map(wks[1]:graphic,in_file[1]:file,opt_args[1]:logical)
> ; function
> wrf_map_overlays(in_file[1]:file,wks:graphic,base[1]:graphic, \
> ;
> plots[*]:graphic,opt_arg[1]:logical,mp_arg[1]:logical)
> ; function wrf_overlays(in_file[1]:file,wks:graphic,
> plots[*]:graphic, \
> ; opt_arg[1]:logical)
> ;
> ; These next three functions are obsolete as of version 5.0.1 of NCL.
> ; Do not use them. Use wrf_overlays and wrf_map_overlays instead.
> ;
> ; function
> wrf_map_zoom(wks[1]:graphic,in_file[1]:file,opt_args[1]:logical, \
> ; y1:integer,y2:integer,x1:integer,x2:integer)
> ; procedure wrf_map_overlay(wks:graphic,base[1]:graphic, \
> ; plots[*]:graphic, \
> ; opt_arg[1]:logical)
> ; procedure wrf_overlay(wks:graphic, plots[*]:graphic, \
> ; opt_arg[1]:logical)
> ;
> ; function add_white_space(str:string,maxlen:integer)
> ; procedure print_opts(opts_name,opts,debug)
> ; procedure print_header(icount:integer,debug)
> ;--------------------------------------------------------------------------------
>
> ;***********************************************************************;
> ; function :
> tointeger ;
> ;
> x:numeric ;
> ; ;
> ; Convert input to
> integer. ;
> ; ;
> ;***********************************************************************;
> undef("tointeger")
> function tointeger(x:numeric)
> local xi
> begin
> if(typeof(x).eq."double")
> xi = doubletointeger(x)
> else
> if(typeof(x).eq."float")
> xi = floattointeger(x)
> else
> if(isatt(x,"_FillValue"))
> xi = new(dimsizes(x),integer,x@_FillValue)
> else
> xi = new(dimsizes(x),integer)
> delete(xi@_FillValue)
> end if
> xi = x
> end if
> end if
> return(xi)
> end
>
> undef("wrf_user_set_xy")
> function wrf_user_set_xy( var:numeric, xp:numeric, yp:numeric,
> x1:numeric, \
> y1:numeric, angle:numeric, opts )
>
> ; mass coordinate version of ncl user routines
>
> local dims,x,y,slope,intercept,distance,dx,dy,dxy,npts,xy
>
> begin
>
> ; find intersection of line and domain boundaries
>
> dims = dimsizes(var)
>
> if (.not. opts) then ; We have a pivot point and location and
> ; need to calculate the start and end points of
> ; the cross section
>
> if ((angle .gt. 315.) .or. (angle .lt. 45.) .or. \
> ((angle .gt. 135.) .and. (angle .lt. 225.)) ) then
>
> ; x = y*slope + intercept
>
> slope = -(360.-angle)/45.
> if( angle .lt. 45. ) then
> slope = angle/45.
> end if
> if( angle .gt. 135.) then
> slope = (angle-180.)/45.
> end if
> intercept = xp - yp*slope
>
> ; find intersections with domain boundaries
>
> y0 = 0.
> x0 = y0*slope + intercept
>
> if( x0 .lt. 0.) then ; intersect outside of left boundary
> x0 = 0.
> y0 = (x0 - intercept)/slope
> end if
> if( x0 .gt. dims(2)-1) then ; intersect outside of right
> boundary
> x0 = dims(2)-1
> y0 = (x0 - intercept)/slope
> end if
>
> y1 = dims(1)-1. ; need to make sure this will be a float?
> x1 = y1*slope + intercept
>
> if( x1 .lt. 0.) then ; intersect outside of left boundary
> x1 = 0.
> y1 = (x1 - intercept)/slope
> end if
> if( x1 .gt. dims(2)-1) then ; intersect outside of right
> boundary
> x1 = dims(2)-1
> y1 = (x1 - intercept)/slope
> end if
>
> else
>
> ; y = x*slope + intercept
>
> slope = (90.-angle)/45.
> if( angle .gt. 225. ) then
> slope = (270.-angle)/45.
> end if
> intercept = yp - xp*slope
>
> ; find intersections with domain boundaries
>
> x0 = 0.
> y0 = x0*slope + intercept
>
> if( y0 .lt. 0.) then ; intersect outside of bottom boundary
> y0 = 0.
> x0 = (y0 - intercept)/slope
> end if
> if( y0 .gt. dims(1)-1) then ; intersect outside of top
> boundary
> y0 = dims(1)-1
> x0 = (y0 - intercept)/slope
> end if
>
> x1 = dims(2)-1. ; need to make sure this will be a float?
> y1 = x1*slope + intercept
>
> if( y1 .lt. 0.) then ; intersect outside of bottom boundary
> y1 = 0.
> x1 = (y1 - intercept)/slope
> end if
> if( y1 .gt. dims(1)-1) then ; intersect outside of top
> boundary
> y1 = dims(1)-1
> x1 = (y1 - intercept)/slope
> end if
>
> end if ; we have beginning and ending points
>
> end if
>
>
> if (opts) then ; We have a specified start and end point
> x0 = xp
> y0 = yp
> if ( x1 .gt. dims(2)-1 ) then
> x1 = dims(2)
> end if
> if ( y1 .gt. dims(1)-1 ) then
> y1 = dims(1)
> end if
> end if
>
> dx = x1 - x0
> dy = y1 - y0
> distance = (dx*dx + dy*dy)^0.5
> npts = tointeger(distance)
> dxy = new(1,typeof(distance))
> dxy = distance/npts
>
> xy = new((/ npts, 2 /),typeof(x1))
>
> dx = dx/npts
> dy = dy/npts
>
> do i=0,npts-1
> xy(i,0) = x0 + i*dx
> xy(i,1) = y0 + i*dy
> end do
>
> ; print(xy)
> return(xy)
>
> end
>
> ;--------------------------------------------------------------------------------
>
> undef("wrf_user_intrp3d")
> function wrf_user_intrp3d( var3d:numeric, z:numeric, \
> plot_type:string, \
> loc_param:numeric, angle:numeric,
> opts:logical )
>
> ; var3d - field to inerpolate (all input fields must be
> unstaggered
> ; z - interpolate to this field (either p/z)
> ; plot_type - interpolate horizontally "h", or vertically "v"
> ; loc_param - level for horizontal plots (eg. 500hPa ; 3000m -
> scalar),
> ; plane for vertical plots (2 values representing an xy
> point
> ; on the model domain through which the vertical plane
> will pass
> ; OR 4 values specifying start and end values
> ; angle - 0.0 for horizontal plots, and
> ; an angle for vertical plots - 90 represent a WE cross
> section
> ; opts Used IF opts is TRUE, else use loc_param and angle to
> determine crosssection
>
> begin
>
> if(plot_type .eq. "h" ) then ; horizontal cross section needed
>
> var2d = wrf_interp_3d_z(var3d,z,loc_param)
>
> if ( z(1,1,1) .gt. 500.) then
> var2d_at_PlotLevelID = loc_param + " hPa"
> else
> var2d_at_PlotLevelID = .001*loc_param + " km"
> end if
>
> end if
>
>
>
>
> if(plot_type .eq. "v" ) then ; vertical cross section needed
>
> ; set vertical cross section
> if (opts) then
> xy = wrf_user_set_xy( z, loc_param(0)-1, loc_param(1)-1,
> \ ; the -1 is for NCL dimensions
> loc_param(2)-1, loc_param(3)-1, \
> angle, opts )
> else
> xy = wrf_user_set_xy( z, loc_param(0), loc_param(1), \
> 0.0, 0.0, angle, opts )
> end if
> xp = dimsizes(xy)
>
>
> ; first we interp the variable, then z
> var2dtmp = wrf_interp_2d_xy( var3d, xy)
> var2dz = wrf_interp_2d_xy( z, xy)
>
>
> ; interp to constant z grid
> z_max = max(var2dz)
> z_min = min(var2dz)
> dz = 0.01*(z_max - z_min)
> nlevels = 101
> z_var2d = new( (/nlevels/), typeof(var2dz))
> z_var2d(0) = z_min
>
> var2d = new( (/nlevels, xp(0)/), typeof(var2dz))
>
> if(var2dz(0,0) .gt. var2dz(1,0) ) then ; monotonically
> decreasing coordinate
> dz = -dz
> z_var2d(0) = z_max
> end if
>
>
> do i=1, nlevels-1
> z_var2d(i) = z_var2d(0)+i*dz
> end do
>
> do i=0,xp(0)-1
> var2d(:,i) = wrf_interp_1d( var2dtmp(:,i), var2dz(:,i),
> z_var2d)
> end do
>
>
> st_x = tointeger(xy(0,0)) + 1
> st_y = tointeger(xy(0,1)) + 1
> ed_x = tointeger(xy(xp(0)-1,0)) + 1
> ed_y = tointeger(xy(xp(0)-1,1)) + 1
> if (opts) then
> var2d_at_Orientation = "Cross-Sesion: (" + \
> st_x + "," + st_y + ") to (" + \
> ed_x + "," + ed_y + ")"
> else
> var2d_at_Orientation = "Cross-Sesion: (" + \
> st_x + "," + st_y + ") to (" + \
> ed_x + "," + ed_y + ") ; center=(" + \
> loc_param(0) + "," + loc_param(1) + \
> ") ; angle=" + angle
> end if
>
> end if
>
>
>
> return(var2d)
>
> end
>
> ;--------------------------------------------------------------------------------
>
> undef("wrf_user_intrp2d")
> function wrf_user_intrp2d( var2d:numeric, \
> loc_param:numeric, angle:numeric,
> opts:logical )
>
> ; var2d - field to inerpolate
> ; loc_param - plane for vertical plots (2 values representing an xy
> point
> ; on the model domain through which the vertical plane
> will pass
> ; OR 4 values specifying start and end values
> ; angle - 0.0 for horizontal plots, and
> ; an angle for vertical plots - 90 represent a WE cross
> section
> ; opts Used IF opts is TRUE, else use loc_param and angle to
> determine crosssection
>
> begin
>
> dims = dimsizes(var2d)
> var2dtmp = new( (/1,dims(0),dims(1)/), typeof(var2d))
> var2dtmp(0,:,:) = var2d(:,:)
>
> ; set vertical cross section
> if (opts) then
> xy = wrf_user_set_xy( var2dtmp, \
> loc_param(0)-1, loc_param(1)-1, \ ;
> the -1 is for NCL dimensions
> loc_param(2)-1, loc_param(3)-1, \
> angle, opts )
> else
> xy = wrf_user_set_xy( var2dtmp, \
> loc_param(0), loc_param(1), \
> 0.0, 0.0, angle, opts )
> end if
> xp = dimsizes(xy)
>
> var1dtmp = wrf_interp_2d_xy( var2dtmp, xy)
>
> dims = dimsizes(var1dtmp)
> var = new( (/dims(1)/), typeof(var1dtmp))
> var(:) = var1dtmp(0,:)
>
> return(var)
>
> end
>
>
> undef("wrf_user_getvar")
> function wrf_user_getvar( nc_file:file, variable:string,
> time:integer )
>
> begin
>
> if( (variable .eq. "uvmet") ) then
> ;; Calculate winds rotated to earth coord.
> ;; Return a 4D array where uvmet(0,:,:,:) is umet and
> ;; uvmet(1,:,:,:) is vmet
>
> pii = 3.14159265
> radians_per_degree = pii/180.
>
> v = nc_file->V(time,:,:,:)
> dimv = dimsizes(v)
> u = nc_file->U(time,:,:,:)
> dimu = dimsizes(u)
>
> uvmet = new( (/ 2, dimu(0), dimu(1), dimv(2)/), typeof(u))
>
> map_projection = nc_file_at_MAP_PROJ
>
>
> if(map_projection .eq. 0) then ; no projection
> uvmet(0,:,:,:) = 0.5*(u(:,:,:dimv(2)-2) + u(:,:,
> 1:dimv(2)-1))
> uvmet(1,:,:,:) = 0.5*(v(:,:dimv(1)-2,:) + v(:,
> 1:dimv(1)-1,:))
> uvmet_at_description = " u,v met velocity"
> uvmet_at_units = "m/s"
> return(uvmet)
> end if
>
>
> cen_lat = nc_file_at_CEN_LAT
> if(isatt(nc_file,"STAND_LON")) then
> cen_long = nc_file_at_STAND_LON
> else
> cen_long = nc_file_at_CEN_LON
> end if
> true_lat1 = nc_file_at_TRUELAT1
> true_lat2 = nc_file_at_TRUELAT2
> if(isfilevar(nc_file,"XLAT"))
> latitude = nc_file->XLAT(0,:,:)
> longitude = nc_file->XLONG(0,:,:)
> else
> latitude = nc_file->XLAT_M(0,:,:)
> longitude = nc_file->XLONG_M(0,:,:)
> end if
>
> cone = 1.
> if( map_projection .eq. 1) then ; Lambert Conformal mapping
> if( (fabs(true_lat1 - true_lat2) .gt. 0.1) .and. \
> (fabs(true_lat2 - 90. ) .gt. 0.1) ) then
> cone = 10^(cos(true_lat1*radians_per_degree)) \
> -10^(cos(true_lat2*radians_per_degree))
> cone = cone/(10^(tan(45. -
> fabs(true_lat1/2.)*radians_per_degree)) - \
> 10^(tan(45. -
> fabs(true_lat2/2.)*radians_per_degree)) )
> else
> cone = sin(fabs(true_lat1)*radians_per_degree)
> end if
> end if
> if(map_projection .eq. 2) then ; polar steraographic
> cone = 1.
> end if
> if(map_projection .eq. 3) then ; Mercator
> cone = 0.
> end if
>
> uvmet = wrf_uvmet( u,v, latitude, longitude, cen_long, cone )
>
> return(uvmet)
>
> end if
>
>
> if( variable .eq. "ua" ) then
> ; U interpolated to mass points
> U = nc_file->U(time,:,:,:)
> dim = dimsizes(U)
> ua = 0.5*(U(:,:,:dim(2)-2) + U(:,:,1:dim(2)-1))
> ua_at_description = "u Velocity"
> ua_at_units = "m/s"
> return(ua)
> end if
>
>
> if( variable .eq. "va" ) then
> ; V interpolated to mass points
> V = nc_file->V(time,:,:,:)
> dim = dimsizes(V)
> va = 0.5*(V(:,:dim(1)-2,:)+V(:,1:dim(1)-1,:))
> va_at_description = "v Velocity"
> va_at_units = "m/s"
> return(va)
> end if
>
>
> if( variable .eq. "wa" ) then
> ; W interpolated to mass points
> W = nc_file->W(time,:,:,:)
> dim = dimsizes(W)
> wa = 0.5*(W(0:dim(0)-2,:,:)+W(1:dim(0)-1,:,:))
> wa_at_description = "Vertical Velocity w"
> wa_at_units = "m/s"
> return(wa)
> end if
>
>
> if( variable .eq. "pressure" ) then
> ; Full model pressure [=base pressure (PB) + pertubation
> pressure (P)]
> P = nc_file->P(time,:,:,:)
> PB = nc_file->PB(time,:,:,:)
> pressure = 0.01 * ( P + PB )
> pressure_at_description = "Pressure"
> pressure_at_units = "hPa"
> return(pressure)
> end if
>
>
> if( variable .eq. "z" ) then
> ; Height [=full geopotentail height / 9.81]
> PH = nc_file->PH(time,:,:,:)
> PHB = nc_file->PHB(time,:,:,:)
> var = ( PH + PHB ) / 9.81
> dim = dimsizes(var)
> z = 0.5*(var(0:dim(0)-2,:,:)+var(1:dim(0)-1,:,:))
> z_at_description = "Height"
> z_at_units = "m"
> return(z)
> end if
>
>
> if( variable .eq. "th" ) then
> ; Potentail Temperature is model output T + 300K
> T = nc_file->T(time,:,:,:)
> th = T + 300.
> th_at_description = "Potential Temperature (theta) "
> th_at_units = "K"
> return(th)
> end if
>
>
> if( variable .eq. "tk" ) then
> ;; function wrf_tk needs theta and pressure (Pa) on input and
> returns temperature in K on return
> T = nc_file->T(time,:,:,:)
> th = T + 300.
> P = nc_file->P(time,:,:,:)
> PB = nc_file->PB(time,:,:,:)
> p = P + PB
> tk = wrf_tk( p , th )
> return(tk)
> end if
>
>
> if( variable .eq. "tc" ) then
> ;; function wrf_tk needs theta and pressure (Pa) on input and
> returns temperature in K on return
> T = nc_file->T(time,:,:,:)
> th = T + 300.
> P = nc_file->P(time,:,:,:)
> PB = nc_file->PB(time,:,:,:)
> p = P + PB
> tc = wrf_tk( p , th )
> tc = tc - 273.16
> tc_at_units = "C" ; Overwrite return units
> return(tc)
> end if
>
>
> if( variable .eq. "td" ) then
> ;; function wrf_td needs qv and pressure (Pa) on input and
> returns dewpoint temperature on return
> QVAPOR = nc_file->QVAPOR(time,:,:,:)
> QVAPOR = QVAPOR > 0.0
> P = nc_file->P(time,:,:,:)
> PB = nc_file->PB(time,:,:,:)
> p = P + PB
> td = wrf_td( p , QVAPOR )
> return(td)
> end if
>
>
> if( variable .eq. "td2" ) then
> ;; function wrf_td needs qv and pressure (Pa) on input and
> returns dewpoint temperature on return
> Q2 = nc_file->Q2(time,:,:)
> Q2 = Q2 > 0.0
> P = nc_file->P(time,0,:,:) ; get only the lowest level,
> to match Q2
> PB = nc_file->PB(time,0,:,:) ; get only the lowest level,
> to match Q2
> p = P + PB
> td2 = wrf_td( p , Q2 )
> td2_at_description = "2m Dewpoint Temperature" ; Overwrite
> return description
> return(td2)
> end if
>
>
> if( variable .eq. "slp" ) then
> ;; first compute theta
> ;; function wrf_tk needs theta and pressure (Pa) on input
> T = nc_file->T(time,:,:,:)
> th = T + 300.
> P = nc_file->P(time,:,:,:)
> PB = nc_file->PB(time,:,:,:)
> p = P + PB
> tk = wrf_tk( p , th )
>
> ;; now compute sea level pressure, from qv, p (Pa), tk, z
> QVAPOR = nc_file->QVAPOR(time,:,:,:)
> QVAPOR = QVAPOR > 0.000
> PH = nc_file->PH(time,:,:,:)
> PHB = nc_file->PHB(time,:,:,:)
> var = ( PH + PHB ) / 9.81
> dim = dimsizes(var)
> z = 0.5 * ( var(0:dim(0)-2,:,:) + var(1:dim(0)-1,:,:) )
>
> ;; get slp
> slp = wrf_slp( z, tk, p, QVAPOR )
> return(slp)
> end if
>
>
> if( variable .eq. "rh" ) then
> ;; first compuet tk
> ;; function wrf_tk needs theta and pressure (Pa) on input
> T = nc_file->T(time,:,:,:)
> th = T + 300.
> P = nc_file->P(time,:,:,:)
> PB = nc_file->PB(time,:,:,:)
> p = P + PB
> tk = wrf_tk( p , th )
>
> ;; now compute rh, using tk, p (Pa), QVAPOR
> QVAPOR = nc_file->QVAPOR(time,:,:,:)
> QVAPOR = QVAPOR > 0.000
> rh = wrf_rh( QVAPOR, p, tk )
> return(rh)
> end if
>
>
>
> ; end of diagnostic variable list - we must want a variable already
> in
> ; the file. check variable dimensionality and pull proper time
> out of file
>
> ndims = dimsizes(filevardimsizes(nc_file,variable))
> if( ndims .eq. 4) then
> var = nc_file->$variable$(time,:,:,:)
> end if
> if( ndims .eq. 3) then
> var = nc_file->$variable$(time,:,:)
> end if
> if( ndims .eq. 2) then
> var = nc_file->$variable$(time,:)
> end if
> if( ndims .eq. 1) then
> var = nc_file->$variable$(time)
> end if
>
> return(var)
>
> end
>
> ;--------------------------------------------------------------------------------
>
> undef("wrf_user_list_times")
> function wrf_user_list_times( nc_file:file )
>
> local times, times_in_file, dims, i
> begin
>
> times_in_file = nc_file->Times
> dims = dimsizes(times_in_file)
> times = new(dims(0),string)
> do i=0,dims(0)-1
> times(i) = chartostring(times_in_file(i,:))
> end do
> times_at_description = "times in file"
> print(times)
> return(times)
>
> end
>
> ;--------------------------------------------------------------------------------
>
> undef("wrf_user_latlon_to_ij")
> function wrf_user_latlon_to_ij( nc_file:file, latitude:numeric, \
> longitude:numeric )
>
> begin
>
> if(isfilevar(nc_file,"XLAT"))
> XLAT = nc_file->XLAT(0,:,:)
> XLONG = nc_file->XLONG(0,:,:)
> else
> XLAT = nc_file->XLAT_M(0,:,:)
> XLONG = nc_file->XLONG_M(0,:,:)
> end if
>
> loc = wrf_latlon_to_ij( XLAT, XLONG, latitude, longitude )
>
> return(loc)
>
> end
>
> ;--------------------------------------------------------------------------------
> ;--------------------------------------------------------------------------------
> ;--------------------------------------------------------------------------------
>
> undef("delete_attrs")
> procedure delete_attrs(opts:logical)
>
> ; This procedure does some cleanup by removing unneeded attributes
> ; so they don't get passed to other routines by accident.
>
> begin
> list_attrs = (/"MainTitle","MainTitlePos","MainTitlePosF", \
> "InitTime","ValidTime","TimePos","TimePosF", \
> "NoHeaderFooter","TimeLabel","LevelLabel", \
> "FieldTitle","UnitLabel","NumVectors","AspectRatio", \
> "SubFieldTitle","PlotOrientation","PlotLevelID", \
> "ContourParameters","FontHeightF","Footer"/)
>
> do i=0,dimsizes(list_attrs)-1
> if(isatt(opts,list_attrs(i))) then
> delete(opts@$list_attrs(i)$)
> end if
> end do
> end
>
> ;--------------------------------------------------------------------------------
> ;--------------------------------------------------------------------------------
> undef("set_cn_resources")
> function set_cn_resources (data[*][*]:numeric, res:logical)
>
> begin
>
> opts = res
>
> ; The ContourParameters resource can either be a scalar that
> ; represents the contour level spacing, or it can be an array
> ; of three elements that represent the minimum level, the maximum
> ; level, and the level spacing.
> ;
> mx = max(data)
> mn = min(data)
>
> if(mn.ne.mx.and.opts.and.isatt(opts,"ContourParameters")) then
> if(dimsizes(opts_at_ContourParameters) .eq. 1) then
>
> ; Only the contour interval is specified.
> nlev = tointeger((mx-mn)/opts_at_ContourParameters)+1
> levels = nice_mnmxintvl(mn,mx,nlev,True)
> if(levels(0) .lt. 0.) then
> ; Set a zero contour.
> nlev = tointeger(levels(0)/opts_at_ContourParameters) - 1
> levels(0) = nlev*opts_at_ContourParameters
> end if
> nlev = tointeger((levels(1)-levels(0))/opts_at_ContourParameters)+1
> levels(1) = levels(0) + nlev*opts_at_ContourParameters
> levels(2) = opts_at_ContourParameters
>
> ; Min level, max level, and level spacing are specified by user.
> else
> if(dimsizes(opts_at_ContourParameters) .eq. 3) then
> levels = opts_at_ContourParameters
> else
> print("wrf_contour: Warning: illegal setting for
> ContourParameters attribute")
> end if
> end if
>
> end if
>
> ; Contour levels
> if(isvar("levels")) then
> opts_at_cnLevelSelectionMode = get_res_value_keep(opts,
> "cnLevelSelectionMode", "ManualLevels")
> opts_at_cnMinLevelValF = get_res_value_keep(opts,
> "cnMinLevelValF", levels(0))
> opts_at_cnMaxLevelValF = get_res_value_keep(opts,
> "cnMaxLevelValF", levels(1))
> opts_at_cnLevelSpacingF = get_res_value_keep(opts,
> "cnLevelSpacingF",levels(2))
> delete(levels)
> end if
>
>
> ; Set the default zero line thickness to 2, and the negative contour
> ; line dash pattern to 1 (0 is solid).
> opts_at_gsnContourZeroLineThicknessF = get_res_value_keep(opts,
> "gsnContourZeroLineThicknessF",2.0)
> opts_at_gsnContourNegLineDashPattern = get_res_value_keep(opts,
> "gsnContourNegLineDashPattern",1)
>
>
> ; Set resources that apply for both filled and line contour plots.
> opts_at_cnFillDrawOrder =
> get_res_value_keep(opts,"cnFillDrawOrder", "PreDraw")
>
>
> opts_at_cnLineLabelAngleF =
> get_res_value_keep(opts,"cnLineLabelAngleF", 0.0)
> opts_at_cnLineLabelFontHeightF =
> get_res_value_keep(opts,"cnLineLabelFontHeightF", 0.015)
> opts_at_cnInfoLabelFontHeightF =
> get_res_value_keep(opts,"cnInfoLabelFontHeightF", 0.015)
> opts_at_cnLineLabelPerimOn =
> get_res_value_keep(opts,"cnLineLabelPerimOn", True)
> opts_at_cnInfoLabelPerimOn =
> get_res_value_keep(opts,"cnInfoLabelPerimOn", False)
> opts_at_cnLineLabelBackgroundColor =
> get_res_value_keep(opts,"cnLineLabelBackgroundColor", -1)
> opts_at_cnHighLabelBackgroundColor =
> get_res_value_keep(opts,"cnHighLabelBackgroundColor", -1)
> opts_at_cnLowLabelBackgroundColor =
> get_res_value_keep(opts,"cnLowLabelBackgroundColor", -1)
> opts_at_cnLineColor =
> get_res_value_keep(opts,"cnLineColor", "Black")
> opts_at_cnLineLabelFontColor = opts_at_cnLineColor
> opts_at_cnLineLabelPerimColor = opts_at_cnLineColor
> opts_at_cnInfoLabelFontColor = opts_at_cnLineColor
> opts_at_cnHighLabelFontColor = opts_at_cnLineColor
> opts_at_cnLowLabelFontColor = opts_at_cnLineColor
>
>
> ; Set field Title and levels if available
> if(isatt(opts,"FieldTitle")) then
> opts_at_cnInfoLabelString = opts_at_FieldTitle + " Contours: $CMN$
> to $CMX$ by $CIU$"
> else
> if(isatt(data,"description")) then
> opts_at_cnInfoLabelString = data_at_description + " Contours: $CMN
> $ to $CMX$ by $CIU$"
> else
> opts_at_cnInfoLabelString = " Contours: $CMN$ to $CMX$ by $CIU$"
> end if
> end if
>
>
> return(opts)
> end
> ;--------------------------------------------------------------------------------
> ;--------------------------------------------------------------------------------
> undef("set_lb_resources")
> function set_lb_resources (data[*][*]:numeric, res:logical)
>
> begin
>
> opts = res
>
>
> ; Somewhat convoluted way to see if a labelbar is not desired.
> if(check_attr(opts,"pmTickMarkDisplayMode","Never",True).or.\
> check_attr(opts,"pmTickMarkDisplayMode",-1,False).or.\
> check_attr(opts,"pmTickMarkDisplayMode",0,False).or. \
> check_attr(opts,"lbLabelBarOn",False,False).or.\
> check_attr(opts,"lbLabelBarOn",0,False)) then
> lbar_on = False
> else
> lbar_on = True
> end if
> atmp = get_res_value(opts,"lbLabelBarOn",True) ; Remove this resource
> delete(atmp) ; just in case.
>
>
> ; Possible title for the labelbar
> if(isatt(opts,"FieldTitle")) then
> lb_desc = opts_at_FieldTitle
> else
> if(isatt(data,"description")) then
> lb_desc = data_at_description
> else
> lb_desc = ""
> end if
> end if
>
> if(isatt(opts,"UnitLabel") ) then
> lb_desc = lb_desc + " (" + opts_at_UnitLabel + ")"
> else
> if(isatt(data,"units") .and. .not.(data_at_units.eq."")) then
> lb_desc = lb_desc + " (" + data_at_units + ")"
> end if
> end if
>
>
> if(.not.isatt(opts,"cnFillColors")) then
> opts_at_gsnSpreadColors = get_res_value_keep(opts,
> "gsnSpreadColors", True)
> end if
> opts_at_cnInfoLabelOn = get_res_value_keep(opts,"cnInfoLabelOn",
> False)
> opts_at_cnLinesOn = get_res_value_keep(opts,"cnLinesOn", False)
> opts_at_cnLineLabelsOn = get_res_value_keep(opts,"cnLineLabelsOn",
> False)
>
> ; Labelbar resources
> if(lbar_on) then
> opts_at_pmLabelBarDisplayMode =
> get_res_value_keep(opts,"pmLabelBarDisplayMode", "Always")
> opts_at_pmLabelBarSide =
> get_res_value_keep(opts,"pmLabelBarSide", "Bottom")
> opts_at_lbAutoManage =
> get_res_value_keep(opts,"lbAutoManage",False)
> opts_at_lbOrientation =
> get_res_value_keep(opts,"lbOrientation", "Horizontal")
> opts_at_lbPerimOn =
> get_res_value_keep(opts,"lbPerimOn", False)
> opts_at_lbLabelJust =
> get_res_value_keep(opts,"lbLabelJust", "BottomCenter")
> opts_at_lbLabelAutoStride =
> get_res_value_keep(opts,"lbLabelAutoStride",True)
> opts_at_lbBoxMinorExtentF =
> get_res_value_keep(opts,"lbBoxMinorExtentF", 0.13)
> opts_at_lbTitleFontHeightF =
> get_res_value_keep(opts,"lbTitleFontHeightF", 0.015)
> opts_at_lbLabelFontHeightF =
> get_res_value_keep(opts,"lbLabelFontHeightF", 0.015)
> opts_at_pmLabelBarOrthogonalPosF =
> get_res_value_keep(opts,"pmLabelBarOrthogonalPosF", -0.1)
>
> opts_at_lbTitleOn =
> get_res_value_keep(opts,"lbTitleOn", True)
> if(lb_desc.ne."" .and. opts_at_lbTitleOn) then
> opts_at_lbTitleOn =
> get_res_value_keep(opts,"lbTitleOn", True)
> opts_at_lbTitleString =
> get_res_value_keep(opts,"lbTitleString", lb_desc)
> opts_at_lbTitleJust =
> get_res_value_keep(opts,"lbTitleJust", "BottomCenter")
> opts_at_lbTitleOffsetF =
> get_res_value_keep(opts,"lbTitleOffsetF", -0.5)
> else
> opts_at_lbTitleOn = False
> end if
> end if
>
>
> return(opts)
> end
> ;--------------------------------------------------------------------------------
> ;--------------------------------------------------------------------------------
> undef("set_title_resources")
> function set_title_resources (data[*][*]:numeric, res:logical)
>
> begin
>
> opts = res
>
> ; Set field Title and levels if available
> if(isatt(opts,"FieldTitle")) then
> SubTitles = opts_at_FieldTitle
> else
> if(isatt(data,"description")) then
> SubTitles = data_at_description
> else
> SubTitles = "UnKnown"
> end if
> end if
>
> if(isatt(opts,"SubFieldTitle")) then
> SubTitles = SubTitles + " " + opts_at_SubFieldTitle
> end if
> if(isatt(opts,"UnitLabel")) then
> SubTitles = SubTitles + " (" + opts_at_UnitLabel + ")"
> else
> if(isatt(data,"units") .and. .not.(data_at_units.eq."")) then
> SubTitles = SubTitles + " (" + data_at_units + ")"
> end if
> end if
> if (isatt(opts,"PlotLevelID")) then
> SubTitles = SubTitles + " at " + opts_at_PlotLevelID
> else
> if (isatt(data,"PlotLevelID")) then
> SubTitles = SubTitles + " at " + data_at_PlotLevelID
> end if
> end if
> opts_at_tiMainString = SubTitles
>
> return(opts)
> end
> ;--------------------------------------------------------------------------------
> ;--------------------------------------------------------------------------------
> undef("set_vc_resources")
> function set_vc_resources (res:logical)
>
> begin
>
> opts = res
>
> if ( isatt(opts,"vpWidthF") ) then
> ; num_vectors is used for vcMinDistanceF and vcRefLengthF
> width = opts_at_vpWidthF
> num_vectors = get_res_value(opts,"NumVectors",25.0)
> opts_at_vcMinDistanceF = get_res_value_keep(opts,"vcMinDistanceF",
> width/num_vectors)
> opts_at_vcRefLengthF = get_res_value_keep(opts,"vcRefLengthF",
> width/num_vectors)
> else
> opts_at_vcMinDistanceF = get_res_value_keep(opts,"vcMinDistanceF",
> 0.02)
> opts_at_vcRefLengthF = get_res_value_keep(opts,"vcRefLengthF",
> 0.02)
> end if
>
>
> opts_at_vcGlyphStyle = get_res_value_keep(opts,"vcGlyphStyle",
> "WindBarb")
> opts_at_vcWindBarbColor = get_res_value_keep(opts,"vcWindBarbColor",
> "Black")
> opts_at_vcRefAnnoOn = get_res_value_keep(opts,"vcRefAnnoOn", False)
> opts_at_vcMinFracLengthF =
> get_res_value_keep(opts,"vcMinFracLengthF", .2)
>
> return(opts)
> end
> ;--------------------------------------------------------------------------------
> ;--------------------------------------------------------------------------------
>
> ;procedure
> _SetMainTitle(nc_file:file,wks[1]:graphic,cn[1]:graphic,opts)
>
> ; This procedure checks the input data for certain attributes, and
> ; based on those, sets MainTitle, InitTime and ValidTime
> ;
> ; Attributes recognized by this procedure:
> ; MainTitle (main title - top left)
> ; (with Init time top right)
> ; TimeLabel (valid time - right under init time)
> ; NoHeaderFooter (switch all headers and footers off
> - mainly for panels)
> ;
> ; If the "NoHeaderFooter" attribute exists and is set True, then
> ; don't create any titles.
>
> begin
> ;
> opts = True
> if(opts.and.isatt(opts,"NoHeaderFooter").and.opts_at_NoHeaderFooter)
> then
> return
> end if
>
> ;
> ; Set basic plot font
> ;
> font_height = get_res_value_keep(opts,"FontHeightF",0.01)
> ;
> ;
> ; If a MainTitle attribute hasn't been set, then set to "WRF"
> ; Also set an Initial time
> ;
> ; MAIN Header of plot
> opts = True
> opts_at_MainTitle = get_res_value_keep(opts,"MainTitle", " ")
> opts_at_MainTitlePos = get_res_value_keep(opts,"MainTitlePos", "Left")
> opts_at_InitTime = get_res_value_keep(opts,"InitTime", True)
> opts_at_ValidTime = get_res_value_keep(opts,"ValidTime", True)
> opts_at_TimePos = get_res_value_keep(opts,"TimePos", "Right")
> opts_at_Footer = get_res_value_keep(opts,"Footer", True)
>
> if (opts_at_MainTitlePos .eq. "Left")
> opts_at_MainTitlePos = "CenterLeft"
> opts_at_MainTitlePosF = 0.0
> end if
> if (opts_at_MainTitlePos .eq. "Center")
> opts_at_MainTitlePos = "CenterCenter"
> opts_at_MainTitlePosF = 0.5
> end if
> if (opts_at_MainTitlePos .eq. "Right")
> opts_at_MainTitlePos = "CenterRight"
> opts_at_MainTitlePosF = 1.0
> end if
>
> if (opts_at_TimePos .eq. "Left")
> MTOPosF = 0.30
> else
> MTOPosF = 0.20
> end if
>
> txt0 = create "MainPlotTilte" textItemClass wks
> "txString" : opts_at_MainTitle
> "txFontHeightF" : font_height*1.5
> end create
> anno = NhlAddAnnotation(cn,txt0)
> setvalues anno
> "amZone" : 3
> "amSide" : "Top"
> "amJust" : opts_at_MainTitlePos
> "amParallelPosF" : opts_at_MainTitlePosF
> "amOrthogonalPosF" : MTOPosF
> "amResizeNotify" : False
> end setvalues
>
> ; Time information on plot
> if (opts_at_TimePos .eq. "Left")
> opts_at_TimePos = "CenterLeft"
> opts_at_TimePosF = 0.0
> if (opts_at_MainTitlePos .eq. "CenterLeft")
> MTOPosF = MTOPosF - 0.05
> end if
> end if
> if (opts_at_TimePos .eq. "Right")
> opts_at_TimePos = "CenterRight"
> opts_at_TimePosF = 1.0
> if (opts_at_MainTitlePos .eq. "CenterRight")
> MTOPosF = MTOPosF - 0.05
> end if
> end if
>
> if( isatt(nc_file,"START_DATE") .and. opts_at_InitTime ) then
> InitTime = "Init: " + nc_file_at_START_DATE
> txt1 = create "InitTime" textItemClass ks
> "txString" : InitTime
> "txFontHeightF" : font_height
> end create
> anno = NhlAddAnnotation(cn,txt1)
> setvalues anno
> "amZone" : 3
> "amSide" : "Top"
> "amJust" : opts_at_TimePos
> "amParallelPosF" : opts_at_TimePosF
> "amOrthogonalPosF" : MTOPosF
> "amResizeNotify" : False
> end setvalues
> end if
>
> plot_narrow = False
> if((opts).and.(isatt(opts,"vpWidthF")).and.
> (isatt(opts,"vpHeightF"))) then
> ph = opts_at_vpHeightF
> pw = opts_at_vpWidthF
> phw = ph/pw
> if ( phw .gt. 1.8 ) then
> plot_narrow = True
> end if
> end if
>
> if( opts_at_ValidTime .and. isatt(opts,"TimeLabel") ) then
>
> ValidTime = "Valid: " + opts_at_TimeLabel
>
> MTOPosF = MTOPosF - 0.03
> txt2 = create "ValidTime" textItemClass wks
> "txString" : ValidTime
> "txFontHeightF" : font_height
> end create
> anno = NhlAddAnnotation(cn,txt2)
> setvalues anno
> "amZone" : 3
> "amSide" : "Top"
> "amJust" : opts_at_TimePos
> "amParallelPosF" : opts_at_TimePosF
> "amOrthogonalPosF" : MTOPosF
> "amResizeNotify" : False
> end setvalues
> end if
>
>
> ; Add Footer if called for
> if( opts_at_Footer ) then
> footer1 = nc_file_at_TITLE
> dis = nc_file_at_DX / 1000.0
> if(isatt(nc_file,"WEST-EAST_PATCH_END_UNSTAG"))
> WE = "WEST-EAST_GRID_DIMENSION"
> SN = "SOUTH-NORTH_GRID_DIMENSION"
> BT = "BOTTOM-TOP_GRID_DIMENSION"
> footer2 = " Phys Opt = " + nc_file_at_MP_PHYSICS + \
> " ; PBL Opt = " + nc_file_at_BL_PBL_PHYSICS + \
> " ; Cu Opt = " + nc_file_at_CU_PHYSICS + \
> " ; WE = " + nc_file@$WE$ + \
> " ; SN = " + nc_file@$SN$ + \
> " ; Levels = " + nc_file@$BT$ + \
> " ; Dis = " + dis + "km"
> else
> WE = "WEST-EAST_GRID_DIMENSION"
> SN = "SOUTH-NORTH_GRID_DIMENSION"
> BT = "BOTTOM-TOP_GRID_DIMENSION"
> footer2 = " WE = " + nc_file@$WE$ + \
> " ; SN = " + nc_file@$SN$ + \
> " ; Levels = " + nc_file@$BT$ + \
> " ; Dis = " + dis + "km"
> end if
> Footer = footer1 + "~C~" + footer2
> else
> Footer = " "
> end if
> txt3 = create "Footer" textItemClass wks
> "txString" : Footer
> "txFontHeightF" : font_height*.9
> end create
> anno = NhlAddAnnotation(cn,txt3)
> setvalues anno
> "amZone" : 1
> ; "amZone" : 7
> "amJust" : "TopLeft"
> "amSide" : "Bottom"
> "amParallelPosF" : 0.0
> "amOrthogonalPosF" : -0.55
> "amResizeNotify" : False
> end setvalues
>
>
> ; Add X-setion information if needed
> if(opts.and.isatt(opts,"PlotOrientation")) then
> ;Xsection = "Cross-Section Orientation : " + opts_at_PlotOrientation
> Xsection = opts_at_PlotOrientation
> txt4 = create "Xsection" textItemClass wks
> "txString" : Xsection
> "txFontHeightF" : font_height*.9
> end create
> anno = NhlAddAnnotation(cn,txt4)
> setvalues anno
> "amZone" : 3
> "amSide" : "Top"
> "amJust" : "CenterRight"
> "amParallelPosF" : 1.0
> "amOrthogonalPosF" : 0.005
> "amResizeNotify" : False
> end setvalues
> end if
>
> end
>
> ;--------------------------------------------------------------------------------
> ;--------------------------------------------------------------------------------
> function wrf_contour_ps(nc_file:file,wks[1]: graphic, data[*]
> [*]:numeric, \
> opt_args[1]:logical)
>
> begin
>
> callFrame = True
> if ( isatt(opt_args,"FrameIT") ) then
> if ( .not.opt_args_at_FrameIT ) then
> callFrame = False
> end if
> delete (opt_args_at_FrameIT)
> end if
>
> lat2U = nc_file->XLAT_U(0,:,:)
> lon2U = nc_file->XLONG_U(0,:,:)
>
> opts = opt_args
> opts_at_sfXArray = lon2U
> opts_at_sfYArray = lat2U
> opts_at_sfDataArray = data
> opts_at_mpProjection = "Stereographic"
> opts_at_mpEllipticalBoundary = True
> opts_at_mpFillOn = False
>
> opts_at_mpGeophysicalLineColor = get_res_value_keep(opts,
> "mpGeophysicalLineColor","Gray")
> opts_at_mpGeophysicalLineThicknessF = get_res_value_keep(opts,
> "mpGeophysicalLineThicknessF",2.0)
>
>
> ; Set the contour resources
> opts = set_cn_resources(data,opts)
> opts_at_cnInfoLabelFontHeightF = 0.012
> opts_at_cnLineLabelPerimOn = False
> opts_at_cnInfoLabelPerimOn = True
>
>
> ; Find out if we are working with a contour or a shaded plot
> ; fill_on = False : line contour plot
> ; fill_on = True : filled contour plot
> fill_on = get_res_value_keep(opts,"cnFillOn",False)
> if(fill_on) then ; set lb resources if needed
> opts_at_pmLabelBarDisplayMode =
> get_res_value_keep(opts,"pmLabelBarDisplayMode", "Never")
> opts = set_lb_resources(data,opts)
> opts_at_pmLabelBarOrthogonalPosF = 0.0
> opts_at_lbTitleJust = "BottomLeft"
> end if
>
>
>
> opts_at_gsnDraw = False
> opts_at_gsnFrame = False
> opts_at_gsnMaximize = False
> delete_attrs(opts)
> cn = gsn_csm_contour_map_polar(wks,data,opts) ; Create the
> plot.
>
> draw(cn)
>
>
> if ( callFrame ) then
> frame(wks)
> end if
>
> return (cn)
> end
> ;--------------------------------------------------------------------------------
> function wrf_vector_ps(nc_file:file,wks[1]: graphic, \
> data_u[*][*]:numeric, data_v[*][*]:numeric, \
> opt_args[1]:logical)
>
> begin
>
> callFrame = True
> if ( isatt(opt_args,"FrameIT") ) then
> if ( .not.opt_args_at_FrameIT ) then
> callFrame = False
> end if
> delete (opt_args_at_FrameIT)
> end if
>
> if(isfilevar(nc_file,"XLAT"))
> lat2T = nc_file->XLAT(0,:,:)
> lon2T = nc_file->XLONG(0,:,:)
> else
> lat2T = nc_file->XLAT_M(0,:,:)
> lon2T = nc_file->XLONG_M(0,:,:)
> end if
>
> opts = opt_args
> opts_at_vfXArray = lon2T
> opts_at_vfYArray = lat2T
> opts_at_vfUDataArray = data_u
> opts_at_vfVDataArray = data_v
> opts_at_mpProjection = "Stereographic"
> opts_at_mpEllipticalBoundary = True
> opts_at_mpFillOn = False
>
> opts_at_mpGeophysicalLineColor = get_res_value_keep(opts,
> "mpGeophysicalLineColor","Gray")
> opts_at_mpGeophysicalLineThicknessF = get_res_value_keep(opts,
> "mpGeophysicalLineThicknessF",2.0)
>
>
> ; Set vector resources
> opts = set_vc_resources(opts)
>
>
> opts_at_gsnDraw = False
> opts_at_gsnFrame = False
> opts_at_gsnMaximize = False
> delete_attrs(opts)
> cn = gsn_csm_vector_map_polar(wks,data_u,data_v,opts) ;
> Create the plot.
>
> draw(cn)
>
>
> if ( callFrame ) then
> frame(wks)
> end if
>
> return (cn)
> end
> ;--------------------------------------------------------------------------------
>
> ;undef("wrf_contour")
> function wrf_contour(nc_file:file,wks[1]: graphic, data[*]
> [*]:numeric, \
> opt_args[1]:logical)
>
> ; This function creates a contour plot and adds some titles to it.
> ;
> ; 1. Determine width to height ratio of plot.
> ;
> ; 2. First determine if this is to be a filled or line
> ; contour plot (fill_on)
> ;
> ; 3. If the ContourParameters attribute is set, then calculate
> ; the contour levels.
> ;
> ; 4. Set two resources for setting the zero contour line to
> ; a larger thickness, and for changing the negative contour
> ; lines to a dashed pattern.
> ;
> ; 5. If doing a filled contour plot, set a title for the labelbar
> ; based on whether a units attribute is set.
> ;
> ; 6. Make a copy of the resource list, and set some additional
> ; resources for filled contour plots.
> ;
> ; 7. Create the contour plot, attach the titles, and draw
> ; and advance the frame (if requested).
>
> local dims
> begin
> opts = opt_args ; Make a copy of the resource list.
>
> if(opts.and.isatt(opts,"mpOutlineBoundarySets")) then
> delete(opts_at_mpOutlineBoundarySets)
> end if
>
>
> ; Calculate ratio of plot width and height. Note that this doesn't
> ; affect the setting of gsnMaximize to True, because gsnMaximize will
> ; retain the aspect ratio of the plot.
>
> if(opts.and.isatt(opts,"AspectRatio")) then
> ratio = opts_at_AspectRatio
> else
> dims = dimsizes(data)
> ratio = 1.*dims(0)/dims(1)
> if(ratio .gt. 1.2) then
> ratio = 1.2
> end if
> if(ratio .lt. 0.6667) then
> ratio = 0.6667
> end if
> end if
>
> if(ratio .gt. 1)
> width = 0.65 * 1.0/ratio
> height = 0.65
> else
> width = 0.85
> height = 0.85 * ratio
> end if
>
> opts_at_vpWidthF = get_res_value_keep(opts,"vpWidthF", width)
> opts_at_vpHeightF = get_res_value_keep(opts,"vpHeightF", height)
>
>
> ; Set some basic contour resources
> opts = set_cn_resources(data,opts)
>
>
> ; Find out if we are working with a contour or a shaded plot
> ; fill_on = False : line contour plot
> ; fill_on = True : filled contour plot
> fill_on = get_res_value_keep(opts,"cnFillOn",False)
> if(fill_on) then ; set lb resources if needed
> opts = set_lb_resources(data,opts)
> atmp = get_res_value(opts,"lbLabelBarOn",True) ; Remove this
> resource
> delete(atmp) ; just in case.
> end if
>
>
> ; Set Title resources
> opts = set_title_resources(data,opts)
>
>
> ; Setting gsnScale to True ensures that the tickmark lengths and
> labels
> ; will be the same size on both axes.
> opts_at_gsnScale = get_res_value_keep(opts,"gsnScale", True)
>
>
> ; The default is not to draw the plot or advance the frame, and
> ; to maximize the plot in the frame.
> opts_at_gsnDraw = False ; Make sure don't draw or frame or,
> opts_at_gsnFrame = False ; maximize, b/c we'll do this later.
> opts_at_gsnMaximize = False
>
>
> opts2 = opts
> delete_attrs(opts2) ; Clean up.
> cn = gsn_contour(wks,data,opts2) ; Create the plot.
> SetMainTitle(nc_file,wks,cn,opts) ; Set some titles
>
>
> opts2_at_gsnDraw = get_res_value_keep(opts2,"gsnDraw", False)
> opts2_at_gsnFrame = get_res_value_keep(opts2,"gsnFrame", False)
> opts2_at_gsnMaximize = get_res_value_keep(opts2,"gsnMaximize", True)
>
> draw_and_frame
> (wks,cn,opts2_at_gsnDraw,opts2_at_gsnFrame,False,opts2_at_gsnMaximize)
>
> return(cn) ; Return
>
> end
>
> ;--------------------------------------------------------------------------------
>
> undef("wrf_vector")
> function wrf_vector(nc_file:file,wks[1]: graphic, data_u[*]
> [*]:numeric, \
> data_v[*][*]:numeric, opt_args[1]:logical)
> ;
> ; This function creates a vector plot and adds some titles to it.
> ;
> ; 1. Determine width to height ratio of plot. Will also be use
> ; to calculate values for vector resources later.
> ;
> ; 2. Make a copy of the resource list, and set some additional
> ; resources.
> ;
> ; 3. Create the vector plot, attach the titles, and draw
> ; and advance the frame (if requested).
>
> local dims
> begin
> opts = opt_args ; Make a copy of the resource list.
>
> if(opts.and.isatt(opts,"mpOutlineBoundarySets")) then
> delete(opts_at_mpOutlineBoundarySets)
> end if
> ;
> ; The ratio is used to determine the width and height of the
> ; plot, and also to determine the value for the vcMinDistanceF
> ; resource.
> ;
> if(opts.and.isatt(opts,"AspectRatio")) then
> ratio = get_res_value(opts,"AspectRatio",0.)
> else
> dims = dimsizes(data_u)
> ratio = 1.*dims(0)/dims(1)
> if(ratio .gt. 1.2) then
> ratio = 1.2
> end if
> if(ratio .lt. 0.6667) then
> ratio = 0.6667
> end if
> end if
>
> if(ratio .gt. 1)
> width = 0.65/ratio
> height = 0.65
> else
> width = 0.95
> height = 0.95 * ratio
> end if
>
> opts_at_vpWidthF = get_res_value_keep(opts,"vpWidthF", width)
> opts_at_vpHeightF = get_res_value_keep(opts,"vpHeightF", height)
>
>
> ; Set Title resources
> opts = set_title_resources(data_u,opts)
>
>
> ; Set vector resources
> opts = set_vc_resources(opts)
>
>
> ; Setting gsnScale to True ensures that the tickmark lengths and
> labels
> ; will be the same size on both axes.
> opts_at_gsnScale = get_res_value_keep(opts,"gsnScale", True)
>
>
> ; The default is not to draw the plot or advance the frame, and
> ; to maximize the plot in the frame.
> opts_at_gsnDraw = False ; Make sure don't draw or frame or,
> opts_at_gsnFrame = False ; maximize, b/c we'll do this later.
> opts_at_gsnMaximize = False
>
> opts2 = opts
> delete_attrs(opts2) ; Clean up.
> vct = gsn_vector(wks,data_u,data_v,opts2) ; Create vector plot.
> _SetMainTitle(nc_file,wks,vct,opts)
>
> opts2_at_gsnDraw = get_res_value_keep(opts,"gsnDraw", False)
> opts2_at_gsnFrame = get_res_value_keep(opts,"gsnFrame", False)
> opts2_at_gsnMaximize = get_res_value_keep(opts,"gsnMaximize", True)
> draw_and_frame(wks,vct,opts2_at_gsnDraw,opts2_at_gsnFrame,False, \
> opts2_at_gsnMaximize)
>
> return(vct) ; Return.
> end
>
>
> ;--------------------------------------------------------------------------------
>
> undef("wrf_map")
> function wrf_map(wks[1]:graphic,in_file[1]:file,opt_args[1]:logical)
>
> begin
> ;
> ; This function creates a map plot, and bases the projection on
> ; the MAP_PROJ attribute in the given file.
> ;
> ; 1. Make a copy of the resource list, and set some resources
> ; common to all map projections.
> ;
> ; 2. Determine the projection being used, and set resources based
> ; on that projection.
> ;
> ; 3. Create the map plot, and draw and advance the frame
> ; (if requested).
> ;
> if(isatt(in_file,"MAP_PROJ"))
> ; if(in_file_at_MAP_PROJ .eq. 0)
> ; print("wrf_map: Error: attribute MAP_PROJ=0.")
> ; print(" Data does not have a valid map projection")
> ; return(new(1,graphic))
> ; end if
> ;
> ; There's a bug with map tickmarks that you can't set the length of
> ; the tickmark to 0.0 in a setvalues call. You have to do it
> ; during the create call. Double-check, though, b/c the bug may have
> ; been fixed by now.
> ;
> opts = opt_args ; Make a copy of the resource list
> opts = True
> ;
> ; Set some resources depending on what kind of map projection is
> ; chosen.
> ;
> ; MAP_PROJ = 1 : "LambertConformal"
> ; MAP_PROJ = 2 : "Stereographic"
> ; MAP_PROJ = 3 : "Mercator"
> ; MAP_PROJ = 6 : "Lat/Lon"
> ;
> ; Set some resources common to all three map projections.
> ;
> if(any(in_file_at_MAP_PROJ.eq.(/1,2,3,6/))) then
> if(isfilevar(in_file,"XLAT"))
> lat = in_file->XLAT(0,:,:)
> lon = in_file->XLONG(0,:,:)
> else
> lat = in_file->XLAT_M(0,:,:)
> lon = in_file->XLONG_M(0,:,:)
> end if
> dims = dimsizes(lat)
> ;
> ; "LowRes" is the default that NCL uses, so you don't need to
> ; set it here. However, if you want a higher resolution, use
> ; "MediumRes". If you want higher resolution for the coastlines,
> ; then set it to "HighRes", but then you also need to download
> ; the RANGS-GSHHS database. Higher resolutions take longer to
> ; draw.
> ;
> opts_at_mpDataBaseVersion = get_res_value_keep(opts,
> "mpDataBaseVersion","MediumRes")
> ;opts_at_mpOutlineBoundarySets = get_res_value_keep(opts,
> "mpOutlineBoundarySets", "AllBoundaries")
> opts_at_mpOutlineBoundarySets = get_res_value_keep(opts,
> "mpOutlineBoundarySets", "GeophysicalAndUSStates")
> opts_at_mpPerimLineThicknessF = get_res_value_keep(opts,
> "mpPerimLineThicknessF", 1.0)
> opts_at_tmXBLabelFontHeightF = get_res_value_keep(opts,
> "tmXBLabelFontHeightF", 0.01)
> opts_at_tmYLLabelFontHeightF = get_res_value_keep(opts,
> "tmYLLabelFontHeightF", 0.01)
> ;
> ; Select portion of the map to view.
> ;
> ;if(any(in_file_at_MAP_PROJ.eq.(/1,2,3/))) then
> opts_at_mpLimitMode = get_res_value_keep(opts,
> "mpLimitMode","Corners")
> opts_at_mpLeftCornerLatF = get_res_value_keep(opts,
> "mpLeftCornerLatF",lat(0,0))
> opts_at_mpLeftCornerLonF = get_res_value_keep(opts,
> "mpLeftCornerLonF",lon(0,0))
> opts_at_mpRightCornerLatF = get_res_value_keep(opts,
> "mpRightCornerLatF", lat(dims(0)-1,dims(1)-1))
> opts_at_mpRightCornerLonF = get_res_value_keep(opts,
> "mpRightCornerLonF", lon(dims(0)-1,dims(1)-1))
> if ( isatt(opts,"ZoomIn") .and. opts_at_ZoomIn ) then
>
> y1 = 0
> x1 = 0
> y2 = dims(0)-1
> x2 = dims(1)-1
> if ( isatt(opts,"Ystart") ) then
> y1 = opts_at_Ystart
> delete(opts_at_Ystart)
> end if
> if ( isatt(opts,"Xstart") ) then
> x1 = opts_at_Xstart
> delete(opts_at_Xstart)
> end if
> if ( isatt(opts,"Yend") ) then
> if ( opts_at_Yend .le. y2 ) then
> y2 = opts_at_Yend
> end if
> delete(opts_at_Yend)
> end if
> if ( isatt(opts,"Xend") ) then
> if ( opts_at_Xend .le. x2 ) then
> x2 = opts_at_Xend
> end if
> delete(opts_at_Xend)
> end if
>
> opts_at_mpLeftCornerLatF = lat(y1,x1)
> opts_at_mpLeftCornerLonF = lon(y1,x1)
> opts_at_mpRightCornerLatF = lat(y2,x2)
> opts_at_mpRightCornerLonF = lon(y2,x2)
>
> delete(opts_at_ZoomIn)
> end if
>
> if ( opts_at_mpRightCornerLonF .lt. 0.0 ) then
> opts_at_mpRightCornerLonF = opts_at_mpRightCornerLonF + 360.0
> ;print(opts_at_mpLeftCornerLonF)
> ;print(opts_at_mpRightCornerLonF)
> end if
> ;end if
> ;
> ; Set some other resources for line colors and grid spacing.
> ;
> opts_at_mpGeophysicalLineColor = get_res_value_keep(opts,
> "mpGeophysicalLineColor","Gray")
> opts_at_mpGeophysicalLineThicknessF = get_res_value_keep(opts,
> "mpGeophysicalLineThicknessF",0.5)
> opts_at_mpGridLineColor = get_res_value_keep(opts,
> "mpGridLineColor","Gray")
> opts_at_mpGridLineThicknessF = get_res_value_keep(opts,
> "mpGridLineThicknessF",0.5)
> opts_at_mpGridMaskMode = get_res_value_keep(opts,
> "mpGridMaskMode",3)
> opts_at_mpGridSpacingF = get_res_value_keep(opts,
> "mpGridSpacingF",10)
> opts_at_mpLimbLineColor = get_res_value_keep(opts,
> "mpLimbLineColor","Gray")
> opts_at_mpLimbLineThicknessF = get_res_value_keep(opts,
> "mpLimbLineThicknessF",0.5)
> opts_at_mpNationalLineColor = get_res_value_keep(opts,
> "mpNationalLineColor","Gray")
> opts_at_mpNationalLineThicknessF = get_res_value_keep(opts,
> "mpNationalLineThicknessF",0.5)
> opts_at_mpPerimLineColor = get_res_value_keep(opts,
> "mpPerimLineColor","Gray")
> opts_at_mpPerimOn = get_res_value_keep(opts,
> "mpPerimOn",True)
> opts_at_mpUSStateLineColor = get_res_value_keep(opts,
> "mpUSStateLineColor","Gray")
> opts_at_mpUSStateLineThicknessF = get_res_value_keep(opts,
> "mpUSStateLineThicknessF",0.5)
> opts_at_pmTickMarkDisplayMode = get_res_value_keep(opts,
> "pmTickMarkDisplayMode","Always")
>
> ;
> ; Tick mark resources
> ;
> ;opts_at_tmXBMajorLengthF = get_res_value(opts,
> "tmXBMajorLengthF",-0.03)
> ;opts_at_tmYLMajorLengthF = get_res_value(opts,
> "tmYLMajorLengthF",-0.03)
> opts_at_tmXTOn = get_res_value(opts,"tmXTOn",False)
> opts_at_tmYROn = get_res_value(opts,"tmYROn",False)
> opts_at_tmYRLabelsOn =
> get_res_value(opts,"tmYRLabelsOn",True)
> opts_at_tmXBBorderOn =
> get_res_value(opts,"tmXBBorderOn",True)
> opts_at_tmXTBorderOn =
> get_res_value(opts,"tmXTBorderOn",True)
> opts_at_tmYLBorderOn =
> get_res_value(opts,"tmYLBorderOn",True)
> opts_at_tmYRBorderOn =
> get_res_value(opts,"tmYRBorderOn",True)
> ;
> end if
>
> ;
> ; LambertConformal projection
> ;
> if(in_file_at_MAP_PROJ .eq. 1)
> projection = "LambertConformal"
>
> opts_at_mpLambertParallel1F = get_res_value_keep(opts,
> "mpLambertParallel1F",in_file_at_TRUELAT1)
> opts_at_mpLambertParallel2F = get_res_value_keep(opts,
> "mpLambertParallel2F",in_file_at_TRUELAT2)
> if(isatt(in_file,"STAND_LON"))
> opts_at_mpLambertMeridianF = get_res_value_keep(opts,
> "mpLambertMeridianF",in_file_at_STAND_LON)
> else
> if(isatt(in_file,"CEN_LON"))
> opts_at_mpLambertMeridianF = get_res_value_keep(opts,
> "mpLambertMeridianF",in_file_at_CEN_LON)
> else
> print("ERROR: Found neither STAND_LON or CEN_LON in file")
> end if
> end if
> end if
>
> ;
> ; Stereographic projection
> ;
> if(in_file_at_MAP_PROJ .eq. 2)
> projection = "Stereographic"
> opts_at_mpCenterLatF = get_res_value_keep(opts, "mpCenterLatF",
> in_file_at_CEN_LAT)
> if(isatt(in_file,"STAND_LON"))
> opts_at_mpCenterLonF = get_res_value_keep(opts,
> "mpCenterLonF",in_file_at_STAND_LON)
> else
> if(isatt(in_file,"CEN_LON"))
> opts_at_mpCenterLonF = get_res_value_keep(opts,
> "mpCenterLonF",in_file_at_CEN_LON)
> else
> print("ERROR: Found neither STAND_LON or CEN_LON in file")
> end if
> end if
> end if
>
> ;
> ; Mercator projection
> ;
> if(in_file_at_MAP_PROJ .eq. 3)
> projection = "Mercator"
> opts_at_mpCenterLatF = get_res_value_keep(opts, "mpCenterLatF",
> 0.0)
> if(isatt(in_file,"STAND_LON"))
> opts_at_mpCenterLonF = get_res_value_keep(opts,
> "mpCenterLonF",in_file_at_STAND_LON)
> else
> if(isatt(in_file,"CEN_LON"))
> opts_at_mpCenterLonF = get_res_value_keep(opts,
> "mpCenterLonF",in_file_at_CEN_LON)
> else
> print("ERROR: Found neither STAND_LON or CEN_LON in file")
> end if
> end if
> end if
>
> ;
> ; global WRF CylindricalEquidistant
> ;
> if(in_file_at_MAP_PROJ .eq. 6)
> projection = "CylindricalEquidistant"
> opts_at_mpGridSpacingF = 45
> opts_at_mpCenterLonF = get_res_value_keep(opts,
> "mpCenterLonF",in_file_at_CEN_LON)
> if( isatt(in_file,"POLE_LAT") ) then
> opts_at_mpCenterRotF = get_res_value_keep(opts, "mpCenterRotF",
> 90.0 - in_file_at_POLE_LAT)
> end if
> end if
>
> ;
> ; CylindricalEquidistant
> ;
> if(in_file_at_MAP_PROJ .eq. 0)
> projection = "CylindricalEquidistant"
> opts_at_mpCenterLatF = get_res_value_keep(opts, "mpCenterLatF",
> 0.0)
> if(isatt(in_file,"STAND_LON"))
> opts_at_mpCenterLonF = get_res_value_keep(opts,
> "mpCenterLonF",in_file_at_STAND_LON)
> else
> if(isatt(in_file,"CEN_LON"))
> opts_at_mpCenterLonF = get_res_value_keep(opts,
> "mpCenterLonF",in_file_at_CEN_LON)
> else
> print("ERROR: Found neither STAND_LON or CEN_LON in file")
> end if
> end if
> end if
> else
> print("wrf_map: Error: no MAP_PROJ attribute in input file")
> return(new(1,graphic))
> end if
>
> ;
> ; The default is not to draw the plot or advance the frame, and
> ; to maximize the plot in the frame.
> ;
> opts_at_gsnDraw = get_res_value_keep(opts,"gsnDraw", False)
> opts_at_gsnFrame = get_res_value_keep(opts,"gsnFrame", False)
> opts_at_gsnMaximize = get_res_value_keep(opts,"gsnMaximize", True)
>
> delete_attrs(opts) ; Clean up.
> mp = gsn_map(wks,projection,opts) ; Create map plot.
> return(mp) ; Return.
>
> end
>
> ;--------------------------------------------------------------------------------
>
> undef("wrf_map_overlays")
> function wrf_map_overlays(in_file[1]:file, \
> wks:graphic, \
> plots[*]:graphic, \
> opt_arg[1]:logical, \
> opt_mp[1]:logical)
>
> ; This procedure takes an array of plots and overlays them on a
> ; base plot - map background.
> ;
> ; It will advance the plot and cleanup, unless you set the
> ; PanelPlot resource to True.
> ;
> ; Attributes recognized by this procedure:
> ; FramePlot
> ; PanelPlot
> ; NoTitles (don't do any titles)
> ; CommonTitle & PlotTile is used to overwrite field titles
> ; CommonTitle will super-seed NoTitles
> ;
> ; If FramePlot False, then Draw the plot but do not Frame.
> ; In this case a user want to add to the drawing, and will
> ; have to advance the Frame manually in the script.
> ;
> ; If the "NoTitles" attribute exists and is set True, then
> ; don't create the top-left titles, and leave the main titles alone.
> ; This resource can be useful if you are planning to panel
> ; the plots.
> ;
> ; If PanelPlot is set to True, then this flags to wrf_map_overlays
> ; that these plots are going to be eventually paneled (likely
> ; by gsn_panel), and hence 1) draw and frame should not be called
> ; (unless gsnDraw and/or gsnFrame are explicitly set to True),
> ; and 2) the overlays and titles should not be removed with
> ; NhlRemoveOverlay and NhlRemoveAnnotation.
> ;
> begin
>
> ; Let's make the map first
> base = wrf_map(wks,in_file,opt_mp)
>
> opts = opt_arg ; Make a copy of the resource list
>
> no_titles = get_res_value(opts,"NoTitles",False) ; Do we want
> field titles?
> com_title = get_res_value(opts,"CommonTitle",False) ; Do we
> have a common title?
> if ( com_title ) then
> plot_title = get_res_value(opts,"PlotTitle"," ")
> no_titles = True
> end if
>
> call_draw = True
> call_frame = get_res_value(opts,"FramePlot",True) ; Do we want
> to frame the plot?
> panel_plot = get_res_value(opts,"PanelPlot",False) ; Are we
> paneling?
> opts_at_gsnMaximize = get_res_value_keep(opts,"gsnMaximize", True)
>
> nplots = dimsizes(plots)
> ; font_color = "Black"
>
> do i=0,nplots-1
> if(.not.ismissing(plots(i))) then
> ; class_name = NhlClassName(plots(i))
> ; print(class_name)
> ; if(class_name.eq."contourPlotClass") then
> ; getvalues plots(i)
> ; "cnFillOn" : fill_on
> ; "cnLineColor" : line_color
> ; end getvalues
> ; if (.not.fill_on) then
> ; font_color = line_color
> ; end if
> ; end if
> if(.not.no_titles) then
> getvalues plots(i)
> "tiMainString" : SubTitle
> end getvalues
> if(i.eq.0) then
> SubTitles = SubTitle
> else
> SubTitles = SubTitles + "~C~" + SubTitle
> end if
> end if
> if(com_title .and. i .eq. nplots-1) then
> getvalues plots(i)
> "tiMainString" : SubTitle
> end getvalues
> SubTitles = plot_title
> end if
> setvalues plots(i)
> "tfDoNDCOverlay" : True
> "tiMainOn" : False
> end setvalues
> overlay(base,plots(i))
> else
> print("wrf_map_overlays: Warning: overlay plot #" + i + " is
> not valid.")
> end if
> end do
>
> if(.not.no_titles .or. com_title) then
> font_height = get_res_value_keep(opts,"FontHeightF",0.01)
> txt = create "map_titles" textItemClass wks
> "txString" : SubTitles
> "txFontHeightF" : font_height
> ;"txFontColor" : font_color
> end create
> anno = NhlAddAnnotation(base,txt)
> setvalues anno
> "amZone" : 3
> "amJust" : "BottomLeft"
> "amSide" : "Top"
> "amParallelPosF" : 0.005
> "amOrthogonalPosF" : 0.03
> "amResizeNotify" : False
> end setvalues
> base_at_map_titles = anno
> end if
> ;
> ; gsnDraw and gsnFrame default to False if panel plot.
> ;
> if(panel_plot) then
> call_draw = False
> call_frame= False
> end if
>
>
> opts_at_gsnDraw = get_res_value_keep(opts,"gsnDraw", call_draw)
> opts_at_gsnFrame = get_res_value_keep(opts,"gsnFrame", call_frame)
>
> draw_and_frame(wks,base,opts_at_gsnDraw,opts_at_gsnFrame,False, \
> opts_at_gsnMaximize)
>
> if(.not.panel_plot) then
> do i=0,nplots-1
> if(.not.ismissing(plots(i))) then
> NhlRemoveOverlay(base,plots(i),False)
> else
> print("wrf_remove_map_overlays: Warning: overlay plot #" + i
> + " is not valid.")
> print(" Nothing to remove.")
> end if
> end do
> end if
>
> if(.not.no_titles.and..not.panel_plot) then
> if(isatt(base,"map_titles")) then
> NhlRemoveAnnotation(base,base_at_map_titles)
> delete(base_at_map_titles)
> end if
> end if
>
> return(base)
> end
>
> ;--------------------------------------------------------------------------------
>
> undef("wrf_overlays")
> function wrf_overlays(in_file[1]:file, \
> wks:graphic, plots[*]:graphic, \
> opt_arg[1]:logical)
>
> ; This procedure takes an array of plots and overlays them.
> ;
> ; It will advance the plot and cleanup, unless you set the
> ; PanelPlot resource to True.
> ;
> ; Attributes recognized by this procedure:
> ; FramePlot
> ; PanelPlot
> ; NoTitles (don't do any titles)
> ; CommonTitle & PlotTile is used to overwrite field titles
> ; CommonTitle will super-seed NoTitles
> ;
> ; If FramePlot False, then Draw the plot but do not Frame.
> ; In this case a user want to add to the drawing, and will
> ; have to advance the Frame manually in the script.
> ;
> ; If the "NoTitles" attribute exists and is set True, then
> ; don't create the top-left titles, and leave the main titles alone.
> ; This resource can be useful if you are planning to panel
> ; the plots.
> ;
> ; If PanelPlot is set to True, then this flags to wrf_overlays
> ; that these plots are going to be eventually paneled (likely
> ; by gsn_panel), and hence 1) draw and frame should not be called
> ; (unless gsnDraw and/or gsnFrame are explicitly set to True),
> ; and 2) the overlays and titles should not be removed with
> ; NhlRemoveOverlay and NhlRemoveAnnotation.
> ;
> begin
> opts = opt_arg ; Make a copy of the resource list.
>
> no_titles = get_res_value(opts,"NoTitles",False) ; Do we want
> field titles?
> com_title = get_res_value(opts,"CommonTitle",False) ; Do we
> have a common title?
> if ( com_title ) then
> plot_title = get_res_value(opts,"PlotTitle"," ")
> no_titles = True
> end if
>
> call_draw = True
> call_frame = get_res_value(opts,"FramePlot",True) ; Do we want
> to frame the plot?
> panel_plot = get_res_value(opts,"PanelPlot",False) ; Are we
> paneling?
> opts_at_gsnMaximize = get_res_value_keep(opts,"gsnMaximize", True)
>
> nplots = dimsizes(plots)
>
> base = plots(0)
> if(.not.no_titles) then
> getvalues plots(0)
> "tiMainString" : SubTitle
> end getvalues
> SubTitles = SubTitle
> setvalues plots(0)
> "tfDoNDCOverlay" : True
> "tiMainOn" : False
> end setvalues
> else
> setvalues plots(0)
> "tfDoNDCOverlay" : True
> end setvalues
> end if
>
> if (nplots.eq.1) then
> blank = create "BlankPlot" logLinPlotClass wks
> ;"cnConstFLabelOn" : False
> end create
> overlay(base,blank)
> end if
>
> do i=1,nplots-1
> if(.not.ismissing(plots(i))) then
> if(.not.no_titles) then
> getvalues plots(i)
> "tiMainString" : SubTitle
> end getvalues
> if(i.eq.0) then
> SubTitles = SubTitle
> else
> SubTitles = SubTitles + "~C~" + SubTitle
> end if
> end if
> if(com_title .and. i .eq. nplots-1) then
> getvalues plots(i)
> "tiMainString" : SubTitle
> end getvalues
> SubTitles = plot_title
> end if
> setvalues plots(i)
> "tfDoNDCOverlay" : True
> "tiMainOn" : False
> end setvalues
> overlay(base,plots(i))
> else
> print("wrf_overlays: Warning: overlay plot #" + i + " is not
> valid.")
> end if
> end do
>
> if(.not.no_titles .or. com_title) then
> font_height = get_res_value_keep(opt_arg,"FontHeightF",0.01)
>
> txt = create "map_titles" textItemClass wks
> "txString" : SubTitles
> "txFontHeightF" : font_height
> end create
> anno = NhlAddAnnotation(base,txt)
> setvalues anno
> "amZone" : 3
> "amJust" : "BottomLeft"
> "amSide" : "Top"
> "amParallelPosF" : 0.005
> "amOrthogonalPosF" : 0.03
> "amResizeNotify" : False
> end setvalues
> base_at_map_titles = anno
> end if
>
> ;
> ; gsnDraw and gsnFrame should default to True if not a panel plot.
> ;
> if(panel_plot) then
> call_draw = False
> call_frame= False
> end if
>
> opts_at_gsnDraw = get_res_value_keep(opts,"gsnDraw", call_draw)
> opts_at_gsnFrame = get_res_value_keep(opts,"gsnFrame", call_frame)
> opts_at_gsnMaximize = get_res_value_keep(opts,"gsnMaximize", True)
>
> draw_and_frame(wks,base,opts_at_gsnDraw,opts_at_gsnFrame,False, \
> opts_at_gsnMaximize)
>
> if(.not.no_titles.and..not.panel_plot) then
> NhlRemoveAnnotation(base,base_at_map_titles)
> delete(base_at_map_titles)
> end if
>
> if(.not.panel_plot) then
> if ( nplots .ge. 2 ) then
> do i=1,nplots-1
> if(.not.ismissing(plots(i))) then
> NhlRemoveOverlay(base,plots(i),False)
> else
> print("wrf_remove_overlays: Warning: overlay plot #" + i +
> " is not valid.")
> print(" Nothing to remove.")
> end if
> end do
> end if
> end if
>
> return(base)
> end
>
> ;--------------------------------------------------------------------------------
>
> undef("wrf_map_zoom")
> function
> wrf_map_zoom(wks[1]:graphic,in_file[1]:file,opt_args[1]:logical, \
> y1:integer,y2:integer,x1:integer,x2:integer)
>
> ; As of version 5.0.1, this routine is redundant. Use the special
> "ZoomIn"
> ; resource in wrf_map to accomplish the same thing. This function is
> ; being kept for backwards capability. There should be no need for it
> ; except to run old WRF-NCL codes. Do not make any changes to it
> except
> ; possibly to fix bugs.
> ;
> begin
> print("wrf_map_zoom: Warning: This function is obsolete.
> Consider using")
> print(" the 'ZoomIn' resource in wrf_map
> instead.")
> if(isfilevar(in_file,"XLAT"))
> lat = in_file->XLAT(0,:,:)
> lon = in_file->XLONG(0,:,:)
> else
> lat = in_file->XLAT_M(0,:,:)
> lon = in_file->XLONG_M(0,:,:)
> end if
> opts = opt_args ; Make a copy of the resource list
> opts = True
> opts_at_mpLeftCornerLatF = lat(y1,x1)
> opts_at_mpLeftCornerLonF = lon(y1,x1)
> opts_at_mpRightCornerLatF = lat(y2,x2)
> opts_at_mpRightCornerLonF = lon(y2,x2)
> mz = wrf_map(wks,in_file,opts)
> return(mz)
> end
>
> ;--------------------------------------------------------------------------------
>
> undef("wrf_map_overlay")
> procedure wrf_map_overlay(wks:graphic,base[1]:graphic, \
> plots[*]:graphic, \
> opt_arg[1]:logical)
>
> ; As of version 5.0.1, this procedure is obsolete. Use
> wrf_map_overlays
> ; instead. It is being kept for backwards capability. Do not make any
> ; changes to it except possibly to fix bugs.
> ;
> ; This procedure takes an array of plots and overlays them on a
> ; base plot - map background.
> ;
> ; It will advance the plot and cleanup, unless you set the
> ; PanelPlot resource to True.
> ;
> ; Attributes recognized by this procedure:
> ; NoTitles (don't do any titles)
> ; PanelPlot
> ;
> ; If the "NoTitles" attribute exists and is set True, then
> ; don't create the top-left titles, and leave the main titles alone.
> ; This resource can be useful if you are planning to panel
> ; the plots.
> ;
> ; If PanelPlot is set to True, then this flags to wrf_map_overlay
> ; that these plots are going to be eventually paneled (likely
> ; by gsn_panel), and hence 1) draw and frame should not be called
> ; (unless gsnDraw and/or gsnFrame are explicitly set to True),
> ; and 2) the overlays and titles should not be removed with
> ; NhlRemoveOverlay and NhlRemoveAnnotation.
> ;
> begin
> print("wrf_map_overlay: Warning: This procedure is obsolete.
> Consider" )
> print(" using wrf_map_overlays instead.")
>
> opts = opt_arg ; Make a copy of the resource list
> opts = True
>
> if(opts.and.isatt(opts,"NoTitles").and.opts_at_NoTitles) then
> no_titles = True
> else
> no_titles = False
> end if
>
> panel_plot = get_res_value(opts,"PanelPlot",False) ; Are we
> paneling?
>
> nplots = dimsizes(plots)
> ; font_color = "Black"
>
> do i=0,nplots-1
> if(.not.ismissing(plots(i))) then
> ; class_name = NhlClassName(plots(i))
> ; print(class_name)
> ; if(class_name.eq."contourPlotClass") then
> ; getvalues plots(i)
> ; "cnFillOn" : fill_on
> ; "cnLineColor" : line_color
> ; end getvalues
> ; if (.not.fill_on) then
> ; font_color = line_color
> ; end if
> ; end if
> if(.not.no_titles) then
> getvalues plots(i)
> "tiMainString" : SubTitle
> end getvalues
> if(i.eq.0) then
> SubTitles = SubTitle
> else
> SubTitles = SubTitles + "~C~" + SubTitle
> end if
> setvalues plots(i)
> "tfDoNDCOverlay" : True
> "tiMainOn" : False
> end setvalues
> else
> setvalues plots(i)
> "tfDoNDCOverlay" : True
> end setvalues
> end if
> overlay(base,plots(i))
> else
> print("wrf_map_overlay: Warning: overlay plot #" + i + " is not
> valid.")
> end if
> end do
>
> if(.not.no_titles) then
> font_height = get_res_value_keep(opt_arg,"FontHeightF",0.01)
> txt = create "map_titles" textItemClass wks
> "txString" : SubTitles
> "txFontHeightF" : font_height
> ;"txFontColor" : font_color
> end create
> anno = NhlAddAnnotation(base,txt)
> setvalues anno
> "amZone" : 3
> "amJust" : "BottomLeft"
> "amSide" : "Top"
> "amParallelPosF" : 0.005
> "amOrthogonalPosF" : 0.03
> "amResizeNotify" : False
> end setvalues
> base_at_map_titles = anno
> end if
> ;
> ; gsnDraw and gsnFrame should default to True if not a panel plot.
> ; gsnFrame will default to False if opt_arg is False.
> ;
> if(panel_plot.or..not.opt_arg) then
> call_frame = False
> else
> call_frame = True
> end if
> if(panel_plot) then
> call_draw = False
> else
> call_draw = True
> end if
> opts_at_gsnDraw = get_res_value_keep(opts,"gsnDraw", call_draw)
> opts_at_gsnFrame = get_res_value_keep(opts,"gsnFrame", call_frame)
> opts_at_gsnMaximize = get_res_value_keep(opts,"gsnMaximize", True)
>
> draw_and_frame(wks,base,opts_at_gsnDraw,opts_at_gsnFrame,False, \
> opts_at_gsnMaximize)
>
> if(.not.panel_plot) then
> do i=0,nplots-1
> if(.not.ismissing(plots(i))) then
> NhlRemoveOverlay(base,plots(i),False)
> else
> print("wrf_remove_map_overlay: Warning: overlay plot #" + i +
> " is not valid.")
> print(" Nothing to remove.")
> end if
> end do
> end if
>
> if(.not.no_titles.and..not.panel_plot) then
> if(isatt(base,"map_titles")) then
> NhlRemoveAnnotation(base,base_at_map_titles)
> delete(base_at_map_titles)
> end if
> end if
> end
>
> ;--------------------------------------------------------------------------------
>
> undef("wrf_overlay")
> procedure wrf_overlay(wks:graphic, plots[*]:graphic, \
> opt_arg[1]:logical)
>
> ; As of version 5.0.1, this procedure is obsolete. Use wrf_overlays
> ; instead. It is being kept for backwards capability. Do not make any
> ; changes to it except possibly to fix bugs.
> ;
> ; This procedure takes an array of plots and overlays them.
> ;
> ; It will advance the plot and cleanup, unless you set the
> ; PanelPlot resource to True.
> ;
> ; Attributes recognized by this procedure:
> ; NoTitles (don't do any titles)
> ; PanelPlot
> ;
> ; If the "NoTitles" attribute exists and is set True, then
> ; don't create the top-left titles, and leave the main titles alone.
> ; This resource can be useful if you are planning to panel
> ; the plots.
> ;
> ; If PanelPlot is set to True, then this flags to wrf_overlay
> ; that these plots are going to be eventually paneled (likely
> ; by gsn_panel), and hence 1) draw and frame should not be called
> ; (unless gsnDraw and/or gsnFrame are explicitly set to True),
> ; and 2) the overlays and titles should not be removed with
> ; NhlRemoveOverlay and NhlRemoveAnnotation.
> ;
> begin
> print("wrf_overlay: Warning: This procedure is obsolete. Consider
> using")
> print(" wrf_overlays instead.")
>
> opts = opt_arg ; Make a copy of the resource list.
> opts = True
>
> if(opts.and.isatt(opts,"NoTitles").and.opts_at_NoTitles) then
> no_titles = True
> else
> no_titles = False
> end if
>
> panel_plot = get_res_value(opts,"PanelPlot",False) ; Are we
> paneling?
>
> nplots = dimsizes(plots)
>
> base = plots(0)
> if(.not.no_titles) then
> getvalues plots(0)
> "tiMainString" : SubTitle
> end getvalues
> SubTitles = SubTitle
> setvalues plots(0)
> "tfDoNDCOverlay" : True
> "tiMainOn" : False
> end setvalues
> else
> setvalues plots(0)
> "tfDoNDCOverlay" : True
> end setvalues
> end if
>
> if (nplots.eq.1) then
> blank = create "BlankPlot" logLinPlotClass wks
> ;"cnConstFLabelOn" : False
> end create
> overlay(base,blank)
> end if
>
> do i=1,nplots-1
> if(.not.ismissing(plots(i))) then
> if(.not.no_titles) then
> getvalues plots(i)
> "tiMainString" : SubTitle
> end getvalues
> SubTitles = SubTitles + "~C~" + SubTitle
> setvalues plots(i)
> "tfDoNDCOverlay" : True
> "tiMainOn" : False
> end setvalues
> else
> setvalues plots(i)
> "tfDoNDCOverlay" : True
> end setvalues
> end if
> overlay(base,plots(i))
> else
> print("wrf_overlay: Warning: overlay plot #" + i + " is not
> valid.")
> end if
> end do
>
> if(.not.no_titles) then
> font_height = get_res_value_keep(opt_arg,"FontHeightF",0.01)
>
> txt = create "map_titles" textItemClass wks
> "txString" : SubTitles
> "txFontHeightF" : font_height
> end create
> anno = NhlAddAnnotation(base,txt)
> setvalues anno
> "amZone" : 3
> "amJust" : "BottomLeft"
> "amSide" : "Top"
> "amParallelPosF" : 0.005
> "amOrthogonalPosF" : 0.03
> "amResizeNotify" : False
> end setvalues
> base_at_map_titles = anno
> end if
>
> ;
> ; gsnDraw and gsnFrame should default to True if not a panel plot.
> ; gsnFrame will default to False if opt_arg is False.
> ;
> if(panel_plot.or..not.opt_arg) then
> call_frame = False
> else
> call_frame = True
> end if
> if(panel_plot) then
> call_draw = False
> else
> call_draw = True
> end if
> opts_at_gsnDraw = get_res_value_keep(opts,"gsnDraw", call_draw)
> opts_at_gsnFrame = get_res_value_keep(opts,"gsnFrame", call_frame)
> opts_at_gsnMaximize = get_res_value_keep(opts,"gsnMaximize", True)
>
> draw_and_frame(wks,base,opts_at_gsnDraw,opts_at_gsnFrame,False, \
> opts_at_gsnMaximize)
>
> if(.not.no_titles.and..not.panel_plot) then
> NhlRemoveAnnotation(base,base_at_map_titles)
> delete(base_at_map_titles)
> end if
>
> if(.not.panel_plot) then
> if ( nplots .ge. 2 ) then
> do i=1,nplots-1
> if(.not.ismissing(plots(i))) then
> NhlRemoveOverlay(base,plots(i),False)
> else
> print("wrf_remove_overlay: Warning: overlay plot #" + i + "
> is not valid.")
> print(" Nothing to remove.")
> end if
> end do
> end if
> end if
> end
>
> ;--------------------------------------------------------------------------------
>
> undef("add_white_space")
> function add_white_space(str:string,maxlen:integer)
>
> begin
> cstr = stringtochar(str)
> len = dimsizes(cstr)-1
> ws = ""
> if(len.lt.maxlen) then
> do i=1,maxlen-len
> ws = ws + " "
> end do
> end if
> return(ws)
>
> end
>
> ;--------------------------------------------------------------------------------
>
> undef("print_opts")
> procedure print_opts(opts_name,opts,debug)
>
> begin
> if(.not.debug) then
> return
> end if
> varatts = getvaratts(opts)
> ;
> ; Sort attributes alphabetically/
> ;
> sqsort(varatts)
> ;
> ; Get length of longest attribute name.
> ;
> cvaratts = stringtochar(varatts)
> cmaxlen = dimsizes(cvaratts(0,:))-1
>
> print("------------------------------------------------------------")
> print(opts_name + "...") ; Print name of option variable.
> ;
> ; Loop through each attribute in the list. If its value is an array,
> ; then print out the array with '(/' and '/)' around it.
> ;
> ; If the value is a string, then put ticks (') around each string.
> ;
> do i=0,dimsizes(varatts)-1
> x = opts@$varatts(i)$
> ;
> ; Use add_white_space to make sure all the equal signs line up.
> ;
> tmp_str = " " + varatts(i) + \
> add_white_space(varatts(i),cmaxlen) + " = "
> ;
> ; Check if attribute is an array.
> ;
> if(dimsizes(x).gt.1) then
> tmp_str = tmp_str + "(/"
> do j=0,dimsizes(x)-1
> if(typeof(x(j)).eq."string") then
> tmp_str = tmp_str + "'" + x(j) + "'"
> else
> tmp_str = tmp_str + x(j)
> end if
> if(j.lt.dimsizes(x)-1) then
> tmp_str = tmp_str + ","
> else
> tmp_str = tmp_str + "/)"
> end if
> end do
> else if(typeof(x).eq."string") then
> tmp_str = tmp_str + "'" + x + "'"
> else
> tmp_str = tmp_str + x
> end if
> end if
> print("" + tmp_str)
> delete(x)
> end do
>
> end
>
> ;--------------------------------------------------------------------------------
>
> undef("print_header")
> procedure print_header(icount:integer,debug)
> begin
> icount = icount + 1
> if(.not.debug) then
> return
> end if
> print("END plot #" + icount)
> print("------------------------------------------------------------")
>
> end
>
> ;--------------------------------------------------------------------------------
> _______________________________________________
> ncl-talk mailing list
> ncl-talk_at_ucar.edu
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk

_______________________________________________
ncl-talk mailing list
ncl-talk_at_ucar.edu
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Sun Sep 28 2008 - 16:52:07 MDT