! revision of Stuart's script for mass conservation<br>
! modernize to use @shf instead of fixed i,j values<br>
! ignores partial bottom cells/BBL<br>
<br>
use subdomain_notides-mean.50s-20n.cdf<br>
use subdomain_notides-mean-w.50s-20n.cdf<br>
<br>
! constants and such<br>
let radius = 6370000.		! may be different in ACOM3?<br>
let pi = 4*atan(1)<br>
let rad = 180/pi<br>
<br>
let uu = u_30d[d=subdomain_notides-mean.50s-20n]/100 ! convert to m/s</br>
let vv = v_30d[d=subdomain_notides-mean.50s-20n]/100<br>
<br>
! gridbox lengths in radians<br>
let dyt = ybox[g=temp_30d[d=subdomain_notides-mean.50s-20n]]/rad<br>
let dxt = xbox[g=temp_30d[d=subdomain_notides-mean.50s-20n]]/rad<br>
<br>
! cos(latitude) at two j's used for differencing<br>
! name yy to avoid transforming pseudo-variables  (should be g=u_30d?)<br>
let yy = y[g=temp_30d[d=subdomain_notides-mean.50s-20n]]<br>
let costheta1 = cos(yy/rad)<br>
let costheta2 = cos(yy[j=@shf:-1]/rad)<br>
<br>
! u_x. Delta-x is 1 gridbox, but average differences over latitudes j and j-1.<br>
let ux = (uu + uu[j=@shf:-1] - uu[i=@shf:-1] - uu[i=@shf:-1,  j=@shf:-1])/(2.*dxt)<br>
! v_y. Delta-y is 1 gridbox, but average differences over longitudes i and i-1.<br>
! multiply each latitude by cos(theta) (will divide by the mean of these later).<br>
let vy = ((vv + vv[i=@shf:-1])*costheta1 - (vv[j=@shf:-1] + vv[i=@shf:-1, j=@shf:-1])*costheta2)/(2.*dyt)<br>
<br>
! horizontal divergence (indefinite integral in z)<br>
! factor of 2 is to take average of the 2 cos(theta)<br>
let horadv = 2.*(ux+vy)/(radius*(costheta1+costheta2))<br>
let horadvin = horadv[z=@iin]<br>
let horgrid = horadvin[gz=w_30d[d=subdomain_notides-mean-w.50s-20n]@asn]<br>
<br>
! magnitude of horizontal divergence<br>
let horma = 2.*(ux*ux + vy*vy)^0.5/(radius*(costheta1+costheta2))<br>
let hormax  =horma[z=@iin]<br>
<br>
! w (m/s)<br>
let vertadv = w_30d[d=subdomain_notides-mean-w.50s-20n]/100<br>
<br>
! compare Int(u_x+v_y)dz with w<br>
let vertdiff = (vertadv - horgrid)<br>