New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
#1715 (Bug in zdfric with simplified equation of state) – NEMO

Opened 8 years ago

Closed 7 years ago

Last modified 7 years ago

#1715 closed Bug (fixed)

Bug in zdfric with simplified equation of state

Reported by: timgraham Owned by: flavoni
Priority: low Milestone:
Component: OCE Version: v3.6
Severity: Keywords:
Cc: flavoni

Description (last modified by nicolasmartin)

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)

Commit History (1)

ChangesetAuthorTimeChangeLog
8595flavoni2017-10-05T11:49:56+02:00

#1715 fix for zdfric if ln_mldw=true

Change History (6)

comment:1 Changed 8 years ago by 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

Last edited 7 years ago by nicolasmartin (previous) (diff)

comment:2 Changed 7 years ago by clevy

  • Owner changed from nemo to mchekki
  • Version changed from trunk to v3.6

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

Last edited 7 years ago by nicolasmartin (previous) (diff)

comment:4 Changed 7 years ago by clevy

  • Cc flavoni added
  • Owner changed from mchekki to flavoni
  • Status changed from new to assigned

comment:5 Changed 7 years ago by nicolasmartin

  • Description modified (diff)

comment:6 Changed 7 years ago by flavoni

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

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

Note: See TracTickets for help on using tickets.