# Noise patterns in QG Leith viscosity field

**URL:** https://nemo-ocean.discourse.group/t/noise-patterns-in-qg-leith-viscosity-field/574
**Category:** v4.0.x
**Created:** [12 October 2023 09:16 UTC](https://nemo-ocean.discourse.group/t/noise-patterns-in-qg-leith-viscosity-field/574 "2023-10-12T09:16:36Z")
**Posts on this page:** 2
**Page:** 1

<div class="post-metadata">

### Author: ![tmwilder](https://yyz2.discourse-cdn.com/free1/user_avatar/nemo-ocean.discourse.group/tmwilder/32/120_2.png) [@tmwilder](https://nemo-ocean.discourse.group/u/tmwilder)
#### Post date: [12 October 2023 09:16 UTC](https://nemo-ocean.discourse.group/t/noise-patterns-in-qg-leith-viscosity-field/574/1 "2023-10-12T09:16:36Z")

</div>

Hello,

We have written the QG Leith viscosity parameterisation into NEMO v4.0.4,

![Screenshot from 2023-10-12 09-51-30](https://global.discourse-cdn.com/free1/uploads/nemo_ocean/original/1X/294f2e2e3b47251691810060065f6c3602bd1d63.png)

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).

 ![visc_qg_64m](https://global.discourse-cdn.com/free1/uploads/nemo_ocean/original/1X/7bf70166ca9eaa2fc11f80178523b719ed03cf6a.png)

The `namelist_cfg` is found [here](https://github.com/thomaswilder/NEMO_activities/blob/main/IDEAL/EXP/namelist_cfg). 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](https://github.com/thomaswilder/NEMO_activities/blob/main/IDEAL/documentation.md).

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](https://github.com/thomaswilder/NEMO_activities/blob/main/QG_Leith/ldfdyn.F90)

```fortran
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.

---

<div class="post-metadata">

### Author: ![acc](https://avatars.discourse-cdn.com/v4/letter/a/b38774/32.png) [@acc](https://nemo-ocean.discourse.group/u/acc)
#### Post date: [17 October 2023 16:46 UTC](https://nemo-ocean.discourse.group/t/noise-patterns-in-qg-leith-viscosity-field/574/2 "2023-10-17T16:46:30Z")

</div>

A couple of possibilities based on a quick scan of the code as shown:

- zbudx and zbudy at jpk are not set but may be used
- rfr2 is not set in the haloes but may be used. Needs to be added to the lbc\_lnk\_multi call.