begin
;************************************************
; Read the file
;************************************************
  ymStrt = 198601                              ; climatology start year/month
  ymLast = 201712                              ;      "      last  year/month  a      = addfile("air.mon.mean.nc","r")
  a      = addfile("air.mon.mean.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
  tmp    = a->air(ntStrt:ntLast,:,:)            ; (time,lat,lon)
  tmp= tmp-273.15
  ntim = dimsizes(time)
 printMinMax(tmp,0)
  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(tmp)     
  ;printVarSummary(prcStd)
  anom = calcMonAnomTLL(tmp,prcClm)
;************************************************
; create colors
;*************************************************
  wks = gsn_open_wks("png","climo_ghcn")             ; send graphics to PNG file

;************************************************
; create panel plots
;*************************************************
  plot   = new (12, graphic)                    ; create graphical array

  res                      = True               ; plot options desired
  res@gsnDraw              = False
  res@gsnFrame             = False
  res@cnFillOn             = True               ; turn on color fill
  res@cnFillPalette        = "amwg256"          ; 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@cnLevelSelectionMode = "ManualLevels"     ; set manual contour levels
  res@cnMinLevelValF       = -10        ; set min contour level
  res@cnMaxLevelValF       = 6  ; set max contour level
  res@cnLevelSpacingF      = 1            ; set contour spacing

  res@mpFillOn             = False              ; turn off gray continents
  res@mpCenterLonF         = 0                  ; Centers the plot at 0
  res@lbLabelBarOn         = False              ; No single label bar


 res@mpMaxLatF            = 42        ; choose subregion
 res@mpMinLatF            = 0
 res@mpMaxLonF            = 67
 res@mpMinLonF            = -60
 res@mpDataBaseVersion    = "MediumRes"      ; better map outlines
 res@gsnRightString      = " mm"




  resP                     = True               ; panel options
  resP@gsnPanelMainString  = " Monthly anomaly temperature : "+(ymStrt/100)+"-"+(ymLast/100)
  resP@gsnMaximize         = True               ; maximize image
  resP@gsnPanelLabelBar    = True               ; Add common label bar


  do nmo=0,2                                  ; loop over the months
     res@gsnLeftString     = months(nmo)
     plot(nmo) = gsn_csm_contour_map(wks,anom(nmo,:,:), res)  ; create plot
  end do


    do nmo=9,11                                  ; loop over the months
     res@gsnLeftString     = months(nmo)
     plott(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,(/2,3/),resP)
end