! revision of Stuart's script for mass conservation
! modernize to use @shf instead of fixed i,j values
! ignores partial bottom cells/BBL

use subdomain_notides-mean.50s-20n.cdf
use subdomain_notides-mean-w.50s-20n.cdf

! constants and such
let radius = 6370000. ! may be different in ACOM3?
let pi = 4*atan(1)

let uu = u_30d[d=subdomain_notides-mean.50s-20n]/100 ! convert to m/s
let vv = v_30d[d=subdomain_notides-mean.50s-20n]/100

! cos(latitude) at two j's used for differencing
! name yy to avoid transforming pseudo-variables (should be g=u_30d?)
let yy = y[g=temp_30d[d=subdomain_notides-mean.50s-20n]]

! u_x. Delta-x is 1 gridbox, but average differences over latitudes j and j-1.
let ux = (uu + uu[j=@shf:-1] - uu[i=@shf:-1] - uu[i=@shf:-1, j=@shf:-1])/(2.*dxt)
! v_y. Delta-y is 1 gridbox, but average differences over longitudes i and i-1.
! multiply each latitude by cos(theta) (will divide by the mean of these later).
let vy = ((vv + vv[i=@shf:-1])*costheta1 - (vv[j=@shf:-1] + vv[i=@shf:-1, j=@shf:-1])*costheta2)/(2.*dyt)

! horizontal divergence (indefinite integral in z)
! factor of 2 is to take average of the 2 cos(theta)