;============================================
; This script reads the output imfs from
; MERRA2 and plots the variable as hovmoller
; ===========================================
;
; These files are loaded by default in NCL V6.2.0 and newer
 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"
; ===========================================
begin
;  dir1  = "/gpfsm/dnb32/mganesha/hht2d/201712.2/output/merra2/20N/v700/"
  dir1  = "/gpfsm/dnb32/mganesha/hht2d/201712.2/output/merra2/20N/v200/"
  f1     = addfile (dir1+"imf.nc","r")
  v_tmp = f1->Decomposed_IMF
  print(dimsizes(v_tmp))
  v     = v_tmp(2,:,:)                                ; get v 
  time_v = f1->time

  print(dimsizes(v))
  print("Min and max meridional wind")
  print(min(v))
  print(max(v))

;  dir2  = "/gpfsm/dnb32/mganesha/hht2d/201712.2/output/merra2/20N/u200/"
  dir2  = "/gpfsm/dnb32/mganesha/hht2d/201712.2/output/merra2/20N/u700/"
  f2     = addfile (dir2+"imf.nc","r")
  u_tmp = f2->Decomposed_IMF
  print(dimsizes(u_tmp))
  u     = u_tmp(0,:,:)                                ; get v
  time_u = f2->time

  print(dimsizes(u))
  print("Min and max zonal wind")
  print(min(u))
  print(max(u))

;Alternatively reading t200
  dir3  = "/gpfsm/dnb32/mganesha/hht2d/201712.2/output/merra2/20N/t200/"
;  dir3  = "/gpfsm/dnb32/mganesha/hht2d/201712.2/output/merra2/20N/t700/"
  f3     = addfile (dir3+"imf.nc","r")
  t     = f3->Decomposed_IMF(4,:,:)                                ; get v
  time_t = f3->time

  t&time = time_t/(60*24)
  u&time = time_u/(60*24)
  v&time = time_v/(60*24)

  print(dimsizes(t))
  print("Min and max temperatures")
  print(min(t))
  print(max(t))
;==============================================
; create color plot
;=============================================
  wks  = gsn_open_wks ("x11", "MERRA2_u700_IMF1_SepStart" )           ; send graphics to PNG file

  res                      = True               ; plot mods desired
  res@cnFillOn             = True               ; turn on color fill
  res@cnFillPalette        = "BlWhRe"           ; set color map

;  res@trYReverse 	   = True

   res@tiYAxisString = "Days"   ; Y-axis title

;Uncomment for v
;  res@tiMainString         = "v700"   ; title
;  res@cnLevelSelectionMode = "ManualLevels"     ; manual contour levels
;  res@cnMinLevelValF       = -10.	        ; min level
;  res@cnMaxLevelValF       =  10.               ; max level
;  res@cnLevelSpacingF      = 2.                 ; contour level spacing

;Uncomment for u
;  res@tiMainString         = "u700"   ; title
;  res@cnLevelSelectionMode = "ManualLevels"     ; manual contour levels
;  res@cnMinLevelValF       = -1.6              ; min level
;  res@cnMaxLevelValF       =  1.6               ; max level
;  res@cnLevelSpacingF      = 0.4                 ; contour level spacing

;Uncomment for t
  res@tiMainString         = "MERRA2 U700 IMF1 @ 20N Sep-Oct"   ; title
;  res@cnLevelSelectionMode = "ManualLevels"     ; manual contour levels
;  res@cnMinLevelValF       = 276.              ; min level
;  res@cnMaxLevelValF       = 286.               ; max level
;  res@cnLevelSpacingF      = 1.                 ; contour level spacing  
  
  res@cnFillColors = (/15,20,25,30,35,40,45,50,50,55,60,65,70,75,80,85,90/)
  res@cnLinesOn = False
  res@cnLevelSelectionMode = "ExplicitLevels"
  res@cnLevels = (/-10,-8,-6,-4,-2,-1,-0.5,0,0.5,1,2,4,6,8,10/)

;***************** CHANGE HERE: Variable and contour levels *********************
;To plot symmetric around 0
; automatically create nice min/max/ci values for blue/red colortable
;  symMinMaxPlt (v,10,False,res)

  plot = gsn_csm_hov(wks, u, res) 
  
end