;*************************************************
; rngpcpampsst_scaterplot.ncl  original- scatter_4.ncl
;;http://www.ncl.ucar.edu/Applications/Scripts/scatter_4.ncl
;
; Concepts illustrated:
;   - Drawing a scatter plot with a regression line
;   - Drawing a time series plot
;   - Calculating the least squared regression for a one dimensional array
;   - Smoothing data so that seasonal cycle is less prominent
;   - Changing the markers in an XY plot
;   - Changing the marker color in an XY plot
;   - Changing the marker size in an XY plot
;
;*************************************************
;
; These files are loaded by default in NCL V6.2.0 and newer
 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
 load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"

begin
;********************************
; Read GPCP monthly rainfall data
;********************************
   ymdStrt = 19790101   ; start yyyymmdd
   ymdLast = 20081201   ; end yyyymmdd
   
   diri    = "/iitm3/erpas-res/skm/data/gpcp/mon/"       ; input directory   
   fili    = "precip.mon.mean.1979jan2017oct.nc" ; input file 
;;;********************************

   in1  = addfile(diri+fili ,"r")
   
   tm    = in1->time                          ; all times on file
   ymd     = cd_calendar(tm, -2)            ; yyyymmdd
   
;;;;
   iStrt   = ind(ymd.eq.ymdStrt)              ; index start
   iLast   = ind(ymd.eq.ymdLast)              ; index last 
   
   delete(tm)
   delete(ymd)

;***********************************************************
    
   rn= in1->precip(iStrt:iLast,{-10:5},{75:100})
   rnio=wgt_areaave(rn,1.0,1.0,0);;;create area averaged rain over equatorial Indian ocean
;;   dimrn=dimsizes(rnio)
;;   print("dimension of rnio=" +dimrn)
;************************************************
; Read monthly AMIP SST  data
;********************************

   ymdSt = 19790116   ; start yyyymmdd t=1309
   ymdLa = 20081216   ; End yyyymmdd t=1668

   vars     = "tos"   ; name of variable in file  

   dirs    = "/iitm3/erpas-res/skm/data/cmp5/gfdl/ampsstseaice/" ; input directory   
   fils    = "tos_input4MIPs_SSTsAndSeaIce_CMIP_PCMDI-AMIP-1-1-3_gn_187001-201706.nc" ; input file 

;***********************************************************
; Read user specified time and create required yyyyddd                    
;***********************************************************
   in2      = addfile (dirs+fils , "r")                          

   tms    = in2->time                 ; all times on file
   ymds   = cd_calendar(tms, -2)    ; yyyymmdd
  
;;;;
   iSts   = ind(ymds.eq.ymdSt)       ; index start
   iLas   = ind(ymds.eq.ymdLa)       ; index last 
      
   delete(tms)
   delete(ymds)

;***********************************************************   
   st= in2->tos(iSts:iLas,{-10:5},{75:100})
   stio=wgt_areaave(st,1.0,1.0,0)   ;;;create area averaged rain
;************************************************
; calculate the regression coefficient for monthly   
; GPCP rainfall and AMIP sst over equatorial Indian Ocean
; Note regline works on one dimensional arrays.   
;************************************************
   rc    = regline(stio,rnio)
;************************************************
; Create an array to hold both the original data
; and the calculated regression line.
;************************************************
 data      = new ( (/2,dimsizes(rnio)/), typeof(rnio))
 data(0,:) = rnio
; y = mx+b 
; m is the slope:       rc      returned from regline
; b is the y intercept: rc@yave attribute of rc returned from regline
 data(1,:) = rc*(stio-rc@xave) + rc@yave
;************************************************
; plotting parameters
;************************************************
 wks  = gsn_open_wks("eps","rngpcpampsst_scater792k8")  ; send graphics to EPS file
 
 res                     = True                   ; plot mods desired
 res@gsnMaximize         = True                   ; maximize plot in frame
 res@xyMarkLineModes     = (/"Markers","Lines"/)  ; choose which have markers
 res@xyMarkers           = 16                     ; choose type of marker 
 res@xyMarkerColor       = "red"                  ; Marker color
 res@xyMarkerSizeF       = 0.005                  ; Marker size (default 0.01)
 res@xyDashPatterns      = 1                      ; solid line 
 res@xyLineThicknesses   = (/1,2/)                ; set second line to 2

 res@tiMainString        = "Output from regline"  ; title

 plot  = gsn_csm_xy (wks,rnio,stio,res)        ; create plot
end