; Calculates the IPO pattern, timeseries, and spectra.
;
; Variables used: ts
;
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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/shea_util.ncl"
load "$CVDP_SCRIPTS/functions.ncl"

begin
  print("Starting: ipo.ncl")
  
  SCALE_TIMESERIES = getenv("SCALE_TIMESERIES")  
  OUTPUT_DATA      = getenv("OUTPUT_DATA")  
  PNG_SCALE        = tofloat(getenv("PNG_SCALE"))
  OPT_CLIMO        = getenv("OPT_CLIMO")
  CLIMO_SYEAR      = toint(getenv("CLIMO_SYEAR"))
  CLIMO_EYEAR      = toint(getenv("CLIMO_EYEAR"))
  OUTPUT_TYPE      = getenv("OUTPUT_TYPE") 
  COLORMAP         = getenv("COLORMAP")  
  
  nsim = numAsciiRow("namelist_byvar/namelist_ts")
  na = asciiread("namelist_byvar/namelist_ts",(/nsim/),"string")
  names = new(nsim,"string")
  paths = new(nsim,"string")
  syear = new(nsim,"integer",-999)
  eyear = new(nsim,"integer",-999)
  delim = "|"

  do gg = 0,nsim-1
     names(gg) = str_strip(str_get_field(na(gg),1,delim))
     paths(gg) = str_strip(str_get_field(na(gg),2,delim))
     syear(gg) = stringtointeger(str_strip(str_get_field(na(gg),3,delim)))
     eyear(gg) = stringtointeger(str_strip(str_get_field(na(gg),4,delim)))
  end do
  nyr = eyear-syear+1
  nyr_max = max(nyr)

  pi=4.*atan(1.0)
  rad=(pi/180.)
         
  wks_type = OUTPUT_TYPE
  if (wks_type.eq."png") then
     wks_type@wkWidth = 1500*PNG_SCALE
     wks_type@wkHeight = 1500*PNG_SCALE
  end if
  wks = gsn_open_wks(wks_type,getenv("OUTDIR")+"ipo")
  wks2 = gsn_open_wks(wks_type,getenv("OUTDIR")+"ipo.powspec")
  wks3 = gsn_open_wks(wks_type,getenv("OUTDIR")+"ipo.timeseries")  
  if (COLORMAP.eq."0") then
     gsn_define_colormap(wks,"ncl_default")   
     gsn_define_colormap(wks2,"cb_9step")
     gsn_define_colormap(wks3,"ncl_default")
  end if
  if (COLORMAP.eq."1") then
     gsn_define_colormap(wks,"BlueDarkRed18")   
     gsn_define_colormap(wks2,"cb_9step")
     gsn_define_colormap(wks3,"ncl_default")
  end if
  map = new(nsim,"graphic")  
  pspec = new(nsim,"graphic")
  xyplot = new(nsim,"graphic")
  xyplot2 = new(nsim,"graphic")
  if (isfilepresent2("obs_ts")) then
     pspec_obs = new(nsim,"graphic")
  end if

  fca = 1./157.       
  fcb = 0.
  ihp = 0
  nsigma = 1
  nwgt = 217         

  wgt = new(nwgt,float)
  wgt = filwgts_lancos(nwgt,ihp,fca,fcb,nsigma)   ; create low pass filter

  do ee = 0,nsim-1
     sst = data_read_in(paths(ee),"TS",syear(ee),eyear(ee))    ; read in data, orient lats/lons correctly, set time coordinate variable up
     if (isatt(sst,"is_all_missing").or.nyr(ee).lt.40) then
        delete(sst)
        continue
     end if
     sst = where(sst.le.-1.8,-1.8,sst)    ; set all values below -1.8 to -1.8
     d = addfile("$NCARG_ROOT/lib/ncarg/data/cdf/landsea.nc","r")   ; mask out land (this is redundant for data that is already masked)
     basemap = d->LSMASK
     lsm = landsea_mask(basemap,sst&lat,sst&lon)
     sst = mask(sst,conform(sst,lsm,(/1,2/)).ge.1,False)
     delete([/lsm,basemap/])
     delete(d)  

     if (OPT_CLIMO.eq."Full") then
        sst = rmMonAnnCycTLL(sst)
     else
        check_custom_climo(names(ee),syear(ee),eyear(ee),CLIMO_SYEAR,CLIMO_EYEAR)
        temp_arr = sst
        delete(temp_arr&time)
        temp_arr&time = cd_calendar(sst&time,-1)
        climo = clmMonTLL(temp_arr({CLIMO_SYEAR*100+1:CLIMO_EYEAR*100+12},:,:))                 
        delete(temp_arr)
        sst   = calcMonAnomTLL(sst,climo) 
        delete(climo)
     end if

     sst = wgt_runave_n_Wrap(sst,wgt,1,0)   ; apply low pass filter

     coswgt=cos(rad*sst&lat)
     coswgt!0 = "lat"
     coswgt&lat= sst&lat    
     do ff = 0,dimsizes(sst&time)-1
        sst(ff,:,:) = (/ sst(ff,:,:) - wgt_areaave(sst(ff,{-60:70},:),coswgt({-60.:70.}),1.0,0) /)
     end do
     delete(coswgt)

     sst_CW = sst
     sst_CW = SqrtCosWeight(sst)  
     evecv = eofunc(sst_CW({lat|-40:60},{lon|110:290},time|:),2,75)
     delete(sst_CW)
     pcts = eofunc_ts(sst({lat|-40:60},{lon|110:290},time|:),evecv,False)
     pctsS = dim_standardize(pcts(0,:),0)
     delete([/pcts/])

     ipo = sst(0,:,:)
     ipo = ipo@_FillValue
     ipo = (/ regCoef(pctsS,sst(lat|:,lon|:,time|:)) /)
     ipo@syear = syear(ee)
     ipo@eyear = eyear(ee)

     pc1 = pctsS
     pc1!0 = "time"
     pc1&time = sst&time
     pc1@units = "1"
     sig_pcv = eofunc_north2(evecv@pcvar,dimsizes(pc1),False)
     if (sig_pcv(0)) then    ; if True then significant
        ipo@pcvar = tofloat(sprintf("%4.1f", evecv@pcvar(0)))+"%*"
     else
        ipo@pcvar = tofloat(sprintf("%4.1f", evecv@pcvar(0)))+"%"
     end if
     delete([/sig_pcv,evecv/])
    
     if (.not.ismissing(ipo({35},{160}))) then
        if (ipo({35},{160}).ge.0) then   ; arbitrary attempt to make all plots have the same sign..
           ipo = ipo*-1.
           pc1 = pc1*-1.
        end if
     end if
     delete([/sst,pctsS/])
