;=======================================================================
; name:
;   NHpolar_tas_sea.ncl
;
; category:
;   NCL script
;
; purpose:
;   This scripts plots seasonal scalar data on the polar stereographic
;     projection.
;
; author:
;   m.yoshimori
;=======================================================================

begin

; specify file names
  pathin = "../monthly/"
  fin = "t2m_y1959-y2021.nc"  ; input file name (netCDF)
  fout = "NHpolar_tas_sea"    ; output file name (png/eps)

; open file
  in = addfile(pathin+fin,"r")

; read data
  tmp = in->var167

; check the reading
  printVarSummary(tmp)

; change unit from (K) to (degC)
  tmp = tmp - 273.15

; extract data for 1979-1993 & 2007-2021
  m0 = (1979-1959)*12
  m1 = m0 + 12*15 - 1
  m2 = (2007-1959)*12
  m3 = m2 + 12*15 - 1
  tmpmon0 = clmMonTLL(tmp(m0:m1,:,:))  ; 1979-1993
  tmpmon1 = clmMonTLL(tmp(m2:m3,:,:))  ; 2007-2021

; difference
  tmpmon0 = tmpmon1 - tmpmon0

; seasonal averages (DJF, MAM, JJA, SON)
  tmp0djf = tmpmon0(0,:,:)
  tmp0mam = tmp0djf
  tmp0jja = tmp0djf
  tmp0son = tmp0djf
  tmp0djf = (tmpmon0(11,:,:)*31.+tmpmon0(0,:,:)*31.+tmpmon0(1,:,:)*28.)/90.
  tmp0mam = (tmpmon0(2,:,:)*31.+tmpmon0(3,:,:)*30.+tmpmon0(4,:,:)*31.)/92.
  tmp0jja = (tmpmon0(5,:,:)*30.+tmpmon0(6,:,:)*31.+tmpmon0(7,:,:)*31.)/92.
  tmp0son = (tmpmon0(8,:,:)*30.+tmpmon0(9,:,:)*31.+tmpmon0(10,:,:)*30.)/91.

;-----
; plot
;-----

  wks = gsn_open_wks("png",fout)  ; open a file

; color table
;  gsn_define_colormap(wks,"temp_19lev")  ; choose a color map
  gsn_define_colormap(wks,"temp_diff_18lev")  ; choose a color map

; create a plot array
  plot = new(4,graphic)

;--------------
; plot map only
;--------------

  mpres = True
  mpres@gsnDraw = False
  mpres@gsnFrame = False
  mpres@mpFillOn = False
  mpres@gsnPolar = "NH"
  mpres@mpMinLatF = 45.
  mpres@mpMaxLatF = 90.
  mpres@tiMainString = ""
  mpres@gsnLeftString = ""
  mpres@gsnCenterString = ""
  mpres@gsnRightString = ""

; fontsize of the main title
  mpres@tiMainFontHeightF = 0.04
  mpres@tiMainOffsetYF = -0.0002

; fontsize of the subtitles
  mpres@txFontHeightF = 0.028

; fontsize of tickmark labels
  mpres@tmXBLabelFontHeightF = 0.018
  mpres@tmYLLabelFontHeightF = 0.018

;---------------
; draw each plot
;---------------

; basic settings
  res = True                  ; plot mods desire
  res@gsnDraw = False         ; draw plot
  res@gsnFrame = False        ; advance frome
  res@cnFillOn = True         ; turn on color fill
  res@cnLinesOn = False       ; turn off contour lines
  res@cnInfoLabelOn = False   ; turn off info label
  res@cnLineLabelsOn = False  ; turn off contour labels

; color bar
  res@lbLabelBarOn = False  ; turn off contour labels

; contour levels
;  res@cnLevelSelectionMode = "ManualLevels"  ; set manual contour levels
;  res@cnMinLevelValF = -5                    ; set min contour level
;  res@cnMaxLevelValF =  5                    ; set max contour level
;  res@cnLevelSpacingF = 1                    ; set contour spacing
  res@cnLevelSelectionMode = "ExplicitLevels"  ; set explicit contour levels
  res@cnLevels             = (/-5, -4, -3, -2, -1, 1, 2, 3, 4, 5/)

; title
  res@tiMainString = ""
  res@gsnLeftString = ""
  res@gsnCenterString = ""
  res@gsnRightString = ""

; change these depending on color or B&W plot
  res@gsnSpreadColors      = True       ; spread out color table
;  res@gsnContourZeroLineThicknessF = 0

  mpres@tiMainString = "(a) DJF"
  plot(0) = gsn_csm_map_polar(wks,mpres)
  map0 = gsn_csm_contour(wks,gsn_add_cyclic_point(tmp0djf),res)
  overlay(plot(0),map0)

  mpres@tiMainString = "(b) MAM"
  plot(1) = gsn_csm_map_polar(wks,mpres)
  map1 = gsn_csm_contour(wks,gsn_add_cyclic_point(tmp0mam),res)
  overlay(plot(1),map1)

  mpres@tiMainString = "(c) JJA"
  plot(2) = gsn_csm_map_polar(wks,mpres)
  map2 = gsn_csm_contour(wks,gsn_add_cyclic_point(tmp0jja),res)
  overlay(plot(2),map2)

  mpres@tiMainString = "(d) SON"
  plot(3) = gsn_csm_map_polar(wks,mpres)
  map3 = gsn_csm_contour(wks,gsn_add_cyclic_point(tmp0son),res)
  overlay(plot(3),map3)

;---------------
; create a panel
;---------------

  resP = True
  resP@gsnFrame         = False                       ; don't advance panel plot
  resP@gsnPanelLabelBar = True                        ; add common colorbar
  resP@lbBoxEndCapStyle = "TriangleBothEnds"          ; triangle edges
  resP@pmLabelBarWidthF = 0.8                         ; label bar width
  resP@pmLabelBarOrthogonalPosF = -0.02               ; shift up-down
  resP@lbLabelStride = 1                              ; label interval

; title
  resP@txFontHeightF = 0.02
  resP@txString = "Surface air temperature change"    ; title

; color bar title
  resP@lbTitleOn          = True
  resP@lbTitleString      = "(K)"
  resP@lbTitleFontHeightF = 0.016
  resP@lbTitlePosition    = "Right"
  resP@lbTitleDirection   = "Across"

  gsn_panel(wks,plot,(/1,4/),resP)  ; now draw as one plot

  frame(wks)

end