Hello,
We have written the QG Leith viscosity parameterisation into NEMO v4.0.4,
and have got idealised and global forced configurations running without errors. However, there is a random pattern of large viscosity coefficients dotted about the domain, which we are struggling to get rid of (see figure below).
The namelist_cfg is found here. The domain is 3000 km deep and has a bump at 3000 km in x that is 1500 m in height.
We have tried many fixes:
- Harmonic diffusion,
- FCT and CEN advection schemes,
- EVD and NPC convection schemes,
- Volume varying vertical grid thickness,
- Increasing horizontal diffusion by an order of magnitude,
- Increasing the viscosity coefficient parameter, \Lambda
Further documentation of the model development is found here.
We know the noise pattern comes from the QG Leith stretching term, which is the term in the equation above with the buoyancy frequency, N^2, and buoyancy gradients.
The code that computes the stretching term is below, and the full ldfdyn.F90 is found here
zqglep1 = 1.e-12_wp
zqglep2 = 1.e-24_wp
!
!== begin calculation of stretching term d/dz[(f/(N**2))*grad(b)] ==!
!== find buoyancy and interpolate onto w-grid ==!
DO jk = 1, jpkm1
DO jj = 1, jpj
DO ji = 1, jpi
IF( jk < 2 ) THEN
!== buoyancy at surface ==!
zbu(ji,jj,jk) = - grav * prd(ji,jj,jk)
ELSE
!== buoyancy below surface ==!
zbuup = - grav * prd(ji,jj,jk-1)
zbulw = - grav * prd(ji,jj,jk )
zbu(ji,jj,jk) = 0.5_wp * ( zbuup + zbulw ) * wmask(ji,jj,jk)
ENDIF
END DO
END DO
END DO
!
!== Calculate horizontal gradients of buoyancy and put on w-grid ==!
DO jk = 1, jpkm1
DO jj = 1, jpjm1
DO ji = 1, jpim1
!== gradients of buoyancy on U and V grid at w point of cell ==!
zbudxup(ji,jj,jk) = r1_e1u(ji,jj) * ( zbu(ji+1,jj,jk) - zbu(ji,jj,jk) ) * umask(ji,jj,jk)
zbudyvp(ji,jj,jk) = r1_e2v(ji,jj) * ( zbu(ji,jj+1,jk) - zbu(ji,jj,jk) ) * vmask(ji,jj,jk)
END DO
END DO
END DO
!
CALL lbc_lnk_multi( 'ldfdyn', zbudxup, 'U', 1., zbudyvp, 'V', 1. )
!
!== gradients of buoyancy on W- points ==!
DO jk = 1, jpkm1
DO jj = 2, jpjm1
DO ji = 2, jpim1
zbudx(ji,jj,jk) = r1_2 * ( zbudxup(ji-1,jj,jk) + zbudxup(ji,jj,jk) ) * wmask(ji,jj,jk)
zbudy(ji,jj,jk) = r1_2 * ( zbudyvp(ji,jj-1,jk) + zbudyvp(ji,jj,jk) ) * wmask(ji,jj,jk)
END DO
END DO
END DO
!
CALL lbc_lnk_multi( 'ldfdyn', zbudx, 'T', 1., zbudy, 'T', 1. )
!
!== take vertical gradient and find stretching d/dz[(f * grad(b))/N^2] (t-point) ==!
DO jk = 1, jpkm1
DO jj = 1, jpj
DO ji = 1, jpi
!== are we below the mixed layer and above the sea floor? ==!
IF( jk > nmlnqg(ji,jj) .AND. jk < mbkt(ji,jj) ) THEN
!== vertical gradient of x component ==!
zker1 = ( ff_t(ji,jj) * zbudx(ji,jj,jk ) ) / MAX( pn2(ji,jj,jk ), zqglep1 )
zker2 = ( ff_t(ji,jj) * zbudx(ji,jj,jk+1) ) / MAX( pn2(ji,jj,jk+1), zqglep1 )
zstx(ji,jj,jk) = ( ( zker1 - zker2 ) / e3t_b(ji,jj,jk) ) * tmask(ji,jj,jk)
!== vertical gradient of y component ==!
zker1 = ( ff_t(ji,jj) * zbudy(ji,jj,jk ) ) / MAX( pn2(ji,jj,jk ), zqglep1 )
zker2 = ( ff_t(ji,jj) * zbudy(ji,jj,jk+1) ) / MAX( pn2(ji,jj,jk+1), zqglep1 )
zsty(ji,jj,jk) = ( ( zker1 - zker2 ) / e3t_b(ji,jj,jk) ) * tmask(ji,jj,jk)
ENDIF
END DO
END DO
END DO
!
!== calculate the Burger number and square of Rossby number on t-point ==!
!== calculate over entire domain for diagnostics ==!
DO jk = 1, jpkm1
DO jj = 2, jpj
DO ji = 2, jpi
!== grid scale velocity squared ==!
zusq = 0.5_wp * ( ( ub(ji-1,jj ,jk) * ub(ji-1,jj ,jk) + ub(ji,jj,jk) * ub(ji,jj,jk) ) + &
& ( vb(ji ,jj-1,jk) * vb(ji ,jj-1,jk) + vb(ji,jj,jk) * vb(ji,jj,jk) ) )
!== square of Rossby number U^2/(f^2 * A) ==!
rro2(ji,jj,jk) = ( zusq / ( MAX( ff_t(ji,jj)**2, zqglep2 ) * esqt(ji,jj) ) ) * tmask(ji,jj,jk)
!== averaging square of buoyancy frequency onto t-grid ==!
IF( jk < mbkt(ji,jj) ) THEN
!== accounting for negative N^2 ==!
znsq = r1_2 * ( MAX( pn2(ji,jj,jk), zqglep1 ) + MAX( pn2(ji,jj,jk+1), zqglep1 ) ) * tmask(ji,jj,jk)
ELSE
!== stratification is continuous at bottom ==!
znsq = MAX( pn2(ji,jj,jk ), zqglep1 )
ENDIF
!== Burger number (N^2 * delta_z^2)/(f^2 * A) ==!
rbu(ji,jj,jk) = ( MAX( znsq , zqglep1 ) * e3t_b(ji,jj,jk)**2 ) / &
& ( MAX( ff_t(ji,jj)**2, zqglep2 ) * esqt(ji,jj) )
!== Froude number squared (Fr^2 = Ro^2/Bu) ==!
rfr2(ji,jj,jk) = rro2(ji,jj,jk)/rbu(ji,jj,jk)
END DO
END DO
END DO
!
CALL lbc_lnk_multi( 'ldfdyn', rro2, 'T', 1., rbu, 'T', 1. )
!
!== are we in the QG limit? Find the stretching value in x and y components ==!
DO jk = 1, jpkm1
DO jj = 1, jpj
DO ji = 1, jpi
!== are we below the mixed layer and above the sea floor? ==!
IF( jk > nmlnqg(ji,jj) .AND. jk < mbkt(ji,jj) ) THEN
!== x component of stretching ==!
zztmpx = MIN( ABS( zstx(ji,jj,jk) ), &
& ABS( ( zwzdx(ji,jj,jk) * rfr2(ji,jj,jk) ) / &
& ( rro2(ji,jj,jk) + rfr2(ji,jj,jk)**2 + zqglep2 ) ) )
zstlimx(ji,jj,jk) = SIGN( zztmpx, zstx(ji,jj,jk) )
!== y component of stretching ==!
zztmpy = MIN( ABS( zsty(ji,jj,jk) ), &
& ABS( ( zwzdy(ji,jj,jk) * rfr2(ji,jj,jk) ) / &
& ( rro2(ji,jj,jk) + rfr2(ji,jj,jk)**2 + zqglep2 ) ) )
zstlimy(ji,jj,jk) = SIGN( zztmpy, zsty(ji,jj,jk) )
ENDIF
END DO
END DO
END DO
Currently we do not know how to fix this. Any suggestions most welcome. Thanks in advance.