With absolute vorticity q_2 = f+\zeta and stretching \partial_z (f/N^2)b, their combined horizontal gradients and magnitude are found on the T point

```
!== calculate gradients of vorticity and stretching, then square of magnitude (t-point) ==!
DO jk = 1, jpkm1
DO jj = 2, jpjm1
DO ji = 2, jpim1
zvsdx = r1_2 * ( r1_e1f(ji,jj) * ( zst(ji,jj-1,jk) + zwz(ji,jj-1,jk) - &
& zst(ji-1,jj-1,jk) - zwz(ji-1,jj-1,jk) ) &
& + r1_e1f(ji,jj) * ( zst(ji,jj,jk) + zwz(ji,jj,jk) - &
& zst(ji-1,jj,jk) - zwz(ji-1,jj,jk) ) )
zvsdy = r1_2 * ( r1_e2f(ji,jj) * ( zst(ji-1,jj,jk) + zwz(ji-1,jj,jk) - &
& zst(ji-1,jj-1,jk) - zwz(ji-1,jj-1,jk) ) &
& + r1_e2f(ji,jj) * ( zst(ji,jj,jk) + zwz(ji,jj,jk) - &
& zst(ji,jj-1,jk) - zwz(ji,jj-1,jk) ) )
dvsmagsq(ji,jj,jk) = ( zvsdx**2 + zvsdy**2 ) * tmask(ji,jj,jk)
END DO
END DO
END DO
!
CALL lbc_lnk_multi( 'ldfdyn', dvsmagsq, 'T', 1. )
```

and `dvsmagsq`

is applied through the lateral boundary condition.

Then, gradients of divergence and its magnitude are found on the F point. Horizontal divergence is calculated in the module `divhor.F90`

as `hdivn`

```
!== calculate gradients of divergence, then square of magnitude (f-point) ==!
DO jk = 1, jpkm1
DO jj = 1, jpjm1
DO ji = 1, jpim1
zdivdx = r1_2 * ( r1_e1t(ji,jj) * ( hdivn(ji+1,jj+1,jk) - hdivn(ji,jj+1,jk) ) &
& + r1_e1t(ji,jj) * ( hdivn(ji+1,jj,jk) - hdivn(ji,jj,jk) ) )
zdivdy = r1_2 * ( r1_e2t(ji,jj) * ( hdivn(ji+1,jj+1,jk) - hdivn(ji+1,jj,jk) ) &
& + r1_e2t(ji,jj) * ( hdivn(ji,jj+1,jk) - hdivn(ji,jj,jk) ) )
ddivmagsq(ji,jj,jk) = ( zdivdx**2 + zdivdy**2 ) * fmask(ji,jj,jk)
END DO
END DO
END DO
```

We now have nearly everything to calculate the viscous coefficient A_{qg}, but first we find (\Lambda/\pi)^6 (n.b. this is later square rooted)

```
zcmqgl = (rn_cqgc/rpi)**6 ! (C_qg/pi)^6
```

The viscous coefficients are found

```
DO jk = 1, jpkm1 !== QG Leith viscosity coefficient on T-point ==!
DO jj = 2, jpjm1
DO ji = fs_2, fs_jpim1 ! vector opt.
zsqQG = r1_4 * ( ddivmag(ji,jj,jk) + ddivmag(ji-1,jj,jk) + ddivmag(ji,jj-1,jk) + &
& ddivmag(ji-1,jj-1,jk) ) + dvsmag(ji,jj,jk)
amht(ji,jj,jk) = SQRT( zcmqgl * esqt(ji,jj)**3 * zsqQG )
END DO
END DO
END DO
!
DO jk = 1, jpkm1 !== QG Leith viscosity coefficient on F-point ==!
DO jj = 1, jpjm1
DO ji = 1, fs_jpim1 ! vector opt.
zsqQG = r1_4 * ( dvsmag(ji,jj,jk) + dvsmag(ji+1,jj,jk) + dvsmag(ji,jj+1,jk) + &
& dvsmag(ji+1,jj+1,jk) ) + ddivmag(ji,jj,jk)
amhf(ji,jj,jk) = SQRT( zcmqgl * esqf(ji,jj)**3 * zsqQG )
END DO
END DO
END DO
```

And then `amht`

and `amhf`

are passed through `lbc_lnk()`

.

I believe this to be almost complete. As always, I would appreciate any further tips if possible e.g. optimisation. For example, is it worth calculating Ro and Fr numbers every n time steps?