;This script computes the trend of a time series using 
;least squares. This example used the GISST SST data set.

;This script will also plot a zonal mean to a .ps file.
;That part of script is optional.

;Author: Sara Rauscher srausche@wisc.edu
;Date: 10.09.02

;load some function scripts
;the following path is for NCAR computers
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

begin

;open file that has GISST data 
f = addfile("gisst.nc","r")

;open file to write computed trend 
system("rm gisst.trend.1949-1994.nc")
out = addfile("gisst.trend.1949-1994.nc","c")

;get needed variables from file
;create dummy time variable if necessary
;X= time = independent var
X = new((/46/),float)
do i=0,45
        X(i)=i+1.0
end do

;get SST data
;SST = Y = dependent variable
Y = f->sst(78:123,:,:)
Y@_FillValue = -999.

;create variables for calculations
;must have same dimensions as Y
SY = new((/180,360/),float)
B = new((/180,360/),float)

;ZERO variables and create some one-dim ones (SX,ST2)
SX = 0.0
SY = 0.0
ST2 = 0.0
B = 0.0

;DO least squares calc
;code modified from http://grads.iges.org/grads/gadoc/udf.html#example

;SUM of X and Y
do i = 0, 45
  SX = SX + X(i)
  SY(:,:) = SY + Y(i,:,:)
end do

SS = 46.    			;number of time steps
SXOSS = SX/SS			;divide sum of x by number of time steps
do i = 0, 45
  T = X(i) - SXOSS
  ST2 = ST2 + T * T     
  B(:,:) = B(:,:) + T * Y(i,:,:) ;calc slope (=B)
end do

B = B/ST2
A = (SY - SX * B)/SS		;calc y-intercept (=A)

slope = B*100			;multiply slope by 100 to get trend/100 years
slope!0 = "lat"			;assign coord vars to slope so grads can read data
slope&lat = f->lat
slope!1 = "lon"
slope&lon = f->lon

;plot zonal avg			;calc & plot zonal average
;plotting code modified from http://www.cgd.ucar.edu/csm/support/Data_P/zonal.shtml

zave = dim_avg(slope)
print(zave)
zave!0 = "lat"
zave&lat = f->lat
zave@long_name = "Trend (degrees/Century)"

wks  = gsn_open_wks("ps","zonal")
r            = True                            ; plot mods desired
r@tiMainString = "Global Zonal Average of GISST SST Trend - 1949-1994"
plot=gsn_csm_xy(wks,zave(:),slope&lat,r)  

out->slope = slope		;write out slope to file for later use/plots
out->lat = f->lat
out->lon = f->lon
end