drawing fronts

Hello,

I try to overlay weather fronts on some wrf output fields. To make a small test, I defined the front as:

  xxlat = (/50., 55.0/) ; Latitudes.
  xxlon = (/-5., 15./) ; Longitudes.

but the front doesn't seem to be located correctly. Some help would be appreciated. Thanks in advance.

PS: Here's the full script:

; Example script to produce plots for a WRF real-data run,
; with the ARW coordinate dynamics option.
; Interpolating to specified pressure levels

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"

begin
;
; The WRF ARW input file.
; This needs to have a ".nc" appended, so just do it.
  a = addfile("./wrfout_d01_2008-08-03_12:30:00.nc","r")

; We generate plots, but what kind do we prefer?
; type = "x11"
; type = "pdf"
 type = "ps"
; type = "ncgm"
  wks = gsn_open_wks(type,"fronts")

; Set some Basic Plot options
  res = True
  res@MainTitle = "REAL-TIME WRF"
  res@Footer = False

  pltres = True
  mpres = True
  mpres@mpGeophysicalLineColor = "Black"
  mpres@mpNationalLineColor = "Black"
  mpres@mpUSStateLineColor = "Black"
  mpres@mpGridLineColor = "Black"
  mpres@mpLimbLineColor = "Black"
  mpres@mpPerimLineColor = "Black"
  mpres@mpGeophysicalLineThicknessF = 2.0
  mpres@mpGridLineThicknessF = 2.0
  mpres@mpLimbLineThicknessF = 2.0
  mpres@mpNationalLineThicknessF = 2.0
  mpres@mpUSStateLineThicknessF = 2.0

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

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

; The specific pressure levels that we want the data interpolated to.
  pressure_levels = (/ 500./) ; pressure levels to plot
  nlevels = dimsizes(pressure_levels) ; number of pressure levels

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  it = 0 ; TIME LOOP

    res@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
    u = wrf_user_getvar(a,"ua",it) ; u averaged to mass points
    v = wrf_user_getvar(a,"va",it) ; v averaged to mass points
    p = wrf_user_getvar(a, "pressure",it) ; pressure is our vertical coordinate
    z = wrf_user_getvar(a, "z",it) ; grid point height
    rh = wrf_user_getvar(a,"rh",it) ; relative humidity

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

    do level = 0,nlevels-1 ; LOOP OVER LEVELS

      pressure = pressure_levels(level)

      z_plane = wrf_user_intrp3d( z,p,"h",pressure,0.,False)
      u_plane = wrf_user_intrp3d( u,p,"h",pressure,0.,False)
      v_plane = wrf_user_intrp3d( v,p,"h",pressure,0.,False)

      spd = (u_plane*u_plane + v_plane*v_plane)^(0.5) ; m/sec
      spd@description = "Wind Speed"
      spd@units = "m/s"
      u_plane = u_plane*1.94386 ; kts
      v_plane = v_plane*1.94386 ; kts
      u_plane@units = "kts"
      v_plane@units = "kts"

      ; Plotting options for Wind Speed
        opts = res
        opts@cnLineColor = "MediumSeaGreen"
        opts@ContourParameters = (/ 10. /)
        opts@cnInfoLabelOrthogonalPosF = 0.07 ; offset second label information
        opts@gsnContourLineThicknessesScale = 3.0
        contour_spd = wrf_contour(a,wks,spd,opts)
        delete(opts)

      ; Plotting options for Wind Vectors
        opts = res
        opts@FieldTitle = "Wind" ; overwrite Field Title
        opts@NumVectors = 47 ; wind barb density
        vector = wrf_vector(a,wks,u_plane,v_plane,opts)
        delete(opts)

      ; Plotting options for Geopotential Heigh
        opts_z = res
        opts_z@cnLineColor = "Black"
        opts_z@gsnContourLineThicknessesScale = 1.0

      ; MAKE PLOTS

; Draw a stationary front.

  wmsetp("ezf",1) ; Tell wmap we are using an existing map projection.
  wmsetp("lin",1.0) ; Line width of front curve.
  wmsetp("fro","stationnary"); Specify stationary front.
  wmsetp("cfc",4) ; Use blue for the triangles.
  wmsetp("wfc",2) ; Use red for the bumps.
  wmsetp("swi",0.02) ; Increase the size of the bumps and triangles.
  xxlat = (/50., 55.0/) ; Latitudes.
  xxlon = (/-5., 15./) ; Longitudes.
  wmdrft(wks, xxlat, xxlon)

          opts_z@ContourParameters = (/ 60.0 /)
          contour_height = wrf_contour(a,wks, z_plane,opts_z)
          plot = wrf_map_overlays(a,wks,(/contour_height, \
                                    vector/),pltres,mpres)
        delete(opts_z)

    end do ; END OF LEVEL LOOP

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

end

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Thu Jun 30 07:35:08 2011