;---------------------------------------------------------------------------------------------     
     if (OUTPUT_DATA.eq."True") then
        modname = str_sub_str(names(ee)," ","_")
        bc = (/"/","'","(",")"/)
        do gg = 0,dimsizes(bc)-1
           modname = str_sub_str(modname,bc(gg),"_")
        end do
        fn = getenv("OUTDIR")+modname+".cvdp_data.ipo."+syear(ee)+"-"+eyear(ee)+".nc"
        if (.not.isfilepresent2(fn)) then
           z = addfile(fn,"c")
           z@source = "NCAR Climate Analysis Section's Climate Variability Diagnostics Package v"+getenv("VERSION")
           z@notes = "Data from "+names(ee)+" from "+syear(ee)+"-"+eyear(ee)
           if (OPT_CLIMO.eq."Full") then
              z@climatology = syear(ee)+"-"+eyear(ee)+" climatology removed prior to all calculations (other than means)"
           else
              z@climatology = CLIMO_SYEAR+"-"+CLIMO_EYEAR+" climatology removed prior to all calculations (other than means)"
           end if
           z@Conventions = "CF-1.6"
        else
           z = addfile(fn,"w")
        end if
        z->ipo_pattern_mon = set_varAtts(ipo,"IPO spatial pattern (monthly)","","")
        z->ipo_timeseries_mon = set_varAtts(pc1,"IPO normalized principal component timeseries (monthly)","1","")
        delete([/modname,fn/])
     end if              
