S3.3-b Restoring variable volume option (key_vvl)

Last edited Timestamp?

This page gives a non exhaustive description of the work done to restore VVL option from NEMO version 3.1
At this stage (april 2009) :

• variable volume option can be used with either z, zps and s coordinate in fitered free surface
• time splitting issue is under investigation

---

1. Getting the code

At this time (04-20-09), the work has been done on a branch starting from svn rev1359.

• Outside modipsl :
  

svn co --username nemo_user http://forge.ipsl.jussieu.fr/nemo/svn/branches/dev_004_VVL/NEMO

• Using modipsl :
  

Edit modipsl/util/mod.def and replace tags/nemo_v3_1 by branches/dev_004_VVL in NEMO section

Note that this branch has been created some time after nemo_v3_1 release (svn rev1332) so it includes intermediate developments and correction not dealing with vvl, including:

• complete work on time origin in outputs (ticket:335) + downward vertical axis (ticket:357)
• Update lib_mpp, see ticket #379
• add s-sigma coordinates option, ticket #378
• Update lib_mpp, see ticket #379 -first implementation of iom_put, see ticket:387
• update of diaptr

In addition, the VVL has been synchronized with the trunk :

• dev_004_VVL:sync: see ticket #417 : corrections for IOM
• dev_004_VVL:sync: synchro with trunk, see ticket #361 : correction in diaptr
• dev_004_VVL:sync: synchronisation with the trunk, see ticket #388 : bound salt exchange

2. Code organisation

2.1 Flow chart

(original from Mathieu Leclair)

2.2 Sea surface height repartition

Repartition on the whole water column, so key_sigma_vll has been suppressed. Reference coordinate is referred as e3t_0 for instance, and pre-processing (in domzgr_substitute.h90) is used to define before, now, after scale factors, for instance :
# define fsdept_n(i,j,k) (fsdept_0(i,j,k)*(1+sshn(i,j)*mut(i,j,k)))

Note: we shouldn't need 3d array mut (muu, muv) , but the model results change with 2D arrays.

3. Modified routines

• domvvl is now used only for initialisation (called from istate) (===> should be better to call it from domain....)
• dynspg_flt, dynspg_ts, dynspg_exp : no more update of ssh, just barotropic contribution for velocities
• trazdf_exp, trazdf_imp : ponderation by the correct scale factors
• wzvmod : suprresion of wzv_mod routine, addition of ssh_wzv to compute ssh after from the ssh equation, update the vertical coordinate and compute vertical velocity, addition of ssh_nxt for ssh time stepping
• dynnxt, tranxt : rewritting, take vvl case into account
• zdfevd : for stability issue, test to apply enhanced diffusion has to be performed on now Brunt-Vaisala frequency

4. Cautions

• use s-coordinate pressure gradient when using variable volume (alternative scheme from Jerome is to be considered)
  
• the time splitting is minimalistic (i.e. Leap-frog + asselin). Predictor-corrector scheme as in ROMS will be incorporated (see jerome's work)
  
• sea ice and vvl : the sea-surface boundary condition is wrong. The salt flux exchanged with sea-ice should be added to the salinity equation
  
• vertical velocity in output is now i.e. computed at the begining of step, not at the end
  

5. To be done

5.1 bug corrections

• correct forcing term in salinity equation:
  

remplace emps by fsalt, a salt flux so that : in vvl case d(e3 S)/dt = ... + fsalt (salt flux only)
without vvl d(e3 S)/dt = ... + fsalt + emp * SSS (salt flux + concentration/dillution effect)
see 2009WP/2009Stream3/LeapFrogC

• correct the left-hand-side of the momentum equation in vector invariant formulation :

time derivative of e3 should not be not included in the left hand side of the momentum equation
solve dU/dt = ... and not d( e3 U)/dt = ... ).
Change to be done in both dynnxt AND dynspg (_ts ; _flt ; _exp)

In dynnxt.F90 perform the thickness weighted time leap-frog and Asselin filter only in ln_dynadv_vec=FALSE and lk_vvl=TRUE., in other word unweighted operation if ln_dynadv_vec=T or NOT lk_vvl. So replace the 2 two

      IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
         part AAA
      ELSE                                            ! applied on thickness weighted velocity
         part BBB
      ENDIF

By

     IF( ln_dynadv_vec .OR. .NOT. lk_vvl ) THEN      ! applied on velocity
        part BBB
     ELSE                                            ! applied on thickness weighted velocity
        part AAA
     ENDIF

