; ==============================================================
; NCL script to compute POP MOC field offline from POP netcdf
; history files.  This routine is designed for the CCSM4 ocean 
; component.
; ==============================================================
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 "/work/yiwen89419/CCSM_ncl/ncl/yeager_util.ncl"

begin
; ==============================================================
; Input/Output
; ==============================================================
; datadir = "/ccsm/ocn/yeager/amoc_ts/b40.1850.track1.1deg.006/"
; fin = datadir+"b40.1850.track1.1deg.006.pop.h.LABrho.0700_cat_1300.nc"
; sfx = get_file_suffix(fin,0)
; fout  = sfx@fBase+".volavg.nc"

  datadir = "/work/yiwen89419/CCSM_ncl/"
; fin = datadir+"Yearly/TSRyr_b_20ct_ctrl_1.pop.h.192101-200012.rho_TS.nc"
  fin = datadir+"f09.B-PI.tn15.cmip6.j01.re/yearly/60yrs/TSRyr_taiesm_501_560.rho_TS.nc"
  sfx = get_file_suffix(fin,0)
  fout  = sfx@fBase+".volavg.nc"

; ==============================================================
; Get variables and reduce and shift
; ==============================================================
  f  = addfile (fin, "r")
  time=f->time
  tvol=f->tvolume
  rho_on=f->rho_on
  rho_off=f->rho_off
  delr_on=f->delr_on
  delr_off=f->delr_off
  delr_tpluss=f->delr_tpluss
  rhot=f->rho_T
  rhos=f->rho_S

  tvolfull = conform(delr_off,tvol,(/1,2,3/))
  delr_on_vavg = dim_sum(dim_sum(dim_sum(delr_on*tvolfull)))/dim_sum(dim_sum(dim_sum(tvolfull)))
  delr_off_vavg = dim_sum(dim_sum(dim_sum(delr_off*tvolfull)))/dim_sum(dim_sum(dim_sum(tvolfull)))
  delr_tpluss_vavg = dim_sum(dim_sum(dim_sum(delr_tpluss*tvolfull)))/dim_sum(dim_sum(dim_sum(tvolfull)))
  rhot_vavg = dim_sum(dim_sum(dim_sum(rhot*tvolfull)))/dim_sum(dim_sum(dim_sum(tvolfull)))
  rhos_vavg = dim_sum(dim_sum(dim_sum(rhos*tvolfull)))/dim_sum(dim_sum(dim_sum(tvolfull)))
  

; ==============================================================
; Save to netcdf
; ==============================================================
  system("/bin/rm -f "+fout)
  outf = addfile(fout,"c")
  filedimdef(outf, "time",-1,True)    ; make "time" UNLIMITED
  outf->time=time
  outf->delr_on_vavg=delr_on_vavg
  outf->delr_off_vavg=delr_off_vavg
  outf->delr_tpluss_vavg=delr_tpluss_vavg
  outf->rhot_vavg=rhot_vavg
  outf->rhos_vavg=rhos_vavg
end