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

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_Surface2")

; Set some Basic Plot options
  ARWres = True
  ARWres@MainTitle                   = "REAL-TIME WRF"

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

; 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

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

    ARWres@TimeLabel = times(0)    ; Set Valid time to use on plots

    mpres = True                   ; Create map background
    map = wrf_map(wks,a,mpres)

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

    ; TERRAIN

      ter = wrf_user_getvar(a,"HGT",0)
      opts = ARWres
      opts@cnFillOn = True
      opts@ContourParameters = (/ 100. /)
      contour = wrf_contour(a,wks,ter,opts)
      delete(opts)
      wrf_map_overlay(wks,map,(/contour/),True)


    ; LANDUSE , SOIL & VEGETATION FRACTION

      opts = ARWres
      opts@cnFillOn = True
      opts@cnRasterModeOn = True
      opts@gsnSpreadColorStart = 32

      if(isfilevar(a,"LU_INDEX"))
        land = wrf_user_getvar(a,"LU_INDEX",0)
        contour = wrf_contour(a,wks,land,opts)
        wrf_map_overlay(wks,map,(/contour/),True)
      end if
  
      if(isfilevar(a,"ISLTYP"))
        isoil = wrf_user_getvar(a,"ISLTYP",0)
        contour = wrf_contour(a,wks,isoil,opts)
        wrf_map_overlay(wks,map,(/contour/),True)
      end if
  
      if(isfilevar(a,"VEGFRA"))
        vgf = wrf_user_getvar(a,"VEGFRA",0)
        opts@gsnSpreadColorStart= 0
        opts@ContourParameters = (/ 10. /);
        contour = wrf_contour(a,wks,vgf,opts)
        wrf_map_overlay(wks,map,(/contour/),True)
      end if
  
      delete(opts)
  
  
    ; SOIL TEMPERATURE

      if(isfilevar(a,"TSLB"))
        tsoil = wrf_user_getvar(a,"TSLB",0)
        dimsl = dimsizes(tsoil)
        do level = 0, dimsl(0)-1
          slevel = level+1
          opts = ARWres
          opts@cnFillOn = True
          opts@gsnSpreadColorEnd = -10
          opts@PlotLevelID = "Level  " + slevel
          contour = wrf_contour(a,wks,tsoil(level,:,:),opts)
          wrf_map_overlay(wks,map,(/contour/),True)
        end do
        delete(opts)
        delete(tsoil)
      end if


      opts = ARWres
      opts@cnFillOn = True
      opts@ContourParameters = (/ 240.,300.,5./)
      opts@gsnSpreadColorEnd = -10

      if(isfilevar(a,"ST000010"))
        tsoil = wrf_user_getvar(a,"ST000010",0)
        contour = wrf_contour(a,wks,tsoil,opts)
        wrf_map_overlay(wks,map,(/contour/),True)
        delete(tsoil)
      end if
      if(isfilevar(a,"ST010040"))
        tsoil = wrf_user_getvar(a,"ST010040",0)
        contour = wrf_contour(a,wks,tsoil,opts)
        wrf_map_overlay(wks,map,(/contour/),True)
        delete(tsoil)
      end if
      if(isfilevar(a,"ST040100"))
        tsoil = wrf_user_getvar(a,"ST040100",0)
        contour = wrf_contour(a,wks,tsoil,opts)
        wrf_map_overlay(wks,map,(/contour/),True)
        delete(tsoil)
      end if
      if(isfilevar(a,"ST100200"))
        tsoil = wrf_user_getvar(a,"ST100200",0)
        contour = wrf_contour(a,wks,tsoil,opts)
        wrf_map_overlay(wks,map,(/contour/),True)
        delete(tsoil)
      end if

      delete(opts)

  
    ; SOIL MOISTURE
    
      if(isfilevar(a,"SMOIS"))
        sf_sfc_phys  = a@SF_SURFACE_PHYSICS
        if ( sf_sfc_phys .eq. 2 ) then         ; only available from LSM option
          soilm = wrf_user_getvar(a,"SMOIS",0)
          dimsl = dimsizes(soilm)
          do level = 0, dimsl(0)-1
            slevel = level+1
            opts = ARWres
            opts@cnFillOn = True
            opts@gsnSpreadColorEnd = -10
            opts@PlotLevelID = "Level  " + slevel
            contour = wrf_contour(a,wks,msoil(level,:,:),opts)
            wrf_map_overlay(wks,map,(/contour/),True)
          end do
          delete(opts)
          delete(soilm)
        end if
      end if


      opts = ARWres
      opts@cnFillOn = True
      opts@gsnSpreadColorEnd = -10

      if(isfilevar(a,"SM000010"))
        soilm = wrf_user_getvar(a,"SM000010",0)
        contour = wrf_contour(a,wks,soilm,opts)
        wrf_map_overlay(wks,map,(/contour/),True)
        delete(soilm)
      end if
      if(isfilevar(a,"SM010040"))
        soilm = wrf_user_getvar(a,"SM010040",0)
        contour = wrf_contour(a,wks,soilm,opts)
        wrf_map_overlay(wks,map,(/contour/),True)
        delete(soilm)
      end if
      if(isfilevar(a,"SM040100"))
        soilm = wrf_user_getvar(a,"SM040100",0)
        contour = wrf_contour(a,wks,soilm,opts)
        wrf_map_overlay(wks,map,(/contour/),True)
        delete(soilm)
      end if
      if(isfilevar(a,"SM100200"))
        soilm = wrf_user_getvar(a,"SM100200",0)
        contour = wrf_contour(a,wks,soilm,opts)
        wrf_map_overlay(wks,map,(/contour/),True)
        delete(soilm)
      end if

      delete(opts)

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

end