;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; This is to plot 4 figures at one page, i.e., panel plot. 
;
;                  mchen, 2012/10
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 load "$NCARG_ROOT/lib/ncarg/nclex/gsun/gsn_code.ncl"
 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
 load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
; load "./WRFUserARW.ncl"

begin
 
  ; define experiment. exp_num=1 plot the first exp. exp_num=2 plot the two exp. If three more experiments needs to
  ; plot here, line_modes, line_markers, line_colors, legend_labels are needed to change to match exp_num.

  titlenum=4
  titlename=(/"Total Precip","Surface T and WIND","300mb Field","500mb Fields"/)
  plot=new(4,graphic)
  wks = gsn_open_wks("x11","2010061200")
  gsn_define_colormap(wks,"wh-bl-gr-ye-re")


; Basic Plot Information
    resB = True
    resB@MainTitle = " "
    resB@Footer = False

    pltres = True      ; 
    pltres@PanelPlot = True      ; 
    pltres@gsnDraw            = False  
    pltres@gsnFrame           = False  
    pltres@gsnMaximize        = True

    mpres = True

  DATADir = "./"
  FILES = systemfunc (" ls -1 " + DATADir + "wrfout_d01* ")
  numFILES = dimsizes(FILES)
  print("numFILES = " + numFILES)
  print(FILES)


  pressure_levels = (/ 500., 300./)   ; pressure levels to plot
  nlevels         = dimsizes(pressure_levels)     ; number of pressure levels


