Opened 4 years ago

Closed 4 years ago

# zdfric: momentum background viscosity propagates into tracer diffusivity

Reported by:reporter: Owned by:owner: The author of the ticket. Person in charge for handling mclaus nemo normal 2015 release-3.6 OCE trunk release-3.6 zdfric

# 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)  )   &
!                                            ! 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.

### Commit History (2)

ChangesetTimeChangeLog
7049emanuelaclementi2016-10-20T11:15:23+02:00

#1768 Bug fixed in zdfric.F90 in NEMO_v3.6_Stable

7048emanuelaclementi2016-10-20T11:12:30+02:00

#1768 Bug fixed in zdfric.F90 in trunk

### Changed 4 years ago by mclaus

Patch file for zdfric.F90

### comment:1 follow-up: ↓ 5 Changed 4 years ago by 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?

### comment:2 Changed 4 years ago by 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) ) &

! ! Add the background coefficient on eddy viscosity

END DO

END DO

### comment:3 Changed 4 years ago by 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

### comment:4 Changed 4 years ago by 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

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

### comment:5 in reply to: ↑ 1 Changed 4 years ago by mclaus

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.

### comment:6 Changed 4 years ago by emanuelaclementi

• Milestone set to 2015 nemo_v3_6_STABLE

### comment:7 Changed 4 years ago by emanuelaclementi

• Resolution set to fixed
• Status changed from new to closed

### comment:8 Changed 4 years ago by emanuelaclementi

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