;*************************************************
; WRF static: panel different variables 
;************************************************
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/WRF_contributed.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
begin
;************************************************
; open file and read in data
;************************************************



h = addfile("./geo_em.d01_india.nc", "r")



  lat2d = h->XLAT_M(0,:,:)
  lon2d = h->XLONG_M(0,:,:)
  lsMask= h->LANDMASK(0,:,:)           ; land (1) water (0) mas

        dirs = "/usr/rmt_share/tgdata/hrldas/workspace/jam2/india_out/" ; change directory 
        filesn = systemfunc("ls " + dirs + "200511*LDASOUT*.nc")
        filesn1 = systemfunc("ls " + dirs + "2006*LDASOUT*.nc")
       f = addfiles(filesn, "r") ; change to filesn1 if 2006 )
       x       = f[:]->HFX(:,:,:)                              ; (Time, south_north, west_east)
        dirk   = "/usr/rmt_share/tgdata/hrldas/workspace/jam2/india_gemout/"
        filesk = systemfunc("ls " + dirk + "200511*LDASOUT*.nc")
        filesk1 = systemfunc("ls " + dirs + "2006*LDASOUT*.nc")
      
       g = addfiles(filesk, "r")
       y = g[:]->HFX(:,:,:)
       ntim = 700 
       nStrt = 8
       nJump = 3
       nSkip = 24

       do n = 8,ntim,nSkip
          n1 = n
          n2 = n1+nJump
       ;print("n1="+n1+" n2="+n2)
           x1 = x(n1:n2,:,:)
           y1 = y(n1:n2,:,:)
       end do  
       ;  print(x1)
 x1@_FillValue = -1e+33 
 y1@_FillValue = -1e+33 
          xavg = dim_avg_n(x1,0)
          yavg = dim_avg_n(y1,0)
    
 xavg@_FillValue = -1e+33 
 yavg@_FillValue = -1e+33 
printVarSummary(yavg)
printVarSummary(xavg)
print(xavg)
   
      
  lat2d@units = "degrees_north"
  lon2d@units = "degrees_east"  
      xavg@lat2d = lat2d
      xavg@lon2d = lon2d
      yavg@lat2d = lat2d
      yavg@lon2d = lon2d

;  xavg@units = "degrees_north"  
; yavg@units = "degrees_north"

;************************************************                            
; The file should be examined via: ncdump -v grid_type static.wrsi          
; This will print the print type. then enter below.                          
;************************************************                            
  projection = "LambertConformal"   
;************************************************
; create plots  
;************************************************
  wks = gsn_open_wks("ps" ,"sh_novavg_noon3")      ; ps,pdf,x11,ncgm,eps
  gsn_define_colormap(wks ,"BlAqGrYeOrReVi200"); choose colormap

  res                       = True             ; plot mods desired
  res@gsnSpreadColors       = True             ; use full range of colormap
  res@cnFillOn              = True             ; color plot desired
   res@mpMinLonF            = 0               ; set min lon
   res@mpMaxLonF            = 200               ; set max lon
   res@mpMinLatF            =  30.               ; set min lat
  res@vpWidthF  = 0.9              ; height and width of plot
  res@vpHeightF = 0.4
		 
  res@cnLinesOn             = False            ; turn off contour lines
  res@cnLineLabelsOn        = False            ; turn off contour labels
  res@lbLabelBarOn          = False            ; turn off individual lb's 
  res@cnFillMode            = "RasterFill"     ; activate raster mode
 res@cnLevelSelectionMode = "ManualLevels"        ; set manual contour levels
 res@cnMinLevelValF       =  -50                 ; set min contour level
 res@cnMaxLevelValF       =  200           ; set max contour level
 res@cnLevelSpacingF      =  20                  ; set contour spacing
;************************************************
   projection =  "Mercator"
;    projection = "LambertConformal"
   
   
  dimll = dimsizes(lat2d)
   nlat  = dimll(0)
   mlon  = dimll(1)

   
   res@mpProjection          = projection
   res@mpLimitMode           = "Corners"
   res@mpLeftCornerLatF      = lat2d(0,0)
   res@mpLeftCornerLonF      = lon2d(0,0)
  res@mpRightCornerLatF     = lat2d(nlat-1,mlon-1)
   res@mpRightCornerLonF     = lon2d(nlat-1,mlon-1)

   res@mpCenterLonF          = 81.45          ; set center logitude

   if (projection.eq."LambertConformal") then
    res@mpLambertParallel1F   = 30
    res@mpLambertParallel2F   = 60
    res@mpLambertMeridianF    = -90
    end if
				    
  
  
  
  
  res@mpFillOn              = False            ; turn off map fill
  res@mpOutlineDrawOrder    = "PostDraw"       ; draw continental outline last
  res@mpOutlineBoundarySets = "National" ; state boundaries

  res@tfDoNDCOverlay        = False             ; True for 'native' grid
  res@gsnAddCyclic          = False            ; data are not cyclic
;  res@lbOrientation         = "vertical"

;************************************************
; allocate array for 3 plots
;************************************************
  plts                      = new (2,"graphic")   

;************************************************
; Tell NCL not to draw or advance frame for individual plots
;************************************************
  res@gsnDraw               = False            ; (a) do not draw 
  res@gsnFrame              = False            ; (b) do not advance 'frame'

  plts(0)                   = gsn_csm_contour_map(wks,x(8,:,:),res)
  plts(1)                   = gsn_csm_contour_map(wks,y(8,:,:),res)
;************************************************
; create panel: panel plots have their own set of resources
;************************************************
  resP                  = True                   ; modify the panel plot
  resP@txString         = "Avg sensible heat flux (W m-2)  Oct 2005"     
  resP@gsnMaximize      = True                   ; maximize panel area
  resP@gsnPanelLabelBar = True                   ; add common colorbar
  resP@lbLabelAutoStride= True                   ; let NCL figure lb stride 
  resP@gsnPanelRowSpec  = True                   ; specify 1 top, 2 lower level  
  gsn_panel(wks,plts,(/1,2/),resP)               ; now draw as one plot

end