Re: NCL not plotting first 200 m in the vertical cross section

From: Song Guiting (Dr) <guitingsong_at_nyahnyahspammersnyahnyah>
Date: Mon Mar 14 2011 - 02:00:43 MDT

Dear Mary,

I also encounter the same problem. My NCL version is 6.00 beta.

Best regards,
Guiting

From: ncl-talk-bounces@ucar.edu [mailto:ncl-talk-bounces@ucar.edu] On Behalf Of Mary Haley
Sent: Monday, February 14, 2011 11:18 PM
To: Basit Khan
Cc: ncl-talk@ucar.edu
Subject: Re: NCL not plotting first 200 m in the vertical cross section

This might be due to a bug we fixed awhile back. What version of NCL (ncl -V) are you running?
If you can provide me with the file, then I can try your script and see if the bug is fixed.

--Mary

On Feb 13, 2011, at 3:15 PM, Basit Khan wrote:

Hi,
I am trying to plot vertical cross section (height x latitude) from a WRF output file. The problem is that NCL not plotting anything in the first 200 m while I have at least 11 vertical levels in the first 200 m asl. I am using terrain following coordinates system and most of the cross section area is either sea or terrain less than 50 m high asl. I will greatly appreciate if somebody could help me fix this problem. The ncl code is copied below and the output of the ncl code is attached.

Regards,

Basit
Centre for Atmospheric Research
Department of Geography,
University of Canterbury, Private Bag 4800,
Christchurch 8140, New Zealand.
Tel: +643 364 2987 ext. 3088
email: basit.khan@canterbury.ac.nz<mailto:basit.khan@canterbury.ac.nz>

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

