[[BR]] == '''S3.3-b Restoring variable volume option (key_vvl)'' == [[BR]] Last edited [[Timestamp]] [[PageOutline]] This page gives a non exhaustive description of the work done to restore VVL option from NEMO version 3.1[[BR]] 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.[[BR]] - Outside modipsl :[[BR]] ''svn co --username nemo_user !http://forge.ipsl.jussieu.fr/nemo/svn/branches/dev_004_VVL/NEMO - Using modipsl :[[BR]] Edit modipsl/util/mod.def and replace t''ags/nemo_v3_1'' by ''branches/dev_004_VVL in NEMO section''[[BR]] 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:[[BR]] - 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 [[BR]] === 2. Code organisation === ==== 2.1 Flow chart ==== (original from Mathieu Leclair)[[BR]] [[Image(NEW_flowchart.png, 50%)]] ==== 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 :[[BR]] # define fsdept_n(i,j,k) (fsdept_0(i,j,k)*(1+sshn(i,j)*mut(i,j,k)))[[BR]] Note: we shouldn't need 3d array mut (muu, muv) , but the model results change with 2D arrays.[[BR]] [[BR]] === 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 [[BR]] === 4. Cautions === - use s-coordinate pressure gradient when using variable volume (alternative scheme from Jerome is to be considered)[[BR]] - the time splitting is minimalistic (i.e. Leap-frog + asselin). Predictor-corrector scheme as in ROMS will be incorporated (see jerome's work)[[BR]] - 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[[BR]] - vertical velocity in output is now i.e. computed at the begining of step, not at the end [[BR]] [[BR]] === 5. To be done === ==== 5.1 bug corrections ==== - correct forcing term in salinity equation: [[BR]] remplace emps by fsalt, a salt flux so that : in vvl case d(e3 S)/dt = ... + fsalt (salt flux only)[[BR]] without vvl d(e3 S)/dt = ... + fsalt + emp * SSS (salt flux + concentration/dillution effect)[[BR]] see [wiki:2009WP/2009Stream3/LeapFrogC#a-trasbc.F90][[BR]] - 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 [[BR]] solve dU/dt = ... and not d( e3 U)/dt = ... ). [[BR]] 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 [[BR]] {{{ 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: [[BR]] {{{ 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.[[BR]] '''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[[BR]] {{{ ! 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). [[BR]] 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. [[BR]] - 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 [[BR]] 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. [[BR]] - 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) [[BR]] 2- remove the forcing from the asselin filter (change in tranxt) [[BR]] - 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 [[BR]] move the time stepping of the momentum equation at the end of dynspg_ts and dynspg_exp [[BR]] 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) [[BR]] === 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) ----