From d5db5702637dec11ef8c1bcbcbd7733873b57470 Mon Sep 17 00:00:00 2001
From: Sibylle TECHENE <techenes@irene193.c-irene.tgcc.ccc.cea.fr>
Date: Wed, 27 Jul 2022 12:10:09 +0200
Subject: [PATCH] reduce local memory usage for 2 halos #80
---
src/OCE/DYN/dynadv.F90 | 25 +-
src/OCE/DYN/dynadv_cen2.F90 | 144 ++++++------
src/OCE/DYN/dynadv_ubs.F90 | 334 +++++++++++++++++++++-----
src/OCE/DYN/dynkeg.F90 | 127 +++++++++-
src/OCE/DYN/dynvor.F90 | 450 ++++++++++++++++++++++++++++++++++--
5 files changed, 924 insertions(+), 156 deletions(-)
diff --git a/src/OCE/DYN/dynadv.F90 b/src/OCE/DYN/dynadv.F90
index bda2c3a4f..6d59583ff 100644
--- a/src/OCE/DYN/dynadv.F90
+++ b/src/OCE/DYN/dynadv.F90
@@ -7,6 +7,7 @@ MODULE dynadv
!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase
!! 3.6 ! 2015-05 (N. Ducousso, G. Madec) add Hollingsworth scheme as an option
!! 4.0 ! 2017-07 (G. Madec) add a linear dynamics option
+ !! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
@@ -50,7 +51,7 @@ MODULE dynadv
!!----------------------------------------------------------------------
CONTAINS
- SUBROUTINE dyn_adv( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )
+ SUBROUTINE dyn_adv( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )
!!---------------------------------------------------------------------
!! *** ROUTINE dyn_adv ***
!!
@@ -64,7 +65,6 @@ CONTAINS
!! (see dynvor module).
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt , Kbb, Kmm, Krhs ! ocean time step and level indices
- INTEGER , OPTIONAL , INTENT(in ) :: no_zad ! no vertical advection compotation
REAL(wp), DIMENSION(:,:,:), OPTIONAL, TARGET, INTENT(in ) :: pau, pav, paw ! advective velocity
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), TARGET, INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum Eq.
!!----------------------------------------------------------------------
@@ -73,12 +73,23 @@ CONTAINS
!
SELECT CASE( n_dynadv ) !== compute advection trend and add it to general trend ==!
CASE( np_VEC_c2 ) != vector form =!
- CALL dyn_keg ( kt, nn_dynkeg , Kmm, puu, pvv, Krhs ) ! horizontal gradient of kinetic energy
- CALL dyn_zad ( kt , Kmm, puu, pvv, Krhs ) ! vertical advection
+ ! !* horizontal gradient of kinetic energy
+ IF (nn_hls==1) THEN ! halo 1 case
+ CALL dyn_keg_hls1( kt, nn_dynkeg , Kmm, puu, pvv, Krhs ) ! lbc needed with Hollingsworth scheme
+ ELSE ! halo 2 case
+ CALL dyn_keg ( kt, nn_dynkeg , Kmm, puu, pvv, Krhs )
+ ENDIF
+ CALL dyn_zad ( kt , Kmm, puu, pvv, Krhs ) !* vertical advection
+ !
CASE( np_FLX_c2 ) != flux form =!
- CALL dyn_adv_cen2( kt , Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad ) ! 2nd order centered scheme
- CASE( np_FLX_ubs )
- CALL dyn_adv_ubs ( kt , Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad ) ! 3rd order UBS scheme (UP3)
+ CALL dyn_adv_cen2( kt , Kmm, puu, pvv, Krhs, pau, pav, paw ) !* 2nd order centered scheme
+ !
+ CASE( np_FLX_ubs ) !* 3rd order UBS scheme (UP3)
+ IF (nn_hls==1) THEN ! halo 1 case
+ CALL dyn_adv_ubs_hls1( kt , Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )
+ ELSE ! halo 2 case
+ CALL dyn_adv_ubs ( kt , Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )
+ ENDIF
END SELECT
!
IF( ln_timing ) CALL timing_stop( 'dyn_adv' )
diff --git a/src/OCE/DYN/dynadv_cen2.F90 b/src/OCE/DYN/dynadv_cen2.F90
index 7fd7f65f5..07dd93112 100644
--- a/src/OCE/DYN/dynadv_cen2.F90
+++ b/src/OCE/DYN/dynadv_cen2.F90
@@ -6,6 +6,7 @@ MODULE dynadv_cen2
!!======================================================================
!! History : 2.0 ! 2006-08 (G. Madec, S. Theetten) Original code
!! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option
+ !! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
@@ -35,7 +36,7 @@ MODULE dynadv_cen2
!!----------------------------------------------------------------------
CONTAINS
- SUBROUTINE dyn_adv_cen2( kt, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )
+ SUBROUTINE dyn_adv_cen2( kt, Kmm, puu, pvv, Krhs, pau, pav, paw )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_adv_cen2 ***
!!
@@ -51,15 +52,17 @@ CONTAINS
!! ** Action : (puu,pvv)(:,:,:,Krhs) updated with the advective trend
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt , Kmm, Krhs ! ocean time-step and level indices
- INTEGER , OPTIONAL , INTENT(in ) :: no_zad ! no vertical advection computation
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), TARGET, INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
REAL(wp), DIMENSION(:,:,:), OPTIONAL, TARGET, INTENT(in ) :: pau, pav, paw ! advective velocity
!
INTEGER :: ji, jj, jk ! dummy loop indices
- REAL(wp) :: zzu, zzv ! local scalars
- REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfu_t, zfu_f, zfu_uw, zfu
- REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zfv_t, zfv_f, zfv_vw, zfv, zfw
+ REAL(wp) :: zzu, zzfu_kp1 ! local scalars
+ REAL(wp) :: zzv, zzfv_kp1 ! - -
+ REAL(wp), DIMENSION(A2D(1)) :: zfu_t, zfu_f, zfu
+ REAL(wp), DIMENSION(A2D(1)) :: zfv_t, zfv_f, zfv
+ REAL(wp), DIMENSION(A2D(1)) :: zfu_uw, zfv_vw, zfw
REAL(wp), DIMENSION(:,:,:) , POINTER :: zpt_u, zpt_v, zpt_w
+ REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: zu_trd, zv_trd
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
@@ -71,8 +74,9 @@ CONTAINS
ENDIF
!
IF( l_trddyn ) THEN ! trends: store the input trends
- zfu_uw(:,:,:) = puu(:,:,:,Krhs)
- zfv_vw(:,:,:) = pvv(:,:,:,Krhs)
+ ALLOCATE( zu_trd(A2D(0),jpkm1), zv_trd(A2D(0),jpkm1) )
+ zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs)
+ zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs)
ENDIF
!
IF( PRESENT( pau ) ) THEN ! RK3: advective velocity (pau,pav,paw) /= advected velocity (puu,pvv,ww)
@@ -89,84 +93,86 @@ CONTAINS
!
DO jk = 1, jpkm1 ! horizontal transport
DO_2D( 1, 1, 1, 1 )
- zfu(ji,jj,jk) = 0.25_wp * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * zpt_u(ji,jj,jk)
- zfv(ji,jj,jk) = 0.25_wp * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * zpt_v(ji,jj,jk)
+ zfu(ji,jj) = 0.25_wp * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * zpt_u(ji,jj,jk)
+ zfv(ji,jj) = 0.25_wp * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * zpt_v(ji,jj,jk)
END_2D
DO_2D( 1, 0, 1, 0 ) ! horizontal momentum fluxes (at T- and F-point)
- zfu_t(ji+1,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj ,jk,Kmm) )
- zfv_f(ji ,jj ,jk) = ( zfv(ji,jj,jk) + zfv(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) )
- zfu_f(ji ,jj ,jk) = ( zfu(ji,jj,jk) + zfu(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) )
- zfv_t(ji ,jj+1,jk) = ( zfv(ji,jj,jk) + zfv(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji ,jj+1,jk,Kmm) )
+ zfu_t(ji+1,jj ) = ( zfu(ji,jj) + zfu(ji+1,jj) ) * ( puu(ji,jj,jk,Kmm) + puu(ji+1,jj ,jk,Kmm) )
+ zfv_f(ji ,jj ) = ( zfv(ji,jj) + zfv(ji+1,jj) ) * ( puu(ji,jj,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) )
+ zfu_f(ji ,jj ) = ( zfu(ji,jj) + zfu(ji,jj+1) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) )
+ zfv_t(ji ,jj+1) = ( zfv(ji,jj) + zfv(ji,jj+1) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji ,jj+1,jk,Kmm) )
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) &
+ puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_t(ji+1,jj) - zfu_t(ji,jj ) &
+ & + zfv_f(ji ,jj) - zfv_f(ji,jj-1) ) * 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) &
+ pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfu_f(ji,jj ) - zfu_f(ji-1,jj) &
+ & + zfv_t(ji,jj+1) - zfv_t(ji ,jj) ) * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,jk,Kmm)
END_2D
END DO
!
IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic
- zfu_uw(:,:,:) = puu(:,:,:,Krhs) - zfu_uw(:,:,:)
- zfv_vw(:,:,:) = pvv(:,:,:,Krhs) - zfv_vw(:,:,:)
- CALL trd_dyn( zfu_uw, zfv_vw, jpdyn_keg, kt, Kmm )
- zfu_t(:,:,:) = puu(:,:,:,Krhs)
- zfv_t(:,:,:) = pvv(:,:,:,Krhs)
+ zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs) - zu_trd(A2D(0),:)
+ zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs) - zv_trd(A2D(0),:)
+ CALL trd_dyn( zu_trd, zv_trd, jpdyn_keg, kt, Kmm )
+ zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs)
+ zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs)
ENDIF
!
- IF( PRESENT( no_zad ) ) THEN !== No vertical advection ==! (except if linear free surface)
- ! ==
- IF( ln_linssh ) THEN ! linear free surface: advection through the surface z=0
- DO_2D( 0, 0, 0, 0 )
- zzu = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
- zzv = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
- puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) - zzu * r1_e1e2u(ji,jj) &
- & / e3u(ji,jj,1,Kmm)
- pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) - zzv * r1_e1e2v(ji,jj) &
- & / e3v(ji,jj,1,Kmm)
- END_2D
- ENDIF
- !
- ELSE !== Vertical advection ==!
+ ! !== Vertical advection ==!
+ !
+ ! ! surface vertical fluxes
+ !
+ 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 * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
+ zfv_vw(ji,jj) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,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 jk = 1, jpk-2 ! divergence of advective fluxes
!
- DO_2D( 0, 0, 0, 0 ) ! surface/bottom advective fluxes set to zero
- zfu_uw(ji,jj,jpk) = 0._wp ; zfv_vw(ji,jj,jpk) = 0._wp
- zfu_uw(ji,jj, 1 ) = 0._wp ; zfv_vw(ji,jj, 1 ) = 0._wp
+ DO_2D( 0, 1, 0, 1 ) ! 1/4 * Vertical transport at level k+1
+ zfw(ji,jj) = 0.25_wp * e1e2t(ji,jj) * zpt_w(ji,jj,jk+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,1) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
- zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
- END_2D
- ENDIF
- DO jk = 2, jpkm1 ! interior advective fluxes
- DO_2D( 0, 1, 0, 1 ) ! 1/4 * Vertical transport
- zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * zpt_w(ji,jj,jk)
- END_2D
- DO_2D( 0, 0, 0, 0 )
- zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji+1,jj ,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji,jj,jk-1,Kmm) )
- zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk) + zfw(ji ,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk-1,Kmm) )
- END_2D
- END DO
- DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! divergence of vertical momentum flux divergence
- puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) &
+ DO_2D( 0, 0, 0, 0 )
+ ! ! vertical flux at level k+1
+ zzfu_kp1 = ( zfw(ji,jj) + zfw(ji+1,jj ) ) * ( puu(ji,jj,jk+1,Kmm) + puu(ji,jj,jk,Kmm) )
+ zzfv_kp1 = ( zfw(ji,jj) + zfw(ji ,jj+1) ) * ( pvv(ji,jj,jk+1,Kmm) + pvv(ji,jj,jk,Kmm) )
+ ! ! divergence of vertical momentum flux
+ puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj) - zzfu_kp1 ) * r1_e1e2u(ji,jj) &
& / e3u(ji,jj,jk,Kmm)
- pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) &
- & / e3v(ji,jj,jk,Kmm)
- END_3D
- !
- IF( l_trddyn ) THEN ! trends: send trend to trddyn 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=' cen2 adv - Ua: ', mask1=umask, &
- & tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
- !
+ pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj) - zzfv_kp1 ) * r1_e1e2v(ji,jj) &
+ & / e3v(ji,jj,jk,Kmm)
+ ! ! store vertical flux for next level calculation
+ zfu_uw(ji,jj) = zzfu_kp1
+ zfv_vw(ji,jj) = zzfv_kp1
+ END_2D
+ END DO
+ !
+ jk = jpkm1
+ 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 ! trends: send trend to trddyn for diagnostic
+ zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs) - zu_trd(A2D(0),:)
+ zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs) - zv_trd(A2D(0),:)
+ CALL trd_dyn( zu_trd, zv_trd, jpdyn_zad, kt, Kmm )
+ DEALLOCATE( zu_trd, zv_trd )
ENDIF
+ ! ! Control print
+ IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' cen2 adv - Ua: ', mask1=umask, &
+ & tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
END SUBROUTINE dyn_adv_cen2
diff --git a/src/OCE/DYN/dynadv_ubs.F90 b/src/OCE/DYN/dynadv_ubs.F90
index 645f6fa80..8ab908124 100644
--- a/src/OCE/DYN/dynadv_ubs.F90
+++ b/src/OCE/DYN/dynadv_ubs.F90
@@ -6,6 +6,7 @@ MODULE dynadv_ubs
!!======================================================================
!! History : 2.0 ! 2006-08 (R. Benshila, L. Debreu) Original code
!! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option
+ !! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
@@ -29,7 +30,8 @@ MODULE dynadv_ubs
REAL(wp), PARAMETER :: gamma1 = 1._wp/3._wp ! =1/4 quick ; =1/3 3rd order UBS
REAL(wp), PARAMETER :: gamma2 = 1._wp/32._wp ! =0 2nd order ; =1/32 4th order centred
- PUBLIC dyn_adv_ubs ! routine called by step.F90
+ PUBLIC dyn_adv_ubs ! routine called by dynadv.F90
+ PUBLIC dyn_adv_ubs_hls1 ! routine called by dynadv.F90
!! * Substitutions
# include "do_loop_substitute.h90"
@@ -41,7 +43,243 @@ MODULE dynadv_ubs
!!----------------------------------------------------------------------
CONTAINS
- SUBROUTINE dyn_adv_ubs( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )
+ SUBROUTINE dyn_adv_ubs( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE dyn_adv_ubs ***
+ !!
+ !! ** 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
+ !! gamma2 = 0 2nd order finite differencing
+ !! = 1/32 4th order finite differencing
+ !! 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.
+ !!
+ !! In RK3 time stepping case, the optional arguments
+ !! (pau,pav,paw) are present. They are used as advective velocity
+ !! while the advected velocity remains (puu,pvv).
+ !!
+ !! ** Action : (puu,pvv)(:,:,:,Krhs) updated with the advective trend
+ !!
+ !! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
+ !!----------------------------------------------------------------------
+ INTEGER , INTENT(in ) :: kt , Kbb, Kmm, Krhs ! ocean time-step and level indices
+ REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), TARGET, INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
+ REAL(wp), DIMENSION(:,:,:), OPTIONAL, TARGET, INTENT(in ) :: pau, pav, paw ! advective velocity
+ !
+ INTEGER :: ji, jj, jk ! dummy loop indices
+ REAL(wp) :: zzu, zui, zfuj, zl_u, zzfu_kp1 ! local scalars
+ REAL(wp) :: zzv, zvj, zfvi, zl_v, zzfv_kp1 ! - -
+ REAL(wp), DIMENSION(A2D(2)) :: zfu_t, zfu_f, zfu
+ REAL(wp), DIMENSION(A2D(2)) :: zfv_t, zfv_f, zfv
+ REAL(wp), DIMENSION(A2D(2),2) :: zlu_uu, zlu_uv
+ REAL(wp), DIMENSION(A2D(2),2) :: zlv_vv, zlv_vu
+ REAL(wp), DIMENSION(:,:,:) , POINTER :: zpt_u, zpt_v, zpt_w
+ REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: zu_trd, zv_trd
+ !!----------------------------------------------------------------------
+ !
+ 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
+ !
+ IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic
+ ALLOCATE( zu_trd(A2D(0),jpkm1), zv_trd(A2D(0),jpkm1) )
+ zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs)
+ zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs)
+ ENDIF
+ !
+ IF( PRESENT( pau ) ) THEN ! RK3: advective velocity (pau,pav,paw) /= advected velocity (puu,pvv,ww)
+ zpt_u => pau(:,:,:)
+ zpt_v => pav(:,:,:)
+ zpt_w => paw(:,:,:)
+ ELSE ! MLF: advective velocity = (puu,pvv,ww)
+ zpt_u => puu(:,:,:,Kmm)
+ zpt_v => pvv(:,:,:,Kmm)
+ zpt_w => ww (:,:,: )
+ ENDIF
+ !
+ ! ! =========================== !
+ DO jk = 1, jpkm1 ! Laplacian of the velocity !
+ ! ! =========================== !
+ ! ! horizontal volume fluxes
+ DO_2D( 2, 2, 2, 2 )
+ zfu(ji,jj) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * zpt_u(ji,jj,jk)
+ zfv(ji,jj) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * zpt_v(ji,jj,jk)
+ END_2D
+ !
+ DO_2D( 1, 1, 1, 1 ) ! laplacian
+ ! round brackets added to fix the order of floating point operations
+ ! needed to ensure halo 1 - halo 2 compatibility (north fold)
+!! zlu_uu(ji,jj,1) = ( ( puu(ji+1,jj ,jk,Kbb) + puu(ji-1,jj ,jk,Kbb) ) - 2._wp * puu(ji,jj,jk,Kbb) ) * umask(ji,jj,jk)
+!! zlv_vv(ji,jj,1) = ( ( pvv(ji ,jj+1,jk,Kbb) + pvv(ji ,jj-1,jk,Kbb) ) - 2._wp * pvv(ji,jj,jk,Kbb) ) * vmask(ji,jj,jk)
+ zlu_uu(ji,jj,1) = ( ( 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) &
+ & ) ) * umask(ji ,jj ,jk)
+ zlv_vv(ji,jj,1) = ( ( 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) &
+ & ) ) * vmask(ji ,jj ,jk)
+ zlu_uv(ji,jj,1) = ( 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,1) = ( 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)
+ !
+!! zlu_uu(ji,jj,2) = ( ( zfu(ji+1,jj ) + zfu(ji-1,jj ) ) - 2._wp * zfu(ji,jj) ) * umask(ji ,jj ,jk)
+!! zlv_vv(ji,jj,2) = ( ( zfv(ji ,jj+1) + zfv(ji ,jj-1) ) - 2._wp * zfv(ji,jj) ) * vmask(ji ,jj ,jk)
+ zlu_uu(ji,jj,2) = ( ( zfu(ji+1,jj ) - zfu(ji ,jj ) &
+ & ) & ! bracket for halo 1 - halo 2 compatibility
+ & + ( zfu(ji-1,jj ) - zfu(ji ,jj ) &
+ & ) ) * umask(ji ,jj ,jk)
+ zlv_vv(ji,jj,2) = ( ( zfv(ji ,jj+1) - zfv(ji ,jj ) &
+ & ) & ! bracket for halo 1 - halo 2 compatibility
+ & + ( zfv(ji ,jj-1) - zfv(ji ,jj ) &
+ & ) ) * vmask(ji ,jj ,jk)
+ zlu_uv(ji,jj,2) = ( zfu(ji ,jj+1) - zfu(ji ,jj ) ) * fmask(ji ,jj ,jk) &
+ & - ( zfu(ji ,jj ) - zfu(ji ,jj-1) ) * fmask(ji ,jj-1,jk)
+ zlv_vu(ji,jj,2) = ( zfv(ji+1,jj ) - zfv(ji ,jj ) ) * fmask(ji ,jj ,jk) &
+ & - ( zfv(ji ,jj ) - zfv(ji-1,jj ) ) * fmask(ji-1,jj ,jk)
+ END_2D
+ !
+ ! ! ====================== !
+ ! ! Horizontal advection !
+ ! ! ====================== !
+ ! ! horizontal volume fluxes
+ DO_2D( 1, 1, 1, 1 )
+ zfu(ji,jj) = 0.25_wp * zfu(ji,jj)
+ zfv(ji,jj) = 0.25_wp * zfv(ji,jj)
+ 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,1)
+ ELSE ; zl_u = zlu_uu(ji+1,jj,1)
+ ENDIF
+ IF( zvj > 0 ) THEN ; zl_v = zlv_vv(ji,jj ,1)
+ ELSE ; zl_v = zlv_vv(ji,jj+1,1)
+ ENDIF
+ !
+ zfu_t(ji+1,jj ) = ( zfu(ji,jj) + zfu(ji+1,jj ) - gamma2 * ( zlu_uu(ji,jj,2) + zlu_uu(ji+1,jj ,2) ) ) &
+ & * ( zui - gamma1 * zl_u )
+ zfv_t(ji ,jj+1) = ( zfv(ji,jj) + zfv(ji ,jj+1) - gamma2 * ( zlv_vv(ji,jj,2) + zlv_vv(ji ,jj+1,2) ) ) &
+ & * ( zvj - gamma1 * zl_v )
+ !
+ zfuj = ( zfu(ji,jj) + zfu(ji ,jj+1) )
+ zfvi = ( zfv(ji,jj) + zfv(ji+1,jj ) )
+ IF( zfuj > 0 ) THEN ; zl_v = zlv_vu(ji ,jj,1)
+ ELSE ; zl_v = zlv_vu(ji+1,jj,1)
+ ENDIF
+ IF( zfvi > 0 ) THEN ; zl_u = zlu_uv(ji,jj ,1)
+ ELSE ; zl_u = zlu_uv(ji,jj+1,1)
+ ENDIF
+ !
+ zfv_f(ji ,jj ) = ( zfvi - gamma2 * ( zlv_vu(ji,jj,2) + zlv_vu(ji+1,jj ,2) ) ) &
+ & * ( puu(ji,jj,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) - gamma1 * zl_u )
+ zfu_f(ji ,jj ) = ( zfuj - gamma2 * ( zlu_uv(ji,jj,2) + zlu_uv(ji ,jj+1,2) ) ) &
+ & * ( 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) - zfu_t(ji,jj ) &
+ & + zfv_f(ji ,jj) - zfv_f(ji,jj-1) ) * r1_e1e2u(ji,jj) &
+ & / e3u(ji,jj,jk,Kmm)
+ pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfu_f(ji,jj ) - zfu_f(ji-1,jj) &
+ & + zfv_t(ji,jj+1) - zfv_t(ji ,jj) ) * r1_e1e2v(ji,jj) &
+ & / e3v(ji,jj,jk,Kmm)
+ END_2D
+ END DO
+ !
+ IF( l_trddyn ) THEN ! trends: send trend to trddyn for diagnostic
+ zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs) - zu_trd(A2D(0),:)
+ zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs) - zv_trd(A2D(0),:)
+ CALL trd_dyn( zu_trd, zv_trd, jpdyn_keg, kt, Kmm )
+ zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs)
+ zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs)
+ ENDIF
+ ! ! ==================== !
+ ! ! Vertical advection !
+ ! ! ==================== !
+ !
+#define zfu_uw zfu_t
+#define zfv_vw zfv_t
+#define zfw zfu
+ !
+ ! ! surface vertical fluxes
+ !
+ 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 * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
+ zfv_vw(ji,jj) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,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 jk = 1, jpk-2 ! divergence of advective fluxes
+ !
+ DO_2D( 0, 1, 0, 1 ) ! 1/4 * Vertical transport at level k+1
+ zfw(ji,jj) = 0.25_wp * e1e2t(ji,jj) * zpt_w(ji,jj,jk+1)
+ END_2D
+ DO_2D( 0, 0, 0, 0 )
+ ! ! vertical flux at level k+1
+ zzfu_kp1 = ( zfw(ji,jj) + zfw(ji+1,jj ) ) * ( puu(ji,jj,jk+1,Kmm) + puu(ji,jj,jk,Kmm) )
+ zzfv_kp1 = ( zfw(ji,jj) + zfw(ji ,jj+1) ) * ( pvv(ji,jj,jk+1,Kmm) + pvv(ji,jj,jk,Kmm) )
+ ! ! divergence of vertical momentum flux
+ puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj) - zzfu_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_kp1 ) * r1_e1e2v(ji,jj) &
+ & / e3v(ji,jj,jk,Kmm)
+ ! ! store vertical flux for next level calculation
+ zfu_uw(ji,jj) = zzfu_kp1
+ zfv_vw(ji,jj) = zzfv_kp1
+ END_2D
+ END DO
+ !
+ jk = jpkm1 ! compute last level (zzfu_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 ! trends: send trend to trddyn for diagnostic
+ zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs) - zu_trd(A2D(0),:)
+ zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs) - zv_trd(A2D(0),:)
+ CALL trd_dyn( zu_trd, zv_trd, jpdyn_zad, kt, Kmm )
+ DEALLOCATE( zu_trd, zv_trd )
+ ENDIF
+ !
+#undef zfu_uw
+#undef zfv_vw
+#undef zfw
+ ! ! 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_ubs
+
+
+ SUBROUTINE dyn_adv_ubs_hls1( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_adv_ubs ***
!!
@@ -73,7 +311,6 @@ CONTAINS
!! Reference : Shchepetkin & McWilliams, 2005, Ocean Modelling.
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt , Kbb, Kmm, Krhs ! ocean time-step and level indices
- INTEGER , OPTIONAL , INTENT(in ) :: no_zad ! no vertical advection compotation
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), TARGET, INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
REAL(wp), DIMENSION(:,:,:), OPTIONAL, TARGET, INTENT(in ) :: pau, pav, paw ! advective velocity
!
@@ -89,7 +326,7 @@ CONTAINS
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,*) 'dyn_adv_ubs_hls1 : UBS flux form momentum advection'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
ENDIF
ENDIF
@@ -230,63 +467,44 @@ CONTAINS
! ! Vertical advection !
! ! ==================== !
!
- ! ! ======================== !
- IF( PRESENT( no_zad ) ) THEN ! No vertical advection ! (except if linear free surface)
- ! ! ======================== ! ------
- !
- IF( ln_linssh ) THEN ! linear free surface: advection through the surface z=0
- DO_2D( 0, 0, 0, 0 )
- zzu = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
- zzv = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
- puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) - zzu * r1_e1e2u(ji,jj) &
- & / e3u(ji,jj,1,Kmm)
- pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) - zzv * r1_e1e2v(ji,jj) &
- & / e3v(ji,jj,1,Kmm)
- END_2D
- ENDIF
- ! ! =================== !
- ELSE ! Vertical advection !
- ! ! =================== !
- DO_2D( 0, 0, 0, 0 ) ! surface/bottom advective fluxes set to zero
- zfu_uw(ji,jj,jpk) = 0._wp
- zfv_vw(ji,jj,jpk) = 0._wp
- zfu_uw(ji,jj, 1 ) = 0._wp
- zfv_vw(ji,jj, 1 ) = 0._wp
+ DO_2D( 0, 0, 0, 0 ) ! surface/bottom advective fluxes set to zero
+ zfu_uw(ji,jj,jpk) = 0._wp
+ zfv_vw(ji,jj,jpk) = 0._wp
+ zfu_uw(ji,jj, 1 ) = 0._wp
+ zfv_vw(ji,jj, 1 ) = 0._wp
+ END_2D
+ IF( ln_linssh ) THEN ! constant volume : advection through the surface
+ DO_2D( 0, 0, 0, 0 )
+ zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
+ zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
END_2D
- IF( ln_linssh ) THEN ! constant volume : advection through the surface
- DO_2D( 0, 0, 0, 0 )
- zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
- zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm)
- END_2D
- ENDIF
- DO jk = 2, jpkm1 ! interior fluxes
- DO_2D( 0, 1, 0, 1 )
- zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * zpt_w(ji,jj,jk)
- END_2D
- DO_2D( 0, 0, 0, 0 )
- zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji,jj,jk-1,Kmm) )
- zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk-1,Kmm) )
- END_2D
- END DO
- DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! divergence of vertical momentum flux divergence
- puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) &
- & / e3u(ji,jj,jk,Kmm)
- pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) &
- & / e3v(ji,jj,jk,Kmm)
- END_3D
- !
- 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' )
- !
ENDIF
+ DO jk = 2, jpkm1 ! interior fluxes
+ DO_2D( 0, 1, 0, 1 )
+ zfw(ji,jj,jk) = 0.25_wp * e1e2t(ji,jj) * zpt_w(ji,jj,jk)
+ END_2D
+ DO_2D( 0, 0, 0, 0 )
+ zfu_uw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji+1,jj,jk) ) * ( puu(ji,jj,jk,Kmm) + puu(ji,jj,jk-1,Kmm) )
+ zfv_vw(ji,jj,jk) = ( zfw(ji,jj,jk)+ zfw(ji,jj+1,jk) ) * ( pvv(ji,jj,jk,Kmm) + pvv(ji,jj,jk-1,Kmm) )
+ END_2D
+ END DO
+ DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! divergence of vertical momentum flux divergence
+ puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_uw(ji,jj,jk) - zfu_uw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) &
+ & / e3u(ji,jj,jk,Kmm)
+ pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfv_vw(ji,jj,jk) - zfv_vw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) &
+ & / e3v(ji,jj,jk,Kmm)
+ END_3D
!
- END SUBROUTINE dyn_adv_ubs
+ 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_ubs_hls1
!!==============================================================================
END MODULE dynadv_ubs
diff --git a/src/OCE/DYN/dynkeg.F90 b/src/OCE/DYN/dynkeg.F90
index d751899d0..10e4134b8 100644
--- a/src/OCE/DYN/dynkeg.F90
+++ b/src/OCE/DYN/dynkeg.F90
@@ -6,7 +6,8 @@ MODULE dynkeg
!! History : 1.0 ! 1987-09 (P. Andrich, M.-A. Foujols) Original code
!! 7.0 ! 1997-05 (G. Madec) Split dynber into dynkeg and dynhpg
!! NEMO 1.0 ! 2002-07 (G. Madec) F90: Free form and module
- !! 3.6 ! 2015-05 (N. Ducousso, G. Madec) add Hollingsworth scheme as an option
+ !! 3.6 ! 2015-05 (N. Ducousso, G. Madec) add Hollingsworth scheme as an option
+ !! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
@@ -27,7 +28,8 @@ MODULE dynkeg
IMPLICIT NONE
PRIVATE
- PUBLIC dyn_keg ! routine called by step module
+ PUBLIC dyn_keg ! routine called by step module
+ PUBLIC dyn_keg_hls1 ! routine called by step module
INTEGER, PARAMETER, PUBLIC :: nkeg_C2 = 0 !: 2nd order centered scheme (standard scheme)
INTEGER, PARAMETER, PUBLIC :: nkeg_HW = 1 !: Hollingsworth et al., QJRMS, 1983
@@ -44,6 +46,123 @@ MODULE dynkeg
CONTAINS
SUBROUTINE dyn_keg( kt, kscheme, Kmm, puu, pvv, Krhs )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE dyn_keg ***
+ !!
+ !! ** Purpose : Compute the now momentum trend due to the horizontal
+ !! gradient of the horizontal kinetic energy and add it to the
+ !! general momentum trend.
+ !!
+ !! ** Method : * kscheme = nkeg_C2 : 2nd order centered scheme that
+ !! conserve kinetic energy. Compute the now horizontal kinetic energy
+ !! zhke = 1/2 [ mi-1( un^2 ) + mj-1( vn^2 ) ]
+ !! * kscheme = nkeg_HW : Hollingsworth correction following
+ !! Arakawa (2001). The now horizontal kinetic energy is given by:
+ !! zhke = 1/6 [ mi-1( 2 * un^2 + ((u(j+1)+u(j-1))/2)^2 )
+ !! + mj-1( 2 * vn^2 + ((v(i+1)+v(i-1))/2)^2 ) ]
+ !!
+ !! Take its horizontal gradient and add it to the general momentum
+ !! trend.
+ !! u(rhs) = u(rhs) - 1/e1u di[ zhke ]
+ !! v(rhs) = v(rhs) - 1/e2v dj[ zhke ]
+ !!
+ !! ** Action : - Update the (puu(:,:,:,Krhs), pvv(:,:,:,Krhs)) with the hor. ke gradient trend
+ !! - send this trends to trd_dyn (l_trddyn=T) for post-processing
+ !!
+ !! ** References : Arakawa, A., International Geophysics 2001.
+ !! Hollingsworth et al., Quart. J. Roy. Meteor. Soc., 1983.
+ !!----------------------------------------------------------------------
+ INTEGER , INTENT(in ) :: kt ! ocean time-step index
+ INTEGER , INTENT(in ) :: kscheme ! =0/1 type of KEG scheme
+ INTEGER , INTENT(in ) :: Kmm, Krhs ! ocean time level indices
+ REAL(wp), 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) :: zu, zv ! local scalars
+ REAL(wp), DIMENSION(:,: ) , ALLOCATABLE :: zhke
+ REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: zu_trd, zv_trd
+ !!----------------------------------------------------------------------
+ !
+ IF( ln_timing ) CALL timing_start('dyn_keg')
+ !
+ 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_keg : kinetic energy gradient trend, scheme number=', kscheme
+ IF(lwp) WRITE(numout,*) '~~~~~~~'
+ ENDIF
+ ENDIF
+ !
+ IF( l_trddyn ) THEN ! Save the input trends
+ ALLOCATE( zu_trd(A2D(0),jpk), zv_trd(A2D(0),jpk) )
+ zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs)
+ zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs)
+ ENDIF
+ !
+ SELECT CASE ( kscheme )
+ !
+ CASE ( nkeg_C2 ) !== Standard scheme ==!
+ ALLOCATE( zhke(A2D(1)) )
+ DO jk = 1, jpkm1
+ DO_2D( 0, 1, 0, 1 ) !* Horizontal kinetic energy at T-point
+ zu = puu(ji-1,jj ,jk,Kmm) * puu(ji-1,jj ,jk,Kmm) &
+ & + puu(ji ,jj ,jk,Kmm) * puu(ji ,jj ,jk,Kmm)
+ zv = pvv(ji ,jj-1,jk,Kmm) * pvv(ji ,jj-1,jk,Kmm) &
+ & + pvv(ji ,jj ,jk,Kmm) * pvv(ji ,jj ,jk,Kmm)
+ zhke(ji,jj) = 0.25_wp * ( zv + zu )
+ END_2D
+ !
+ DO_2D( 0, 0, 0, 0 ) !* grad( KE ) added to the general momentum trends
+ puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zhke(ji+1,jj ) - zhke(ji,jj) ) * r1_e1u(ji,jj)
+ pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zhke(ji ,jj+1) - zhke(ji,jj) ) * r1_e2v(ji,jj)
+ END_2D
+ END DO
+ DEALLOCATE( zhke )
+ !
+ CASE ( nkeg_HW ) !* Hollingsworth scheme
+ ALLOCATE( zhke(A2D(1)) )
+ DO jk = 1, jpkm1
+ DO_2D( 0, 1, 0, 1 )
+ ! round brackets added to fix the order of floating point operations
+ ! needed to ensure halo 1 - halo 2 compatibility
+ zu = ( puu(ji-1,jj ,jk,Kmm) * puu(ji-1,jj ,jk,Kmm) &
+ & + puu(ji ,jj ,jk,Kmm) * puu(ji ,jj ,jk,Kmm) ) * 8._wp &
+ & + ( ( puu(ji-1,jj-1,jk,Kmm) + puu(ji-1,jj+1,jk,Kmm) ) * ( puu(ji-1,jj-1,jk,Kmm) + puu(ji-1,jj+1,jk,Kmm) ) &
+ & + ( puu(ji ,jj-1,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) ) * ( puu(ji ,jj-1,jk,Kmm) + puu(ji ,jj+1,jk,Kmm) ) &
+ & ) ! bracket for halo 1 - halo 2 compatibility
+ zv = ( pvv(ji ,jj-1,jk,Kmm) * pvv(ji ,jj-1,jk,Kmm) &
+ & + pvv(ji ,jj ,jk,Kmm) * pvv(ji ,jj ,jk,Kmm) ) * 8._wp &
+ & + ( ( pvv(ji-1,jj-1,jk,Kmm) + pvv(ji+1,jj-1,jk,Kmm) ) * ( pvv(ji-1,jj-1,jk,Kmm) + pvv(ji+1,jj-1,jk,Kmm) ) &
+ & + ( pvv(ji-1,jj ,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) ) * ( pvv(ji-1,jj ,jk,Kmm) + pvv(ji+1,jj ,jk,Kmm) ) &
+ & ) ! bracket for halo 1 - halo 2 compatibility
+ zhke(ji,jj) = r1_48 * ( zv + zu )
+ END_2D
+ !
+ DO_2D( 0, 0, 0, 0 ) !* grad( KE ) added to the general momentum trends
+ puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zhke(ji+1,jj ) - zhke(ji,jj) ) * r1_e1u(ji,jj)
+ pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zhke(ji ,jj+1) - zhke(ji,jj) ) * r1_e2v(ji,jj)
+ END_2D
+ END DO
+ DEALLOCATE( zhke )
+ !
+ END SELECT
+ !
+ IF( l_trddyn ) THEN ! save the Kinetic Energy trends for diagnostic
+ zu_trd(A2D(0),:) = puu(A2D(0),:,Krhs) - zu_trd(A2D(0),:)
+ zv_trd(A2D(0),:) = pvv(A2D(0),:,Krhs) - zv_trd(A2D(0),:)
+ CALL trd_dyn( zu_trd, zv_trd, jpdyn_keg, kt, Kmm )
+ DEALLOCATE( zu_trd, zv_trd )
+ ENDIF
+ !
+ IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' keg - Ua: ', mask1=umask, &
+ & tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
+ !
+ IF( ln_timing ) CALL timing_stop('dyn_keg')
+ !
+ END SUBROUTINE dyn_keg
+
+
+ SUBROUTINE dyn_keg_hls1( kt, kscheme, Kmm, puu, pvv, Krhs )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_keg ***
!!
@@ -86,7 +205,7 @@ CONTAINS
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_keg : kinetic energy gradient trend, scheme number=', kscheme
+ IF(lwp) WRITE(numout,*) 'dyn_keg_hls1 : kinetic energy gradient trend, scheme number=', kscheme
IF(lwp) WRITE(numout,*) '~~~~~~~'
ENDIF
ENDIF
@@ -147,7 +266,7 @@ CONTAINS
!
IF( ln_timing ) CALL timing_stop('dyn_keg')
!
- END SUBROUTINE dyn_keg
+ END SUBROUTINE dyn_keg_hls1
!!======================================================================
END MODULE dynkeg
diff --git a/src/OCE/DYN/dynvor.F90 b/src/OCE/DYN/dynvor.F90
index 03dda78ce..b7a84124a 100644
--- a/src/OCE/DYN/dynvor.F90
+++ b/src/OCE/DYN/dynvor.F90
@@ -22,6 +22,7 @@ MODULE dynvor
!! - ! 2018-04 (G. Madec) add pre-computed gradient for metric term calculation
!! 4.x ! 2020-03 (G. Madec, A. Nasser) make ln_dynvor_msk truly efficient on relative vorticity
!! 4.2 ! 2020-12 (G. Madec, E. Clementi) add vortex force trends (ln_vortex_force=T)
+ !! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
@@ -172,11 +173,20 @@ CONTAINS
CALL vor_enT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force
ENDIF
CASE( np_EET ) !* energy conserving scheme (een scheme using e3t)
+ IF( nn_hls == 1 ) THEN
+ CALL vor_eeT_hls1( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
+ IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN
+ CALL vor_eeT_hls1( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
+ ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN
+ CALL vor_eeT_hls1( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force
+ ENDIF
+ ELSE
CALL vor_eeT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
- IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN
+ IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN
CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
- ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN
+ ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN
CALL vor_eeT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force
+ ENDIF
ENDIF
CASE( np_ENE ) !* energy conserving scheme
CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
@@ -199,11 +209,20 @@ CONTAINS
IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
IF( ln_vortex_force ) CALL vor_ens( kt, Kmm, nrvm, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add vortex force
CASE( np_EEN ) !* energy and enstrophy conserving scheme
+ IF( nn_hls == 1 ) THEN
+ CALL vor_een_hls1( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
+ IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN
+ CALL vor_een_hls1( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
+ ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN
+ CALL vor_een_hls1( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force
+ ENDIF
+ ELSE
CALL vor_een( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
- IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN
+ IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN
CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend
- ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN
+ ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN
CALL vor_een( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force
+ ENDIF
ENDIF
END SELECT
!
@@ -320,7 +339,122 @@ CONTAINS
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zx1, zy1, zx2, zy2 ! local scalars
- REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx, zwy, zwt ! 2D workspace
+ REAL(wp), DIMENSION(A2D(1)) :: zwx, zwy, zwt ! 2D workspace
+ REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zwz ! 3D workspace, jpkm1 -> avoid lbc_lnk on jpk that is not defined
+ !!----------------------------------------------------------------------
+ !
+ 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:vor_enT : vorticity term: t-point energy conserving scheme'
+ IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+ ENDIF
+ ENDIF
+ !
+ !
+ ! ! ===============
+ DO jk = 1, jpkm1 ! Horizontal slab
+ ! ! ===============
+ !
+ SELECT CASE( kvor ) !== relative vorticity considered ==!
+ !
+ CASE ( np_RVO , np_CRV ) !* relative vorticity at f-point is used
+ ALLOCATE( zwz(A2D(1)) )
+ DO_2D( 1, 1, 1, 1 )
+ zwz(ji,jj) = ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
+ & - e1u(ji,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)
+ END_2D
+ IF( ln_dynvor_msk ) THEN ! mask relative vorticity
+ DO_2D( 1, 1, 1, 1 )
+ zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
+ END_2D
+ ENDIF
+ !
+ END SELECT
+ !
+ SELECT CASE( kvor ) !== volume weighted vorticity considered ==!
+ !
+ CASE ( np_COR ) !* Coriolis (planetary vorticity)
+ DO_2D( 0, 1, 0, 1 )
+ zwt(ji,jj) = ff_t(ji,jj) * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
+ END_2D
+ CASE ( np_RVO ) !* relative vorticity
+ DO_2D( 0, 1, 0, 1 )
+ zwt(ji,jj) = r1_4 * ( zwz(ji-1,jj ) + zwz(ji,jj ) &
+ & + zwz(ji-1,jj-1) + zwz(ji,jj-1) ) &
+ & * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
+ END_2D
+ CASE ( np_MET ) !* metric term
+ DO_2D( 0, 1, 0, 1 )
+ zwt(ji,jj) = ( ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) &
+ & - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) &
+ & * e3t(ji,jj,jk,Kmm)
+ END_2D
+ CASE ( np_CRV ) !* Coriolis + relative vorticity
+ DO_2D( 0, 1, 0, 1 )
+ zwt(ji,jj) = ( ff_t(ji,jj) + r1_4 * ( zwz(ji-1,jj ) + zwz(ji,jj ) &
+ & + zwz(ji-1,jj-1) + zwz(ji,jj-1) ) ) &
+ & * e1e2t(ji,jj)*e3t(ji,jj,jk,Kmm)
+ END_2D
+ CASE ( np_CME ) !* Coriolis + metric
+ DO_2D( 0, 1, 0, 1 )
+ zwt(ji,jj) = ( ff_t(ji,jj) * e1e2t(ji,jj) &
+ & + ( pv(ji,jj,jk) + pv(ji,jj-1,jk) ) * di_e2u_2(ji,jj) &
+ & - ( pu(ji,jj,jk) + pu(ji-1,jj,jk) ) * dj_e1v_2(ji,jj) ) &
+ & * e3t(ji,jj,jk,Kmm)
+ END_2D
+ CASE DEFAULT ! error
+ CALL ctl_stop('STOP','dyn_vor: wrong value for kvor')
+ END SELECT
+ !
+ ! !== compute and add the vorticity term trend =!
+ DO_2D( 0, 0, 0, 0 )
+ pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + r1_4 * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) &
+ & * ( zwt(ji+1,jj) * ( pv(ji+1,jj,jk) + pv(ji+1,jj-1,jk) ) &
+ & + zwt(ji ,jj) * ( pv(ji ,jj,jk) + pv(ji ,jj-1,jk) ) )
+ !
+ pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) - r1_4 * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) &
+ & * ( zwt(ji,jj+1) * ( pu(ji,jj+1,jk) + pu(ji-1,jj+1,jk) ) &
+ & + zwt(ji,jj ) * ( pu(ji,jj ,jk) + pu(ji-1,jj ,jk) ) )
+ END_2D
+ ! ! ===============
+ END DO ! End of slab
+ ! ! ===============
+ !
+ SELECT CASE( kvor ) ! deallocate zwz if necessary
+ CASE ( np_RVO , np_CRV ) ; DEALLOCATE( zwz )
+ END SELECT
+ !
+ END SUBROUTINE vor_enT
+
+
+ SUBROUTINE vor_enT_hls1( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE vor_enT ***
+ !!
+ !! ** Purpose : Compute the now total vorticity trend and add it to
+ !! the general trend of the momentum equation.
+ !!
+ !! ** Method : Trend evaluated using now fields (centered in time)
+ !! and t-point evaluation of vorticity (planetary and relative).
+ !! conserves the horizontal kinetic energy.
+ !! The general trend of momentum is increased due to the vorticity
+ !! term which is given by:
+ !! voru = 1/bu mj[ ( mi(mj(bf*rvor))+bt*f_t)/e3t mj[vn] ]
+ !! vorv = 1/bv mi[ ( mi(mj(bf*rvor))+bt*f_t)/e3f mj[un] ]
+ !! where rvor is the relative vorticity at f-point
+ !!
+ !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
+ !!----------------------------------------------------------------------
+ INTEGER , INTENT(in ) :: kt ! ocean time-step index
+ INTEGER , INTENT(in ) :: Kmm ! ocean time level index
+ INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
+ REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
+ REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
+ !
+ INTEGER :: ji, jj, jk ! dummy loop indices
+ REAL(wp) :: zx1, zy1, zx2, zy2 ! local scalars
+ REAL(wp), DIMENSION(A2D(1)) :: zwx, zwy, zwt ! 2D workspace
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zwz ! 3D workspace, jpkm1 -> avoid lbc_lnk on jpk that is not defined
!!----------------------------------------------------------------------
!
@@ -336,7 +470,7 @@ CONTAINS
SELECT CASE( kvor ) !== relative vorticity considered ==!
!
CASE ( np_RVO , np_CRV ) !* relative vorticity at f-point is used
- ALLOCATE( zwz(A2D(nn_hls),jpk) )
+ ALLOCATE( zwz(A2D(1),jpk) )
DO jk = 1, jpkm1 ! Horizontal slab
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
zwz(ji,jj,jk) = ( e2v(ji+1,jj) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
@@ -409,7 +543,7 @@ CONTAINS
CASE ( np_RVO , np_CRV ) ; DEALLOCATE( zwz )
END SELECT
!
- END SUBROUTINE vor_enT
+ END SUBROUTINE vor_enT_hls1
SUBROUTINE vor_ene( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
@@ -440,7 +574,7 @@ CONTAINS
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zx1, zy1, zx2, zy2, ze3f, zmsk ! local scalars
- REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx, zwy, zwz ! 2D workspace
+ REAL(wp), DIMENSION(A2D(1)) :: zwx, zwy, zwz ! 2D workspace
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
@@ -573,7 +707,7 @@ CONTAINS
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zuav, zvau, ze3f, zmsk ! local scalars
- REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx, zwy, zwz ! 2D workspace
+ REAL(wp), DIMENSION(A2D(1)) :: zwx, zwy, zwz ! 2D workspace
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
@@ -705,16 +839,176 @@ CONTAINS
INTEGER :: ierr ! local integer
REAL(wp) :: zua, zva ! local scalars
REAL(wp) :: zmsk, ze3f ! local scalars
- REAL(wp), DIMENSION(A2D(nn_hls)) :: z1_e3f
+ REAL(wp), DIMENSION(A2D(1)) :: z1_e3f
#if defined key_loop_fusion
REAL(wp) :: ztne, ztnw, ztnw_ip1, ztse, ztse_jp1, ztsw_jp1, ztsw_ip1
REAL(wp) :: zwx, zwx_im1, zwx_jp1, zwx_im1_jp1
REAL(wp) :: zwy, zwy_ip1, zwy_jm1, zwy_ip1_jm1
#else
- REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx , zwy
- REAL(wp), DIMENSION(A2D(nn_hls)) :: ztnw, ztne, ztsw, ztse
+ REAL(wp), DIMENSION(A2D(1)) :: zwx , zwy
+ REAL(wp), DIMENSION(A2D(1)) :: ztnw, ztne, ztsw, ztse
#endif
- REAL(wp), DIMENSION(A2D(nn_hls),jpkm1) :: zwz ! 3D workspace, jpkm1 -> jpkm1 -> avoid lbc_lnk on jpk that is not defined
+ REAL(wp), DIMENSION(A2D(1)) :: zwz ! 3D workspace, jpkm1 -> jpkm1 -> avoid lbc_lnk on jpk that is not defined
+ !!----------------------------------------------------------------------
+ !
+ 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:vor_een : vorticity term: energy and enstrophy conserving scheme'
+ IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+ ENDIF
+ ENDIF
+ !
+ ! ! ===============
+ DO jk = 1, jpkm1 ! Horizontal slab
+ ! ! ===============
+ !
+#if defined key_qco || defined key_linssh
+ DO_2D( 1, 1, 1, 1 ) ! == reciprocal of e3 at F-point (key_qco)
+ z1_e3f(ji,jj) = 1._wp / e3f_vor(ji,jj,jk)
+ END_2D
+#else
+ SELECT CASE( nn_e3f_typ ) ! == reciprocal of e3 at F-point
+ CASE ( 0 ) ! original formulation (masked averaging of e3t divided by 4)
+ DO_2D( 1, 1, 1, 1 )
+ ! round brackets added to fix the order of floating point operations
+ ! needed to ensure halo 1 - halo 2 compatibility
+ ze3f = ( (e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) &
+ & + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)) &
+ & + (e3t(ji ,jj ,jk,Kmm)*tmask(ji ,jj ,jk) &
+ & + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk)) )
+ IF( ze3f /= 0._wp ) THEN ; z1_e3f(ji,jj) = 4._wp / ze3f
+ ELSE ; z1_e3f(ji,jj) = 0._wp
+ ENDIF
+ END_2D
+ CASE ( 1 ) ! new formulation (masked averaging of e3t divided by the sum of mask)
+ DO_2D( 1, 1, 1, 1 )
+ ! round brackets added to fix the order of floating point operations
+ ! needed to ensure halo 1 - halo 2 compatibility
+ ze3f = ( (e3t(ji ,jj+1,jk,Kmm)*tmask(ji ,jj+1,jk) &
+ & + e3t(ji+1,jj+1,jk,Kmm)*tmask(ji+1,jj+1,jk)) &
+ & + (e3t(ji ,jj ,jk,Kmm)*tmask(ji ,jj ,jk) &
+ & + e3t(ji+1,jj ,jk,Kmm)*tmask(ji+1,jj ,jk)) )
+ zmsk = ( tmask(ji,jj+1,jk) + tmask(ji+1,jj+1,jk) &
+ & + tmask(ji,jj ,jk) + tmask(ji+1,jj ,jk) )
+ IF( ze3f /= 0._wp ) THEN ; z1_e3f(ji,jj) = zmsk / ze3f
+ ELSE ; z1_e3f(ji,jj) = 0._wp
+ ENDIF
+ END_2D
+ END SELECT
+#endif
+ !
+ SELECT CASE( kvor ) !== vorticity considered ==!
+ !
+ CASE ( np_COR ) !* Coriolis (planetary vorticity)
+ DO_2D( 1, 1, 1, 1 )
+ zwz(ji,jj) = ff_f(ji,jj) * z1_e3f(ji,jj)
+ END_2D
+ CASE ( np_RVO ) !* relative vorticity
+ DO_2D( 1, 1, 1, 1 )
+ zwz(ji,jj) = ( e2v(ji+1,jj ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
+ & - e1u(ji ,jj+1) * pu(ji,jj+1,jk) + e1u(ji,jj) * pu(ji,jj,jk) ) * r1_e1e2f(ji,jj)*z1_e3f(ji,jj)
+ END_2D
+ IF( ln_dynvor_msk ) THEN ! mask the relative vorticity
+ DO_2D( 1, 1, 1, 1 )
+ zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
+ END_2D
+ ENDIF
+ CASE ( np_MET ) !* metric term
+ DO_2D( 1, 1, 1, 1 )
+ zwz(ji,jj) = ( ( pv(ji+1,jj,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
+ & - ( pu(ji,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj)
+ END_2D
+ CASE ( np_CRV ) !* Coriolis + relative vorticity
+ DO_2D( 1, 1, 1, 1 )
+ ! round brackets added to fix the order of floating point operations
+ ! needed to ensure halo 1 - halo 2 compatibility
+ zwz(ji,jj) = ( ff_f(ji,jj) + ( ( e2v(ji+1,jj ) * pv(ji+1,jj,jk) - e2v(ji,jj) * pv(ji,jj,jk) &
+ & ) & ! bracket for halo 1 - halo 2 compatibility
+ & - ( e1u(ji ,jj+1) * pu(ji,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk) &
+ & ) & ! bracket for halo 1 - halo 2 compatibility
+ & ) * r1_e1e2f(ji,jj) ) * z1_e3f(ji,jj)
+ END_2D
+ IF( ln_dynvor_msk ) THEN ! mask the relative vorticity
+ DO_2D( 1, 1, 1, 1 )
+ zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
+ END_2D
+ ENDIF
+ CASE ( np_CME ) !* Coriolis + metric
+ DO_2D( 1, 1, 1, 1 )
+ zwz(ji,jj) = ( ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
+ & - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj) ) * z1_e3f(ji,jj)
+ END_2D
+ CASE DEFAULT ! error
+ CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' )
+ END SELECT
+ !
+ ! !== horizontal fluxes ==!
+ DO_2D( 1, 1, 1, 1 )
+ zwx(ji,jj) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * pu(ji,jj,jk)
+ zwy(ji,jj) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pv(ji,jj,jk)
+ END_2D
+ !
+ ! !== compute and add the vorticity term trend =!
+ DO_2D( 0, 1, 0, 1 )
+ ztne(ji,jj) = zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1)
+ ztnw(ji,jj) = zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj )
+ ztse(ji,jj) = zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1)
+ ztsw(ji,jj) = zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj )
+ END_2D
+ !
+ DO_2D( 0, 0, 0, 0 )
+ zua = + r1_12 * r1_e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) &
+ & + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
+ zva = - r1_12 * r1_e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) &
+ & + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) )
+ pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
+ pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
+ END_2D
+ END DO
+ ! ! ===============
+ ! ! End of slab
+ ! ! ===============
+ END SUBROUTINE vor_een
+
+
+ SUBROUTINE vor_een_hls1( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE vor_een ***
+ !!
+ !! ** Purpose : Compute the now total vorticity trend and add it to
+ !! the general trend of the momentum equation.
+ !!
+ !! ** Method : Trend evaluated using now fields (centered in time)
+ !! and the Arakawa and Lamb (1980) flux form formulation : conserves
+ !! both the horizontal kinetic energy and the potential enstrophy
+ !! when horizontal divergence is zero (see the NEMO documentation)
+ !! Add this trend to the general momentum trend (pu_rhs,pv_rhs).
+ !!
+ !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
+ !!
+ !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
+ !!----------------------------------------------------------------------
+ INTEGER , INTENT(in ) :: kt ! ocean time-step index
+ INTEGER , INTENT(in ) :: Kmm ! ocean time level index
+ INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
+ REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
+ REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
+ !
+ INTEGER :: ji, jj, jk ! dummy loop indices
+ INTEGER :: ierr ! local integer
+ REAL(wp) :: zua, zva ! local scalars
+ REAL(wp) :: zmsk, ze3f ! local scalars
+ REAL(wp), DIMENSION(A2D(1)) :: z1_e3f
+#if defined key_loop_fusion
+ REAL(wp) :: ztne, ztnw, ztnw_ip1, ztse, ztse_jp1, ztsw_jp1, ztsw_ip1
+ REAL(wp) :: zwx, zwx_im1, zwx_jp1, zwx_im1_jp1
+ REAL(wp) :: zwy, zwy_ip1, zwy_jm1, zwy_ip1_jm1
+#else
+ REAL(wp), DIMENSION(A2D(1)) :: zwx , zwy
+ REAL(wp), DIMENSION(A2D(1)) :: ztnw, ztne, ztsw, ztse
+#endif
+ REAL(wp), DIMENSION(A2D(1),jpkm1) :: zwz ! 3D workspace, jpkm1 -> jpkm1 -> avoid lbc_lnk on jpk that is not defined
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
@@ -874,7 +1168,7 @@ CONTAINS
! ! ===============
! ! End of slab
! ! ===============
- END SUBROUTINE vor_een
+ END SUBROUTINE vor_een_hls1
SUBROUTINE vor_eeT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
@@ -904,9 +1198,129 @@ CONTAINS
INTEGER :: ierr ! local integer
REAL(wp) :: zua, zva ! local scalars
REAL(wp) :: zmsk, z1_e3t ! local scalars
- REAL(wp), DIMENSION(A2D(nn_hls)) :: zwx , zwy
- REAL(wp), DIMENSION(A2D(nn_hls)) :: ztnw, ztne, ztsw, ztse
- REAL(wp), DIMENSION(A2D(nn_hls),jpkm1) :: zwz ! 3D workspace, avoid lbc_lnk on jpk that is not defined
+ REAL(wp), DIMENSION(A2D(1)) :: zwx , zwy
+ REAL(wp), DIMENSION(A2D(1)) :: ztnw, ztne, ztsw, ztse
+ REAL(wp), DIMENSION(A2D(1)) :: zwz ! 3D workspace, avoid lbc_lnk on jpk that is not defined
+ !!----------------------------------------------------------------------
+ !
+ 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:vor_eeT : vorticity term: energy and enstrophy conserving scheme'
+ IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
+ ENDIF
+ ENDIF
+ !
+ ! ! ===============
+ DO jk = 1, jpkm1 ! Horizontal slab
+ ! ! ===============
+ !
+ !
+ SELECT CASE( kvor ) !== vorticity considered ==!
+ CASE ( np_COR ) !* Coriolis (planetary vorticity)
+ DO_2D( 1, 1, 1, 1 )
+ zwz(ji,jj) = ff_f(ji,jj)
+ END_2D
+ CASE ( np_RVO ) !* relative vorticity
+ DO_2D( 1, 1, 1, 1 )
+ ! round brackets added to fix the order of floating point operations
+ ! needed to ensure halo 1 - halo 2 compatibility
+ zwz(ji,jj) = ( (e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk)) &
+ & - (e1u(ji ,jj+1) * pu(ji ,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk)) ) &
+ & * r1_e1e2f(ji,jj)
+ END_2D
+ IF( ln_dynvor_msk ) THEN ! mask the relative vorticity
+ DO_2D( 1, 1, 1, 1 )
+ zwz(ji,jj) = zwz(ji,jj) * fmask(ji,jj,jk)
+ END_2D
+ ENDIF
+ CASE ( np_MET ) !* metric term
+ DO_2D( 1, 1, 1, 1 )
+ zwz(ji,jj) = ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
+ & - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
+ END_2D
+ CASE ( np_CRV ) !* Coriolis + relative vorticity
+ DO_2D( 1, 1, 1, 1 )
+ ! round brackets added to fix the order of floating point operations
+ ! needed to ensure halo 1 - halo 2 compatibility
+ zwz(ji,jj) = ( ff_f(ji,jj) + ( (e2v(ji+1,jj ) * pv(ji+1,jj ,jk) - e2v(ji,jj) * pv(ji,jj,jk)) &
+ & - (e1u(ji ,jj+1) * pu(ji ,jj+1,jk) - e1u(ji,jj) * pu(ji,jj,jk)) ) &
+ & * r1_e1e2f(ji,jj) )
+ END_2D
+ IF( ln_dynvor_msk ) THEN ! mask the relative vorticity
+ DO_2D( 1, 1, 1, 1 )
+ zwz(ji,jj) = ( zwz(ji,jj) - ff_f(ji,jj) ) * fmask(ji,jj,jk) + ff_f(ji,jj)
+ END_2D
+ ENDIF
+ CASE ( np_CME ) !* Coriolis + metric
+ DO_2D( 1, 1, 1, 1 )
+ zwz(ji,jj) = ff_f(ji,jj) + ( pv(ji+1,jj ,jk) + pv(ji,jj,jk) ) * di_e2v_2e1e2f(ji,jj) &
+ & - ( pu(ji ,jj+1,jk) + pu(ji,jj,jk) ) * dj_e1u_2e1e2f(ji,jj)
+ END_2D
+ CASE DEFAULT ! error
+ CALL ctl_stop('STOP','dyn_vor: wrong value for kvor' )
+ END SELECT
+ !
+ !
+ ! !== horizontal fluxes ==!
+ DO_2D( 1, 1, 1, 1 )
+ zwx(ji,jj) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * pu(ji,jj,jk)
+ zwy(ji,jj) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pv(ji,jj,jk)
+ END_2D
+ !
+ ! !== compute and add the vorticity term trend =!
+ DO_2D( 0, 1, 0, 1 )
+ z1_e3t = 1._wp / e3t(ji,jj,jk,Kmm)
+ ztne(ji,jj) = ( zwz(ji-1,jj ) + zwz(ji ,jj ) + zwz(ji ,jj-1) ) * z1_e3t
+ ztnw(ji,jj) = ( zwz(ji-1,jj-1) + zwz(ji-1,jj ) + zwz(ji ,jj ) ) * z1_e3t
+ ztse(ji,jj) = ( zwz(ji ,jj ) + zwz(ji ,jj-1) + zwz(ji-1,jj-1) ) * z1_e3t
+ ztsw(ji,jj) = ( zwz(ji ,jj-1) + zwz(ji-1,jj-1) + zwz(ji-1,jj ) ) * z1_e3t
+ END_2D
+ !
+ DO_2D( 0, 0, 0, 0 )
+ zua = + r1_12 * r1_e1u(ji,jj) * ( ztne(ji,jj ) * zwy(ji ,jj ) + ztnw(ji+1,jj) * zwy(ji+1,jj ) &
+ & + ztse(ji,jj ) * zwy(ji ,jj-1) + ztsw(ji+1,jj) * zwy(ji+1,jj-1) )
+ zva = - r1_12 * r1_e2v(ji,jj) * ( ztsw(ji,jj+1) * zwx(ji-1,jj+1) + ztse(ji,jj+1) * zwx(ji ,jj+1) &
+ & + ztnw(ji,jj ) * zwx(ji-1,jj ) + ztne(ji,jj ) * zwx(ji ,jj ) )
+ pu_rhs(ji,jj,jk) = pu_rhs(ji,jj,jk) + zua
+ pv_rhs(ji,jj,jk) = pv_rhs(ji,jj,jk) + zva
+ END_2D
+ ! ! ===============
+ END DO ! End of slab
+ ! ! ===============
+ END SUBROUTINE vor_eeT
+
+
+ SUBROUTINE vor_eeT_hls1( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE vor_eeT ***
+ !!
+ !! ** Purpose : Compute the now total vorticity trend and add it to
+ !! the general trend of the momentum equation.
+ !!
+ !! ** Method : Trend evaluated using now fields (centered in time)
+ !! and the Arakawa and Lamb (1980) vector form formulation using
+ !! a modified version of Arakawa and Lamb (1980) scheme (see vor_een).
+ !! The change consists in
+ !! Add this trend to the general momentum trend (pu_rhs,pv_rhs).
+ !!
+ !! ** Action : - Update (pu_rhs,pv_rhs) with the now vorticity term trend
+ !!
+ !! References : Arakawa and Lamb 1980, Mon. Wea. Rev., 109, 18-36
+ !!----------------------------------------------------------------------
+ INTEGER , INTENT(in ) :: kt ! ocean time-step index
+ INTEGER , INTENT(in ) :: Kmm ! ocean time level index
+ INTEGER , INTENT(in ) :: kvor ! total, planetary, relative, or metric
+ REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu, pv ! now velocities
+ REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pu_rhs, pv_rhs ! total v-trend
+ !
+ INTEGER :: ji, jj, jk ! dummy loop indices
+ INTEGER :: ierr ! local integer
+ REAL(wp) :: zua, zva ! local scalars
+ REAL(wp) :: zmsk, z1_e3t ! local scalars
+ REAL(wp), DIMENSION(A2D(1)) :: zwx , zwy
+ REAL(wp), DIMENSION(A2D(1)) :: ztnw, ztne, ztsw, ztse
+ REAL(wp), DIMENSION(A2D(1),jpkm1) :: zwz ! 3D workspace, avoid lbc_lnk on jpk that is not defined
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
@@ -1003,7 +1417,7 @@ CONTAINS
! ! ===============
END DO ! End of slab
! ! ===============
- END SUBROUTINE vor_eeT
+ END SUBROUTINE vor_eeT_hls1
SUBROUTINE dyn_vor_init
--
GitLab