begin
;
  a = addfile("/hpc/bluefern/WRF/basit/basit-archive/b-haze/bh30/wrfout_d03_2006-05-30_12:00: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,"plt_CrossSection4")

; Set some basic resources
  res = True
  res@MainTitle = "REAL-TIME WRF"
  res@Footer = False

  pltres = True

;;;;;;;;;;;;;;;; HIGHRES RANGS/GSHHS COASTLINE PLOT RESOURCES ;;;;;;;;;;;;;;;;;;;;;;;;
  mpres =True
  mpres@gsnMaximize = True
  mpres@mpOutlineOn = True
  mpres@mpFillOn = False
  mpres@mpProjection = "Mercator"
  mpres@mpLimitMode = "LatLon"
  mpres@mpMinLonF = 174.4061
  mpres@mpMaxLonF = 175.2502
  mpres@mpMinLatF = -37.27655
  mpres@mpMaxLatF = -36.60195
  mpres@mpOutlineOn = True
  mpres@mpGridLineColor = -1
  mpres@gsnDraw = False ; don't draw the plots
  mpres@gsnFrame = False ; don't advance the frame
  mpres@tmYROn = False
  mpres@tmXTOn = False
  mpres@mpGeophysicalLineColor= "Red"
  mpres@mpLimbLineColor = "Blue"
  mpres@mpGeophysicalLineThicknessF = 1.5
  mpres@mpDataBaseVersion = "HighRes" ; Use the high-res database
  mpres@pmTickMarkDisplayMode = "Always" ; Turn on map tickmarks.
; mpres@tiMainString = "Using the RANGS-GSHHS coastal database"

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  FirstTime = True
  FirstTimeMap = True
  times = wrf_user_list_times(a) ; get times in the file
  ntimes = dimsizes(times) ; number of times in the file

  mdims = getfilevardimsizes(a,"P") ; get some dimension sizes for the file; print(mdims)=97x50x150x150
  nd = dimsizes(mdims) ; incase of 'P' the print(nd)=4 cause P has four dims

  xlat = wrf_user_getvar(a, "XLAT",0)
  xlon = wrf_user_getvar(a, "XLONG",0)
  ter = wrf_user_getvar(a, "HGT",0)

;---------------------------------------------------------------

  do it = 58,74,1 ; TIME LOOP

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

    tc = wrf_user_getvar(a,"tc",it) ; T in C
    rh = wrf_user_getvar(a,"rh",it) ; relative humidity
    z = wrf_user_getvar(a, "height",it) ; grid point height

    if ( FirstTime ) then ; get height info for labels
      zmin = 0
      zmax = 1. ; We are only interested in the first 1000 m
      nz = floattoint(zmax + 10)
    end if

;---------------------------------------------------------------

        plane = new(2,float)
        plane = (/ mdims(nd-1)/2, mdims(nd-2)/2 /) ; pivot point is center of domain (x,y)
        opts = False

         angle = 0.
         X_plane = wrf_user_intrp2d(xlat,plane,angle,opts)
         X_desc = "latitude"

        rh_plane = wrf_user_intrp3d(rh,z,"v",plane,angle,opts) ; "v" is for vertical and "h" is for horinzontal interporlation
        tc_plane = wrf_user_intrp3d(tc,z,"v",plane,angle,opts) ; Array of vertical coordninates to interpolate

      ; Find the index where 6km is - only need to do this once
        if ( FirstTime ) then

                  zz = wrf_user_intrp3d(z,z,"v",plane,angle,opts) ; printVarSummary(zz) = dimsizesV x Dimsize H = 100 x 149

                  b = ind(zz(:,0) .gt. zmax*1000. ) ; (zmax *1000 --> conver km into metres

                  zmax_pos = b(0) - 1
          if ( abs(zz(zmax_pos,0)-zmax*1000.) .lt. abs(zz(zmax_pos+1,0)-zmax*1000.) ) then
            zspan = b(0) - 1
          else
            zspan = b(0)
          end if
          delete(zz)
          delete(b)
          FirstTime = False
        end if

      ; X-axis lables
      dimsX = dimsizes(X_plane)
      xmin = X_plane(0)
      xmax = X_plane(dimsX(0)-1)
      xspan = dimsX(0)-1
      nx = floattoint( (xmax-xmin)/2 + 10)
; nx = (xmax-xmin)/2
      ;---------------------------------------------------------------
   ; print("the xmin value is "+xmin)
    ; print("the xmax value is "+xmax)

      ; Options for XY Plots
        opts_xy = res
        opts_xy@tiXAxisString = X_desc
        opts_xy@tiYAxisString = "Height (km)"
        opts_xy@cnMissingValPerimOn = True
        opts_xy@cnMissingValFillColor = 0
        opts_xy@cnMissingValFillPattern = 11
        opts_xy@tmXTOn = False
        opts_xy@tmYROn = False
        opts_xy@tmXBMode = "Explicit"
        opts_xy@tmXBValues = fspan(0,xspan,nx) ; Create tick marks
; opts_xy@tmXBValues = fspan(0,xspan,nx) ; Create tick marks
        opts_xy@tmXBLabels = sprintf("%6.2f",fspan(xmin,xmax,nx)) ; Create labels
        opts_xy@tmXBLabelFontHeightF = 0.015
        opts_xy@tmYLMode = "Explicit"
        opts_xy@tmYLValues = fspan(0,zspan,nz) ; Create tick marks
        opts_xy@tmYLLabels = sprintf("%3.1f",fspan(zmin,zmax,nz)) ; Create labels
        opts_xy@tiXAxisFontHeightF = 0.020
        opts_xy@tiYAxisFontHeightF = 0.020
        opts_xy@tmXBMajorLengthF = 0.02
        opts_xy@tmYLMajorLengthF = 0.02
        opts_xy@tmYLLabelFontHeightF = 0.015
        opts_xy@PlotOrientation = tc_plane@Orientation

      ; Plotting options for RH
        opts_rh = opts_xy
        opts_rh@ContourParameters = (/ 10., 90., 10. /)
        opts_rh@pmLabelBarOrthogonalPosF = -0.1
        opts_rh@cnFillOn = True
        opts_rh@cnFillColors = (/"White","White","White", \
                                            "White","Chartreuse","Green", \
                                            "Green3","Green4", \
                                            "ForestGreen","PaleGreen4"/)

      ; Plotting options for Temperature
        opts_tc = opts_xy
        opts_tc@cnInfoLabelZone = 1
        opts_tc@cnInfoLabelSide = "Top"
        opts_tc@cnInfoLabelPerimOn = True
        opts_tc@cnInfoLabelOrthogonalPosF = -0.00005
        opts_tc@ContourParameters = (/ 1. /)

      ; 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

        if (FirstTimeMap) then
          lat_plane = wrf_user_intrp2d(xlat,plane,angle,opts)
          lon_plane = wrf_user_intrp2d(xlon,plane,angle,opts)
          mpres = True
          pltres = True
          pltres@FramePlot = False
          optsM = res
          optsM@NoHeaderFooter = True
          optsM@cnFillOn = True
          optsM@lbTitleOn = False
          contour = wrf_contour(a,wks,ter,optsM)
          plot = wrf_map_overlays(a,wks,(/contour/),pltres,mpres)
          lnres = True
          lnres@gsLineThicknessF = 3.0
          lnres@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@FramePlot = True
       end if

       plot = wrf_overlays(a,wks,(/contour_rh,contour_tc/),pltres) ; plot x-section

  ; 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

This email may be confidential and subject to legal privilege, it may

not reflect the views of the University of Canterbury, and it is not

guaranteed to be virus free. If you are not an intended recipient,

please notify the sender immediately and erase all copies of the message

and any attachments.

Please refer to http://www.canterbury.ac.nz/emaildisclaimer for more

information.

<vcs_ncl.eps>_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk

________________________________
CONFIDENTIALITY: This email is intended solely for the person(s) named and may be confidential and/or privileged. If you are not the intended recipient, please delete it, notify us and do not copy, use, or disclose its content. Thank you.

Towards A Sustainable Earth: Print Only When Necessary

_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Mon Mar 14 02:01:08 2011

This archive was generated by hypermail 2.1.8 : Wed Mar 16 2011 - 09:22:37 MDT