;***********************
;SOI 20CRV3 ensemble
;***********************

;***********************
begin

  tot_ens=80
  debug=0
  nvar=10

  optEOF = True       
  optEOF@jopt = 0   ; This is the default; most commonly used; no need to specify.
  optETS = False

    iy1     = 1981
    iy2     = 2010
    iystart = 1806
    iyend   = 2015


  stdtotal=new( (/(iyend-(iystart-1))*12 /),float)
  stdtotal=0
  do nens=0,tot_ens-1
    cens=sprinti("%.03i",nens+1)
    filein = "/Projects/reanl20/ncfiles_20CRv3/mnmean/PRMSL/PRMSL.1806-2015.mnmean_mem"+cens+".nc"
    print(filein+" ")

    fm   = addfile (filein,"r")
    pres = fm->PRMSL(:,:,:)
;   lat  = fm->lat
    latD = -12.4667
    lonD = 130.8333

    latT = -17.6667
    lonT = 360-149.45

    xdar = new( (/30 /), float)
    xtah = new( (/30 /), float)

    cl1            = (iy1 - iystart)*12
    cl2            = 11+(iy2 - iystart)*12
    slpt           = pres(lat|:,lon|:,time|:)
    slpclim        = clmMonLLT( slpt(:,:,cl1:cl2 ))
    slpta          = calcMonAnomLLT(slpt(:,:,:),slpclim)
    slptadarwin    = slpta({latD},{lonD},:)
    slptatahiti    = slpta({latT},{lonT},:)
    slptdarwin_std = slptadarwin
    slpttahiti_std = slptatahiti


; loop through months
; get std for each month

do im = 0, 11
  k = im + (iy1 - iystart)*12
  do iy = iy1,iy2
    xdar(iy-iy1) = slptadarwin(k)
    xtah(iy-iy1) = slptatahiti(k)
    k = k + 12
  end do
; get std
  stddar = stddev(xdar)
  stdtah = stddev(xtah)
  if(debug.eq.1)then
	print("mon "+im+" "+stddar)
	print("mon "+im+" "+stdtah)
  end if
  k = im
  do iy = iystart,iyend
   slptdarwin_std(k) = slptadarwin(k)/stddar
   slpttahiti_std(k) = slptatahiti(k)/stdtah
   k = k+12
  end do

end do

; get difference

  slpsoi     = slptdarwin_std-slpttahiti_std
  slpsoi_std = -1*(slpsoi/stddev(slpsoi(cl1:cl2)))
  stdtotal=stdtotal+slpsoi_std
end do

stdtotal=stdtotal/tot_ens

soiout     = "soi.1951-1980.txt0"
darout     = "dar.1951-1980.txt0"
tahout     = "tah.1951-1980.txt0"
nlines     = (2014-1805)+7
datasoi    = new( nlines, "string")
datadar    = new( nlines, "string")
datatah    = new( nlines, "string")
datasoi(0) = iystart+" "+iyend
datadar(0) = iystart+" "+iyend
datatah(0) = iystart+" "+iyend
print (iystart+" "+iyend)

do n = iystart,iyend
    nn         = (n-iystart)*12
    m          = n-(iystart-1)
    datasoi(m) = (n+" "+sprintf("%7.3f",stdtotal(nn))+" "+sprintf("%7.3f",stdtotal(nn+1))+" "+sprintf("%7.3f",stdtotal(nn+2))+" "+sprintf("%7.3f",stdtotal(nn+3))+" "+sprintf("%7.3f",stdtotal(nn+4))+" "+sprintf("%7.3f",stdtotal(nn+5))+" "+sprintf("%7.3f",stdtotal(nn+6))+" "+sprintf("%7.3f",stdtotal(nn+7))+" "+sprintf("%7.3f",stdtotal(nn+8))+" "+sprintf("%7.3f",stdtotal(nn+9))+" "+sprintf("%7.3f",stdtotal(nn+10))+" "+sprintf("%7.3f",stdtotal(nn+11)))


datadar(m) = (n+" "+sprintf("%7.3f",slptdarwin_std(nn))+" "+sprintf("%7.3f",slptdarwin_std(nn+1))+" "+sprintf("%7.3f",slptdarwin_std(nn+2))+" "+sprintf("%7.3f",slptdarwin_std(nn+3))+" "+sprintf("%7.3f",slptdarwin_std(nn+4))+" "+sprintf("%7.3f",slptdarwin_std(nn+5))+" "+sprintf("%7.3f",slptdarwin_std(nn+6))+" "+sprintf("%7.3f",slptdarwin_std(nn+7))+" "+sprintf("%7.3f",slptdarwin_std(nn+8))+" "+sprintf("%7.3f",slptdarwin_std(nn+9))+" "+sprintf("%7.3f",slptdarwin_std(nn+10))+" "+sprintf("%7.3f",slptdarwin_std(nn+11)))

    datatah(m) = (n+" "+sprintf("%7.3f",slpttahiti_std(nn))+" "+sprintf("%7.3f",slpttahiti_std(nn+1))+" "+sprintf("%7.3f",slpttahiti_std(nn+2))+" "+sprintf("%7.3f",slpttahiti_std(nn+3))+" "+sprintf("%7.3f",slpttahiti_std(nn+4))+" "+sprintf("%7.3f",slpttahiti_std(nn+5))+" "+sprintf("%7.3f",slpttahiti_std(nn+6))+" "+sprintf("%7.3f",slpttahiti_std(nn+7))+" "+sprintf("%7.3f",slpttahiti_std(nn+8))+" "+sprintf("%7.3f",slpttahiti_std(nn+9))+" "+sprintf("%7.3f",slpttahiti_std(nn+10))+" "+sprintf("%7.3f",slpttahiti_std(nn+11)))
end do

print(m+" m nlines "+nlines)
datasoi(m+1) = " -999"
datasoi(m+2) = "SOI from the 20CRv3"
datasoi(m+3) = " noaa/psl"
datasoi(m+4) = " https://psl.noaa.gov/data/20thC_Rean/timeseries/"
datasoi(m+5) = " created using individual ensemble members"
asciiwrite (soiout , datasoi)

datadar(m+1)=" -999"
datadar(m+2)="DARWIN SLP (std) from the 20CRv3"
datadar(m+3)=" noaa/psl"
datasoi(m+4) = " https://psl.noaa.gov/data/20thC_Rean/timeseries/"
datasoi(m+5) = " created using individual ensemble members"
asciiwrite (darout , datadar)

datatah(m+1)=" -999"
datatah(m+2)="TAHITI SLP (std)from the 20CRv3"
datatah(m+3)=" noaa/psl"
datasoi(m+4) = " https://psl.noaa.gov/data/20thC_Rean/timeseries/"
datasoi(m+5) = " created using individual ensemble members"
asciiwrite (tahout , datatah)


 end