S3.3b 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 (042009), 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 ssigma 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 preprocessing (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 BruntVaisala frequency
4. Cautions
 use scoordinate pressure gradient when using variable volume (alternative scheme from Jerome is to be considered)
 the time splitting is minimalistic (i.e. Leapfrog + asselin). Predictorcorrector scheme as in ROMS will be incorporated (see jerome's work)
 sea ice and vvl : the seasurface boundary condition is wrong. The salt flux exchanged with seaice 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/Leap FrogC
 correct the lefthandside 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 leapfrog 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 : Leapfrog 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: LeapFrog 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: LeapFrog 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 ! LeapFrog 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 timesplitting algorithm when starting with an Euler time step. Indeed, even when an Euler time step, the timesplitting 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 timestepping 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 LeclairMadec conservative leapFrog in order to enable strict conservation tests
1 shift all ocean forcing term to kt+½ (change in sbcmod and trasbc)
2 remove the forcing from the asselin filter (change in tranxt)
 replace the leapfrog by a predictorcorrector shceme for time splitting (as Jérome does)
 develop a pressure gradient algorithm specific and more accurate for zpsvvl (not only the use of the standard scoord 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 seaice (Gurvan and Yevgeni) : introduce 2 others ways of nandling seaice (in both vvl and constant cases):
1 standard case: levitating seaice with virtual salt flux : iceocean exchange are salt flux (true + C/D effct) only no volume flux 2 first new case: levitating seaice with salt conservation : iceocean exchange true salt flux + volume flux, but not suface pressure associated to seaice 3 2nd new case: full seaice embeddement : iceocean exchange true salt flux + volume flux + surface pressure gradient associated with the mass of ice+snow
 introduce a optional Roulletmadec 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)
Attachments (5)

NEW_flowchart.png
(216.1 KB) 
added by rblod 11 years ago.
NEMO/OPA flowchart
 dynnxt.F90 (12.4 KB)  added by gm 11 years ago.
 dynspg_flt.F90 (20.3 KB)  added by gm 11 years ago.
 dynspg_exp.F90 (5.2 KB)  added by gm 11 years ago.
 dynspg.F90 (8.8 KB)  added by gm 11 years ago.
Download all attachments as: .zip