load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

 begin
  diri  = "/project/cas/shea/"
  fili  = "co_tmax1895-1997.nc"
 
  a     = addfile(diri+fili, "r")
;;tmax  = a->tmax                ; get data
  tmax  = dble2flt( a->tmax )    ; This works!
  tmax  = tmax*0.01              ; make units degrees C
  printVarSummary(tmax)          ; Note: NCL imports all meta data
  printMinMax(tmax, True)

  dimtx = dimsizes( tmax )
  ntim  = dimtx(0)
  nlat  = dimtx(1)
  mlon  = dimtx(2)

  nMsg  = num(ismissing(tmax))
  print("     ")
  print("nMsg="+nMsg)

  lat   = tmax&lat
  lon   = tmax&lon
  
  yrs   = ispan(1,ntim,1)        
  TMAX  = tmax(lat|:,lon|:,time|:)  ; temporary work array
;*********************************************************
; regression coefficient: add meta data for plotting
;*********************************************************
  rc    = regCoef(yrs, TMAX )
  
  rc@long_name   = "regression coefficient"
  rc!0   = "lat"      ; name dimensions
  rc!1   = "lon"
  rc&lat =  lat       ; assign coordinate values to named dimensions
  rc&lon =  lon

  printVarSummary(rc)           ; variable overview
  printMinMax(rc, True)

;*********************************************************
; statistics for regression coefficient: add meta for plot
;*********************************************************
  tval = onedtond(rc@tval , dimsizes(rc))           ; tval(nlat,mlon)
  printVarSummary(tval)
  printMinMax(tval, True)

  df   = onedtond(rc@nptxy, dimsizes(rc)) - 2       ; df(nlat,mlon)

  prob = cdft_p(tval, df)                           ; prob(nlat,nlon)
  prob@long_name = "probability: reg coef"
  copy_VarCoords(rc, prob)
  printVarSummary(prob)
  printMinMax(prob, True)

;*********************************************************
; detrend 
;*********************************************************
 
  yDtrend   = dtrend_msg(yrs, TMAX,False,True)      ; (lat,lon,time)  
  yDtrend@long_name = tmax@long_name
  yDtrend@units     = tmax@lunits
  copy_VarCoords(TMAX, yDtrend)
  printVarSummary(yDtrend)

  slope     = onedtond(yDtrend@slope, (/nlat,mlon/) )
  slope@long_name = "linear slope"
  slope@units     = "C/year"
  slope!0   = "lat"
  slope!1   = "lon"
  slope&lat =  lat
  slope&lon =  lon

  printVarSummary(slope)
  printMinMax(slope, True)

;************************************************
; create plots
;************************************************
  wks = gsn_open_wks("ps" ,"baker")                ; open a ps file
  gsn_define_colormap(wks,"BlAqGrYeOrRe")        ; choose colormap

  res                       = True     ; plot mods desired

  res@cnFillOn              = True     ; turn on color fill
  res@cnLinesOn             = False    ; turn of contour lines
 ;res@cnLevelSpacingF       = 0.5      ; contour spacing
  res@gsnSpreadColors       = True     ; use full range of color map
  res@lbLabelAutoStride     = True

  res@pmTickMarkDisplayMode = "Always"; use NCL default lat/lon labels
  res@gsnAddCyclic          = False   

  res@mpMinLatF             = min(lat)
  res@mpMaxLatF             = max(lat)
  res@mpMinLonF             = min(lon)
  res@mpMaxLonF             = max(lon) 

  plot = gsn_csm_contour_map_ce(wks, rc, res)
  plot = gsn_csm_contour_map_ce(wks, prob, res)

  plot = gsn_csm_contour_map_ce(wks, slope, res)

end