;   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 from a a given point A to point B
;   Vertical coordinate is pressure

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

begin
;
; The WRF ARW input file.  
; This needs to have a ".nc" appended, so just do it.
  a = addfile("../wrfout_d01_2000-01-24_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_CrossSection3")

; Set some Basic Plot options
  ARWres = True
  ARWres@MainTitle          = "REAL-TIME WRF"
  ARWres@vpWidthF           = .9               ; overwrite basic plot size
  ARWres@vpHeightF          = 1.0

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

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

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

  do it = 0,ntimes-1               ; TIME LOOP

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

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need        

    th = wrf_user_getvar(a,"th",it)        ; theta
    rh = wrf_user_getvar(a,"rh",it)        ; relative humidity
    z  = wrf_user_getvar(a, "z",it)        ; grid point height
    p  = wrf_user_getvar(a, "pressure",it) ; grid point height

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

    if (FirstTime) then
      ; THIS IS NEEDED FOR Y LABLES - ALWAYS DO
      zz = z/100.
      dims = dimsizes(zz) 
      zmax = -999.
      do imax = 0,dims(0)-1
        if ( .not.ismissing(zz(imax,1,1)) .and. zz(imax,1,1) .gt. zmax ) then
          zmax = zz(imax,1,1)
        end if
      end do
      imax = floattoint(zmax)
      zmax = imax/10.
      print("Y labels set to: 0 - " + zmax + " km")
      FirstTime = False
      ; END OF ALWAYS DO
    end if



    do ip = 1, 1      	; we are just going to do one crossection for now

  ;   Specify the start and end points of the cross section              
  ;   For this case, set opts in wrf_user_intrp3d to True
        opts = True
        plane = new(4,float)

        if(ip .eq. 1) then
          plane = (/  10,30 ,  70,50  /) ; approx. start x;y and end x;y point
        end if
        if(ip .eq. 2) then
          plane = (/  130,1 , 130,162  /) ; approx. start x;y and end x;y point
        end if
        if(ip .eq. 3) then
          plane = (/  49,1 , 210,162  /) ; approx. start x;y and end x;y point
        end if

        ; Interpolate data vertically (in z)

        rh_plane = wrf_user_intrp3d(rh,p,"v",plane,0.,opts)
        th_plane = wrf_user_intrp3d(th,p,"v",plane,0.,opts)
        
  ; Options for XY Plots
        opts_xy                         = ARWres
        opts_xy@tiXAxisString           = "Number of Grid Points"
        opts_xy@tiYAxisString           = "Pressure (mb)"
        opts_xy@AspectRatio             = 0.75
        opts_xy@cnMissingValPerimOn     = True
        opts_xy@cnMissingValFillColor   = 0
        opts_xy@cnMissingValFillPattern = 11
        opts_xy@tmYLMode                = "Explicit"
        opts_xy@tmYLValues              = fspan(0,100,10) ; over the 0-100 span of the Y label
							  ; we are creating 10+1 tick marks
							  ; 100/10 + 1
        opts_xy@tmYLLabels              = fspan(1000,100,10); Corresponding 11 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         = th_plane@Orientation


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

  ; Plotting options for Temperature
        opts_th = opts_xy
        opts_th@ContourParameters  = (/ 3. /)



  ; Get the contour info for the rh and temp
        contour_th = wrf_contour(a,wks,th_plane,opts_th)
        contour_rh = wrf_contour(a,wks,rh_plane,opts_rh)


  ; MAKE PLOTS         
        wrf_overlay(wks,(/contour_rh,contour_th/),True)


  ; Delete options and fields, so we don't have carry over
        delete(opts_th)
        delete(opts_rh)
        delete(th_plane)
        delete(rh_plane)

    end do  ; make next cross section

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

  end do        ; END OF TIME LOOP

end