begin
; CHANGE Aug by the intital month (e.g. Jan, Feb, Mar, Apr, May, Jun
;                                  (Jul, Aug, Sep, Oct, Nov, Dec) 
iMon = OMONS
YYYY = CURRENTYR  ; forecast  year
SEASON = SEASONNAME 
yrS = FIRSTYR    
yrE = LASTYR     
Leadi = LEADI
Leadf = LEADF
opt  = 0  ; use sample standard deviation

;gcmnames <- c("CanSIP","NASA","GFDLS","CCSM4", "CFSv2" /)
cansip  = True
ccsm4   = True 
cfsv2   = True 
nasa    = True 
gfdlS   = True 

monthdef     = (/"", "Jan","Feb","Mar","Apr","May","Jun", \
               "Jul","Aug","Sep","Oct","Nov","Dec" /)

XXX=monthdef(iMon) ; Month in 3 letters


latwnpS = 13.   ; WNP -- From Funk et al BAMS 2018   
latwnpE = 30.   ;        Western V sst = (z(WNP)+z(WEV))/2  
lonwnpS = 160.  ; start  lon
lonwnpE = 220.  ; end   lon for western north Pacific 

latwnpS = 20.   ; WNP -- As defined in WC draft on 17 Aug 2022
latwnpE = 35.   ; end  lat
lonwnpS = 160.  ; start  lon
lonwnpE = 210.  ; end   lon for western north Pacific 

latwepS = -20.  ; WEP -- From Funk et al. BAMS 2018        
latwepE = 20.   ; end
lonwepS = 120.  ; start
lonwepE = 170.  ; end  lon for western equatorial Pacific 

latwepS = -15.  ; WEP -- As defined in WC draft 17Aug2022
latwepE = 20.   ; end
lonwepS = 120.  ; start
lonwepE = 160.  ; end  lon for western equatorial Pacific 

latwspS = -30.  ; WSP -- As defined in WC draft 17Aug2022
latwspE = -15.  ; end
lonwspS = 155.  ; start
lonwspE = 210.  ; end  lon for western equatorial Pacific 

latni4S = -5.    ; start lat for Nino 4
latni4E = 5.    ; end
lonni4S = 160.  ; start
lonni4E = 210.  ; end  lon for Nino 4 

latn34S = -5.    ; start lat for Nino34
latn34E = 5.    ; end
lonn34S = 190.  ; start
lonn34E = 240.  ; end  lon for Nino34

dir0 = "./analdat/"
system("if ! test -d " + dir0 +" ; then mkdir " + dir0 + " ; fi")

dir0 = "./epsfiles/"
system("if ! test -d " + dir0 +" ; then mkdir " + dir0 + " ; fi")

dir0 = "./sstdat/"
system("if ! test -d " + dir0 +" ; then mkdir " + dir0 + " ; fi")
dirsst = dir0 

;;;;;;;;;;;;;;;;;;;;;;
;;;;; CanSIP
;;;;;;;;;;;;;;;;;;;;;;

if (cansip) then

  print("... downloading canSIP ")
   
a = addfile("http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.CanSIPS-IC3/.HINDCAST/.MONTHLY/.sst/X/%28115%29/%28290%29/RANGEEDGES/L/%28"+Leadi+"%29/%28"+Leadf+"%29/RANGEEDGES/M/%281.0%29/%2820.0%29/RANGEEDGES/Y/%2835S%29/%2835N%29/RANGEEDGES/S/%280000%201%20"+XXX+"%201982-2020%20%29/VALUES/%5BM/%5Daverage/%5BL/%5Daverage/-9999/setmissing_value/dods","r")

b = addfile("http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.CanSIPS-IC3/.FORECAST/.MONTHLY/.sst/X/%28115E%29/%28290E%29/RANGEEDGES/L/%28"+Leadi+"%29/%28"+Leadf+"%29/RANGEEDGES/M/%281.0%29/%2820.0%29/RANGEEDGES/Y/%2835S%29/%2835N%29/RANGEEDGES/S/%280000%201%20"+XXX+"%202021-"+YYYY+"%29/VALUES/%5BM/%5Daverage/%5BL/%5Daverage/-9999/setmissing_value/dods","r")
 
  print("... downloaded data for canSIP ")

  x := a->sst
  y := b->sst
  if (any(isnan_ieee(x))) then
     if(.not.isatt(x,"_FillValue")) then
         x@_FillValue = -9999.  ;default_fillvalue(typeof(prec))
         x@missing_value = -9999.
     end if
     replace_ieeenan (x, -9999., 0)
     x@missing_value = -9999.
  end if

  times = a->S
  lat   = a->Y
  lon   = a->X
  lat!0 = "lat"
  lon!0 = "lon"
  lat&lat = lat
  lon&lon = lon

  x!0 = "time"
  x!1 = "lat"
  x!2 = "lon"

  x&time = times
  x&lat  = lat
  x&lon  = lon
  x@_FillValue = -9999
  x@description = "Pacific SST "


  time2  = b->S
  y!0    = "time"
  y!1    = "lat"
  y!2    = "lon"
  y&time = time2
  y&lat  = lat
  y&lon  = lon

  z =  array_append_record (x, y, 0)
  time = z&time
  
  yyyymm = cd_calendar(time,-1)
  allyrs = yyyymm/100
  ii     = ind(allyrs .ge. yrS .and. allyrs .le. toint(YYYY))
  
  sst = dim_standardize_n_Wrap(z(ii,:,:),0,0)
  printVarSummary(z)
  delete([/a,b,x,lat,lon,y,z /])

  wnp =  sst(:,{latwnpS:latwnpE},{lonwnpS:lonwnpE})
  wnp := dim_avg_n_Wrap(wnp,(/1,2/))

  wep =  sst(:,{latwepS:latwepE},{lonwepS:lonwepE})
  wep := dim_avg_n_Wrap(wep,(/1,2/))

  ni4 =  sst(:,{latni4S:latni4E},{lonni4S:lonni4E})
  ni4 := dim_avg_n_Wrap(ni4,(/1,2/))

  n34 =  sst(:,{latn34S:latn34E},{lonn34S:lonni4E})
  n34 := dim_avg_n_Wrap(n34,(/1,2/))

;  wvg  = ni4
;  wvg  = (/ ni4 - (wnp + wep)/2 /)

  wsp =  sst(:,{latwspS:latwspE},{lonwspS:lonwspE})
  wsp := dim_avg_n_Wrap(wsp,(/1,2/))

  wv  = wnp
  wv  = (/ (wnp + wep + wsp)/3 /)

  wvg  = n34
  wvg  = (/ n34 - wv /)

; Save data on local disk

  fname = dirsst+"sst"+SEASON+"WPG"+"CanSIP.nc"
  system("/bin/rm -f "+fname)
  ncdf = addfile(fname, "c")

       fAtt               = True            ; assign file attributes
       fAtt@title         = "Regridded Global Data from IRIDL"
       fAtt@source_file   =  "IRIDL-CanSIP IC3   "
       fAtt@Conventions   = "None"
       fAtt@creation_date = systemfunc ("date")
       fileattdef( ncdf, fAtt )            ; copy file attributes

       filedimdef(ncdf,"time",-1,True)

       ncdf->wnp  = wnp
       ncdf->wep  = wep
       ncdf->ni4  = ni4
       ncdf->n34  = n34
       ncdf->wvg  = wvg

   delete([/ yyyymm, allyrs, ii, time, sst,times,time2, ncdf,fname/])
   delete([/wv, wsp, wvg, wnp, wep, ni4, n34/])
 
end if

;;;;;;;;;;;;;;;;;;;;;;
;;;;; CCSM4
;;;;;;;;;;;;;;;;;;;;;;

if (ccsm4) then
   print("... downloading CCM4 ")

    a = addfile("http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.COLA-RSMAS-CCSM4/.MONTHLY/.sst/X/%28115E%29/%28290%29/RANGEEDGES/L/%28"+Leadi+"%29/%28"+Leadf+"%29/RANGEEDGES/M/%281.0%29/%2810.0%29/RANGEEDGES/Y/%2835S%29/%2835N%29/RANGEEDGES/S/%280000%201%20"+XXX+"%201982-2020%29/VALUES/%5BM/%5Daverage/%5BL/%5Daverage/-9999/setmissing_value/dods","r")

    b = addfile("http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.COLA-RSMAS-CCSM4/.MONTHLY/.sst/X/%28115E%29/%28290E%29/RANGEEDGES/L/%28"+Leadi+"%29/%28"+Leadf+"%29/RANGEEDGES/M/%281.0%29/%2828.0%29/RANGEEDGES/Y/%2835S%29/%2835N%29/RANGEEDGES/S/%280000%201%20"+XXX+"%202021-"+YYYY+"%29/VALUES/%5BL/%5Daverage/%5BM/%5Daverage/-9999/setmissing_value/dods","r")


   print("... downloaded data for CCSM4 ...")

;;;;;;;; use logical to get only available model data

  x = a->sst
  y = b->sst
  if (any(isnan_ieee(x))) then
     if(.not.isatt(x,"_FillValue")) then
         x@_FillValue = -9999.  ;default_fillvalue(typeof(prec))
         x@missing_value = -9999.
     end if
     replace_ieeenan (x, -9999., 0)
     x@missing_value = -9999.
  end if

  times = a->S
  lat   = a->Y
  lon   = a->X
  lat!0 = "lat"
  lon!0 = "lon"
  lat&lat = lat
  lon&lon = lon

  x!0 = "time"
  x!1 = "lat"
  x!2 = "lon"

  x&time = times
  x&lat  = lat  
  x&lon  = lon
  x@_FillValue = -9999
  x@description = "Pacific SST "

  time2  = b->S
  y!0    = "time"
  y!1    = "lat"
  y!2    = "lon"
  y&time = time2
  y&lat  = lat  
  y&lon  = lon

  z =  array_append_record (x, y, 0)
  time = z&time

  yyyymm = cd_calendar(time,-1)
  allyrs = yyyymm/100
  ii     = ind(allyrs .ge. yrS .and. allyrs .le. toint(YYYY))

  sst = dim_standardize_n_Wrap(z(ii,:,:),0,0)
  printVarSummary(z)
  delete([/a,x,lat,lon,y,z /])
 
  wnp =  sst(:,{latwnpS:latwnpE},{lonwnpS:lonwnpE})
  wnp := dim_avg_n_Wrap(wnp,(/1,2/))

  wep =  sst(:,{latwepS:latwepE},{lonwepS:lonwepE})
  wep := dim_avg_n_Wrap(wep,(/1,2/))

  ni4 =  sst(:,{latni4S:latni4E},{lonni4S:lonni4E})
  ni4 := dim_avg_n_Wrap(ni4,(/1,2/))

  n34 =  sst(:,{latn34S:latn34E},{lonn34S:lonni4E})
  n34 := dim_avg_n_Wrap(n34,(/1,2/))

;  wvg  = ni4
;  wvg  = (/ ni4 - (wnp + wep)/2 /)

  wsp =  sst(:,{latwspS:latwspE},{lonwspS:lonwspE})
  wsp := dim_avg_n_Wrap(wsp,(/1,2/))
  
  wv   = wnp 
  wv   = (/ (wnp + wep + wsp)/3 /)
  
  wvg  = n34
  wvg  = (/ n34 - wv /)

;  wnpZ = dim_standardize_Wrap(wnp,opt)
;  wnpAve = avg(wnp)
;  wnpStd = stddev(wnp)
  
; Save data on local disk

  fname = dirsst+"sst"+SEASON+"WPG"+"CCSM4.nc"
  system("/bin/rm -f "+fname)
  ncdf = addfile(fname, "c")

       fAtt               = True            ; assign file attributes
       fAtt@title         = "Regridded Global Data from IRIDL"
       fAtt@source_file   =  "IRIDL-CCM4  "
       fAtt@Conventions   = "None"
       fAtt@creation_date = systemfunc ("date")
       fileattdef( ncdf, fAtt )            ; copy file attributes

       filedimdef(ncdf,"time",-1,True)

       ncdf->wnp  = wnp
       ncdf->wep  = wep
       ncdf->ni4  = ni4
       ncdf->n34  = n34
       ncdf->wvg  = wvg

   delete([/ yyyymm, allyrs, ii, time, sst,times,time2, ncdf,fname/])
   delete([/wv, wsp, wvg, wnp, wep, ni4, n34/])

end if ; ccsm4

;;;;;;;;;;;;;;;;;;;;;;
;;;;; CFSv2
;;;;;;;;;;;;;;;;;;;;;;

if(cfsv2) then
;
; SOURCES .Models .NMME .NCEP-CFSv2
;  a: .HINDCAST .PENTAD_SAMPLES .MONTHLY .sst
;    :a: .FORECAST .PENTAD_SAMPLES .MONTHLY .sst :a
;   appendstream
;  M (1.0) (28.0) RANGEEDGES
;  Y (35S) (35N) RANGEEDGES
;  S (0000 1 Feb 1982-2020) VALUES
;  X (115E) (290E) RANGEEDGES
;  L (1.5) (3.5) RANGEEDGES
;  [M ]average
;  [L ]average
;  -9999 setmissing_value

    a = addfile("http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.NCEP-CFSv2/a:/.HINDCAST/.PENTAD_SAMPLES/.MONTHLY/.sst/:a:/.FORECAST/.PENTAD_SAMPLES/.MONTHLY/.sst/:a/appendstream/M/%281.0%29/%2828.0%29/RANGEEDGES/Y/%2835S%29/%2835N%29/RANGEEDGES/S/%280000%201%20"+XXX+"%201982-2020%29/VALUES/X/%28115E%29/%28290E%29/RANGEEDGES/L/%28"+Leadi+"%29/%28"+Leadf+"%29/RANGEEDGES/-999/setmissing_value/%5BM/%5Daverage/%5BL/%5Daverage/-9999/setmissing_value/dods","r")

   b = addfile("http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.NCEP-CFSv2/.FORECAST/.PENTAD_SAMPLES/.MONTHLY/.sst/X/%28115E%29/%28290E%29/RANGEEDGES/L/%28"+Leadi+"%29/%28"+Leadf+"%29/RANGEEDGES/M/%281.0%29/%2828.0%29/RANGEEDGES/Y/%2835S%29/%2835N%29/RANGEEDGES/S/%280000%201%20"+XXX+"%202021-"+YYYY+"%29/VALUES/%5BL/%5Daverage/%5BM/%5Daverage/-9999/setmissing_value/dods","r")

  print("... downloaded data for CVSv2 ...")

;;;;;;;; use logical to get only available model data

  x = a->sst
  y = b->sst
  if (any(isnan_ieee(x))) then
     if(.not.isatt(x,"_FillValue")) then
         x@_FillValue = -9999.  ;default_fillvalue(typeof(prec))
         x@missing_value = -9999.
     end if
     replace_ieeenan (x, -9999., 0)
     x@missing_value = -9999.
  end if

  times = a->S
  lat   = a->Y
  lon   = a->X
  lat!0 = "lat"
  lon!0 = "lon"
  lat&lat = lat
  lon&lon = lon

  x!0 = "time"
  x!1 = "lat"
  x!2 = "lon"

  x&time = times
  x&lat  = lat
  x&lon  = lon
  x@_FillValue = -9999
  x@description = "Pacific SST "

  time2  = b->S
  y!0    = "time"
  y!1    = "lat"
  y!2    = "lon"
  y&time = time2
  y&lat  = lat
  y&lon  = lon

  z =  array_append_record (x, y, 0)
  time = z&time

  yyyymm = cd_calendar(time,-1)
  allyrs = yyyymm/100
  ii     = ind(allyrs .ge. yrS .and. allyrs .le. toint(YYYY))

  sst = dim_standardize_n_Wrap(z(ii,:,:),0,0)
  printVarSummary(z)
  delete([/a,x,lat,lon,y,z /])

  wnp :=  sst(:,{latwnpS:latwnpE},{lonwnpS:lonwnpE})
  wnp := dim_avg_n_Wrap(wnp,(/1,2/))

  wep :=  sst(:,{latwepS:latwepE},{lonwepS:lonwepE})
  wep := dim_avg_n_Wrap(wep,(/1,2/))

  ni4 :=  sst(:,{latni4S:latni4E},{lonni4S:lonni4E})
  ni4 := dim_avg_n_Wrap(ni4,(/1,2/))

  n34 :=  sst(:,{latn34S:latn34E},{lonn34S:lonn34E})
  n34 := dim_avg_n_Wrap(n34,(/1,2/))

;  wvg  = ni4
;  wvg  = (/ ni4 - (wnp + wep)/2 /)

  wsp =  sst(:,{latwspS:latwspE},{lonwspS:lonwspE})
  wsp := dim_avg_n_Wrap(wsp,(/1,2/))
  
  wv  = wnp 
  wv  = (/ (wnp + wep + wsp)/3 /)
  
  wvg  = n34
  wvg  = (/ n34 - wv /)

; wnpZ = dim_standardize_Wrap(wnp,opt)
; wnpAve = avg(wnp)
; wnpStd = stddev(wnp)

; Save data on local disk

  fname = dirsst+"sst"+SEASON+"WPG"+"CFSv2.nc"
  system("/bin/rm -f "+fname)
  ncdf = addfile(fname, "c")

       fAtt               = True            ; assign file attributes
       fAtt@title         = "Regridded Global Data from IRIDL"
       fAtt@source_file   =  "IRIDL-CFSv2 "
       fAtt@Conventions   = "None"
       fAtt@creation_date = systemfunc ("date")
       fileattdef( ncdf, fAtt )            ; copy file attributes

       filedimdef(ncdf,"time",-1,True)

       ncdf->wnp  = wnp
       ncdf->wep  = wep
       ncdf->ni4  = ni4
       ncdf->n34  = n34
       ncdf->wvg  = wvg

   delete([/ yyyymm, allyrs, ii, time, sst,times,time2, ncdf,fname/])
   delete([/wv, wsp, wvg, wnp, wep, ni4, n34/])

end if ;cfsv2 

;;;;;;;;;;;;;;;;;;;;;;
;;;;; GFDL-SPEAR
;;;;;;;;;;;;;;;;;;;;;;

if (gfdlS) then

a = addfile("http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.GFDL-SPEAR/.HINDCAST/.MONTHLY/.sst_regridded/X/%28115E%29/%28290E%29/RANGEEDGES/L/%28"+Leadi+"%29/%28"+Leadf+"%29/RANGEEDGES/M/%281.0%29/%2815.0%29/RANGEEDGES/Y/%2835S%29/%2835N%29/RANGEEDGES/S/%280000%201%20"+XXX+"%201991-2020%29/VALUES/%5BM/%5Daverage/%5BL/%5Daverage/-9999/setmissing_value/dods","r")

b = addfile("http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.GFDL-SPEAR/.FORECAST/.MONTHLY/.sst_regridded/X/%28115E%29/%28290E%29/RANGEEDGES/L/%28"+Leadi+"%29/%28"+Leadf+"%29/RANGEEDGES/M/%281.0%29/%2830.0%29/RANGEEDGES/Y/%2835S%29/%2835N%29/RANGEEDGES/S/%280000%201%20"+XXX+"%202021-"+YYYY+"%29/VALUES/%5BM/%5Daverage/%5BL/%5Daverage/-9999/setmissing_value/dods","r")

  print("... downloaded data for GFDLS ...")

;;;;;;;; use logical to get only available model data

  x = a->sst_regridded
  y = b->sst_regridded
  if (any(isnan_ieee(x))) then
     if(.not.isatt(x,"_FillValue")) then
         x@_FillValue = -9999.  ;default_fillvalue(typeof(prec))
         x@missing_value = -9999.
     end if
     replace_ieeenan (x, -9999., 0)
     x@missing_value = -9999.
  end if

  times = a->S
  lat   = a->Y
  lon   = a->X
  lat!0 = "lat"
  lon!0 = "lon"
  lat&lat = lat
  lon&lon = lon

  x!0 = "time"
  x!1 = "lat"
  x!2 = "lon"

  x&time = times
  x&lat  = lat
  x&lon  = lon
  x@_FillValue = -9999
  x@description = "Pacific SST "

  time2  = b->S
  y!0    = "time"
  y!1    = "lat"
  y!2    = "lon"
  y&time = time2
  y&lat  = lat
  y&lon  = lon

  z =  array_append_record (x, y, 0)
  time := z&time

  yyyymm = cd_calendar(time,-1)
  allyrs = yyyymm/100
  ii     = ind(allyrs .ge. yrS .and. allyrs .le. toint(YYYY))

  sst = dim_standardize_n_Wrap(z(ii,:,:),0,0)
  printVarSummary(z)
  delete([/a,x,lat,lon,y,z /])

  wnp :=  sst(:,{latwnpS:latwnpE},{lonwnpS:lonwnpE})
  wnp := dim_avg_n_Wrap(wnp,(/1,2/))

  wep :=  sst(:,{latwepS:latwepE},{lonwepS:lonwepE})
  wep := dim_avg_n_Wrap(wep,(/1,2/))

  ni4 :=  sst(:,{latni4S:latni4E},{lonni4S:lonni4E})
  ni4 := dim_avg_n_Wrap(ni4,(/1,2/))

  n34 :=  sst(:,{latn34S:latn34E},{lonn34S:lonn34E})
  n34 := dim_avg_n_Wrap(n34,(/1,2/))
 
;  wvg  = ni4
;  wvg  = (/ ni4 - (wnp + wep)/2 /)

  wsp =  sst(:,{latwspS:latwspE},{lonwspS:lonwspE})
  wsp := dim_avg_n_Wrap(wsp,(/1,2/))
  
  wv  = wnp 
  wv  = (/ (wnp + wep + wsp)/3 /)
  
  wvg  = n34
  wvg  = (/ n34 - wv /)

; Save data on local disk

  fname = dirsst+"sst"+SEASON+"WPG"+"GFDLS.nc"
  system("/bin/rm -f "+fname)
  ncdf = addfile(fname, "c")

       fAtt               = True            ; assign file attributes
       fAtt@title         = "Regridded Global Data from IRIDL"
       fAtt@source_file   =  "IRIDL-GFDL SPEAR"
       fAtt@Conventions   = "None"
       fAtt@creation_date = systemfunc ("date")
       fileattdef( ncdf, fAtt )            ; copy file attributes

       filedimdef(ncdf,"time",-1,True)

       ncdf->wnp  = wnp
       ncdf->wep  = wep
       ncdf->ni4  = ni4
       ncdf->n34  = n34
       ncdf->wvg  = wvg

   delete([/ yyyymm, allyrs, ii, time, sst,times,time2, ncdf,fname/])
   delete([/wv, wsp, wvg, wnp, wep, ni4, n34/])

end if  ; gfdlS

;;;;;;;;;;;;;;;;;;;;;;
;;;;; NASA-GEOS
;;;;;;;;;;;;;;;;;;;;;;

if (nasa) then

 a = addfile("http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.NASA-GEOSS2S/.HINDCAST/.MONTHLY/.sst/X/%28115E%29/%28290E%29/RANGEEDGES/L/%28"+Leadi+"%29/%28"+Leadf+"%29/RANGEEDGES/M/%281.0%29/%284.0%29/RANGEEDGES/Y/%2835S%29/%2835N%29/RANGEEDGES/S/%280000%201%20"+XXX+"%201982-2016%29/VALUES/%5BM/%5Daverage/%5BL/%5Daverage/-9999/setmissing_value/dods","r")


 b = addfile("http://iridl.ldeo.columbia.edu/SOURCES/.Models/.NMME/.NASA-GEOSS2S/.FORECAST/.MONTHLY/.sst/X/%28115E%29/%28290E%29/RANGEEDGES/L/%28"+Leadi+"%29/%28"+Leadf+"%29/RANGEEDGES/M/%281.0%29/%2810.0%29/RANGEEDGES/Y/%2835S%29/%2835N%29/RANGEEDGES/S/%280000%201%20"+XXX+"%202017-"+YYYY+"%29/VALUES/%5BM/%5Daverage/%5BL/%5Daverage/-9999/setmissing_value/dods","r")


  print("... downloaded data for NASA GEOS2S  ...")

;;;;;;;; use logical to get only available model data

  x = a->sst
  y = b->sst
  if (any(isnan_ieee(x))) then
     if(.not.isatt(x,"_FillValue")) then
         x@_FillValue = -9999.  ;default_fillvalue(typeof(prec))
         x@missing_value = -9999.
     end if
     replace_ieeenan (x, -9999., 0)
     x@missing_value = -9999.
  end if

  times = a->S
  lat   = a->Y
  lon   = a->X
  lat!0 = "lat"
  lon!0 = "lon"
  lat&lat = lat
  lon&lon = lon

  x!0 = "time"
  x!1 = "lat"
  x!2 = "lon"

  x&time = times
  x&lat  = lat
  x&lon  = lon
  x@_FillValue = -9999
  x@description = "Pacific SST "

  time2  := b->S
  y!0    = "time"
  y!1    = "lat"
  y!2    = "lon"
  y&time = time2
  y&lat  = lat
  y&lon  = lon

  z =  array_append_record (x, y, 0)
  time = z&time

  yyyymm = cd_calendar(time,-1)
  allyrs = yyyymm/100
  ii     = ind(allyrs .ge. yrS .and. allyrs .le. toint(YYYY))

  sst = dim_standardize_n_Wrap(z(ii,:,:),0,0)
  printVarSummary(z)
  delete([/a,x,lat,lon,y,z /])

  wnp :=  sst(:,{latwnpS:latwnpE},{lonwnpS:lonwnpE})
  wnp := dim_avg_n_Wrap(wnp,(/1,2/))

  wep :=  sst(:,{latwepS:latwepE},{lonwepS:lonwepE})
  wep := dim_avg_n_Wrap(wep,(/1,2/))

  ni4 :=  sst(:,{latni4S:latni4E},{lonni4S:lonni4E})
  ni4 := dim_avg_n_Wrap(ni4,(/1,2/))

  n34 :=  sst(:,{latn34S:latn34E},{lonn34S:lonni4E})
  n34 := dim_avg_n_Wrap(n34,(/1,2/))

;  wvg  = ni4
;  wvg  = (/ ni4 - (wnp + wep)/2 /)

  wsp =  sst(:,{latwspS:latwspE},{lonwspS:lonwspE})
  wsp := dim_avg_n_Wrap(wsp,(/1,2/))
  
  wv  = wnp 
  wv  = (/ (wnp + wep + wsp)/3 /)
  
  wvg  = n34
  wvg  = (/ n34 - wv /)

; Save data on local disk

  fname = dirsst+"sst"+SEASON+"WPG"+"NASA.nc"
  system("/bin/rm -f "+fname)
  ncdf = addfile(fname, "c")

       fAtt               = True            ; assign file attributes
       fAtt@title         = "Regridded Global Data from IRIDL"
       fAtt@source_file   =  "IRIDL-NASA  "
       fAtt@Conventions   = "None"
       fAtt@creation_date = systemfunc ("date")
       fileattdef( ncdf, fAtt )            ; copy file attributes

       filedimdef(ncdf,"time",-1,True)

       ncdf->wnp  = wnp
       ncdf->wep  = wep
       ncdf->ni4  = ni4
       ncdf->n34  = n34
       ncdf->wvg  = wvg

   delete([/ yyyymm, allyrs, ii, time, sst,times,time2, ncdf,fname/])
   delete([/wv, wsp, wvg, wnp, wep, ni4, n34/])

end if  ; nasa


end