NCL not plotting first 200 m in the vertical cross section

From: Basit Khan <basit.khan_at_nyahnyahspammersnyahnyah>
Date: Sun Feb 13 2011 - 15:15:52 MST

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

 

------------------------------------------------------------------------
-------

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_20
06-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.

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

Received on Sun Feb 13 15:16:36 2011

This archive was generated by hypermail 2.1.8 : Tue Feb 15 2011 - 09:43:19 MST