In addition : (1) the masking of ua, va can be moved inside the key_bdy case ; (2) all the comment have been updated. The resulting subroutine can be found in attachment (https://forge.ipsl.jussieu.fr/nemo/attachment/wiki/2009WP/2009Stream3/VVL/dynnxt.F90).

In dynspg_flt.F90, the time stepping of (ua,va) has to be modified as follows:

      IF( lk_vvl ) THEN          ! variable volume 

         ! Evaluate the masked next velocity (effect of the additional force not included)
         ! -------------------   (thickness weighted velocity, surface pressure gradient already included in dyn_hpg)
         DO jk = 1, jpkm1
            DO jj = 2, jpjm1
               DO ji = fs_2, fs_jpim1   ! vector opt.
                  ua(ji,jj,jk) = (        ub(ji,jj,jk) * fse3u_b(ji,jj,jk)      &
                     &           + z2dt * ua(ji,jj,jk) * fse3u_n(ji,jj,jk)  )   &
                     &         / fse3u_a(ji,jj,jk) * umask(ji,jj,jk)
                  va(ji,jj,jk) = (        vb(ji,jj,jk) * fse3v_b(ji,jj,jk)      &
                     &           + z2dt * va(ji,jj,jk) * fse3v_n(ji,jj,jk)  )   &
                     &         / fse3v_a(ji,jj,jk) * vmask(ji,jj,jk)
               END DO
            END DO
         END DO

      ELSE                       ! fixed volume 

         ! Surface pressure gradient (now)
         DO jj = 2, jpjm1
            DO ji = fs_2, fs_jpim1   ! vector opt.
               spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj)
               spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)
            END DO 
         END DO 
         !
         ! add the surface pressure trend to the general trend and
         ! evaluate the masked next velocity (effect of the additional force not included)
         DO jk = 1, jpkm1
            DO jj = 2, jpjm1
               DO ji = fs_2, fs_jpim1   ! vector opt.
                  ua(ji,jj,jk) = (  ub(ji,jj,jk) + z2dt * ( ua(ji,jj,jk) + spgu(ji,jj) )  ) * umask(ji,jj,jk)
                  va(ji,jj,jk) = (  vb(ji,jj,jk) + z2dt * ( va(ji,jj,jk) + spgv(ji,jj) )  ) * vmask(ji,jj,jk)
               END DO
            END DO
         END DO
         !
      ENDIF

becomes

      ! Next velocity :   Leap-frog time stepping  (effect of the additional force not included)
      ! -------------
      IF( lk_vvl ) THEN          !* variable volume : surface pressure gradient already included in dyn_hpg
         !
         IF( ln_dynadv_vec ) THEN       ! vector invariant advection form: Leap-Frog applied on velocity
            DO jk = 1, jpkm1
               ua(:,:,jk) = ( ub(:,:,jk) + z2dt * ua(:,:,jk) ) * umask(:,:,jk)
               va(:,:,jk) = ( vb(:,:,jk) + z2dt * va(:,:,jk) ) * vmask(:,:,jk)
            END DO
         ELSE                           ! flux advection form: Leap-Frog applied on thickness weighted velocity
            DO jk = 1, jpkm1
               ua(:,:,jk) = (          ub(:,:,jk) * fse3u_b(:,:,jk)     &
                  &           + z2dt * ua(:,:,jk) * fse3u_n(:,:,jk) )   &
                  &         / fse3u_a(:,:,jk) * umask(:,:,jk)
               va(:,:,jk) = (          vb(:,:,jk) * fse3v_b(:,:,jk)     &
                  &           + z2dt * va(:,:,jk) * fse3v_n(:,:,jk) )   &
                  &         / fse3v_a(:,:,jk) * vmask(:,:,jk)            
            END DO
         ENDIF
         !
      ELSE                       !* fixed volume : add the surface pressure gradient trend
         !
         DO jj = 2, jpjm1               ! now surface pressure gradient
            DO ji = fs_2, fs_jpim1   ! vector opt.
               spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj)
               spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)
            END DO 
         END DO 
         DO jk = 1, jpkm1               ! Leap-Frog applied on velocity
            DO jj = 2, jpjm1
               DO ji = fs_2, fs_jpim1   ! vector opt.
                  ua(ji,jj,jk) = (  ub(ji,jj,jk) + z2dt * ( ua(ji,jj,jk) + spgu(ji,jj) )  ) * umask(ji,jj,jk)
                  va(ji,jj,jk) = (  vb(ji,jj,jk) + z2dt * ( va(ji,jj,jk) + spgv(ji,jj) )  ) * vmask(ji,jj,jk)
               END DO
            END DO
         END DO
         !
      ENDIF

