;******************************************************************
; Creator: Sho Kawazoe
; 3019 Agronomy Hall
; Department of Geological and Atmospheric Sciences
; Iowa State University, Ames IA
;
; Computes the SPI (Standardozed Precipitaion Index) for CNH
; analysis to quantify wet and dry season with association to
; summer results 
;
; Reference: NCL webpage: http://www.ncl.ucar.edu/
;*******************************************************************

;*******************************************************************
;Load NCL library files. Always load "Shea_util" AFTER contributed
;*******************************************************************

  load "/usr/local/ncl/lib/ncarg/nclscripts/csm/gsn_code.ncl"
  load "/usr/local/ncl/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
  load "/usr/local/ncl/lib/ncarg/nclscripts/csm/contributed.ncl"
  load "/usr/local/ncl/lib/ncarg/nclscripts/csm/shea_util.ncl"

begin

;  variable = "pr"
;  diri = "/chinook/shomtm62/Observational/"+variable+"/"   ; Obtain input data
;  prfx       = "UW"
;  fin        = addfile(diri + "Mon" + variable + "_" + prfx + ".nc", "r") 
;  fin =addfile("/chinook/shomtm62/Observational/pr/Monpr_UW.nc","r")  
  fin = addfile("/chinook/shomtm62/CMIP5/pr/test.nc","r")
;  fin = addfile("/tera0/shomtm62/RCMs/MonPR_CRCM.nc","r")

  prc1 = fin->pr 		; Define input file 
  prc = short2flt(prc1)
  time = fin->time
;  time = flt2dble(time1)
  lat = fin->lat 		; and predefine dimensions
  lon = fin->lon
 
  printVarSummary(prc) 
  printVarSummary(time)

;===============================================
; Time conversion (not needed)
;===============================================
;  dimx3 = dimsizes(prc)
;  NTIM  = dimx3(0)
;  time@units = "months since 1970-1-1" 
   
;  ntadd = 0 
;  do nt=0,NTIM-1
;    time(nt) = time(nt)*0 + ntadd
;    ntadd = ntadd + 1
;  end do

;  prc&time = time
;  time&time = time
;==============================================

  runlen = 12
  spi = dim_spi_n(prc,runlen,False,0) 
  spi!0 = "time"
  spi!1 = "lat"
  spi!2 = "lon"

  spi&lat = prc&lat
  spi&lon = prc&lon  

  printVarSummary(spi)
  printMinMax(spi,True)

  wks          = gsn_open_wks ("ps","spi")
  gsn_define_colormap(wks,"BlueWhiteOrangeRed") ; suggested colormap for SPI
  gsn_reverse_colormap(wks)           ; Reverse the color map.

  res          = True                   

  res@gsnMaximize      = True             ; uncomment to maximize size
  res@gsnSpreadColors  = True              
  res@cnFillOn         = True               ; color Fill 
;  res@cnFillMode       = "RasterFill"       ; Raster Mode (box-y plot)
  res@cnLinesOn        =  False             ; Turn off contour lines
  res@cnLevelSelectionMode = "ManualLevels" ; set manual contour levels
  res@cnMinLevelValF   = -3.0               ; set min contour level
  res@cnMaxLevelValF   =  3.0               ; set max contour level
  res@cnLevelSpacingF  =  0.5               ; set contour spacing
  res@mpOutlineBoundarySets = "GeophysicalAndUSStates"
  res@gsnAddCyclic = False

  res@mpLimitMode           = "Corners"      ; Corners method of zooming
  res@mpLeftCornerLatF      = 22.5             ; Need to be this method for
  res@mpLeftCornerLonF      = -120.5           ; ALL map projections
  res@mpRightCornerLatF     = 48.5
  res@mpRightCornerLonF     = -61

  res@mpProjection        = "LambertConformal"
  res@mpLambertParallel1F = 30.
  res@mpLambertParallel2F = 60.
  res@mpLambertMeridianF  = -97.0

  res@tiMainString  = "SPI: 1969-1999 (UW)" 
  res@pmTickMarkDisplayMode = "Always"
  res@tmXTOn = False
  res@tmYROn = False

  plot = gsn_csm_contour_map(wks,spi(3,:,:), res) 
end