;   Example script to produce plots for a WRF real-data run,
;   with the ARW coordinate dynamics option.
;
;   Terrain ONLY to test
;
; Try to add coastlines
;

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-30_00:00:00_32layers.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_Oct302016")
  wks1 = gsn_open_wks(type,"zgrad_gradient_Oct302016")


; Set some basic resources
  res = True
  res@MainTitle = "REAL-TIME WRF"

  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,3             ; 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

    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((/101,101/),"float")
;    grdx = new((/101,101/),"float")
;    do j = 50,150
;    do i = 50,150
;    dzdx(j-50,i-50) = ter(j,i) - ter(j,i-1)
;    dzdx(j-50,i-50) = dzdx(j-50,i-50) / 4000.
;    grdx(j-50,i-50) = dzdx(j-50,i-50) * ua(6,j-50,i-50)
;;refine this! use avg wind
;    end do
;    end do
;; print(max(ua(6,:,:)))
;; print(max(dzdx))
;; print(max(grdx))

    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) * ua(6,j-50,i-120)
    grdx(j-50,i-120) = dzdx(j-50,i-120) * uav
;refine this! use avg wind
    end do
    end do

;    dzdy = new((/101,101/),"float")
;    grdy = new((/101,101/),"float")
;    do i = 50,150
;    do j = 50,150
;    dzdy(j-50,i-50) = ter(j,i) - ter(j-1,i)
;    dzdy(j-50,i-50) = dzdy(j-50,i-50) / 4000.
;    grdy(j-50,i-50) = dzdy(j-50,i-50) * va(6,j-50,i-50)
;    end do
;    end do
;; print(max(va(6,:,:)))
;; print(max(dzdy))
;; print(max(grdy))
;    grad = new((/101,101/),"float")
;    grad = grdx + grdy
;; print(max(grad))

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

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

    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)
;    v6 = va(6,:,:)

; Try to compute V dot grad(Z).


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

      ; 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@cnLevelSpacingF  = 0.1    ; spacing
  res@cnMonoLineColor  = True
  res@cnLineColor  = "Black"    ; spacing
  res@cnLineThicknesses  = 3    ; spacing

  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
  res@cnLevelSpacingF  = 25    ; spacing
  res@cnLineColor  = "Black"    ; spacing
  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 do        ; END OF TIME LOOP

end