; --------------  LOAD FUNCTIONS AND PROCEDURES ----------------

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

; --------------  BEGINING OF NCL SCRIPT ----------------
begin

;********************************************************
; read in netCDF file and make a loop for all time steps
;********************************************************

path ="/wrf/test/em_real/wrfout_d01_2018-06-01_00:00:00.nc"
fili = systemfunc("cd + path")
nfili = dimsizes(fili)

do nf=0,nfili-1

in = addfile(path + fili(nf), "r")

;--------------------------------------

times = wrf_user_list_times(in)            ; get times in the file
ntimes = dimsizes(times)            ; number of times in the file
printVarSummary(times)

dy = 399
dx = 399

;************************************************
;  - Select lon & lat of the point of interest -
;************************************************

 PH=in->PH(0,0:30,:,:)
 PHB=in->PHB(0,0:30,:,:)
 g=9.81
 z=(PH+PHB)/g

 lvl=7;

 terrain=in->HGT(0,:,:)
 H_above_grnd = z(lvl,:,:)-terrain(:,:)

 print("height at level " + (lvl+1) + " and ocean gridcell is " + H_above_grnd(190,260))

;************************************************
;  - extract wind components and unstagger them -
;************************************************
; ----- Hence U and V are staggered, we need to  average them to get the point on the mass grid using wrf_user_unstagger
 
  U     = wrf_user_getvar(in, "U", -1)
  ua_in = wrf_user_unstagger(U,U@stagger)
  ua    = ua_in(:,lvl,:,:)             ;ua(bottom_up,grid_lon,grid_lat)
	;ua = U(y,x)
;printVarSummary(ua)
  V     = wrf_user_getvar(in, "V", -1)
  va_in = wrf_user_unstagger(V,V@stagger)
  va    = va_in(:,lvl,:,:)            ;va(bottom_up,grid_lon,grid_lat)
	;va = V(y,x)
	

   wind_speed_var=new(dimsizes(ua),typeof(ua))
   wind_speed_var(:,:,:) = sqrt(ua^2+va^2)                          ; wind speed at 10 meters
;   copy_VarCoords(ua,wind_speed_var(it,:,:))                         ; copy coord vars to speed

LAT2D=in->XLAT(0,:,:)
LON2D=in->XLONG(0,:,:)
;ua@lat2d=LAT2D(:,:)
;ua@lon2d=LON2D(:,:)
;va@lat2d=LAT2D(:,:)
;va@lon2d=LON2D(:,:)
wind_speed_var@lat2d=LAT2D(:,:)
wind_speed_var@lon2d=LON2D(:,:)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;Plot the wind speeds
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

wks = gsn_open_wks("png","animate")

res1                 = True
res1@gsnMaximize     = True                ; maximize size
res1@cnFillOn        = True                ; color plot desired

;cmap = read_colormap_file("cmocean_thermal")
;res1@cnFillPalette = cmap(::-1:) ; reverse color map
;res1@cnFillPalette = cmap(45:230,:)   ; subset color map
res1@cnFillPalette  = "WhBlGrYeRe"

res1@cnLinesOn       = False               ; turn off contour lines
res1@cnLineLabelsOn  = False               ; turn off contour labels
res1@gsnAddCyclic    = False               ; turn off longitude cyclic point

res1@cnLevelSelectionMode="ExplicitLevels"
res1@cnLevels = (/0,1,2,3,4,5,6,7,8,10,12,14/)

res1@lbLabelFontHeightF = 0.02
res1@lbOrientation =   "vertical"

res1@mpDataBaseVersion = "MediumRes"
res1@mpOutlineBoundarySets = "USStates"

res1@mpProjection        = "LambertConformal"
res1@mpLambertParallel1F = 33.
res1@mpLambertParallel2F = 45.
res1@mpLambertMeridianF  = -97


res1@mpLimitMode        = "Corners"
res1@mpLeftCornerLatF   =LAT2D(0,0)
res1@mpLeftCornerLonF   =LON2D(0,0)
res1@mpRightCornerLatF  =LAT2D(dy-1,dx-1)
res1@mpRightCornerLonF  =LON2D(dy-1,dx-1)


ntimes = dimsizes(wind_speed_var(:,0,0))
do i=0,ntimes-1

    ;printVarSummary(ntimes)

    res1@tiMainString =  "Wind speed at " + H_above_grnd(190,260)+"(m), june 1 at hr"+ (i) +" (UTC)"
    res1@tiMainFontHeightF  = 0.02
    res1@tiXAxisFontHeightF = 0.016

    res1@lbTitleString    = "m/s"
    res1@lbTitleFontHeightF = 0.015


    plot = gsn_csm_contour_map(wks,wind_speed_var(i,:,:),res1)

    setvalues plot
    end setvalues

    ;Draw the plot and advance the frame.
    draw(plot)
    frame(wks)
    maximize_output(wks,False)

end do

cmd = "convert -delay 12 animate*.png windPlot.gif"
   system(cmd)


end do


; --------------  END OF NCL SCRIPT ----------------

end