;=======================================================================
; name
;   diff_signif.ncl
;
; category
;   NCL script
;
; description
;   This script computes statistical significance (p-value) using the
;   two-tailed t-test.
;
; usage
;   ncl diff_signif.ncl
;
; input
;   two time series files to be compared (e.g., control vs. perturbed)
;
; output
;   diff: difference of the two data (e.g., control vs. perturbed)
;   pval: p-value (<0.05 -> significant at 5% level of falsely rejecting
;         the null hypothesis, i.e., control=perturbed;
;         <0.1 for single-tailed t-test)
;
; refrence
;   NCL web page (ttest)
;
; author
;   m.yoshimori
;=======================================================================

begin

; specify file names (input&output, netCDF)
  pathi = "./inp/"                    ; input directory
  fin1 = "ts_t2m_y1979-y1993_djf.nc"  ; input file name #1
  fin2 = "ts_t2m_y2001-y2015_djf.nc"  ; input file name #2
  patho = "./out/"                    ; output directory
  fout = "diff_eraint_t2m_DJF"      ; output file name

; open input files
  in1 = addfile(pathi+fin1,"r")
  in2 = addfile(pathi+fin2,"r")

; read data
  tmp1 = in1->t2m
  tmp2 = in2->t2m
  printVarSummary(tmp1)
  printVarSummary(tmp2)

; t-test
  xtmp = tmp1(latitude|:,longitude|:,time|:)
  ytmp = tmp2(latitude|:,longitude|:,time|:)
  Xave = dim_avg_Wrap(xtmp)
  Yave = dim_avg_Wrap(ytmp)
  Xvar = dim_variance(xtmp)
  Yvar = dim_variance(ytmp)
  Xs   = dimsizes(xtmp(0,0,:))
  Ys   = dimsizes(ytmp(0,0,:))
  print(Xs)
  print(Ys)
  iflag    = False  ; population variance similar
  tval_opt = False  ; p-value only
  prob = ttest(Xave,Xvar,Xs,Yave,Yvar,Ys,iflag,tval_opt)

; calc. the difference
  diff = Xave
  diff = Yave - Xave
  pval = Xave
  pval = (/prob/)
  diff@long_name = "difference of the means"
  pval@long_name = "probability"

; write output file
  system("rm " + patho+fout+".nc")     ; remove any pre-existing file
  out = addfile(patho+fout+".nc","c")  ; open netCDF file
  out->diff = diff
  out->pval = pval
 
;***********************************************************************
; remove the following part unless the graphical output is needed.
;***********************************************************************

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

; 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 = -4.
;  res@cnMaxLevelValF =  4.
;  res@cnLevelSpacingF = 0.4
  res@cnLevelSelectionMode = "ExplicitLevels"  ; set explicit contour levels
  res@cnLevels = (/-4.4,-3.6,-2.8,-2.0,-1.2,-0.4,0.4,1.2,2.0,2.8,3.6,4.4/)

; 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 = "Difference (2001-2015) - (1979-1993)"      ; title
  res@gsnCenterString = ""	                    ; subtitle
  res@gsnLeftString = "DJF"                         ; upper-left subtitle
  res@gsnRightString = ""  ; upper-right subtitle

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

; plot
  plot = gsn_csm_contour_map_ce(wks,diff,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