;Scripts using the external ".shp" files.
;Target plot is the two separated plots.

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

;=============================================================
;Regridding section
;Open source file
;=============================================================
  setfileoption("h5","FileStructure","Advanced")
  
  files = systemfunc("ls /home/rshydro/heqing/AMSR2_L3/GW1AM2_201*_01D_EQMD_L3SGSMCLF2210210.h5")
  printVarSummary(files)
  ;print(files)
  
  nfiles = dimsizes(files)
  
  sm_t =new((/720,1440,368/),float,-9999.)
 
  lat  = fspan(90.00, -90.00, 720)
  lon  = fspan(0.00,360.00,1440);must be matched to the res@mpMinLonF=0. res@mpMaxLonF=360.
  ;lon = fspan(-180.00,180.00,1440) ;why res@mpMinLonF=-180. res@mpMaxLonF=180. invalid?
  time = ispan(0,367,1) 

  do nf  = 0,nfiles-1,1
    setfileoption("h5","FileStructure","Advanced")
    f    := addfile(files(nf),"r") 
    sm   := f->Geophysical_Data
    sm1  := sm*0.1

    sm1!0 = "xlat"
    sm1!1 = "xlon"
    sm1&xlat  = lat
    sm1&xlon  = lon
    sm1&xlat@units  = "degrees_north"
    sm1&xlon@units  = "degrees_east"
    sm1@_FillValue  = -9999.
    sm_t(:,:,nf) = sm(:,:,0)
    
    print("working on "+files(nf))

  end do 

  printVarSummary(sm_t)

  sm_t!0 = "xlat"
  sm_t!1 = "xlon"
  sm_t!2 = "xtime"
  sm_t&xlat  = lat
  sm_t&xlon  = lon
  sm_t@xtime = time
  sm_t&xlat@units  = "degrees_north"
  sm_t&xlon@units  = "degrees_east"
  sm_t&xtime@units = "day"
  sm_t@_FillValue  = -9999. 
 
  sm_t = where(sm_t.gt.0,sm_t,sm_t@_FillValue)
  printVarSummary(sm_t)

  sm_avg = dim_avg_n_Wrap(sm_t,2)
  printVarSummary(sm_avg)
   

  wks = gsn_open_wks("pdf","./plot/amsr2_sm_annual")

  gsn_define_colormap(wks,"GMT_seis")

  
;=============================================================
;Set some common plotting resources
;=============================================================
  res                            = True

  res@gsnDraw                    = False
  res@gsnFrame                   = False 
   
  res@gsnAddCyclic               = False

  res@lbLabelBarOn               = True
  res@tiMainString               = "Annual mean AMSR2 soil moisture"
  res@gsnRightString             = "%"

  res@mpLimitMode                = "LatLon"
  res@mpLandFillColor            = "transparent"
  ;res@mpMinLonF                  =    0.
  ;res@mpMaxLonF                  =  360.
  ;res@mpCenterLonF               =    0.
  
  res@pmTickMarkDisplayMode      = "Always"
  res@tmXTOn                     = False  ;Turn off the top tick marks(coordinate X)
  res@tmYROn                     = False  ;Turn off the right tick marks(coordinate Y)
  res@tmYLMode                   = "Manual"  ;Setting the range of the left tick marks in "Manual" mode ,default is "Automatic".
  
  res@gsnSpreadColors            = True
  res@cnLinesOn                  = False
  
   res@cnFillOn                   = True
  res@cnLevelSelectionMode       = "ManualLevels"
  res@cnMinLevelValF             =  0.
  res@cnMaxLevelValF             =  100.
  res@cnLevelSpacingF            =  5.
  res@cnInfoLabelOn              = False
  res@cnConstFLabelOn            = False

  res@mpLandFillColor            = "transparent"
  res@mpGridAndLimbOn            = True
  res@mpGridLineColor            = "gray"

  plot = gsn_csm_contour_map_ce(wks,sm1(:,:,0),res) 

  draw(plot)
  frame(wks)

end