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

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/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")
  a = addfile("wrfout_d01_2016-10-28_00:00:00.nc","r")

; We generate plots, but what kind do we prefer?
;  type = "x11"
; type = "pdf"
 type = "ps"
; type = "ncgm"
  wks2 = gsn_open_wks(type,"zgrad_terrain_Oct282016")
  wks1 = gsn_open_wks(type,"zgrad_gradient_Oct282016")


; Set some basic resources
  res = True
  res@MainTitle = "RAIN SHADOW"

  pltres = True
  mpres = True

  mpres@ZoomIn = True        ; Tell wrf_map_resources we want to zoom in.
  mpres@Xstart = 120     ; Set these four special WRF resources
  mpres@Xend   = 160       ; required for zooming.
  mpres@Ystart = 50
  mpres@Yend   = 120

  mpres@mpGeophysicalLineColor = "Black"
  mpres@mpNationalLineColor    = "Black"
  mpres@mpUSStateLineColor     = "Black"
  mpres@mpGridLineColor        = "Black"
  mpres@mpLimbLineColor        = "Black"
  mpres@mpPerimLineColor       = "Black"

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

; What times and how many time steps are in the data set?
  FirstTime = True
  times = wrf_user_getvar(a,"times",-1)  ; get all times in the file
  ntimes = dimsizes(times)         ; number of times in the file

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

;  do it = 0,ntimes-1,2             ; TIME LOOP
  do it = 0,36             ; TIME LOOP

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

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
; First get the variables we will need      
; Terrain Z(x,y) static 2D
; u and v winds on sigma = 0.9 surface (variable 2D)  
;
; variable "ter" is terrain in meters!

    ter = wrf_user_getvar(a,"ter",it)  ; terrain

; variables "ua" and "va" are u- and v-winds on sigma surfaces

; currently on level 6 --- p(6) about 912 mb

    ua = wrf_user_getvar(a,"ua",it)  ; u-wind
    va = wrf_user_getvar(a,"va",it)  ; v-wind

    u6 = ua(6,50:100,50:100)
    v6 = va(6,50:100,50:100)

; Compute grad(Z).

    dzdx = new((/71,41/),"float")
    grdx = new((/71,41/),"float")
    do j = 50,120
    do i = 120,160
    dzdx(j-50,i-120) = ter(j,i) - ter(j,i-1)
    dzdx(j-50,i-120) = dzdx(j-50,i-120) / 4000.
    uav = ua(6,j-50,i-120) + ua(6,j-50,i-119)
    grdx(j-50,i-120) = dzdx(j-50,i-120) * uav
;refine this! use avg wind
    end do
    end do

    dzdy = new((/71,41/),"float")
    grdy = new((/71,41/),"float")
    do j = 50,120
    do i = 120,160
    dzdy(j-50,i-120) = ter(j,i) - ter(j,i-1)
    dzdy(j-50,i-120) = dzdy(j-50,i-120) / 4000.
    grdy(j-50,i-120) = dzdy(j-50,i-120) * va(6,j-50,i-120)
;refine this! use avg wind
    end do
    end do

    grad = new((/71,41/),"float")
    grad = grdx + grdy

; How to normalize?
; Try divide by max(dxdx or dzdy)

    maxy = new((/2/),"float")
    maxy(0) = max(dzdx)
    maxy(1) = max(dzdy)
    maxzz = max(maxy)
;print(max(dzdx))
;print(max(dzdy))
;print(maxzz)

;    grad = grad / maxzz

; Try instead divide by max(grad)

    grad = grad / max(grad)

print(max(grad))
print(min(grad))

; variable "pres" is pressure

    pres = wrf_user_getvar(a,"pressure",it)  ; u-wind
    p6 = pres(6,:,:)

; p(level 6 is about 912 hPa). Use for now.
 print(max(p6))


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

      ; Plotting options for gradient
  res = True
  res@cnFillOn  = True
  res@gsnDraw   = True    ; Forces the plot to be drawn
  res@gsnFrame  = True    ; Frame advance
  res@cnMonoLineColor = True
  res@cnLinesOn = True
  res@cnLevelSpacingF  = 0.1    ; spacing
  res@cnLevelSelectionMode = "ExplicitLevels"   ; set explicit contour levels
  res@cnLevels    = (/-1.1,-1.,-.9,-.8,-.7,-.6,-.5,-.4,-.3,-.2,-.1,0,.1,.2,.3,.4,.5,.6,.7,.8,.9,1.,1.1/)
  res@cnMonoLineColor  = True
  res@cnLineColor  = "Black"    ; spacing
  mpres@cnLineColor  = "Black"    ; spacing
  res@cnLineThicknessF = 0.5 
  mpres@cnLineThicknessF = 0.5
  res@FieldTitle = "V . grad(terrain), normalized"    ; Change variable description
  res@lbTitleOn  = False                   ; Turn off title on labelbar

  contoura = wrf_contour(a,wks1, grad(:,:), res)
;  plot  = wrf_map_overlays(a,wks1,contoura,pltres,res)
  plotz = wrf_map_overlays(a,wks1,contoura,pltres,mpres)

      ; Plotting options for terrain in 1st loop
  if (it.eq.0)
  res@cnLinesOn = True
  res@cnLevelSelectionMode = "AutomaticLevels"
  res@cnLevelSpacingF  = 50    ; spacing
  res@cnLineColor  = "Black"    ; spacing
  res@FieldTitle = "TERRAIN"    ; Change variable description
  res@UnitLabel  = "m"                    ; Change variable units
  res@lbTitleOn  = False                   ; Turn off title on labelbar
  contourb = wrf_contour(a,wks2,ter(50:120,120:160),res)
; contourb = wrf_contour(a,wks2,ter(   Y    X ),res)
  plotz = wrf_map_overlays(a,wks2,contourb,pltres,mpres)
  end if



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

  end do        ; END OF TIME LOOP

end