;*************************************************
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
;************************************************
; create pointer to file and read in temperature
;************************************************
   diri = "./"
   ;diri = "/media/srathore/binny/phd1/phd/argo/ensemble/modify/modify/modify_soda/ncl/" 
   fili = "heat_gd1.nc"
   in  = addfile(diri+fili,"r") 
  
   tmp  = in->HEAT_GD1
   printVarSummary(tmp)
   ts   = tmp(DEPTH_AXIS |:, TMONTHLY |:)               ; reorder variable
   end
   delete(tmp)                                  ; no longer needed

;************************************************
; create x and calculate the regression coefficients (slopes, trends)    
;************************************************
   ;time         = ts&TMONTHLY                     ; days since 1850-01-01
   ;time=ispan(1, 1800, 1)
   ;rc           = regCoef(time,ts)        
   ;rc            = regline_stats(time,ts)
   ;print(rc)

;************************************************
; Please check this step, it is showing error
;************************************************
   s = 58
   rc = new ((/58/),double)
   up = new(58, double)
   time         = ts&TMONTHLY                       ; days since 1850-01-01
   ;time=ispan(1, 1800, 1)
   i = 0
   do i = 0,s-1
      temp = regline_stats(time,ts(i,:))
      rc(i)= temp
      b95_temp(i)= temp&b95
      up(i) = b95_temp(0) 
   end do

   rc@long_name = "regression coefficient (trend)" 
   copy_VarCoords(ts(:,0), rc)          ; copy lat,lon coords
   rc           = rc                  ; (Pa/day)*(365_day/year)
   rc@units     = "j/m^2/year"
   

   fout=addfile("ensm_gd_trend.nc","c")
   fout->trend=rc
   fout->pvalue=up

;************************************************
; for illustration: change units to j/m^2/year
;***********************************************

   ;rc@long_name = "regression coefficient (trend)" 
   ;copy_VarCoords(ts(:,:,0), rc)          ; copy lat,lon coords
   ;rc           = rc                  ; (Pa/day)*(365_day/year)
   ;rc@units     = "C/year"