In addition, all the comment have been updated. The resulting subroutine can be found in attachment (https://forge.ipsl.jussieu.fr/nemo/attachment/wiki/2009WP/2009Stream3/VVL/dynspg_flt.F90).[[BR]]

NB: There is a bug in the time-splitting algorithm when starting with an Euler time step. Indeed, even when an Euler time step, the time-splitting algorithm have to be performed over a 2*rdt period! This can be easily understood with a sketch of the time algorithm. I (Gurvan) may add this sketch here soon.

In dynspg_exp.F90, the surface pressure gradient is already taken into account in lk_vvl=TRUE case. Therefore it must not be added in dynspg_exp. So replace

      ! Surface pressure gradient (now)
      DO jj = 2, jpjm1
         DO ji = fs_2, fs_jpim1   ! vector opt.
            spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj)
            spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)
         END DO
      END DO

      ! Add the surface pressure trend to the general trend
      DO jk = 1, jpkm1
         DO jj = 2, jpjm1
            DO ji = fs_2, fs_jpim1   ! vector opt.
               ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)
               va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)
            END DO
         END DO
      END DO

by

      IF( .NOT. lk_vvl ) THEN          !* fixed volume : add the surface pressure gradient trend
         !
         DO jj = 2, jpjm1                    ! now surface pressure gradient
            DO ji = fs_2, fs_jpim1   ! vector opt.
               spgu(ji,jj) = - grav * ( sshn(ji+1,jj) - sshn(ji,jj) ) / e1u(ji,jj)
               spgv(ji,jj) = - grav * ( sshn(ji,jj+1) - sshn(ji,jj) ) / e2v(ji,jj)
            END DO 
         END DO 
         DO jk = 1, jpkm1                    ! Add it to the general trend
            DO jj = 2, jpjm1
               DO ji = fs_2, fs_jpim1   ! vector opt.
                  ua(ji,jj,jk) = ua(ji,jj,jk) + spgu(ji,jj)
                  va(ji,jj,jk) = va(ji,jj,jk) + spgv(ji,jj)
               END DO
            END DO
         END DO
         !
      ENDIF

In addition, all the comment have been updated. The resulting subroutine can be found in attachment (https://forge.ipsl.jussieu.fr/nemo/attachment/wiki/2009WP/2009Stream3/VVL/dynspg_exp.F90).
Note that dynspg.F90 must also be modifed in order to extract the correct dynspg trends when explicit free surface is used. See the attachment: https://forge.ipsl.jussieu.fr/nemo/attachment/wiki/2009WP/2009Stream3/VVL/dynspg.F90.[[BR]]

In dynspg_ts.F90 the whole module have to be modified to solve a velocity equation, not a thickness weighted velocity equation. Then in the specific case of ln_dynadv_vec=FALSE and lk_vvl=TRUE, an term associated with the time derivative of the ssh must be added to the equation solved.

• ensure that BDY and OBC are working with the new time-stepping strategy

1- in dynnxt, obc_dyn_bt seems to modify the sshn field. This should be wrong with the new time strategy
2- in dynspg_exp, a call to obc_dta_bt appears. I don't understand why here! If it is required, it should appear in ssh_wzv routine.

• label all fse3 with its adequate time label (before, now, after) everywhere in the code. Some may be wrong today.
• suppress mu* 3D arrays, and anly use instead 2D arrays

5.2 improvments

• Add the Leclair-Madec conservative leap-Frog in order to enable strict conservation tests

1- shift all ocean forcing term to kt+1/2 (change in sbcmod and trasbc)
2- remove the forcing from the asselin filter (change in tranxt)

• replace the leap-frog by a predictor-corrector shceme for time splitting (as Jérome does)
• develop a pressure gradient algorithm specific and more accurate for zps-vvl (not only the use of the standard s-coord pressure gradient) (see Jérome work)
• Simplification for code :

work on the tracer content tendancy everywhere, not on the tracer tendency
move the time stepping of the momentum equation at the end of dynspg_ts and dynspg_exp
put the call to domvvl in domain.F90 (not in istate)

• embedded sea-ice (Gurvan and Yevgeni) : introduce 2 others ways of nandling sea-ice (in both vvl and constant cases):

1- standard case: levitating sea-ice with virtual salt flux : ice-ocean exchange are salt flux (true + C/D effct) only no volume flux 2- first new case: levitating sea-ice with salt conservation : ice-ocean exchange true salt flux + volume flux, but not suface pressure associated to sea-ice 3- 2nd new case: full sea-ice embeddement : ice-ocean exchange true salt flux + volume flux + surface pressure gradient associated with the mass of ice+snow

• introduce a optional Roullet-madec damping force in the time splitting routine (no more use of elliptic solver)

6. Tests

6.1 Paris

• ORCA2_LIM zps fixed volume and filtred free surface
• ORCA2_LIM zps variable volume and filtred free surface
• ORCA2_LIM zps fixed volume and time splitting free surface
• ORCA2_LIM zps variable volume and time splitting free surface
• restartability
• Agrif test case
• parallelisation
• conservation (with Asselin filter to 0 and ice tracer fluxes to 0)

---