FirstTime = True
do ifil = 0, numFILES-1      ; BIG FILES LOOP
   a = addfile(FILES(ifil)+".nc","r")
     times = wrf_user_getvar(a,"times",-1)  ; get all times in the file
     ntimes = dimsizes(times)         ; number of times in the file

    print("Working on time: " + times(0) )
    resB@TimeLabel = times(0)      ; Set Valid time to use on plots
    if (FirstTime) then
      times_sav = times(0)
    end if

    res = resB
    res@cnFillOn = False

    it=0
    slp = wrf_user_getvar(a,"slp",it)       ; psl
    wrf_smooth_2d(slp, 3)                 ; smooth psl          
    td2 =  wrf_user_getvar(a,"td2",it)   ; Td2 in C
    tc2 = wrf_user_getvar(a,"T2",it)     ; T2 in Kelvin
       tc2 = tc2-273.16                  ; T2 in C
    u10 = wrf_user_getvar(a,"U10",it)    ; u at 10 m, mass point
    v10 = wrf_user_getvar(a,"V10",it)    ; v at 10 m, mass point

    tf2 = 1.8*tc2+32.                    ; Turn temperature into Fahrenheit
      tf2@description = "Surface Temperature"
      tf2@units = "F"
    td_f = 1.8*td2+32.                   ; Turn temperature into Fahrenheit
      td_f@description = "Surface Dew Point Temp"
      td_f@units = "F"
    u10 = u10*1.94386                    ; Turn wind into knots
    v10 = v10*1.94386
      u10@units = "kts"
      v10@units = "kts"
    tc = wrf_user_getvar(a,"tc",it)      ; 3D tc
    u  = wrf_user_getvar(a,"ua",it)      ; 3D U at mass points
    v  = wrf_user_getvar(a,"va",it)      ; 3D V at 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

    rain_exp = wrf_user_getvar(a,"RAINNC",it)
    rain_con = wrf_user_getvar(a,"RAINC",it)
    snow_exp = wrf_user_getvar(a,"SNOWNC",it)
    rain_tot = rain_exp + rain_con
    if( FirstTime ) then
      rain_exp_save = rain_exp
      rain_con_save = rain_con
      rain_tot_save = rain_tot
    end if
    rain_exp_tend = rain_exp - rain_exp_save
    rain_con_tend = rain_con - rain_con_save
    rain_tot_tend = rain_tot - rain_tot_save
    rain_exp_save = rain_exp
    rain_con_save = rain_con
    rain_tot_save = rain_tot

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

      pressure = pressure_levels(level)

      tc_plane = wrf_user_intrp3d(tc,p,"h",pressure,0.,False)
      z_plane  = wrf_user_intrp3d( z,p,"h",pressure,0.,False)
      rh_plane = wrf_user_intrp3d(rh,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 T                
      opts = res
      opts@cnFillOn = True
      opts@ContourParameters = (/ -65., 5., 5./)
      opts@gsnSpreadColorEnd = -3  ; End third from the last color in color map
      contour_tc = wrf_contour(a,wks,tc_plane,opts)
      delete(opts)

      ; Plotting options for RH                
        opts = res                          
        opts@cnFillOn = True  
        opts@pmLabelBarOrthogonalPosF = -0.1
        opts@ContourParameters = (/ 10., 90., 10./)
        opts@cnFillColors = (/"White","White","White", \
                              "White","Chartreuse","Green",\
                              "Green3","Green4", \
                              "ForestGreen","PaleGreen4"/)
        contour_rh = wrf_contour(a,wks,rh_plane,opts)
        delete(opts)


      ; 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
        opts@vcRefAnnoArrowLineColor   = "white"         ; change ref vector color
        opts@vcLineArrowColor        = "white"           ; change vector color
        vector = wrf_vector(a,wks,u_plane,v_plane,opts)
        delete(opts)


      ; Plotting options for Geopotential Heigh
        opts_z = res                          
        opts_z@cnLineColor = "Red"
        opts_z@gsnContourLineThicknessesScale = 3.0


      ; MAKE PLOTS                                       

        if ( pressure .eq. 850 ) then   ; plot temp, rh, height, wind barbs
          opts_z@ContourParameters = (/ 20.0 /)
          contour_height = wrf_contour(a,wks,z_plane,opts_z)
          plot(level) = wrf_map_overlays(a,wks,(/contour_rh,contour_tc,contour_height, \
                                    vector/),pltres,mpres)
        end if

        if ( pressure .eq. 700 ) then   ; plot temp, height, wind barbs
          opts_z@ContourParameters = (/ 30.0 /)
          contour_height = wrf_contour(a,wks, z_plane,opts_z)
          plot(level) = wrf_map_overlays(a,wks,(/contour_tc,contour_height, \
                                    vector/),pltres,mpres)
        end if

        if ( pressure .eq. 500 ) then   ; plot temp, height, wind barbs
          opts_z@ContourParameters = (/ 60.0 /)
          contour_height = wrf_contour(a,wks, z_plane,opts_z)
          plot(level) = wrf_map_overlays(a,wks,(/contour_tc,contour_height, \
                                    vector/),pltres,mpres)
        end if

        if ( pressure .eq. 300 ) then   ; plot windspeed, height, wind barbs
          opts_z@ContourParameters = (/ 60.0 /)
          contour_height = wrf_contour(a,wks, z_plane,opts_z)
          plot(level) = wrf_map_overlays(a,wks,(/contour_tc,contour_height, \
                                    vector/),pltres,mpres)
        end if
        delete(opts_z)

    end do      ; END OF LEVEL LOOP
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ; Plotting options for Sea Level Pressure
      opts = res          
      opts@ContourParameters = (/ 900., 1100., 10. /)
      opts@cnLineColor = "Blue"
      opts@cnInfoLabelOn = False
      ;opts@cnHighLabelsOn = True
      ;opts@cnLowLabelsOn = True
      opts@gsnContourLineThicknessesScale = 1.0
      contour_psl = wrf_contour(a,wks,slp,opts)
      delete(opts)
  
    ; Plotting options for Precipitation
      opts = res                        
      opts@cnLevelSelectionMode = "ExplicitLevels"
      opts@cnLevels = (/ .1, .2, .4, .8, 1.6, 3.2, 6.4, \
                         12.8, 25.6, 51.2, 102.4/)
      opts@cnFillColors = (/"White","White","DarkOliveGreen1", \
                            "DarkOliveGreen3","Chartreuse", \
                            "Chartreuse3","Green","ForestGreen", \
                            "Yellow","Orange","Red","Violet"/)
      opts@cnInfoLabelOn = False
      opts@cnConstFLabelOn = False
  
    ; Precipitation (color fill)
      opts@cnFillOn             = True
      opts@FieldTitle = "Total Precipitation"
      contour_tot = wrf_contour(a,wks, rain_tot, opts)
      delete(opts)
  
    ; MAKE PLOTS                                       
      plot(2) = wrf_map_overlays(a,wks,contour_tot,pltres,mpres)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    ; Plotting options for T                
      opts = res                         
      opts@cnFillOn = True  
      opts@ContourParameters = (/ -20., 90., 5./)
      opts@gsnSpreadColorEnd = -3  ; End third from the last color in color map
      contour_tc = wrf_contour(a,wks,tf2,opts)
      delete(opts)

    ; Plotting options for SLP                     
      opts = res         
      opts@cnLineColor = "Blue"
      opts@cnHighLabelsOn = True
      opts@cnLowLabelsOn = True
      opts@ContourParameters = (/ 900., 1100., 4. /)
      opts@cnLineLabelBackgroundColor = -1
      opts@gsnContourLineThicknessesScale = 2.0
      contour_psl = wrf_contour(a,wks,slp,opts)
      delete(opts)

    ; Plotting options for Wind Vectors                 
      opts = res         
      opts@FieldTitle = "Wind"       ; overwrite Field Title
      opts@NumVectors = 47           ; density of wind barbs
      vector = wrf_vector(a,wks,u10,v10,opts)
      delete(opts)

    ; MAKE PLOTS                                       
      plot(3) = wrf_map_overlays(a,wks,(/contour_tc,contour_psl,vector/),pltres,mpres)

end do  ;       END IF BIG FILE LOOP
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; panel plot
    pnlres                    =  True              ; plot mods desired
    pnlres@gsnPanelYWhiteSpacePercent = 5
    pnlres@gsnPanelXWhiteSpacePercent = 5
    gsn_panel(wks,(/plot/),(/2,2/),pnlres)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


end