Bug in zdfric with simplified equation of state

[timgraham]

Description

I think there is a bug in the following section of code but I don't fully understand what it is trying to do.

At line 177 of zdfric.F90 in the trunk:

      zflageos = ( 0.5 + SIGN( 0.5, nn_eos - 1. ) ) * rau0
      DO jj = 2, jpjm1
            DO ji = fs_2, fs_jpim1
            zrhos          = rhop(ji,jj,1) + zflageos * ( 1. - tmask(ji,jj,1) )

Analysis

zflageos = rau0 in the case of the simplified equation of state only. Otherwise it is zero so this bug doesn't affect most model runs

In the last line zflageos will be added to rhop over land (where 1 - tmask = 1) but not over ocean points

I think this has been in the code since NEMO 3.4

Recommendation

Should the last line be changed to

zrhos          = rhop(ji,jj,1) + zflageos * tmask(ji,jj,1)

[gm]

Two things :

1st - rhop is an obsolescent feature of the code that should disappear from v4.0

here using the in situ density at the 1st ocean level is equivalent to consider potential density (i.e. in the code rdn = (rho - rau0) / rau0

2nd here, I guess the issue was not to divide by zero so add a non zero value

There is a more clever way to do so:

zrhos = rau0 * ( 1._wp + rdn(ji,jj,1) )

no need of mask, rdn is masked over land, so zrhos=rau0 over land

then simply divide the calculation in the following lines by zrhos (no need of adding small)

so the following lines 175-186 :

      !  Compute Ekman depth from wind stress forcing.
      ! -------------------------------------------------------
      zflageos = ( 0.5 + SIGN( 0.5, nn_eos - 1. ) ) * rau0
      DO jj = 2, jpjm1
         DO ji = fs_2, fs_jpim1
            zrhos          = rhop(ji,jj,1) + zflageos * ( 1. - tmask(ji,jj,1) )
            zustar         = SQRT( taum(ji,jj) / ( zrhos +  rsmall ) )
            ekm_dep(ji,jj) = rn_ekmfc * zustar / ( ABS( ff(ji,jj) ) + rsmall )
            ekm_dep(ji,jj) = MAX(ekm_dep(ji,jj),rn_mldmin) ! Minimun allowed
            ekm_dep(ji,jj) = MIN(ekm_dep(ji,jj),rn_mldmax) ! Maximum allowed
         END DO
      END DO

becomes:

      !  Compute Ekman depth from wind stress forcing.
      ! -------------------------------------------------------
      DO jj = 2, jpjm1
         DO ji = fs_2, fs_jpim1
            zustar = SQRT( taum(ji,jj) * r1_rau0 / ( 1._wp +rdn(ji,jj,1) )
            ekm_dep(ji,jj) = rn_ekmfc * zustar / ( ABS( ff(ji,jj) ) + rsmall )
            ekm_dep(ji,jj) = MAX( ekm_dep(ji,jj) , rn_mldmin )   ! Minimun allowed
            ekm_dep(ji,jj) = MIN( ekm_dep(ji,jj) , rn_mldmax )   ! Maximum allowed
         END DO
      END DO

therefore no need a reference to the EOS used, no more need to used zrhos scalar.

Gurvan

[gm]

Some additional remarks concerning zdfric and the minimum Kz in the Ekman layer:

• First :

The Ekman layer depth is an approximation, there is no need to use the local density here. The reference density is really sufficient.
Therefore, I suggest to remove the " USE eosbn2 " in the header of the module and use the following piece of code :

      INTEGER  ::   ji, jj, jk                 ! dummy loop indices
      REAL(wp) ::   zcfRi, zav, zustar, zhek   ! local scalars
      REAL(wp), DIMENSION(jpi,jpj)     ::   zh_ekm   ! 2D workspace
      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zsh2     ! 3D     -

      ...

      !
      IF( ln_mldw ) THEN      !==  set a minimum value in the Ekman layer  ==!
         !
         DO jj = 2, jpjm1        !* Ekman depth
            DO ji = 2, jpim1
               zustar = SQRT( taum(ji,jj) * r1_rau0 )
               zhek   = rn_ekmfc * zustar / ( ABS( ff_t(ji,jj) ) + rsmall )   ! Ekman depth
               zh_ekm(ji,jj) = MAX(  rn_mldmin , MIN( zhek , rn_mldmax )  )   ! set allowed rang
            END DO
         END DO
         DO jk = 2, jpkm1        !* minimum mixing coeff. within the Ekman layer
            DO jj = 2, jpjm1
               DO ji = 2, jpim1
                  IF( gdept_n(ji,jj,jk) < zh_ekm(ji,jj) ) THEN
                     avm(ji,jj,jk) = MAX(  avm(ji,jj,jk), rn_wvmix  ) * wmask(ji,jj,jk)
                     avt(ji,jj,jk) = MAX(  avt(ji,jj,jk), rn_wtmix  ) * wmask(ji,jj,jk)
                  ENDIF
               END DO
            END DO
         END DO
      ENDIF

• Second :

The option ln_mldw =T in zdfric can't work at low latitude where Coriolis parameter becomes too low.
it provides there much to thick mixed layer. ( summer MLD of 150m over most of the GYRE configuration domain !!! )
This option needs to be more general in order to remain is the reference NEMO code.

Gurvan

[flavoni]

Done for 3.6 stable, following the comment of Gurvan in this ticket, in changeset [8595].

Nothing done for trunk because this part has be rewritten in ZDF development branch (in workplan 2017).