From d68030c40eb36c40e34e7e440557f13a9b75d760 Mon Sep 17 00:00:00 2001
From: Sibylle Techene <sibylle.techene@locean-ipsl.upmc.fr>
Date: Fri, 24 Nov 2023 14:57:44 +0000
Subject: [PATCH] Resolve "add up3 advection scheme for momentum"
---
cfgs/SHARED/namelist_ref | 3 +-
src/OCE/DYN/dynadv.F90 | 22 ++-
src/OCE/DYN/dynadv_up3.F90 | 274 +++++++++++++++++++++++++++++++++++++
src/OCE/DYN/dynvor.F90 | 2 +-
4 files changed, 292 insertions(+), 9 deletions(-)
create mode 100644 src/OCE/DYN/dynadv_up3.F90
diff --git a/cfgs/SHARED/namelist_ref b/cfgs/SHARED/namelist_ref
index 79ace60c..76946a98 100644
--- a/cfgs/SHARED/namelist_ref
+++ b/cfgs/SHARED/namelist_ref
@@ -1020,7 +1020,8 @@
ln_dynadv_vec = .false. ! vector form - 2nd centered scheme
nn_dynkeg = 0 ! grad(KE) scheme: =0 C2 ; =1 Hollingsworth correction
ln_dynadv_cen2 = .false. ! flux form - 2nd order centered scheme
- ln_dynadv_ubs = .false. ! flux form - 3rd order UBS scheme
+ ln_dynadv_ubs = .false. ! flux form - 3rd order UBS OLD scheme
+ ln_dynadv_up3 = .false. ! flux form - 3rd order UBS NEW scheme
/
!-----------------------------------------------------------------------
&namdyn_vor ! Vorticity / Coriolis scheme (default: NO selection)
diff --git a/src/OCE/DYN/dynadv.F90 b/src/OCE/DYN/dynadv.F90
index c1500fd4..7271ec7b 100644
--- a/src/OCE/DYN/dynadv.F90
+++ b/src/OCE/DYN/dynadv.F90
@@ -16,6 +16,7 @@ MODULE dynadv
USE dom_oce ! ocean space and time domain
USE dynadv_cen2 ! centred flux form advection (dyn_adv_cen2 routine)
USE dynadv_ubs ! UBS flux form advection (dyn_adv_ubs routine)
+ USE dynadv_up3 ! UBS flux form advection (NEW) (dyn_adv_up3 routine)
USE dynkeg ! kinetic energy gradient (dyn_keg routine)
USE dynzad ! vertical advection (dyn_zad routine)
!
@@ -34,14 +35,16 @@ MODULE dynadv
LOGICAL, PUBLIC :: ln_dynadv_vec !: vector form
INTEGER, PUBLIC :: nn_dynkeg !: scheme of grad(KE): =0 C2 ; =1 Hollingsworth
LOGICAL, PUBLIC :: ln_dynadv_cen2 !: flux form - 2nd order centered scheme flag
- LOGICAL, PUBLIC :: ln_dynadv_ubs !: flux form - 3rd order UBS scheme flag
-
+ LOGICAL, PUBLIC :: ln_dynadv_ubs !: flux form - 3rd order UBS scheme flag (OLD)
+ LOGICAL, PUBLIC :: ln_dynadv_up3 !: flux form - 3rd order UBS scheme flag (NEW)
+
INTEGER, PUBLIC :: n_dynadv !: choice of the formulation and scheme for momentum advection
! ! associated indices:
INTEGER, PUBLIC, PARAMETER :: np_LIN_dyn = 0 ! no advection: linear dynamics
INTEGER, PUBLIC, PARAMETER :: np_VEC_c2 = 1 ! vector form : 2nd order centered scheme
INTEGER, PUBLIC, PARAMETER :: np_FLX_c2 = 2 ! flux form : 2nd order centered scheme
- INTEGER, PUBLIC, PARAMETER :: np_FLX_ubs = 3 ! flux form : 3rd order Upstream Biased Scheme
+ INTEGER, PUBLIC, PARAMETER :: np_FLX_ubs = 3 ! flux form : 3rd order Upstream Biased Scheme (OLD)
+ INTEGER, PUBLIC, PARAMETER :: np_FLX_up3 = 4 ! flux form : 3rd order Upstream Biased Scheme (NEW)
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -77,7 +80,9 @@ CONTAINS
CASE( np_FLX_c2 )
CALL dyn_adv_cen2( kt, Kmm, puu, pvv, Krhs ) ! 2nd order centered scheme
CASE( np_FLX_ubs )
- CALL dyn_adv_ubs ( kt, Kbb, Kmm, puu, pvv, Krhs ) ! 3rd order UBS scheme (UP3)
+ CALL dyn_adv_ubs ( kt, Kbb, Kmm, puu, pvv, Krhs ) ! 3rd order OLD UBS scheme (UP3)
+ CASE( np_FLX_up3 )
+ CALL dyn_adv_up3 ( kt, Kbb, Kmm, puu, pvv, Krhs ) ! 3rd order NEW UBS scheme (UP3)
END SELECT
!
IF( ln_timing ) CALL timing_stop( 'dyn_adv' )
@@ -94,7 +99,7 @@ CONTAINS
!!----------------------------------------------------------------------
INTEGER :: ioptio, ios ! Local integer
!
- NAMELIST/namdyn_adv/ ln_dynadv_OFF, ln_dynadv_vec, nn_dynkeg, ln_dynadv_cen2, ln_dynadv_ubs
+ NAMELIST/namdyn_adv/ ln_dynadv_OFF, ln_dynadv_vec, nn_dynkeg, ln_dynadv_cen2, ln_dynadv_ubs, ln_dynadv_up3
!!----------------------------------------------------------------------
!
IF(lwp) THEN
@@ -115,7 +120,8 @@ CONTAINS
WRITE(numout,*) ' Vector form: 2nd order centered scheme ln_dynadv_vec = ', ln_dynadv_vec
WRITE(numout,*) ' with Hollingsworth scheme (=1) or not (=0) nn_dynkeg = ', nn_dynkeg
WRITE(numout,*) ' flux form: 2nd order centred scheme ln_dynadv_cen2 = ', ln_dynadv_cen2
- WRITE(numout,*) ' 3rd order UBS scheme ln_dynadv_ubs = ', ln_dynadv_ubs
+ WRITE(numout,*) ' 3rd order UBS scheme (OLD) ln_dynadv_ubs = ', ln_dynadv_ubs
+ WRITE(numout,*) ' 3rd order UBS scheme (NEW) ln_dynadv_up3 = ', ln_dynadv_up3
ENDIF
ioptio = 0 ! parameter control and set n_dynadv
@@ -123,6 +129,7 @@ CONTAINS
IF( ln_dynadv_vec ) THEN ; ioptio = ioptio + 1 ; n_dynadv = np_VEC_c2 ; ENDIF
IF( ln_dynadv_cen2 ) THEN ; ioptio = ioptio + 1 ; n_dynadv = np_FLX_c2 ; ENDIF
IF( ln_dynadv_ubs ) THEN ; ioptio = ioptio + 1 ; n_dynadv = np_FLX_ubs ; ENDIF
+ IF( ln_dynadv_up3 ) THEN ; ioptio = ioptio + 1 ; n_dynadv = np_FLX_up3 ; ENDIF
IF( ioptio /= 1 ) CALL ctl_stop( 'choose ONE and only ONE advection scheme' )
IF( nn_dynkeg /= nkeg_C2 .AND. nn_dynkeg /= nkeg_HW ) CALL ctl_stop( 'KEG scheme wrong value of nn_dynkeg' )
@@ -138,7 +145,8 @@ CONTAINS
IF( nn_dynkeg == nkeg_C2 ) WRITE(numout,*) ' with Centered standard keg scheme'
IF( nn_dynkeg == nkeg_HW ) WRITE(numout,*) ' with Hollingsworth keg scheme'
CASE( np_FLX_c2 ) ; WRITE(numout,*) ' ==>>> flux form : 2nd order scheme is used'
- CASE( np_FLX_ubs ) ; WRITE(numout,*) ' ==>>> flux form : UBS scheme is used'
+ CASE( np_FLX_ubs ) ; WRITE(numout,*) ' ==>>> flux form : OLD UBS scheme is used'
+ CASE( np_FLX_up3 ) ; WRITE(numout,*) ' ==>>> flux form : NEW UBS scheme is used'
END SELECT
ENDIF
!
diff --git a/src/OCE/DYN/dynadv_up3.F90 b/src/OCE/DYN/dynadv_up3.F90
new file mode 100644
index 00000000..043eb027
--- /dev/null
+++ b/src/OCE/DYN/dynadv_up3.F90
@@ -0,0 +1,274 @@
+MODULE dynadv_up3
+ !!======================================================================
+ !! *** MODULE dynadv_up3 ***
+ !! Ocean dynamics: Update the momentum trend with the flux form advection
+ !! trend using a 3rd order upstream biased scheme
+ !!======================================================================
+ !! History : 2.0 ! 2006-08 (R. Benshila, L. Debreu) Original code
+ !! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option
+ !!----------------------------------------------------------------------
+
+ !!----------------------------------------------------------------------
+ !! dyn_adv_up3 : flux form momentum advection using (ln_dynadv=T)
+ !! an 3rd order Upstream Biased Scheme or Quick scheme
+ !! combined with 2nd or 4th order finite differences
+ !!----------------------------------------------------------------------
+ USE oce ! ocean dynamics and tracers
+ USE dom_oce ! ocean space and time domain
+ USE trd_oce ! trends: ocean variables
+ USE trddyn ! trend manager: dynamics
+ !
+ USE in_out_manager ! I/O manager
+ USE prtctl ! Print control
+ USE lbclnk ! ocean lateral boundary conditions (or mpp link)
+ USE lib_mpp ! MPP library
+
+ IMPLICIT NONE
+ PRIVATE
+
+ REAL(wp), PARAMETER :: gamma1 = 1._wp/3._wp ! =1/4 quick ; =1/3 3rd order UBS
+
+ PUBLIC dyn_adv_up3 ! routine called by step.F90
+
+ !! * Substitutions
+# include "do_loop_substitute.h90"
+# include "domzgr_substitute.h90"
+ !!----------------------------------------------------------------------
+ !! NEMO/OCE 4.0 , NEMO Consortium (2018)
+ !! $Id: dynadv_up3.F90 14834 2021-05-11 09:24:44Z hadcv $
+ !! Software governed by the CeCILL license (see ./LICENSE)
+ !!----------------------------------------------------------------------
+CONTAINS
+
+ SUBROUTINE dyn_adv_up3( kt, Kbb, Kmm, puu, pvv, Krhs )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE dyn_adv_up3 ***
+ !!
+ !! ** Purpose : Compute the now momentum advection trend in flux form
+ !! and the general trend of the momentum equation.
+ !!
+ !! ** Method : The scheme is the one implemeted in ROMS. It depends
+ !! on two parameter gamma1 and gamma2. The former control the
+ !! upstream baised part of the scheme and the later the centred
+ !! part: gamma1 = 0 pure centered (no diffusive part)
+ !! = 1/4 Quick scheme
+ !! = 1/3 3rd order Upstream biased scheme
+ !! For stability reasons, the first term of the fluxes which cor-
+ !! responds to a second order centered scheme is evaluated using
+ !! the now velocity (centered in time) while the second term which
+ !! is the diffusive part of the scheme, is evaluated using the
+ !! before velocity (forward in time).
+ !! Default value (hard coded in the begining of the module) are
+ !! gamma1=1/3 and gamma2=1/32.
+ !!
+ !! ** Action : - (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) updated with the 3D advective momentum trends
+ !!
+ !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
+ !!----------------------------------------------------------------------
+ INTEGER , INTENT( in ) :: kt ! ocean time-step index
+ INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs ! ocean time level indices
+ REAL(dp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
+ !
+ INTEGER :: ji, jj, jk ! dummy loop indices
+ REAL(wp) :: zui, zvj, zfuj, zfvi, zl_u, zl_v ! local scalars
+ REAL(wp) :: zfwi, zzfu_uw_kp1, zlu_uw_kp1 ! local scalars
+ REAL(wp) :: zfwj, zzfv_vw_kp1, zlv_vw_kp1 ! local scalars
+ REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfu_f, zfu
+ REAL(dp), DIMENSION(A2D(nn_hls),jpk) :: zfu_t
+ REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfv_f, zfv
+ REAL(dp), DIMENSION(A2D(nn_hls),jpk) :: zfv_t
+ REAL(dp), DIMENSION(A2D(nn_hls),jpk) :: ztrdu, ztrdv
+ REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zlu_uu, zlu_uv
+ REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zlv_vv, zlv_vu
+ REAL(dp), DIMENSION(A2D(nn_hls)) :: zfw, zfu_uw, zfv_vw, zlu_uw, zlv_vw
+ !!----------------------------------------------------------------------
+ !
+ IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
+ IF( kt == nit000 ) THEN
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) 'dyn_adv_ubs : UBS flux form momentum advection'
+ IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+ ENDIF
+ ENDIF
+ !
+ zfu_t(:,:,:) = 0._wp
+ zfv_t(:,:,:) = 0._wp
+ zfu_f(:,:,:) = 0._wp
+ zfv_f(:,:,:) = 0._wp
+ !
+ zlu_uu(:,:,:) = 0._wp
+ zlv_vv(:,:,:) = 0._wp
+ zlu_uv(:,:,:) = 0._wp
+ zlv_vu(:,:,:) = 0._wp
+ !
+ IF( l_trddyn ) THEN ! trends: store the input trends
+ ztrdu(:,:,:) = puu(:,:,:,Krhs)
+ ztrdv(:,:,:) = pvv(:,:,:,Krhs)
+ ENDIF
+ ! ! =========================== !
+ DO jk = 1, jpkm1 ! Laplacian of the velocity !
+ ! ! =========================== !
+ ! ! horizontal volume fluxes
+ DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+ zfu(ji,jj,jk) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Kmm)
+ zfv(ji,jj,jk) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Kmm)
+ END_2D
+ !
+ DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! laplacian
+ ! round brackets added to fix the order of floating point operations
+ ! needed to ensure halo 1 - halo 2 compatibility
+ zlu_uu(ji,jj,jk) = ( ( puu (ji+1,jj ,jk,Kbb) - puu (ji ,jj ,jk,Kbb) &
+ & ) & ! bracket for halo 1 - halo 2 compatibility
+ & + ( puu (ji-1,jj ,jk,Kbb) - puu (ji ,jj ,jk,Kbb) &
+ & ) & ! bracket for halo 1 - halo 2 compatibility
+ & ) * umask(ji ,jj ,jk)
+ zlv_vv(ji,jj,jk) = ( ( pvv (ji ,jj+1,jk,Kbb) - pvv (ji ,jj ,jk,Kbb) &
+ & ) & ! bracket for halo 1 - halo 2 compatibility
+ & + ( pvv (ji ,jj-1,jk,Kbb) - pvv (ji ,jj ,jk,Kbb) &
+ & ) & ! bracket for halo 1 - halo 2 compatibility
+ & ) * vmask(ji ,jj ,jk)
+ zlu_uv(ji,jj,jk) = ( puu (ji ,jj+1,jk,Kbb) - puu (ji ,jj ,jk,Kbb) ) * fmask(ji ,jj ,jk) &
+ & - ( puu (ji ,jj ,jk,Kbb) - puu (ji ,jj-1,jk,Kbb) ) * fmask(ji ,jj-1,jk)
+ zlv_vu(ji,jj,jk) = ( pvv (ji+1,jj ,jk,Kbb) - pvv (ji ,jj ,jk,Kbb) ) * fmask(ji ,jj ,jk) &
+ & - ( pvv (ji ,jj ,jk,Kbb) - pvv (ji-1,jj ,jk,Kbb) ) * fmask(ji-1,jj ,jk)
+ END_2D
+ END DO
+ !
+ IF( nn_hls == 1 ) CALL lbc_lnk( 'dynadv_ubs', zlu_uu(:,:,:), 'U', -1.0_wp , zlu_uv(:,:,:), 'U', -1.0_wp, &
+ & zlv_vv(:,:,:), 'V', -1.0_wp , zlv_vu(:,:,:), 'V', -1.0_wp )
+ !
+ ! ! ====================== !
+ ! ! Horizontal advection !
+ DO jk = 1, jpkm1 ! ====================== !
+ ! ! horizontal volume fluxes
+ DO_2D( 1, 1, 1, 1 )
+ zfu(ji,jj,jk) = 0.25_wp * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Kmm)
+ zfv(ji,jj,jk) = 0.25_wp * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Kmm)
+ END_2D
+ !
+ DO_2D( 1, 0, 1, 0 ) ! horizontal momentum fluxes at T- and F-point
+ zui = ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj ,jk,Kmm) )
+ zvj = ( pvv(ji,jj,jk,Kmm) + pvv(ji ,jj+1,jk,Kmm) )
+ !
+ IF( zui > 0 ) THEN ; zl_u = zlu_uu(ji ,jj,jk)
+ ELSE ; zl_u = zlu_uu(ji+1,jj,jk)
+ ENDIF
+ IF( zvj > 0 ) THEN ; zl_v = zlv_vv(ji,jj ,jk)
+ ELSE ; zl_v = zlv_vv(ji,jj+1,jk)
+ ENDIF
+ !
+ zfu_t(ji+1,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj ,jk) ) &
+ & * ( zui - gamma1 * zl_u)
+ zfv_t(ji ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji ,jj+1,jk) ) &
+ & * ( zvj - gamma1 * zl_v)
+ !
+ zfuj = ( zfu(ji,jj,jk) + zfu(ji ,jj+1,jk) )
+ zfvi = ( zfv(ji,jj,jk) + zfv(ji+1,jj ,jk) )
+ IF( zfuj > 0 ) THEN ; zl_v = zlv_vu( ji ,jj ,jk )
+ ELSE ; zl_v = zlv_vu( ji+1,jj ,jk )
+ ENDIF
+ IF( zfvi > 0 ) THEN ; zl_u = zlu_uv( ji,jj ,jk )
+ ELSE ; zl_u = zlu_uv( ji,jj+1,jk )
+ ENDIF
+ !
+ zfv_f(ji ,jj ,jk) = ( zfvi ) &
+ & * ( puu(ji,jj,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) - gamma1 * zl_u )
+ zfu_f(ji ,jj ,jk) = ( zfuj ) &
+ & * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) - gamma1 * zl_v )
+ END_2D
+ DO_2D( 0, 0, 0, 0 ) ! divergence of horizontal momentum fluxes
+ puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_t(ji+1,jj,jk) - zfu_t(ji,jj ,jk) &
+ & + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * r1_e1e2u(ji,jj) &
+ & / e3u(ji,jj,jk,Kmm)
+ pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfu_f(ji,jj ,jk) - zfu_f(ji-1,jj,jk) &
+ & + zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) * r1_e1e2v(ji,jj) &
+ & / e3v(ji,jj,jk,Kmm)
+ END_2D
+ END DO
+ IF( l_trddyn ) THEN ! trends: send trends to trddyn for diagnostic
+ ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
+ ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
+ CALL trd_dyn( ztrdu, ztrdv, jpdyn_keg, kt, Kmm )
+ zfu_t(:,:,:) = puu(:,:,:,Krhs)
+ zfv_t(:,:,:) = pvv(:,:,:,Krhs)
+ ENDIF
+ ! ! ==================== !
+ ! ! Vertical advection !
+ ! ! ==================== !
+ !
+ ! != surface vertical fluxes =! (jk = 1)
+ !
+ DO_2D( 0, 1, 0, 1 )
+ zfw(ji,jj) = e1e2t(ji,jj) * ww(ji,jj,1)
+ END_2D
+ !
+ IF( ln_linssh ) THEN ! linear free surface: advection through the surface z=0
+ DO_2D( 0, 0, 0, 0 )
+ zfu_uw(ji,jj) = 0.5_wp * ( zfw(ji,jj) + zfw(ji+1,jj ) ) * puu(ji,jj,1,Kmm)
+ zfv_vw(ji,jj) = 0.5_wp * ( zfw(ji,jj) + zfw(ji ,jj+1) ) * pvv(ji,jj,1,Kmm)
+ END_2D
+ ELSE ! non linear free: surface advective fluxes set to zero
+ DO_2D( 0, 0, 0, 0 )
+ zfu_uw(ji,jj) = 0._wp
+ zfv_vw(ji,jj) = 0._wp
+ END_2D
+ ENDIF
+ !
+ DO_2D( 0, 0, 0, 0 ) ! uniform shear in the 1st layer : dz(u(k=1)) = dz(u(k=2)) ==> zlu = 0
+ zlu_uw(ji,jj) = 0._wp
+ zlv_vw(ji,jj) = 0._wp
+ END_2D
+ !
+ DO jk = 1, jpk-2 != divergence of advective fluxes =! (jk = 1 to jpk-2)
+ DO_2D( 0, 1, 0, 1 ) ! vertical transport at level k+1
+ zfw(ji,jj) = e1e2t(ji,jj) * ww(ji,jj,jk+1)
+ END_2D
+ !
+ DO_2D( 0, 0, 0, 0 )
+ zlu_uw_kp1 = ( puu(ji,jj,jk ,Kbb) - puu(ji,jj,jk+1,Kbb) ) * wumask(ji,jj,jk+1) &
+ & - ( puu(ji,jj,jk+1,Kbb) - puu(ji,jj,jk+2,Kbb) ) * wumask(ji,jj,jk+2)
+ zlv_vw_kp1 = ( pvv(ji,jj,jk ,Kbb) - pvv(ji,jj,jk+1,Kbb) ) * wvmask(ji,jj,jk+1) &
+ & - ( pvv(ji,jj,jk+1,Kbb) - pvv(ji,jj,jk+2,Kbb) ) * wvmask(ji,jj,jk+2)
+ !
+ zfwi = zfw(ji,jj) + zfw(ji+1,jj)
+ IF( zfwi > 0 ) THEN ; zl_u = zlu_uw_kp1
+ ELSE ; zl_u = zlu_uw(ji,jj)
+ ENDIF
+ zfwj = zfw(ji,jj) + zfw(ji,jj+1)
+ IF( zfwj > 0 ) THEN ; zl_v = zlv_vw_kp1
+ ELSE ; zl_v = zlv_vw(ji,jj)
+ ENDIF
+ ! ! vertical flux at level k+1
+ zzfu_uw_kp1 = 0.25_wp * ( zfw(ji,jj) + zfw(ji+1,jj ) ) * ( puu(ji,jj,jk+1,Kmm) + puu(ji,jj,jk,Kmm) - gamma1 * zl_u )
+ zzfv_vw_kp1 = 0.25_wp * ( zfw(ji,jj) + zfw(ji ,jj+1) ) * ( pvv(ji,jj,jk+1,Kmm) + pvv(ji,jj,jk,Kmm) - gamma1 * zl_v )
+ ! ! divergence of vertical momentum flux
+ puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj) - zzfu_uw_kp1 ) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
+ pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj) - zzfv_vw_kp1 ) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
+ ! ! store vertical flux for next level calculation
+ zfu_uw(ji,jj) = zzfu_uw_kp1
+ zfv_vw(ji,jj) = zzfv_vw_kp1
+ !
+ zlu_uw(ji,jj) = zlu_uw_kp1
+ zlv_vw(ji,jj) = zlv_vw_kp1
+ END_2D
+ END DO
+ !
+ jk = jpkm1 != compute last level =! (zzfu_uw_kp1 = zzfu_uw_kp1 = 0)
+ DO_2D( 0, 0, 0, 0 )
+ puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - zfu_uw(ji,jj) * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
+ pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - zfv_vw(ji,jj) * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
+ END_2D
+ !
+ IF( l_trddyn ) THEN ! save the vertical advection trend for diagnostic
+ zfu_t(:,:,:) = puu(:,:,:,Krhs) - zfu_t(:,:,:)
+ zfv_t(:,:,:) = pvv(:,:,:,Krhs) - zfv_t(:,:,:)
+ CALL trd_dyn( zfu_t, zfv_t, jpdyn_zad, kt, Kmm )
+ ENDIF
+ ! ! Control print
+ IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' ubs2 adv - Ua: ', mask1=umask, &
+ & tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
+ !
+ END SUBROUTINE dyn_adv_up3
+
+ !!==============================================================================
+END MODULE dynadv_up3
diff --git a/src/OCE/DYN/dynvor.F90 b/src/OCE/DYN/dynvor.F90
index bf350a0f..ee767e91 100644
--- a/src/OCE/DYN/dynvor.F90
+++ b/src/OCE/DYN/dynvor.F90
@@ -1001,7 +1001,7 @@ CONTAINS
IF(lwp) WRITE(numout,*) ' ==>>> vector form dynamics : total vorticity = Coriolis + relative vorticity'
nrvm = np_RVO ! relative vorticity
ntot = np_CRV ! relative + planetary vorticity
- CASE( np_FLX_c2 , np_FLX_ubs )
+ CASE( np_FLX_c2 , np_FLX_ubs , np_FLX_up3 )
IF(lwp) WRITE(numout,*) ' ==>>> flux form dynamics : total vorticity = Coriolis + metric term'
nrvm = np_MET ! metric term
ntot = np_CME ! Coriolis + metric term
--
GitLab