! implement an idea of Stuart's: compare (ux+vy)@iin with w
! use doubled, subsampled grid to get divergence on w grid
! properly accounts for partial bottom cells but not BBL (separate script)

use subdomain_notides-mean.50s-20n.cdf
use subdomain_notides-mean-w.50s-20n.cdf
use acom3-cell-thickness.50s-20n.cdf ! from topog_acom3.1.dta.nc

! make layer-integral transports:
! first step just defines the first 36 values of e435 (u layer thickness in meters)
let dhu = e435[d=acom3-cell-thickness.50s-20n, k=1:36]

! second step makes transport in each layer (factor of 100 converts u from cm/s)
! using these thicknesses takes the partial bottom cell into account
let vdhu = v_30d[d=subdomain_notides-mean.50s-20n]*dhu/100
let udhu = u_30d[d=subdomain_notides-mean.50s-20n]*dhu/100

! double the grid for (u,v)
define axis/x=90e:180:.25/unit=lon xhf
define axis/y=50s:20n:`1/6`/unit=lat yhf
define grid/like=v_30d[d=subdomain_notides-mean.50s-20n]/x=xhf/y=yhf ghf

! transport components on the doubled grid
let udhuhf = udhu[g=ghf]
let vdhuhf = vdhu[g=ghf]

! divergence doubled, then subsampled
let divhf = udhuhf[x=@ddc] + vdhuhf[y=@ddc]
let divhfrg = divhf[g=eta_30d[d=subdomain_notides-mean.50s-20n]]

! integrate divergence down
let divhfrgint = divhfrg[k=@rsum]

! implicitly regrid div to w z-grid (thus lowering the value half a gridbox,
! which is where the sum should accumulate to w).
let diff = w_30d[d=subdomain_notides-mean-w.50s-20n]/100 - divhfrgint