;------------------------------------------------------------------------
     iopt = 0
     jave = (3*nyr(ee))/100
     val1 = .95
     val2 = .99
     if (jave.eq.0) then
       jave = 1
     end if
     pct = 0.1    
     spectra_mvf = False        ; missing value flag for nino3.4
     spectra_mvf_obs = True     ; missing value flag for observational nino3.4 default
     if (any(ismissing(pc1))) then
        print("Missing data detected for "+names(ee)+", power spectra function does not allow missing data, not creating IPO spectra")
        spectra_mvf = True
        if (isfilepresent2("obs_ts").and.ee.eq.0) then
           spectra_mvf_obs = True    ; missing value flag for obs nino3.4
        end if
     else
        if (isfilepresent2("obs_ts").and.ee.eq.0) then
           spectra_mvf_obs = False   ; missing value flag for obs nino3.4
        end if 
        sdof = specx_anal(pc1,iopt,jave,pct)   ; pc1 already standardized
        splt1 = specx_ci(sdof,val1,val2)
        if (OUTPUT_DATA.eq."True") then
           splt1!0 = "ncurves"
           splt1&ncurves = ispan(0,3,1)
           splt1&ncurves@long_name = "power spectra curves"
           splt1&ncurves@units = "1"
           splt1!1 = "frequency"
           splt1&frequency = sdof@frq
           splt1&frequency@units = "1"
           splt1@units_info = "df refers to frequency interval; data are standardized so there are no physical units"
           splt1@units = "1/df"
           splt1@info = "(0,:)=spectrum,(1,:)=Markov red noise spectrum, (2,:)="+val1+"% confidence bound for Markhov, (3,:)="+val2+"% confidence bound for Markhov"
           z->ipo_spectra = set_varAtts(splt1,"IPO (monthly) power spectra, Markov spectrum and confidence curves","","") 
        end if
        if (isfilepresent2("obs_ts").and.ee.eq.0) then
           sdof_obs = sdof
        end if
        delete([/iopt,jave,pct/])
     end if
     if (isvar("z")) then
        delete(z)
     end if
;========================================================================
     res = True
     res@mpProjection = "WinkelTripel"
     res@mpGeophysicalLineColor = "gray42"
     res@mpPerimOn    = False
     res@mpGridLatSpacingF =  90            ; change latitude  line spacing
     res@mpGridLonSpacingF = 180.           ; change longitude line spacing
     res@mpGridLineColor   = "transparent"  ; trick ncl into drawing perimeter
     res@mpGridAndLimbOn   = True           ; turn on lat/lon lines  
     res@mpFillOn = False
     res@mpCenterLonF = 210.
     res@mpOutlineOn = True  
     res@gsnDraw      = False
     res@gsnFrame     = False
     res@vpYF = 0.95
     res@vpHeightF = 0.3
     res@vpXF = 0.2
     res@vpWidthF = 0.6
  
