;============================================
; This script reads the input MERRA2 data for 
; EEMD/HHT 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
  diri  = "/gpfsm/dnb32/mganesha/hht2d/201712.2/input/merra2/"
;  f1     = addfile (diri+"MERRA2.v700_12N_60W_0W_201409_201410.nc","r")
  f1     = addfile (diri+"MERRA2.v200_12N_60W_0W_201409_201410.nc","r")
;  v     = f1->v700                                ; get v 
  v = f1->v200
  time_v = f1->time

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

  f2     = addfile (diri+"MERRA2.u700_12N_60W_0W_201409_201410.nc","r")
  u     = f2->u700                                ; get u
;  f2 = addfile (diri+"MERRA2.u200_12N_60W_0W_201409_201410.nc","r")
;  u = f2->u200
  time_u = f2->time

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

;Alternatively, reading t200
;  f3     = addfile (diri+"MERRA2.t700_12N_60W_0W_201409_201410.nc","r")
;  t     = f3->t700                                ; get v
  f3     = addfile (diri+"MERRA2.t200_12N_60W_0W_201409_201410.nc","r")
  t     = f3->t200 
  time_t = f3->time

  print(dimsizes(t))
  print("Min and max temperatures")
  print(min(t))
  print(max(t))

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

  v_anom = new(dimsizes(v),typeof(v),v@_FillValue)
  u_anom = new(dimsizes(u),typeof(u),u@_FillValue)
  t_anom = new(dimsizes(t),typeof(t),t@_FillValue)

  do i = 0,64
  	v_anom(:,i) = v(:,i) - dim_avg_n(v(:,i),0)
        u_anom(:,i) = u(:,i) - dim_avg_n(u(:,i),0)
 	t_anom(:,i) = t(:,i) - dim_avg_n(t(:,i),0)
  end do

  print(dimsizes(v_anom))
  print("Min and max meridional wind anomaly")
  print(min(v_anom))
  print(max(v_anom))

  print(dimsizes(u_anom))
  print("Min and max zonal wind anomaly")
  print(min(u_anom))
  print(max(u_anom))

  print(dimsizes(t_anom))
  print("Min and max temperature anomaly")
  print(min(t_anom))
  print(max(t_anom))
;==============================================
; create color plot
;=============================================
  wks  = gsn_open_wks ("x11", "MERRA2_U700_input_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@cnFillPalette        = "BlRe"           ; 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       = -24.              ; min level
;  res@cnMaxLevelValF       =  16.               ; max level
;  res@cnLevelSpacingF      = 4.                 ; contour level spacing

;Uncomment for t
  res@tiMainString         = "MERRA2 U700 @ 12N Sep-Oct"   ; title
;  res@cnLevelSelectionMode = "ManualLevels"     ; manual contour levels
;  res@cnMinLevelValF       = 278.              ; min level
;  res@cnMaxLevelValF       = 287.               ; max level
;  res@cnLevelSpacingF      = 0.5                 ; contour level spacing  

;  cmap = read_colormap_file("BlWhRe") 
;  do i = 46,56
;  	cmap(46:56,0) = 1
;	cmap(46:56,1) = 0.9215686
;	cmap(46:56,2) = 0.9215686
;  end do

  res@cnFillColors = (/5,15,25,30,35,40,45,50,50,55,60,65,70,75,80,90,100/)  
  res@cnLinesOn = False
  res@cnLevelSelectionMode = "ExplicitLevels"
;  res@cnLevels = (/-30,-25,-20,-15,-10,-5,-1,0,1,5,10,15,20,25,30/)
  res@cnLevels = (/-24,-20,-16,-12,-8,-4,-1,0,1,4,8,12,16,20,24/)

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

  plot = gsn_csm_hov(wks, u, res) 
  
end