;script to plot the u,v,geopotential height plot of MERRA-2 data at location BBY.
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
   dir = "/home/MERRA_DATA_ARV/"
   fname = "MERRA2_400.inst3_3d_asm_Np.20160728.nc"
   a = addfile(dir + fname,"r")

   type = "ps"
   wks = gsn_open_wks(type, "plt_uvz_epv_MERRA2_asm_Jul28_loop") 

   z  = a->H
   z  = z/1000.0 ; geopotential height , convert unit to km. 
   u  = a->U
   v  = a->V
   epv = a->EPV
   epv = epv*1e6

   do kk= 250,500,250  ; level to be plotted
   do it= 6,7  ;MERRA2 data, every 3hr.
   z@long_name = "geopotential height @" +kk+ "hPa"
   z@units     = "km"

   z_kk = z(time|it,{lev|kk},lat|:,lon|:)
   u_kk = u(time|it,{lev|kk},lat|:,lon|:)
   v_kk = v(time|it,{lev|kk},lat|:,lon|:)
   epv_kk = epv(time|it,{lev|kk},lat|:,lon|:)

   printVarSummary(v_kk)

  max_z = floattointeger(max(z_kk))
  min_z = floattointeger(min(z_kk))

   printMinMax(z_kk,0)

   res = True
   res@gsnDraw        = False
   res@gsnFrame       = False
;   res@gsnMaximize    = True   ; maximize size of plot
   res@mpGeophysicalLineThicknessF = 2
   res@mpNationalLineThicknessF    = 2

   res@cnFillOn       = True
   res@cnLinesOn      = False
   res@cnLineLabelsOn = False
   res@gsnAddCyclic   = True    ; add longitude cyclic point, mostly 
                                ; only useful if plotting global data
   res@mpMinLatF      = 20 ;30
   res@mpMaxLatF      = 60 ;45
   res@mpMinLonF      = -150 ;-130
   res@mpMaxLonF      = -90  ;-110

   res@tiMainString   = "Synoptic Condition, MERRA2 on July 28, 2016" 
   res@pmTitleZone    = 4   ; move title down
   res@pmTickMarkDisplayMode = "Always"    ; nicer tickmarks

   res@cnLevelSelectionMode = "ManualLevels"        ; manual contour levels
   res@cnLevelSpacingF      = 0.05                   ; contour interval
   res@cnMinLevelValF       = min_z                ; min level
   res@cnMaxLevelValF       = max_z                 ; max level
   res@cnLineLabelsOn       = True                  ; turn on line labels


   plot = gsn_csm_contour_map(wks,z_kk,res) ; use gsn_csm_contour_map to examine the data variable for lat/lon coordinate arrays,helpful for plotting over a map.

  ;;set vector;;
        res_vc                = True
        res_vc@gsnDraw        = False
        res_vc@gsnFrame       = False
        res_vc@vcRefLengthF            = 0.045            ; define length of vec ref
        res_vc@vcGlyphStyle            = "CurlyVector"    ; turn on curly vectors
        res_vc@vcMinDistanceF          = 0.017

        res_vc@gsnLeftString        = "Wind"
        res_vc@vcRefAnnoOn               = True
        res_vc@vcRefMagnitudeF           = 30 
        res_vc@vcRefAnnoString1          = "30"

         vector  = gsn_csm_vector(wks,u_kk,v_kk,res_vc) 

  ;;set epv ;;
        res_epv    = True
        res_epv@gsnDraw        = False
        res_epv@gsnFrame       = False
        
        res_epv@cnLevelSelectionMode = "ExplicitLevels"  ; use explicit levels
        res_epv@cnLevels             = ispan(1,2,1 )   ; set the contour levels overlay(plot,vector)
        res_epv@cnFillOn              = False
        res_epv@cnLinesOn             = True       
        res_epv@cnLineLabelsOn       = False
        res_epv@cnMonoLineColor = False 

;        cmap       = read_colormap_file("BlueDarkOrange18")
;        cmap(0,3)         = 0.0    ; first color fully transparent
        res_epv@cnLineColors  = (/"black","blue"/)
 
        res_epv@gsnLeftString        = ""
        res_epv@gsnRightString        = ""
        plot_epv = gsn_csm_contour(wks,epv_kk,res_epv) 
      
   overlay(plot, vector)
   overlay(plot, plot_epv)
   draw(plot)
   frame(wks)   

  end do ; time loop
  end do ; lev loop

end