;     res@cnFillMode = "RasterFill"
     res@cnLevelSelectionMode = "ExplicitLevels"
     
     if (COLORMAP.eq."0") then
        res@cnLevels = fspan(-.65,.65,27)
     end if
     if (COLORMAP.eq."1") then
        res@cnLevels = fspan(-.8,.8,17)
     end if

     res@cnLineLabelsOn = False
     res@cnFillOn        = True
     res@cnLinesOn       = False
     res@lbLabelBarOn    = False

     res@gsnLeftStringOrthogonalPosF = -0.05
     res@gsnLeftStringParallelPosF = .005
     res@gsnRightStringOrthogonalPosF = -0.05
     res@gsnRightStringParallelPosF = 0.96
     res@gsnRightString = ""
     res@gsnLeftString = ""
     res@gsnLeftStringFontHeightF = 0.014
     res@gsnCenterStringFontHeightF = 0.018
     res@gsnRightStringFontHeightF = 0.014
     res@gsnLeftString = syear(ee)+"-"+eyear(ee)
 
     res@gsnRightString = ipo@pcvar
     res@gsnCenterString = names(ee)
     map(ee) = gsn_csm_contour_map(wks,ipo,res)
     delete([/ipo/])
     
     pres = True
     pres@vpXF = 0.07
     pres@trYMinF = 0.
     pres@trXMinF = 0.0
     if (isfilepresent2("obs_ts").and.ee.ge.1.and.spectra_mvf_obs.eq.False) then   
        if (max((/max(splt1(0,:)),max(sdof_obs@spcx)/)).lt.235) then
           pres@trYMaxF = 250.
        else
           pres@trYMaxF = max((/max(splt1(0,:)),max(sdof_obs@spcx)/))+(.1*max((/max(splt1(0,:)),max(sdof_obs@spcx)/)))
        end if
     else
        if (max(splt1(0,:)).lt.235) then
           pres@trYMaxF = 250.
        else
           pres@trYMaxF = max(splt1(0,:))+(.1*max(splt1(0,:)))
        end if
     end if
     pres@trXMaxF = 0.0832
     pres@tiYAxisString = "Power"              ; yaxis
     pres@xyLineColor   =  "black"         
     pres@gsnFrame      = False
     pres@gsnDraw       = False
  
     pres@tmXBLabelDeltaF = -.8
     pres@tmXTLabelDeltaF = -.8
     pres@pmLegendDisplayMode    = "Never"
     pres@xyLineThicknesses   = (/3.5,2.,1.,1./) 
     pres@xyDashPatterns      = (/0,0,0,0/)
     pres@xyLineColors        = (/"foreground","red","blue","green"/)      
     pres@xyLabelMode = "custom"
     pres@xyLineLabelFontColors = pres@xyLineColors 
     pres@xyExplicitLabels = (/"","",val1*100+"%",val2*100+"%"/)
     pres@tmXTOn = True
     pres@tmYROn = False
     pres@tmXTLabelsOn = True
     pres@tmXUseBottom = False
     pres@tmXTMode   = "Explicit"  
     pres@tmXBMode   = "Explicit"            
     pres@tmXTValues = (/".00167",".00833",".01667",".02778",".0416",".0556",".0832"/)
     pres@tmXTLabels = (/"50","10","5","3","2","1.5","1"/)           
     pres@tmXBValues = (/".0",".01",".02",".03",".042",".056",".083"/)
     pres@tmXBLabels =     pres@tmXBValues
     pres@tmXTLabelFontHeightF = 0.018
     pres@tmXBLabelFontHeightF = 0.018
     pres@tmYLLabelFontHeightF = 0.018
     pres@tiYAxisString = "Variance" ;"Power (~S~o~N~C~S~2~N~ / cycles mo~S~-1~N~)"              ; yaxis
     pres@tiXAxisString = "Frequency (cycles mo~S~-1~N~)"
     pres@tiMainString = ""
     pres@txFontHeightF = 0.015
     pres@xyLineLabelFontHeightF = 0.022
     pres@tiXAxisFontHeightF = 0.025
     pres@tiYAxisFontHeightF = 0.025
     pres@tiMainFontHeightF = 0.03
     pres@gsnRightStringOrthogonalPosF = -0.115
     
     pres@tiMainOn = False
     pres@gsnCenterString = "Period (years)"
     pres@gsnCenterStringFontHeightF =    pres@tiYAxisFontHeightF
     pres@gsnRightStringFontHeightF  =    pres@tiYAxisFontHeightF - 0.005
     pres@gsnRightString = syear(ee)+"-"+eyear(ee)+"  "
     pres@gsnLeftString = "" 
     if (wks_type.eq."png") then
        pres@xyLineThicknessF   = 3.5
        res@mpGeophysicalLineThicknessF = 2.  
     else
        pres@xyLineThicknessF   = 1.5
        res@mpGeophysicalLineThicknessF = 1.  
     end if
     pres@gsnCenterString = names(ee)
     if (spectra_mvf.eq.False) then
        pspec(ee) = gsn_csm_xy(wks2,sdof@frq,splt1,pres) 
     
        if (isfilepresent2("obs_ts").and.ee.ge.1.and.spectra_mvf_obs.eq.False) then
           pres@xyLineColors        = (/"gray70","black","black","black"/)
           pres@xyCurveDrawOrder = "PreDraw"
           pres@gsnCenterString = ""
           pres@gsnRightString = ""
           pspec_obs(ee) = gsn_csm_xy(wks2,sdof_obs@frq,sdof_obs@spcx,pres)
           overlay(pspec(ee),pspec_obs(ee)) 
           delete(pres@xyCurveDrawOrder)
        end if     
        delete([/sdof,splt1/])
     end if 
     
     xyres = True
     xyres@gsnDraw = False
     xyres@gsnFrame = False
     xyres@gsnRightString = ""
     xyres@gsnLeftString = ""
     xyres@gsnFrame = False
     xyres@gsnYRefLine = 0.0
     xyres@gsnYRefLineColor = "gray42"
     xyres@gsnXYBarChart = False
     xyres@gsnAboveYRefLineColor = 185
     xyres@gsnBelowYRefLineColor = 35
     xyres@xyLineThicknessF = 0.1
     xyres@xyLineColor = "gray70"
     xyres@tiYAxisString = ""
     if (nsim.le.5) then
        xyres@tmXBLabelFontHeightF = 0.0125
        xyres@tmYLLabelFontHeightF = 0.0125
        xyres@gsnStringFontHeightF = 0.017     
     else
        xyres@tmXBLabelFontHeightF = 0.018
        xyres@tmYLLabelFontHeightF = 0.018
        xyres@gsnStringFontHeightF = 0.024
     end if
     xyres@vpXF = 0.05
     xyres@vpHeightF = 0.15
     if (SCALE_TIMESERIES.eq."True") then
        xyres@vpWidthF = 0.9*((nyr(ee)*1.)/nyr_max)
     else
        xyres@vpWidthF = 0.9
     end if
     xyres@gsnCenterString = ""
     
     xyres@trXMinF = syear(ee)-.5
     xyres@trXMaxF = eyear(ee)+1.5
     
