zdfric: momentum background viscosity propagates into tracer diffusivity

[mclaus]

Context

I think, the following code of zdfric.F90 (line 158, trunk@6489) produces to large values of vertical eddy viscosity for the tracers.

    DO jj = 2, jpjm1                                   ! Eddy diffusivity coefficients (avt)
       DO ji = fs_2, fs_jpim1
          avt(ji,jj,jk) = tmric(ji,jj,jk) / ( 1._wp + rn_alp * zwx(ji,jj) )           &
             &                            * (  avmu(ji,jj,jk) + avmu(ji-1,jj,jk)      &
             &                               + avmv(ji,jj,jk) + avmv(ji,jj-1,jk)  )   &
             &          + avtb(jk) * tmask(ji,jj,jk)
          !                                            ! Add the background coefficient on eddy viscosity
          avmu(ji,jj,jk) = avmu(ji,jj,jk) + avmb(jk) * umask(ji,jj,jk)
          avmv(ji,jj,jk) = avmv(ji,jj,jk) + avmb(jk) * vmask(ji,jj,jk)
       END DO
    END DO

Analysis

The problem is that within the same loop body, the vertical eddy viscosity for the tracer is computed using a
four point average of the eddy viscosity for the zonal and meridional momentum and the background eddy viscosity
is added to the viscosity of the momentum. This way, avmu(ji-1,jj,jk) has the background viscosity added but
avmu(ji,jj,jk) has not. As a result, the four point average is larger by about 50% of the momentums' background viscosity.

As a side effect, there is a memory dependence in the loop preventing the compiler to optimize the loop.

Affected versions are 3.4, 3.6 and trunk and possibly older version of NEMO

Fix

I suggest to split the loop into two nested loops. A patch file is attached.

[mclaus]

Patch file for zdfric.F90

[timgraham]

I agree that this is incorrect and should be fixed.

Luckily very few configurations use this option (I don't know of any). Are you using it in your own configurations?

[emanuelaclementi]

I'm using the Richardson parameterization for the Mediterranean configuration.
I also agree that it must be corrected.
But according to Pacanowski & Philander JPO 1981 eq. 1 & 2, the background eddy viscosity should be added to the eddy viscosity for the zonal and meridional momentum before the loop on vertical eddy viscosity for the tracers.

z05alp = 0.5_wp * rn_alp
DO jj = 1, jpjm1 ! Eddy viscosity coefficients (avm)

DO ji = 1, jpim1

avmu(ji,jj,jk) = (rn_avmri / ( 1. + z05alp*(zwx(ji+1,jj)+zwx(ji,jj) ) )nn_ric &

& + avmb(jk) ) * umask(ji,jj,jk)

avmv(ji,jj,jk) = (rn_avmri / ( 1. + z05alp*(zwx(ji+1,jj)+zwx(ji,jj) ) )nn_ric &

& + avmb(jk) ) * vmask(ji,jj,jk)

END DO

END DO

DO jj = 2, jpjm1 ! Eddy diffusivity coefficients (avt)

DO ji = 2, jpim1

avt(ji,jj,jk) = tmric(ji,jj,jk) / ( 1._wp + rn_alp * zwx(ji,jj) ) &

& * ( avmu(ji,jj,jk) + avmu(ji-1,jj,jk) &
& + avmv(ji,jj,jk) + avmv(ji,jj-1,jk) ) &
& + avtb(jk) * tmask(ji,jj,jk)

! ! Add the background coefficient on eddy viscosity

END DO

END DO

[emanuelaclementi]

I've made further investigations and it seems that eq 2 in Pacanowski & Philander JPO 1981 is not correct.
Following Lermusiaux 2001 (Evolving the subspace of the three-dimensional multiscale ocean
variability: Massachusetts Bay, Journal of Marine Systems) it can be deduced that:

avm = avm0 + avmb (where avm0 = rn_avmri/(1+rn_alp*ri)^nn_ric

avt = (avm0/(1+rn_alp*ri)) + avtb

Where:
rn_avmri is the maximum value reaches by avm and avt
avmb and avtb are the background values
rn_alp and nn_ric are adjustable parameters from namelist

So the first proposition from mclaus should be the correct one

[gm]

Clearly the correct solution is just to split the loop in two parts, as it was the case 7 years ago in revision 1537

(https://forge.ipsl.jussieu.fr/nemo/browser/trunk/NEMO/OPA_SRC/ZDF/zdfric.F90?rev=1537)

The error was introduced in rev 1601, probably while thinking to improve the code efficiency, regardless of the physics.

[mclaus]

Replying to timgraham:

I agree that this is incorrect and should be fixed.

Luckily very few configurations use this option (I don't know of any). Are you using it in your own configurations?

I am sorry for the late reply. Yes, I am using the scheme in idealized equatorial basin configurations. My proposed patch does affect the vertical viscosity, however, it does not have a significant impact on my results.

[emanuelaclementi]

Bug fixed
Committed in rev 7048 of the trunk and NEMO_v3.6_table