! 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