;     xyres2 =    xyres
;     delete(xyres2@gsnXYBarChart)
;     delete(xyres2@gsnAboveYRefLineColor)
;     delete(xyres2@gsnBelowYRefLineColor)
;     xyres2@xyLineColor = "black"
;     if (wks_type.eq."png") then
;        xyres2@xyLineThicknessF = 3.5
;     else
;        xyres2@xyLineThicknessF = 2.5
;     end if

     xyres@gsnCenterString = names(ee)
     xyplot(ee) = gsn_csm_xy(wks3,fspan(syear(ee),eyear(ee)+.91667,dimsizes(pc1)),pc1,xyres)   ; use standardized timeseries  
;     xyplot2(ee) = gsn_csm_xy(wks3,fspan(syear(ee),eyear(ee)+.91667,dimsizes(pc1)),runave(pc1,61,0),xyres2)
;     overlay(xyplot(ee),xyplot2(ee))
     
     delete([/val1,val2,pc1/])
  end do

  panres = True
  panres@gsnMaximize = True
  panres@gsnPaperOrientation = "portrait"
  panres@gsnPanelLabelBar = True
  panres@gsnPanelYWhiteSpacePercent = 3.0
  panres@pmLabelBarHeightF = 0.05
  panres@pmLabelBarWidthF = 0.55
  panres@lbTitleOn = False
  panres@lbBoxLineColor = "gray70"
  if (nsim.le.4) then
     if (nsim.eq.1) then
        panres@txFontHeightF = 0.022
        panres@gsnPanelBottom = 0.50
     else
        panres@txFontHeightF = 0.0145
        panres@gsnPanelBottom = 0.50
     end if
  else
     panres@txFontHeightF = 0.016
     panres@gsnPanelBottom = 0.05
  end if
  panres@txString = "IPO (Monthly)"
  ncol = floattointeger(sqrt(nsim))
  nrow = (nsim/ncol)+mod(nsim,ncol)  
  gsn_panel2(wks,map,(/nrow,ncol/),panres)
  delete(wks)

  delete(panres@gsnPanelLabelBar)
  panres@txString = "IPO (Monthly)"
  gsn_panel2(wks2,pspec,(/nrow,ncol/),panres)  
  delete(wks2)
  
  if (SCALE_TIMESERIES.eq."True") then
     tt = ind(nyr.eq.nyr_max)
     panres@gsnPanelScalePlotIndex = tt(0)
     delete(tt)
  end if
  if (nsim.le.12) then
     lp = (/nsim,1/)
  else
     lp = (/nrow,ncol/)  ;(/nsim/2+1,nsim/8+1/)  
  end if
  panres@txString = "IPO (Monthly)"
  gsn_panel2(wks3,xyplot,lp,panres)  
  delete(wks3)
  delete([/map,pspec,syear,eyear,nyr,nyr_max,lp/])
  print("Finished: ipo.ncl")
end