begin
;************************************************
; Read the file
;************************************************
  ymStrt = 198601                              ; climatology start year/month
  ymLast = 202112                              ;      "      last  year/month

  a    = addfile("soilw.mon.mean.v2.nc","r")
  time = a->time
  yyyymm = cd_calendar(time, -1)
  ntStrt = ind(yyyymm.eq.ymStrt)                ; start time index
  ntLast = ind(yyyymm.eq.ymLast)                ; last  time index
  prc    = a->soilw(ntStrt:ntLast,:,:)            ; (time,lat,lon)
  ntim = dimsizes(time)

  months = (/"January", "February", "March", "April" \
            ,"May", "June", "July", "August"         \
            ,"September", "October", "November"      \
            ,"December" /)
;************************************************
; Compute the climatology using a function in contributed.ncl
;************************************************
  ;prcStd = stdMonTLL(tmp)            ; monthly interannual variability
  prcClm   = clmMonTLL(prc)
  ;printVarSummary(prcStd)
  anom = calcMonAnomTLL(prc,prcClm)
;************************************************
; create colors
;*************************************************


 wks = gsn_open_wks("png","soil-moist")             ; send graphics to PNG file

  ;colors = (/ "azure1","beige","lavender" \
            ;,"PaleGreen","SeaGreen3","LightYellow" \
            ;,"Yellow","HotPink","Red"/)         ; choose colors for color map
;************************************************
; create panel plots
;*************************************************
  plot   = new (12 , graphic)                    ; create graphical array

  res                      = True               ; plot options desired

  res@cnFillOn             = True               ; turn on color fill
  res@cnFillPalette        = "perc2_9lev"              ; set color map
  res@cnInfoLabelOn        = False              ; turn off contour info label
  res@cnLinesOn            = False              ; turn off contour lines
  res@cnLineLabelsOn       = False              ; turn off line labels
; res@cnFillMode = "RasterFill"
 res@cnLevelSelectionMode = "ManualLevels"     ; set manual contour levels
  res@cnMinLevelValF       = -15        ; set min contour level
  res@cnMaxLevelValF       = 15 ; set max contour level
  res@cnLevelSpacingF      = 1           ; set contour spacing






  res@mpMaxLatF            = 70        ; choose subregion
  res@mpMinLatF            = -10
  res@mpMaxLonF            = 70
  res@mpMinLonF            = -60

;  res@mpDataBaseVersion    = "MediumRes"      ; better map outlines
  res@mpFillOn             = False              ; turn off gray continents
  res@mpCenterLonF         = 0.5*(res@mpMinLonF + res@mpMaxLonF)         ; Centers the plot at 180
   res@lbLabelBarOn         = False            ; No single label bar
  res@gsnDraw              = False
  res@gsnFrame             = False
  ;res@lbBoxLineThicknessF   = 3.0
  ;res@gsnPaperOrientation = "auto"

;-----------------------------------------------------------------------------------
  resP                     = True               ; panel options
  resP@gsnPanelMainString  = "Soil moisture anomaly"   ; common title
  resP@gsnMaximize         = True               ; maximize image
  resP@gsnPanelLabelBar    = True              ; Add common label bar
  ;resP@lbOrientation        = "vertical"  ; vertical label bar
  ;resP@pmLabelBarHeightF     = 0.1               ; default is taller
  resP@lbLabelFontHeightF= 0.01
;-----------------------------------------------------------------------------------

printVarSummary(anom)
printMinMax(anom(0,:,:),0)
printMinMax(anom(1,:,:),0)
replace_ieeenan (anom, anom@_FillValue, 0)
printMinMax(anom(0,:,:),0)
  do nmo=0,2
      res@tmXBLabelFontHeightF = 0.02
      res@tmYLLabelFontHeightF = 0.02
     ;res@tmYRLabelFontQuality = High
     ;res@tmYRLabelFontThicknessF= 120
     ;res@tmYMajorGridThicknessF = 4
     res@tmYLMajorThicknessF = 3.5
     res@tmBorderThicknessF=3

;     res@gsnMaximize         = True
     res@mpShapeMode  = "FreeAspect"
     res@vpWidthF      = .7
     res@vpHeightF     = .35

     res@gsnRightString      = "(mm)"                                 ; loop over the months
     res@gsnStringFontHeightF = 0.028         ; make larger than default
     res@gsnLeftString     = months(nmo)
     ;res@gsnCenterString   = months(nmo)+":"+"1986"+"-"+"2021"
     plot(nmo) = gsn_csm_contour_map(wks,anom(nmo,:,:), res)  ; create plot
  end do
;------------------------------------------------------------------------------------------------
do nmo=9,11
        res@tmXBLabelFontHeightF = 0.02
        res@tmYLLabelFontHeightF = 0.02
        res@tmYLMajorThicknessF = 3.5
        res@tmBorderThicknessF=3
       ;res@tmXBLabelsOn         = False              ; do not draw bottom labels
       ;res@tmXBOn               = False              ; no bottom tickmarks

;        res@gsnMaximize         = True                                 ; loop over the months
        res@mpShapeMode  = "FreeAspect"
        res@vpWidthF      = .7
        res@vpHeightF     = .35


        res@gsnRightString      = "(mm)"
        res@gsnStringFontHeightF = 0.028         ; make larger than default
        res@gsnLeftString     = months(nmo)
       ;res@gsnCenterString   = months(nmo)+":"+time(0)/100 +"-"+ time(ntim-1)/100
 plot(nmo) = gsn_csm_contour_map(wks,anom(nmo,:,:), res)  ; create plot
  end do
;--------------------------------------------------------------------------------------------------



  ;shapefile_name1 = "clip.basin3.shp"
  shapefile_name2 = "IRN_adm0.shp"
  lnres = True
  lnres@gsLineThicknessF = 5
  lnres@gsLineColor = "Black"

  ;shape1 = gsn_add_shapefile_polylines(wks,plot,shapefile_name1,lnres)
;  shape2 = gsn_add_shapefile_polylines(wks,plot,shapefile_name2,lnres)
;---------------------------------------------------------------------------------------------------


 gsn_panel(wks,(/plot(0),plot(1),plot(2),plot(9),plot(10),plot(11)/),(/2,3/),resP)


end