;=======================================================================
; name
;   trend_signif.ncl
;
; category
;   NCL script
;
; description
;   This script creates a linear regression map (and netCDF file)
;   with a statistical significance.
;
; usage
;   ncl trend_signif.ncl
;
; output
;   rcoef: regression coefficient (per year)
;   pval: p-value (<0.05 -> significant at 5% level of falsely rejecting
;         the null hypothesis, i.e., rcoef=0)
;
; refrence
;   NCL web page (rcoef)
;
; author
;   m.yoshimori
;=======================================================================

begin

; specify file names (input&output, netCDF)
  pathi = "./inp/"                   ; input directory
  fin = "ts_t2m_y1979-y2015_djf.nc"  ; input file name
  patho = "./out/"                   ; output directory
  fout = "trend_eraint_t2m_DJF"      ; output file name

; open input files
  in = addfile(pathi+fin,"r")

; read data
  tmp = in->t2m(latitude|:,longitude|:,time|:)
  printVarSummary(tmp)

; regression analysis
  x = ispan(0,dimsizes(tmp&time)-1,1)*1.
  rc           = regCoef(x,tmp)             
  rc@long_name = "regression coefficient"
  rc@units     = "unitless"
  copyatt(rc,tmp)  ; copy other atts and cv's
  tval = onedtond(rc@tval , dimsizes(rc))
  df   = onedtond(rc@nptxy, dimsizes(rc)) - 2
  b    = tval
  b    = 0.5        ; b must be same size as tval (and df)
  prob = betainc(df/(df+tval^2),df/2.0,b)

; output arrays
  rcoef = tmp(:,:,0)
  rcoef = rc
  pval = rcoef
  pval = (/prob/)  ; copy data
  rcoef@long_name = "regression coefficient"
  pval@long_name = "probability"

; write output netCDF file
  system("rm " + patho+fout+".nc")     ; remove any pre-existing file
  out = addfile(patho+fout+".nc","c")  ; open netCDF file
  out->rcoef = rcoef
  out->pval = pval

;***********************************************************************
; remove the following part unless the graphical output is needed.
;***********************************************************************

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

; change unit for display purpose
  rcoef = rcoef*10.  ; change (K/yr) to (K/10yr)

; judge significance
  siglvl = 0.05
  sig = pval
  ndim = dimsizes(pval)
  nlat = ndim(0)
  nlon = ndim(1)
  do j=0,nlat-1
    do i=0,nlon-1
      if (pval(j,i).lt.siglvl) then
        sig(j,i) = 1.
      else
        sig(j,i) = 0.
      end if
    end do
  end do

; create a plot
  wks = gsn_open_wks("png",patho+fout)	; open an eps file
  gsn_define_colormap(wks,"temp_diff_18lev")	; choose a color map

; basic settings
  res = True		; plot mods desired
  res@gsnDraw = False	; draw plot
  res@gsnFrame = False	; advance frame
  res@cnFillOn = True	; turn on color fill
  res@cnLinesOn = False	; turn off contour lines
  res@cnLineLabelsOn = False    ; turn off contour labels

; contour levels
;  res@cnLevelSelectionMode = "ManualLevels"
;  res@cnMinLevelValF = -1.6
;  res@cnMaxLevelValF =  1.6
;  res@cnLevelSpacingF = 0.2
  res@cnLevelSelectionMode = "ExplicitLevels"  ; set explicit contour levels
  res@cnLevels = (/-1.8,-1.4,-1.0,-0.6,-0.2,0.2,0.6,1.0,1.4,1.8/)

; reshape color bars
  res@lbOrientation = "vertical"        ; vertical label bars
  res@pmLabelBarOrthogonalPosF = -0.02  ; shift left-right
  res@lbBoxEndCapStyle = "TriangleBothEnds"  ; triangle edges
  res@lbLabelStride = 1                 ; label interval
  res@lbBoxMinorExtentF = 0.2           ; slim width

; adjust aspect ratio
  res@mpShapeMode = "FreeAspect"
  res@vpWidthF = 0.75
  res@vpHeightF = 0.45

; map
  res@mpFillOn = False         ; turn off gray continents
  res@mpCenterLonF = 180       ; set centers the plot at 180E
  res@gsnMajorLonSpacing = 60  ; tick mark interval for longitude

; colar table
  res@gsnSpreadColors = True   ; spread out color table

; font size
  res@tiMainFontHeightF = 0.022     ; fontsize of the main title
  res@txFontHeightF = 0.02          ; fontsize of the subtitles
  res@tmXBLabelFontHeightF = 0.018  ; fontsize of tickmark labels (x-axis)
  res@tmYLLabelFontHeightF = 0.018  ; fontsize of tickmark labels (x-axis)

; title
  res@tiMainString = "Linear trend (1979-2015)"  ; title
  res@gsnCenterString = ""	                 ; subtitle
  res@gsnLeftString = "DJF"                      ; upper-left subtitle
  res@gsnRightString = ""                        ; upper-right subtitle

; color bar title
  res@lbTitleOn          = True
  res@lbTitleString      = "(K/10yr)"
  res@lbTitleFontHeightF = 0.018
  res@lbTitlePosition    = "Bottom"

; plot
  plot = gsn_csm_contour_map_ce(wks,rcoef,res)

;------------------------------
; plot statistical significance
;------------------------------

  sgres                      = True		; significance
  sgres@gsnDraw              = False		; draw plot
  sgres@gsnFrame             = False		; advance frome
  sgres@cnInfoLabelOn        = False		; turn off info label
  sgres@cnLinesOn            = False		; draw contour lines
  sgres@cnLineLabelsOn       = False		; draw contour labels
;  sgres@cnFillScaleF         = 0.6		; add extra density
  sgres@cnFillDotSizeF       = 0.003

; activate if gray shading for B&W plot  
  sgres@cnFillOn = True
  sgres@cnFillColors = (/"transparent","transparent"/)	; choose one color for our single cn level
  sgres@cnLevelSelectionMode = "ExplicitLevels"	; set explicit contour levels
  sgres@cnLevels = 0.5	; only set one level
  sgres@lbLabelBarOn = False

  sgres@tiMainString = ""     ; title
  sgres@gsnCenterString = ""  ; subtitle
  sgres@gsnLeftString = ""    ; upper-left subtitle
  sgres@gsnRightString = ""   ; upper-right subtitle

  sig_plot = gsn_csm_contour(wks,sig,sgres)

  opt = True
  opt@gsnShadeFillType = "pattern"
  opt@gsnShadeHigh     = 17
  sig_plot = gsn_contour_shade(sig_plot,-999.,0.5,opt)

  overlay(plot,sig_plot)

  draw(plot)
  frame(wks)

;***********************************************************************
; remove the above part unless the graphical output is needed.
;***********************************************************************

end