diff --git a/src/OCE/DYN/dynadv.F90 b/src/OCE/DYN/dynadv.F90
index bda2c3a4f806cd4fab1bcc5a6b3bc1ff89f77dd1..6d59583ffb868f71fdaf9900ca231e9419030b8d 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 7fd7f65f5c19ca959bd2437cb7352cf5bb1fd1b6..07dd931129ed0b4474092915e3b6c83abcaf2ff5 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 645f6fa800606d41375281e38c4bc077f1a848e8..8ab908124d5b85beb120d5726d2d6ec12b8d28f7 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 d751899d0de9266477ed7f9ecc26941906cba84d..10e4134b8a2d1efd9ab751ce09b5b8b82c23a6f2 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 03dda78ceb4f6ce993419062208a9240e3ae7d75..b7a84124a27df69f6ea82b1622683c7a2be8470b 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
diff --git a/src/OCE/DYN/dynzdf.F90 b/src/OCE/DYN/dynzdf.F90
index 98c65e4604e951f0119fbd4767381487ec3596a0..fed10f23182fea2131de19cd81d668d71ed3aec5 100644
--- a/src/OCE/DYN/dynzdf.F90
+++ b/src/OCE/DYN/dynzdf.F90
@@ -6,6 +6,7 @@ MODULE dynzdf
!! History : 1.0 ! 2005-11 (G. Madec) Original code
!! 3.3 ! 2010-10 (C. Ethe, G. Madec) reorganisation of initialisation phase
!! 4.0 ! 2017-06 (G. Madec) remove the explicit time-stepping option + avm at t-point
+ !! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
@@ -79,7 +80,7 @@ CONTAINS
REAL(wp) :: zWu , zWv ! - -
REAL(wp) :: zWui, zWvi ! - -
REAL(wp) :: zWus, zWvs ! - -
- REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwi, zwd, zws ! 3D workspace
+ REAL(wp), DIMENSION(A1Di(0),jpk) :: zwi, zwd, zws ! 2D workspace
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdu, ztrdv ! - -
!!---------------------------------------------------------------------
!
@@ -105,315 +106,329 @@ CONTAINS
ztrdv(:,:,:) = pvv(:,:,:,Krhs)
ENDIF
!
- ! !== RHS: Leap-Frog time stepping on all trends but the vertical mixing ==! (put in puu(:,:,:,Kaa),pvv(:,:,:,Kaa))
- !
- ! ! time stepping except vertical diffusion
- IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity
- DO_3D( 0, 0, 0, 0, 1, jpkm1 )
- puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kbb) + rDt * puu(ji,jj,jk,Krhs) ) * umask(ji,jj,jk)
- pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kbb) + rDt * pvv(ji,jj,jk,Krhs) ) * vmask(ji,jj,jk)
- END_3D
- ELSE ! applied on thickness weighted velocity
- DO_3D( 0, 0, 0, 0, 1, jpkm1 )
- puu(ji,jj,jk,Kaa) = ( e3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb ) &
- & + rDt * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Krhs) ) &
- & / e3u(ji,jj,jk,Kaa) * umask(ji,jj,jk)
- pvv(ji,jj,jk,Kaa) = ( e3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb ) &
- & + rDt * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Krhs) ) &
- & / e3v(ji,jj,jk,Kaa) * vmask(ji,jj,jk)
- END_3D
- ENDIF
- ! ! add top/bottom friction
- ! With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom.
- ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does
- ! not lead to the effective stress seen over the whole barotropic loop.
- ! G. Madec : in linear free surface, e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) = e3u_0, so systematic use of e3u(:,:,:,Kaa)
- IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
- DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! remove barotropic velocities
- puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - uu_b(ji,jj,Kaa) ) * umask(ji,jj,jk)
- pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - vv_b(ji,jj,Kaa) ) * vmask(ji,jj,jk)
- END_3D
- DO_2D( 0, 0, 0, 0 ) ! Add bottom/top stress due to barotropic component only
- iku = mbku(ji,jj) ! ocean bottom level at u- and v-points
- ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
- puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + zDt_2 *( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) &
- & / e3u(ji,jj,iku,Kaa)
- pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + zDt_2 *( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) &
- & / e3v(ji,jj,ikv,Kaa)
- END_2D
- IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! Ocean cavities (ISF)
- DO_2D( 0, 0, 0, 0 )
- iku = miku(ji,jj) ! top ocean level at u- and v-points
- ikv = mikv(ji,jj) ! (first wet ocean u- and v-points)
- puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + zDt_2 *( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) &
+ ! ! ================= !
+ DO_1Dj( 0, 0 ) ! i-k slices loop !
+ ! ! ================= !
+ !
+ ! !== RHS: Leap-Frog time stepping on all trends but the vertical mixing ==! (put in puu(:,:,:,Kaa),pvv(:,:,:,Kaa))
+ !
+ ! ! time stepping except vertical diffusion
+ IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! applied on velocity
+ DO_2Dik( 0, 0, 1, jpkm1, 1 )
+ puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kbb) + rDt * puu(ji,jj,jk,Krhs) ) * umask(ji,jj,jk)
+ pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kbb) + rDt * pvv(ji,jj,jk,Krhs) ) * vmask(ji,jj,jk)
+ END_2D
+ ELSE ! applied on thickness weighted velocity
+ DO_2Dik( 0, 0, 1, jpkm1, 1 )
+ puu(ji,jj,jk,Kaa) = ( e3u(ji,jj,jk,Kbb) * puu(ji,jj,jk,Kbb ) &
+ & + rDt * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Krhs) ) &
+ & / e3u(ji,jj,jk,Kaa) * umask(ji,jj,jk)
+ pvv(ji,jj,jk,Kaa) = ( e3v(ji,jj,jk,Kbb) * pvv(ji,jj,jk,Kbb ) &
+ & + rDt * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Krhs) ) &
+ & / e3v(ji,jj,jk,Kaa) * vmask(ji,jj,jk)
+ END_2D
+ ENDIF
+ ! ! add top/bottom friction
+ ! With split-explicit free surface, barotropic stress is treated explicitly Update velocities at the bottom.
+ ! J. Chanut: The bottom stress is computed considering after barotropic velocities, which does
+ ! not lead to the effective stress seen over the whole barotropic loop.
+ ! G. Madec : in linear free surface, e3u(:,:,:,Kaa) = e3u(:,:,:,Kmm) = e3u_0, so systematic use of e3u(:,:,:,Kaa)
+ IF( ln_drgimp .AND. ln_dynspg_ts ) THEN
+ DO_2Dik( 0, 0, 1, jpkm1, 1 ) ! remove barotropic velocities
+ puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - uu_b(ji,jj,Kaa) ) * umask(ji,jj,jk)
+ pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - vv_b(ji,jj,Kaa) ) * vmask(ji,jj,jk)
+ END_2D
+ DO_1Di( 0, 0 ) ! Add bottom/top stress due to barotropic component only
+ iku = mbku(ji,jj) ! ocean bottom level at u- and v-points
+ ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
+ puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + zDt_2 *( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) &
& / e3u(ji,jj,iku,Kaa)
- pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + zDt_2 *( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) &
+ pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + zDt_2 *( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) &
& / e3v(ji,jj,ikv,Kaa)
- END_2D
- END IF
- ENDIF
- !
- ! !== Vertical diffusion on u ==!
- !
- ! !* Matrix construction
- IF( ln_zad_Aimp ) THEN !!
- SELECT CASE( nldf_dyn )
- CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu)
- DO_3D( 0, 0, 0, 0, 1, jpkm1 )
- z1_e3ua = 1._wp / e3u(ji,jj,jk,Kaa) ! after scale factor at U-point
- zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) &
- & / e3uw(ji,jj,jk ,Kmm) * z1_e3ua * wumask(ji,jj,jk )
- zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) &
- & / e3uw(ji,jj,jk+1,Kmm) * z1_e3ua * wumask(ji,jj,jk+1)
- zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) * z1_e3ua
- zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) * z1_e3ua
- zwi(ji,jj,jk) = zzwi + zDt_2 * MIN( zWui, 0._wp )
- zws(ji,jj,jk) = zzws - zDt_2 * MAX( zWus, 0._wp )
- zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zDt_2 * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
- END_3D
- CASE DEFAULT ! iso-level lateral mixing
- DO_3D( 0, 0, 0, 0, 1, jpkm1 )
- z1_e3ua = 1._wp / e3u(ji,jj,jk,Kaa) ! after scale factor at U-point
- zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) &
- & / e3uw(ji,jj,jk ,Kmm) * z1_e3ua * wumask(ji,jj,jk )
- zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) &
- & / e3uw(ji,jj,jk+1,Kmm) * z1_e3ua * wumask(ji,jj,jk+1)
- zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) * z1_e3ua
- zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) * z1_e3ua
- zwi(ji,jj,jk) = zzwi + zDt_2 * MIN( zWui, 0._wp )
- zws(ji,jj,jk) = zzws - zDt_2 * MAX( zWus, 0._wp )
- zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zDt_2 * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
- END_3D
- END SELECT
- DO_2D( 0, 0, 0, 0 ) !* Surface boundary conditions
- zwi(ji,jj,1) = 0._wp
- zzws = - zDt_2 * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) &
- & / ( e3u(ji,jj,1,Kaa) * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2)
- zWus = ( wi(ji ,jj,2) + wi(ji+1,jj,2) ) / e3u(ji,jj,1,Kaa)
- zws(ji,jj,1 ) = zzws - zDt_2 * MAX( zWus, 0._wp )
- zwd(ji,jj,1 ) = 1._wp - zzws - zDt_2 * ( MIN( zWus, 0._wp ) )
- END_2D
- ELSE
- SELECT CASE( nldf_dyn )
- CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu)
- DO_3D( 0, 0, 0, 0, 1, jpkm1 )
- zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) &
- & / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk )
- zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) &
- & / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
- zwi(ji,jj,jk) = zzwi
- zws(ji,jj,jk) = zzws
- zwd(ji,jj,jk) = 1._wp - zzwi - zzws
- END_3D
- CASE DEFAULT ! iso-level lateral mixing
- DO_3D( 0, 0, 0, 0, 1, jpkm1 )
- zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) &
- & / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk )
- zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) &
- & / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
- zwi(ji,jj,jk) = zzwi
- zws(ji,jj,jk) = zzws
- zwd(ji,jj,jk) = 1._wp - zzwi - zzws
- END_3D
- END SELECT
- DO_2D( 0, 0, 0, 0 ) !* Surface boundary conditions
- zwi(ji,jj,1) = 0._wp
- zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
- END_2D
- ENDIF
- !
- !
- ! !== Apply semi-implicit bottom friction ==!
- !
- ! Only needed for semi-implicit bottom friction setup. The explicit
- ! bottom friction has been included in "u(v)a" which act as the R.H.S
- ! column vector of the tri-diagonal matrix equation
- !
- IF ( ln_drgimp ) THEN ! implicit bottom friction
- DO_2D( 0, 0, 0, 0 )
- iku = mbku(ji,jj) ! ocean bottom level at u- and v-points
- zwd(ji,jj,iku) = zwd(ji,jj,iku) - zDt_2 *( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) &
- & / e3u(ji,jj,iku,Kaa)
- END_2D
- IF ( ln_isfcav.OR.ln_drgice_imp ) THEN ! top friction (always implicit)
- DO_2D( 0, 0, 0, 0 )
- !!gm top Cd is masked (=0 outside cavities) no need of test on mik>=2 ==>> it has been suppressed
- iku = miku(ji,jj) ! ocean top level at u- and v-points
- zwd(ji,jj,iku) = zwd(ji,jj,iku) - zDt_2 *( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) &
+ END_1D
+ IF( ln_isfcav.OR.ln_drgice_imp ) THEN ! Ocean cavities (ISF)
+ DO_1Di( 0, 0 )
+ iku = miku(ji,jj) ! top ocean level at u- and v-points
+ ikv = mikv(ji,jj) ! (first wet ocean u- and v-points)
+ puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + zDt_2 *( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) &
+ & / e3u(ji,jj,iku,Kaa)
+ pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + zDt_2 *( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) &
+ & / e3v(ji,jj,ikv,Kaa)
+ END_1D
+ END IF
+ ENDIF
+ !
+ ! !== Vertical diffusion on u ==!
+ !
+ !
+ ! !* Matrix construction
+ IF( ln_zad_Aimp ) THEN !- including terms associated with partly implicit vertical advection
+ SELECT CASE( nldf_dyn )
+ CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu)
+ DO_2Dik( 0, 0, 1, jpkm1, 1 )
+ z1_e3ua = 1._wp / e3u(ji,jj,jk,Kaa) ! after scale factor at U-point
+ zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) &
+ & / e3uw(ji,jj,jk ,Kmm) * z1_e3ua * wumask(ji,jj,jk )
+ zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) &
+ & / e3uw(ji,jj,jk+1,Kmm) * z1_e3ua * wumask(ji,jj,jk+1)
+ zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) * z1_e3ua
+ zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) * z1_e3ua
+ zwi(ji,jk) = zzwi + zDt_2 * MIN( zWui, 0._wp )
+ zws(ji,jk) = zzws - zDt_2 * MAX( zWus, 0._wp )
+ zwd(ji,jk) = 1._wp - zzwi - zzws + zDt_2 * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
+ END_2D
+ CASE DEFAULT ! iso-level lateral mixing
+ DO_2Dik( 0, 0, 1, jpkm1, 1 )
+ z1_e3ua = 1._wp / e3u(ji,jj,jk,Kaa) ! after scale factor at U-point
+ zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) &
+ & / e3uw(ji,jj,jk ,Kmm) * z1_e3ua * wumask(ji,jj,jk )
+ zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) &
+ & / e3uw(ji,jj,jk+1,Kmm) * z1_e3ua * wumask(ji,jj,jk+1)
+ zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) * z1_e3ua
+ zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) * z1_e3ua
+ zwi(ji,jk) = zzwi + zDt_2 * MIN( zWui, 0._wp )
+ zws(ji,jk) = zzws - zDt_2 * MAX( zWus, 0._wp )
+ zwd(ji,jk) = 1._wp - zzwi - zzws + zDt_2 * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
+ END_2D
+ END SELECT
+ !
+ zwi(:,1) = 0._wp
+ DO_1Di( 0, 0 ) !* Surface boundary conditions
+ zwi(ji,1) = 0._wp
+ zzws = - zDt_2 * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) &
+ & / ( e3u(ji,jj,1,Kaa) * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2)
+ zWus = ( wi(ji ,jj,2) + wi(ji+1,jj,2) ) / e3u(ji,jj,1,Kaa)
+ zws(ji,1) = zzws - zDt_2 * MAX( zWus, 0._wp )
+ zwd(ji,1) = 1._wp - zzws - zDt_2 * ( MIN( zWus, 0._wp ) )
+ END_1D
+ ELSE !- only vertical diffusive terms
+ SELECT CASE( nldf_dyn )
+ CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu)
+ DO_2Dik( 0, 0, 1, jpkm1, 1 )
+ zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) &
+ & / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk )
+ zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) &
+ & / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
+ zwi(ji,jk) = zzwi
+ zws(ji,jk) = zzws
+ zwd(ji,jk) = 1._wp - zzwi - zzws
+ END_2D
+ CASE DEFAULT ! iso-level lateral mixing
+ DO_2Dik( 0, 0, 1, jpkm1, 1 )
+ zzwi = - zDt_2 * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) &
+ & / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk )
+ zzws = - zDt_2 * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) &
+ & / ( e3u(ji,jj,jk,Kaa) * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
+ zwi(ji,jk) = zzwi
+ zws(ji,jk) = zzws
+ zwd(ji,jk) = 1._wp - zzwi - zzws
+ END_2D
+ END SELECT
+ !
+ zwi(:,1) = 0._wp
+ DO_1Di( 0, 0 ) !* Surface boundary conditions
+ zwd(ji,1) = 1._wp - zws(ji,1)
+ END_1D
+ ENDIF
+ !
+ !
+ ! !== Apply semi-implicit bottom friction ==!
+ !
+ ! Only needed for semi-implicit bottom friction setup. The explicit
+ ! bottom friction has been included in "u(v)a" which act as the R.H.S
+ ! column vector of the tri-diagonal matrix equation
+ !
+ IF ( ln_drgimp ) THEN ! implicit bottom friction
+ DO_1Di( 0, 0 )
+ iku = mbku(ji,jj) ! ocean bottom level at u- and v-points
+ zwd(ji,iku) = zwd(ji,iku) - zDt_2 *( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) &
& / e3u(ji,jj,iku,Kaa)
- END_2D
- END IF
- ENDIF
- !
- ! Matrix inversion starting from the first level
- !-----------------------------------------------------------------------
- ! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
- !
- ! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
- ! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
- ! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
- ! ( ... )( ... ) ( ... )
- ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
- !
- ! m is decomposed in the product of an upper and a lower triangular matrix
- ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
- ! The solution (the after velocity) is in puu(:,:,:,Kaa)
- !-----------------------------------------------------------------------
- !
- DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
- zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
- END_3D
- !
- DO_2D( 0, 0, 0, 0 ) !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==!
+ END_1D
+ IF ( ln_isfcav.OR.ln_drgice_imp ) THEN ! top friction (always implicit)
+ DO_1Di( 0, 0 )
+ !!gm top Cd is masked (=0 outside cavities) no need of test on mik>=2 ==>> it has been suppressed
+ iku = miku(ji,jj) ! ocean top level at u- and v-points
+ zwd(ji,iku) = zwd(ji,iku) - zDt_2 *( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) &
+ & / e3u(ji,jj,iku,Kaa)
+ END_1D
+ ENDIF
+ ENDIF
+ !
+ ! Matrix inversion starting from the first level
+ !-----------------------------------------------------------------------
+ ! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
+ !
+ ! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
+ ! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
+ ! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
+ ! ( ... )( ... ) ( ... )
+ ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
+ !
+ ! m is decomposed in the product of an upper and a lower triangular matrix
+ ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
+ ! The solution (the after velocity) is in puu(:,:,:,Kaa)
+ !-----------------------------------------------------------------------
+ !
+ DO_2Dik( 0, 0, 2, jpkm1, 1 ) !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
+ zwd(ji,jk) = zwd(ji,jk) - zwi(ji,jk) * zws(ji,jk-1) / zwd(ji,jk-1)
+ END_2D
+ !
+ DO_1Di( 0, 0 ) !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==!
#if defined key_RK3
- ! ! RK3: use only utau (not utau_b)
- puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * utauU(ji,jj) &
- & / ( e3u(ji,jj,1,Kaa) * rho0 ) * umask(ji,jj,1)
+ ! ! RK3: use only utau (not utau_b)
+ puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * utauU(ji,jj) &
+ & / ( e3u(ji,jj,1,Kaa) * rho0 ) * umask(ji,jj,1)
#else
- puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + zDt_2 * ( utau_b(ji,jj) + utauU(ji,jj) ) &
- & / ( e3u(ji,jj,1,Kaa) * rho0 ) * umask(ji,jj,1)
+ ! ! MLF: average of utau and utau_b
+ puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + zDt_2 * ( utau_b(ji,jj) + utauU(ji,jj) ) &
+ & / ( e3u(ji,jj,1,Kaa) * rho0 ) * umask(ji,jj,1)
#endif
- END_2D
- DO_3D( 0, 0, 0, 0, 2, jpkm1 )
- puu(ji,jj,jk,Kaa) = puu(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * puu(ji,jj,jk-1,Kaa)
- END_3D
- !
- DO_2D( 0, 0, 0, 0 ) !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==!
- puu(ji,jj,jpkm1,Kaa) = puu(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1)
- END_2D
- DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 )
- puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - zws(ji,jj,jk) * puu(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk)
- END_3D
- !
- ! !== Vertical diffusion on v ==!
- !
- ! !* Matrix construction
- IF( ln_zad_Aimp ) THEN !!
- SELECT CASE( nldf_dyn )
- CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzv)
- DO_3D( 0, 0, 0, 0, 1, jpkm1 )
- z1_e3va = 1._wp / e3v(ji,jj,jk,Kaa) ! after scale factor at V-point
- zzwi = - zDt_2 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) &
- & / e3vw(ji,jj,jk ,Kmm) * z1_e3va * wvmask(ji,jj,jk )
- zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) &
- & / e3vw(ji,jj,jk+1,Kmm) * z1_e3va * wvmask(ji,jj,jk+1)
- zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) * z1_e3va
- zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * z1_e3va
- zwi(ji,jj,jk) = zzwi + zDt_2 * MIN( zWvi, 0._wp )
- zws(ji,jj,jk) = zzws - zDt_2 * MAX( zWvs, 0._wp )
- zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zDt_2 * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
- END_3D
- CASE DEFAULT ! iso-level lateral mixing
- DO_3D( 0, 0, 0, 0, 1, jpkm1 )
- z1_e3va = 1._wp / e3v(ji,jj,jk,Kaa) ! after scale factor at V-point
- zzwi = - zDt_2 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) &
- & / e3vw(ji,jj,jk ,Kmm) * z1_e3va * wvmask(ji,jj,jk )
- zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) &
- & / e3vw(ji,jj,jk+1,Kmm) * z1_e3va * wvmask(ji,jj,jk+1)
- zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) * z1_e3va
- zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * z1_e3va
- zwi(ji,jj,jk) = zzwi + zDt_2 * MIN( zWvi, 0._wp )
- zws(ji,jj,jk) = zzws - zDt_2 * MAX( zWvs, 0._wp )
- zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zDt_2 * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
- END_3D
- END SELECT
- DO_2D( 0, 0, 0, 0 ) !* Surface boundary conditions
- zwi(ji,jj,1) = 0._wp
- zzws = - zDt_2 * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) &
- & / ( e3v(ji,jj,1,Kaa) * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2)
- zWvs = ( wi(ji,jj ,2) + wi(ji,jj+1,2) ) / e3v(ji,jj,1,Kaa)
- zws(ji,jj,1 ) = zzws - zDt_2 * MAX( zWvs, 0._wp )
- zwd(ji,jj,1 ) = 1._wp - zzws - zDt_2 * ( MIN( zWvs, 0._wp ) )
- END_2D
- ELSE
- SELECT CASE( nldf_dyn )
- CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu)
- DO_3D( 0, 0, 0, 0, 1, jpkm1 )
- zzwi = - zDt_2 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) &
- & / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk )
- zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) &
- & / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
- zwi(ji,jj,jk) = zzwi
- zws(ji,jj,jk) = zzws
- zwd(ji,jj,jk) = 1._wp - zzwi - zzws
- END_3D
- CASE DEFAULT ! iso-level lateral mixing
- DO_3D( 0, 0, 0, 0, 1, jpkm1 )
- zzwi = - zDt_2 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) &
- & / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk )
- zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) &
- & / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
- zwi(ji,jj,jk) = zzwi
- zws(ji,jj,jk) = zzws
- zwd(ji,jj,jk) = 1._wp - zzwi - zzws
- END_3D
- END SELECT
- DO_2D( 0, 0, 0, 0 ) !* Surface boundary conditions
- zwi(ji,jj,1) = 0._wp
- zwd(ji,jj,1) = 1._wp - zws(ji,jj,1)
+ END_1D
+ DO_2Dik( 0, 0, 2, jpkm1, 1 )
+ puu(ji,jj,jk,Kaa) = puu(ji,jj,jk,Kaa) - zwi(ji,jk) / zwd(ji,jk-1) * puu(ji,jj,jk-1,Kaa)
END_2D
- ENDIF
- !
- ! !== Apply semi-implicit top/bottom friction ==!
- !
- ! Only needed for semi-implicit bottom friction setup. The explicit
- ! bottom friction has been included in "u(v)a" which act as the R.H.S
- ! column vector of the tri-diagonal matrix equation
- !
- IF( ln_drgimp ) THEN
- DO_2D( 0, 0, 0, 0 )
- ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
- zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - zDt_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) &
- & / e3v(ji,jj,ikv,Kaa)
+ !
+ DO_1Di( 0, 0 ) !== thrid recurrence : SOLk = ( Lk - Uk * Ek+1 ) / Dk ==!
+ puu(ji,jj,jpkm1,Kaa) = puu(ji,jj,jpkm1,Kaa) / zwd(ji,jpkm1)
+ END_1D
+ DO_2Dik( 0, 0, jpk-2, 1, -1 )
+ puu(ji,jj,jk,Kaa) = ( puu(ji,jj,jk,Kaa) - zws(ji,jk) * puu(ji,jj,jk+1,Kaa) ) / zwd(ji,jk)
END_2D
- IF ( ln_isfcav.OR.ln_drgice_imp ) THEN
- DO_2D( 0, 0, 0, 0 )
- ikv = mikv(ji,jj) ! (first wet ocean u- and v-points)
- zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - zDt_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) &
+ !
+ !
+ ! !== Vertical diffusion on v ==!
+ !
+ ! !* Matrix construction
+ IF( ln_zad_Aimp ) THEN !!
+ SELECT CASE( nldf_dyn )
+ CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzv)
+ DO_2Dik( 0, 0, 1, jpkm1, 1 )
+ z1_e3va = 1._wp / e3v(ji,jj,jk,Kaa) ! after scale factor at V-point
+ zzwi = - zDt_2 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) &
+ & / e3vw(ji,jj,jk ,Kmm) * z1_e3va * wvmask(ji,jj,jk )
+ zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) &
+ & / e3vw(ji,jj,jk+1,Kmm) * z1_e3va * wvmask(ji,jj,jk+1)
+ zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) * z1_e3va
+ zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * z1_e3va
+ zwi(ji,jk) = zzwi + zDt_2 * MIN( zWvi, 0._wp )
+ zws(ji,jk) = zzws - zDt_2 * MAX( zWvs, 0._wp )
+ zwd(ji,jk) = 1._wp - zzwi - zzws - zDt_2 * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
+ END_2D
+ CASE DEFAULT ! iso-level lateral mixing
+ DO_2Dik( 0, 0, 1, jpkm1, 1 )
+ z1_e3va = 1._wp / e3v(ji,jj,jk,Kaa) ! after scale factor at V-point
+ zzwi = - zDt_2 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) &
+ & / e3vw(ji,jj,jk ,Kmm) * z1_e3va * wvmask(ji,jj,jk )
+ zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) &
+ & / e3vw(ji,jj,jk+1,Kmm) * z1_e3va * wvmask(ji,jj,jk+1)
+ zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) * z1_e3va
+ zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) * z1_e3va
+ zwi(ji,jk) = zzwi + zDt_2 * MIN( zWvi, 0._wp )
+ zws(ji,jk) = zzws - zDt_2 * MAX( zWvs, 0._wp )
+ zwd(ji,jk) = 1._wp - zzwi - zzws - zDt_2 * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
+ END_2D
+ END SELECT
+ DO_1Di( 0, 0 ) !* Surface boundary conditions
+ zwi(ji,1) = 0._wp
+ zzws = - zDt_2 * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) &
+ & / ( e3v(ji,jj,1,Kaa) * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2)
+ zWvs = ( wi(ji,jj ,2) + wi(ji,jj+1,2) ) / e3v(ji,jj,1,Kaa)
+ zws(ji,1 ) = zzws - zDt_2 * MAX( zWvs, 0._wp )
+ zwd(ji,1 ) = 1._wp - zzws - zDt_2 * ( MIN( zWvs, 0._wp ) )
+ END_1D
+ ELSE
+ SELECT CASE( nldf_dyn )
+ CASE( np_lap_i ) ! rotated lateral mixing: add its vertical mixing (akzu)
+ DO_2Dik( 0, 0, 1, jpkm1, 1 )
+ zzwi = - zDt_2 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) &
+ & / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk )
+ zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) &
+ & / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
+ zwi(ji,jk) = zzwi
+ zws(ji,jk) = zzws
+ zwd(ji,jk) = 1._wp - zzwi - zzws
+ END_2D
+ CASE DEFAULT ! iso-level lateral mixing
+ DO_2Dik( 0, 0, 1, jpkm1, 1 )
+ zzwi = - zDt_2 * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) &
+ & / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk )
+ zzws = - zDt_2 * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) &
+ & / ( e3v(ji,jj,jk,Kaa) * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
+ zwi(ji,jk) = zzwi
+ zws(ji,jk) = zzws
+ zwd(ji,jk) = 1._wp - zzwi - zzws
+ END_2D
+ END SELECT
+ DO_1Di( 0, 0 ) !* Surface boundary conditions
+ zwi(ji,1) = 0._wp
+ zwd(ji,1) = 1._wp - zws(ji,1)
+ END_1D
+ ENDIF
+ !
+ ! !== Apply semi-implicit top/bottom friction ==!
+ !
+ ! Only needed for semi-implicit bottom friction setup. The explicit
+ ! bottom friction has been included in "u(v)a" which act as the R.H.S
+ ! column vector of the tri-diagonal matrix equation
+ !
+ IF( ln_drgimp ) THEN
+ DO_1Di( 0, 0 )
+ ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
+ zwd(ji,ikv) = zwd(ji,ikv) - zDt_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) &
& / e3v(ji,jj,ikv,Kaa)
- END_2D
+ END_1D
+ IF ( ln_isfcav.OR.ln_drgice_imp ) THEN
+ DO_1Di( 0, 0 )
+ ikv = mikv(ji,jj) ! (first wet ocean u- and v-points)
+ zwd(ji,ikv) = zwd(ji,ikv) - zDt_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) &
+ & / e3v(ji,jj,ikv,Kaa)
+ END_1D
+ ENDIF
ENDIF
- ENDIF
-
- ! Matrix inversion
- !-----------------------------------------------------------------------
- ! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
- !
- ! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
- ! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
- ! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
- ! ( ... )( ... ) ( ... )
- ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
- !
- ! m is decomposed in the product of an upper and lower triangular matrix
- ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
- ! The solution (after velocity) is in 2d array va
- !-----------------------------------------------------------------------
- !
- DO_3D( 0, 0, 0, 0, 2, jpkm1 ) !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
- zwd(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwd(ji,jj,jk-1)
- END_3D
- !
- DO_2D( 0, 0, 0, 0 ) !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==!
+
+ ! Matrix inversion
+ !-----------------------------------------------------------------------
+ ! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
+ !
+ ! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
+ ! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
+ ! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
+ ! ( ... )( ... ) ( ... )
+ ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
+ !
+ ! m is decomposed in the product of an upper and lower triangular matrix
+ ! The 3 diagonal terms are in 2d arrays: zwd, zws, zwi
+ ! The solution (after velocity) is in 2d array va
+ !-----------------------------------------------------------------------
+ !
+ DO_2Dik( 0, 0, 2, jpkm1, 1 ) !== First recurrence : Dk = Dk - Lk * Uk-1 / Dk-1 (increasing k) ==
+ zwd(ji,jk) = zwd(ji,jk) - zwi(ji,jk) * zws(ji,jk-1) / zwd(ji,jk-1)
+ END_2D
+ !
+ DO_1Di( 0, 0 ) !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==!
#if defined key_RK3
- ! ! RK3: use only vtau (not vtau_b)
- pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * vtauV(ji,jj) &
- & / ( e3v(ji,jj,1,Kaa) * rho0 ) * vmask(ji,jj,1)
+ ! ! RK3: use only vtau (not vtau_b)
+ pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * vtauV(ji,jj) &
+ & / ( e3v(ji,jj,1,Kaa) * rho0 ) * vmask(ji,jj,1)
#else
- pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + zDt_2 * ( vtau_b(ji,jj) + vtauV(ji,jj) ) &
- & / ( e3v(ji,jj,1,Kaa) * rho0 ) * vmask(ji,jj,1)
+ ! ! MLF: average of vtau and vtau_b
+ pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + zDt_2*( vtau_b(ji,jj) + vtauV(ji,jj) ) &
+ & / ( e3v(ji,jj,1,Kaa) * rho0 ) * vmask(ji,jj,1)
#endif
- END_2D
- DO_3D( 0, 0, 0, 0, 2, jpkm1 )
- pvv(ji,jj,jk,Kaa) = pvv(ji,jj,jk,Kaa) - zwi(ji,jj,jk) / zwd(ji,jj,jk-1) * pvv(ji,jj,jk-1,Kaa)
- END_3D
- !
- DO_2D( 0, 0, 0, 0 ) !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==!
- pvv(ji,jj,jpkm1,Kaa) = pvv(ji,jj,jpkm1,Kaa) / zwd(ji,jj,jpkm1)
- END_2D
- DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 )
- pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - zws(ji,jj,jk) * pvv(ji,jj,jk+1,Kaa) ) / zwd(ji,jj,jk)
- END_3D
+ END_1D
+ DO_2Dik( 0, 0, 2, jpkm1, 1 )
+ pvv(ji,jj,jk,Kaa) = pvv(ji,jj,jk,Kaa) - zwi(ji,jk) / zwd(ji,jk-1) * pvv(ji,jj,jk-1,Kaa)
+ END_2D
+ !
+ DO_1Di( 0, 0 ) !== third recurrence : SOLk = ( Lk - Uk * SOLk+1 ) / Dk ==!
+ pvv(ji,jj,jpkm1,Kaa) = pvv(ji,jj,jpkm1,Kaa) / zwd(ji,jpkm1)
+ END_1D
+ DO_2Dik( 0, 0, jpk-2, 1, -1 )
+ pvv(ji,jj,jk,Kaa) = ( pvv(ji,jj,jk,Kaa) - zws(ji,jk) * pvv(ji,jj,jk+1,Kaa) ) / zwd(ji,jk)
+ END_2D
+ ! ! ================= !
+ END_1D ! i-k slices loop !
+ ! ! ================= !
!
IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics
ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) )*r1_Dt - ztrdu(:,:,:)
diff --git a/src/OCE/TRA/traadv.F90 b/src/OCE/TRA/traadv.F90
index 190a6a7ad3cc55d62bed5e3023620e3177faf32e..924ee9331d45c6343a8ee7ab12ad8041909ee983 100644
--- a/src/OCE/TRA/traadv.F90
+++ b/src/OCE/TRA/traadv.F90
@@ -10,6 +10,7 @@ MODULE traadv
!! - ! 2014-12 (G. Madec) suppression of cross land advection option
!! 3.6 ! 2015-06 (E. Clementi) Addition of Stokes drift in case of wave coupling
!! 4.5 ! 2021-04 (G. Madec, S. Techene) add advective velocities as optional arguments
+ !! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
@@ -178,7 +179,7 @@ CONTAINS
!
!!gm ???
! TEMP: [tiling] This copy-in not necessary after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
- IF( l_diaptr ) CALL dia_ptr( kt, Kmm, zvv(A2D(0),:) ) ! diagnose the effective MSF
+ IF( l_diaptr ) CALL dia_ptr( kt, Kmm, zvv(A2D(nn_hls),:) ) ! diagnose the effective MSF
!!gm ???
!
@@ -191,11 +192,19 @@ CONTAINS
SELECT CASE ( nadv ) !== compute advection trend and add it to general trend ==!
!
CASE ( np_CEN ) ! Centered scheme : 2nd / 4th order
- CALL tra_adv_cen ( kt, nit000, 'TRA', zuu, zvv, zww, Kmm, pts, jpts, Krhs, nn_cen_h, nn_cen_v )
+ IF( nn_hls == 1 ) THEN
+ CALL tra_adv_cen_hls1( kt, nit000, 'TRA', zuu, zvv, zww, Kmm, pts, jpts, Krhs, nn_cen_h, nn_cen_v )
+ ELSE
+ CALL tra_adv_cen ( kt, nit000, 'TRA', zuu, zvv, zww, Kmm, pts, jpts, Krhs, nn_cen_h, nn_cen_v )
+ ENDIF
CASE ( np_FCT ) ! FCT scheme : 2nd / 4th order
CALL tra_adv_fct ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_fct_h, nn_fct_v )
CASE ( np_MUS ) ! MUSCL
+ IF( nn_hls == 1 ) THEN
+ CALL tra_adv_mus_hls1( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups )
+ ELSE
CALL tra_adv_mus( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, ln_mus_ups )
+ END IF
CASE ( np_UBS ) ! UBS
CALL tra_adv_ubs ( kt, nit000, 'TRA', rDt, zuu, zvv, zww, Kbb, Kmm, pts, jpts, Krhs, nn_ubs_v )
CASE ( np_QCK ) ! QUICKEST
diff --git a/src/OCE/TRA/traadv_cen.F90 b/src/OCE/TRA/traadv_cen.F90
index 3351e1935a5503648b2e9e0ca590bb5a7813c037..1666eee41455854145befb571038452abd560ceb 100644
--- a/src/OCE/TRA/traadv_cen.F90
+++ b/src/OCE/TRA/traadv_cen.F90
@@ -4,6 +4,7 @@ MODULE traadv_cen
!! Ocean tracers: advective trend (2nd/4th order centered)
!!======================================================================
!! History : 3.7 ! 2014-05 (G. Madec) original code
+ !! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
@@ -29,7 +30,8 @@ MODULE traadv_cen
IMPLICIT NONE
PRIVATE
- PUBLIC tra_adv_cen ! called by traadv.F90
+ PUBLIC tra_adv_cen ! called by traadv.F90
+ PUBLIC tra_adv_cen_hls1 ! called by traadv.F90
REAL(wp) :: r1_6 = 1._wp / 6._wp ! =1/6
@@ -47,7 +49,321 @@ MODULE traadv_cen
!!----------------------------------------------------------------------
CONTAINS
+ SUBROUTINE tra_adv_cen_test( kt, kit000, cdtype, pU, pV, pW, &
+ & Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE tra_adv_cen ***
+ !!
+ !! ** Purpose : Compute the now trend due to the advection of tracers
+ !! and add it to the general trend of passive tracer equations.
+ !!
+ !! ** Method : The advection is evaluated by a 2nd or 4th order scheme
+ !! using now fields (leap-frog scheme).
+ !! kn_cen_h = 2 ==>> 2nd order centered scheme on the horizontal
+ !! = 4 ==>> 4th order - - - -
+ !! kn_cen_v = 2 ==>> 2nd order centered scheme on the vertical
+ !! = 4 ==>> 4th order COMPACT scheme - -
+ !!
+ !! ** Action : - update pt(:,:,:,:,Krhs) with the now advective tracer trends
+ !! - send trends to trdtra module for further diagnostcs (l_trdtra=T)
+ !! - poleward advective heat and salt transport (l_diaptr=T)
+ !!----------------------------------------------------------------------
+ INTEGER , INTENT(in ) :: kt ! ocean time-step index
+ INTEGER , INTENT(in ) :: Kmm, Krhs ! ocean time level indices
+ INTEGER , INTENT(in ) :: kit000 ! first time step index
+ CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
+ INTEGER , INTENT(in ) :: kjpt ! number of tracers
+ INTEGER , INTENT(in ) :: kn_cen_h ! =2/4 (2nd or 4th order scheme)
+ INTEGER , INTENT(in ) :: kn_cen_v ! =2/4 (2nd or 4th order scheme)
+ ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
+ REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume flux components
+ REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation
+ !
+ INTEGER :: ji, jj, jk, jn ! dummy loop indices
+ INTEGER :: ierr ! local integer
+ REAL(wp) :: zC2t_u, zC4t_u ! local scalars
+ REAL(wp) :: zC2t_v, zC4t_v ! - -
+ REAL(wp) :: zftw_kp1
+ REAL(wp), DIMENSION(A2D(1)) :: zft_u, zft_v !, zft_w
+ REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: zdt_u, zdt_v
+ REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: ztw
+ !!----------------------------------------------------------------------
+ !
+#if defined key_loop_fusion
+ CALL tra_adv_cen_lf ( kt, nit000, cdtype, pU, pV, pW, Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v )
+#else
+ IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
+ IF( kt == kit000 ) THEN
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) 'tra_adv_cen : centered advection scheme on ', cdtype, ' order h/v =', kn_cen_h,'/', kn_cen_v
+ IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~ '
+ ENDIF
+ ! ! set local switches
+ l_trd = .FALSE.
+ l_hst = .FALSE.
+ l_ptr = .FALSE.
+ IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
+ IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.
+ IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
+ & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE.
+ ENDIF
+ !
+ ! !* 2nd order centered
+ DO jn = 1, kjpt !== loop over the tracers ==!
+ !
+ DO jk = 1, jpkm1
+ !
+ DO_2D( 1, 0, 1, 0 ) ! Horizontal fluxes at layer jk
+ zft_u(ji,jj) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj ,jk,jn,Kmm) )
+ zft_v(ji,jj) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) )
+ END_2D
+ !
+ DO_2D( 0, 0, 0, 0 ) ! Horizontal divergence of advective fluxes
+ pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zft_u(ji,jj) - zft_u(ji-1,jj ) &
+ & + zft_v(ji,jj) - zft_v(ji ,jj-1) ) * r1_e1e2t(ji,jj) &
+ & / e3t(ji,jj,jk,Kmm)
+ END_2D
+ END DO
+ !
+#define zft_w zft_u
+ IF( ln_linssh ) THEN !* top value (linear free surf. only as zwz is multiplied by wmask)
+ DO_2D( 0, 0, 0, 0 )
+ zft_w(ji,jj) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kmm)
+ END_2D
+ ELSE
+ DO_2D( 0, 0, 0, 0 )
+ zft_w(ji,jj) = 0._wp
+ END_2D
+ ENDIF
+ DO jk = 1, jpk-2
+ DO_2D( 0, 0, 0, 0 ) ! Vertical fluxes
+ zftw_kp1 = 0.5 * pW(ji,jj,jk+1) * ( pt(ji,jj,jk+1,jn,Kmm) + pt(ji,jj,jk,jn,Kmm) ) * wmask(ji,jj,jk+1)
+ !
+ pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zft_w(ji,jj) - zftw_kp1 ) * r1_e1e2t(ji,jj) &
+ & / e3t(ji,jj,jk,Kmm)
+ zft_w(ji,jj) = zftw_kp1
+ END_2D
+ END DO
+ jk = jpkm1 ! bottom vertical flux set to zero for all tracers
+ DO_2D( 0, 0, 0, 0 )
+ pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - zft_w(ji,jj) * r1_e1e2t(ji,jj) &
+ & / e3t(ji,jj,jk,Kmm)
+ END_2D
+ !
+ END DO
+ !
+
+ !
+#undef zft_w
+ ! ! trend diagnostics
+!!gm + !!st to be done with the whole rewritting of trd
+!! trd routine arguments MUST be changed adding jk and zwx, zwy in 2D
+!! IF( l_trd ) THEN
+!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kmm) )
+!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) )
+!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) )
+!! ENDIF
+!! ! ! "Poleward" heat and salt transports
+!! IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
+!! ! ! heat and salt transport
+!! IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
+ !
+ !
+ !
+#endif
+ END SUBROUTINE tra_adv_cen_test
+
+
SUBROUTINE tra_adv_cen( kt, kit000, cdtype, pU, pV, pW, &
+ & Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE tra_adv_cen ***
+ !!
+ !! ** Purpose : Compute the now trend due to the advection of tracers
+ !! and add it to the general trend of passive tracer equations.
+ !!
+ !! ** Method : The advection is evaluated by a 2nd or 4th order scheme
+ !! using now fields (leap-frog scheme).
+ !! kn_cen_h = 2 ==>> 2nd order centered scheme on the horizontal
+ !! = 4 ==>> 4th order - - - -
+ !! kn_cen_v = 2 ==>> 2nd order centered scheme on the vertical
+ !! = 4 ==>> 4th order COMPACT scheme - -
+ !!
+ !! ** Action : - update pt(:,:,:,:,Krhs) with the now advective tracer trends
+ !! - send trends to trdtra module for further diagnostcs (l_trdtra=T)
+ !! - poleward advective heat and salt transport (l_diaptr=T)
+ !!----------------------------------------------------------------------
+ INTEGER , INTENT(in ) :: kt ! ocean time-step index
+ INTEGER , INTENT(in ) :: Kmm, Krhs ! ocean time level indices
+ INTEGER , INTENT(in ) :: kit000 ! first time step index
+ CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
+ INTEGER , INTENT(in ) :: kjpt ! number of tracers
+ INTEGER , INTENT(in ) :: kn_cen_h ! =2/4 (2nd or 4th order scheme)
+ INTEGER , INTENT(in ) :: kn_cen_v ! =2/4 (2nd or 4th order scheme)
+ ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
+ REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume flux components
+ REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation
+ !
+ INTEGER :: ji, jj, jk, jn ! dummy loop indices
+ INTEGER :: ierr ! local integer
+ REAL(wp) :: zC2t_u, zC4t_u ! local scalars
+ REAL(wp) :: zC2t_v, zC4t_v ! - -
+ REAL(wp) :: zftw_kp1
+ REAL(wp), DIMENSION(A2D(1)) :: zft_u, zft_v
+ REAL(wp), DIMENSION(:,:) , ALLOCATABLE :: zdt_u, zdt_v
+ REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: ztw
+ !!----------------------------------------------------------------------
+ !
+#if defined key_loop_fusion
+ CALL tra_adv_cen_lf ( kt, nit000, cdtype, pU, pV, pW, Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v )
+#else
+ IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
+ IF( kt == kit000 ) THEN
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) 'tra_adv_cen : centered advection scheme on ', cdtype, ' order h/v =', kn_cen_h,'/', kn_cen_v
+ IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~ '
+ ENDIF
+ ! ! set local switches
+ l_trd = .FALSE.
+ l_hst = .FALSE.
+ l_ptr = .FALSE.
+ IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
+ IF( cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.
+ IF( cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
+ & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE.
+ ENDIF
+ !
+ IF( kn_cen_h == 4 ) ALLOCATE( zdt_u(A2D(2)) , zdt_v(A2D(2)) ) ! horizontal 4th order only
+ IF( kn_cen_v == 4 ) ALLOCATE( ztw(A2D(nn_hls),jpk) ) ! vertical 4th order only
+ !
+ DO jn = 1, kjpt !== loop over the tracers ==!
+ !
+ SELECT CASE( kn_cen_h ) !-- Horizontal divergence of advective fluxes --!
+ !
+!!st limitation : does not take into acccount iceshelf specificity
+!! in case of linssh
+ CASE( 2 ) !* 2nd order centered
+ DO jk = 1, jpkm1
+ !
+ DO_2D( 1, 0, 1, 0 ) ! Horizontal fluxes at layer jk
+ zft_u(ji,jj) = 0.5_wp * pU(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj ,jk,jn,Kmm) )
+ zft_v(ji,jj) = 0.5_wp * pV(ji,jj,jk) * ( pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm) )
+ END_2D
+ !
+ DO_2D( 0, 0, 0, 0 ) ! Horizontal divergence of advective fluxes
+ pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zft_u(ji,jj) - zft_u(ji-1,jj ) &
+ & + zft_v(ji,jj) - zft_v(ji ,jj-1) ) * r1_e1e2t(ji,jj) &
+ & / e3t(ji,jj,jk,Kmm)
+ END_2D
+ END DO
+ !
+ CASE( 4 ) !* 4th order centered
+ DO jk = 1, jpkm1
+ DO_2D( 2, 1, 2, 1 ) ! masked gradient
+ zdt_u(ji,jj) = ( pt(ji+1,jj ,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * umask(ji,jj,jk)
+ zdt_v(ji,jj) = ( pt(ji ,jj+1,jk,jn,Kmm) - pt(ji,jj,jk,jn,Kmm) ) * vmask(ji,jj,jk)
+ END_2D
+ !
+ DO_2D( 1, 0, 1, 0 ) ! Horizontal advective fluxes
+ zC2t_u = pt(ji,jj,jk,jn,Kmm) + pt(ji+1,jj ,jk,jn,Kmm) ! C2 interpolation of T at u- & v-points (x2)
+ zC2t_v = pt(ji,jj,jk,jn,Kmm) + pt(ji ,jj+1,jk,jn,Kmm)
+ ! ! C4 interpolation of T at u- & v-points (x2)
+ zC4t_u = zC2t_u + r1_6 * ( zdt_u(ji-1,jj ) - zdt_u(ji+1,jj ) )
+ zC4t_v = zC2t_v + r1_6 * ( zdt_v(ji ,jj-1) - zdt_v(ji ,jj+1) )
+ ! ! C4 fluxes
+ zft_u(ji,jj) = 0.5_wp * pU(ji,jj,jk) * zC4t_u
+ zft_v(ji,jj) = 0.5_wp * pV(ji,jj,jk) * zC4t_v
+ END_2D
+ !
+ DO_2D( 0, 0, 0, 0 ) ! Horizontal divergence of advective fluxes
+ pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zft_u(ji,jj) - zft_u(ji-1,jj ) &
+ & + zft_v(ji,jj) - zft_v(ji ,jj-1) ) * r1_e1e2t(ji,jj) &
+ & / e3t(ji,jj,jk,Kmm)
+ END_2D
+ END DO
+ !
+ CASE DEFAULT
+ CALL ctl_stop( 'traadv_cen: wrong value for nn_cen' )
+ END SELECT
+ !
+#define zft_w zft_u
+ !
+ IF( ln_linssh ) THEN !* top value (linear free surf. only as zwz is multiplied by wmask)
+ DO_2D( 0, 0, 0, 0 )
+ zft_w(ji,jj) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kmm)
+ END_2D
+ ELSE
+ DO_2D( 0, 0, 0, 0 )
+ zft_w(ji,jj) = 0._wp
+ END_2D
+ ENDIF
+ !
+ SELECT CASE( kn_cen_v ) !-- Vertical divergence of advective fluxes --! (interior)
+ !
+ CASE( 2 ) !* 2nd order centered
+ DO jk = 1, jpk-2
+ DO_2D( 0, 0, 0, 0 ) ! Vertical fluxes
+ zftw_kp1 = 0.5 * pW(ji,jj,jk+1) * ( pt(ji,jj,jk+1,jn,Kmm) + pt(ji,jj,jk,jn,Kmm) ) * wmask(ji,jj,jk+1)
+ !
+ pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zft_w(ji,jj) - zftw_kp1 ) * r1_e1e2t(ji,jj) &
+ & / e3t(ji,jj,jk,Kmm)
+ zft_w(ji,jj) = zftw_kp1
+ END_2D
+ END DO
+ jk = jpkm1 ! bottom vertical flux set to zero for all tracers
+ DO_2D( 0, 0, 0, 0 )
+ pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - zft_w(ji,jj) * r1_e1e2t(ji,jj) &
+ & / e3t(ji,jj,jk,Kmm)
+ END_2D
+ !
+ CASE( 4 ) !* 4th order compact
+ CALL interp_4th_cpt( pt(:,:,:,jn,Kmm) , ztw ) ! ztw = interpolated value of T at w-point
+ !
+ DO jk = 1, jpk-2
+ !
+ DO_2D( 0, 0, 0, 0 )
+ zftw_kp1 = pW(ji,jj,jk+1) * ztw(ji,jj,jk+1) * wmask(ji,jj,jk+1)
+ ! ! Divergence of advective fluxes
+ pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zft_w(ji,jj) - zftw_kp1 ) * r1_e1e2t(ji,jj) &
+ & / e3t(ji,jj,jk,Kmm)
+ ! ! update
+ zft_w(ji,jj) = zftw_kp1
+ END_2D
+ !
+ END DO
+ !
+ jk = jpkm1 ! bottom vertical flux set to zero for all tracers
+ DO_2D( 0, 0, 0, 0 )
+ pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - zft_w(ji,jj) * r1_e1e2t(ji,jj) &
+ & / e3t(ji,jj,jk,Kmm)
+ END_2D
+ !
+ END SELECT
+ !
+#undef zft_w
+ ! ! trend diagnostics
+!!gm + !!st to be done with the whole rewritting of trd
+!! trd routine arguments MUST be changed adding jk and zwx, zwy in 2D
+!! IF( l_trd ) THEN
+!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kmm) )
+!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_yad, zwy, pV, pt(:,:,:,jn,Kmm) )
+!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwz, pW, pt(:,:,:,jn,Kmm) )
+!! ENDIF
+!! ! ! "Poleward" heat and salt transports
+!! IF( l_ptr ) CALL dia_ptr_hst( jn, 'adv', zwy(:,:,:) )
+!! ! ! heat and salt transport
+!! IF( l_hst ) CALL dia_ar5_hst( jn, 'adv', zwx(:,:,:), zwy(:,:,:) )
+ !
+ END DO
+ !
+ IF( kn_cen_h == 4 ) DEALLOCATE( zdt_u , zdt_v ) ! horizontal 4th order only
+ IF( kn_cen_v == 4 ) DEALLOCATE( ztw ) ! vertical 4th order only
+ !
+#endif
+ END SUBROUTINE tra_adv_cen
+
+
+ SUBROUTINE tra_adv_cen_hls1( kt, kit000, cdtype, pU, pV, pW, &
& Kmm, pt, kjpt, Krhs, kn_cen_h, kn_cen_v )
!!----------------------------------------------------------------------
!! *** ROUTINE tra_adv_cen ***
@@ -141,6 +457,12 @@ CONTAINS
CASE DEFAULT
CALL ctl_stop( 'traadv_cen: wrong value for nn_cen' )
END SELECT
+ DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !-- Divergence of advective fluxes --!
+ pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) &
+ & - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
+ & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) &
+ & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+ END_3D
!
SELECT CASE( kn_cen_v ) !-- Vertical fluxes --! (interior)
!
@@ -171,11 +493,17 @@ CONTAINS
!
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !-- Divergence of advective fluxes --!
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) &
- & - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
- & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
- & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) &
+ & - ( zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
END_3D
+!
+!!st DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !-- Divergence of advective fluxes --!
+!!st pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) &
+!!st & - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
+!!st & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) &
+!!st & + zwz(ji,jj,jk) - zwz(ji ,jj ,jk+1) ) &
+!!st & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+!!st END_3D
! ! trend diagnostics
IF( l_trd ) THEN
CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_xad, zwx, pU, pt(:,:,:,jn,Kmm) )
@@ -190,7 +518,7 @@ CONTAINS
END DO
!
#endif
- END SUBROUTINE tra_adv_cen
+ END SUBROUTINE tra_adv_cen_hls1
!!======================================================================
END MODULE traadv_cen
diff --git a/src/OCE/TRA/traadv_mus.F90 b/src/OCE/TRA/traadv_mus.F90
index 4c4b0437ca01cb250cf6f560527e8b0b5a632803..c22fcfd01f44ed4ad18071d829c32b6853371ade 100644
--- a/src/OCE/TRA/traadv_mus.F90
+++ b/src/OCE/TRA/traadv_mus.F90
@@ -9,6 +9,7 @@ MODULE traadv_mus
!! 3.2 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA + switch from velocity to transport
!! 3.4 ! 2012-06 (P. Oddo, M. Vichi) include the upstream where needed
!! 3.7 ! 2015-09 (G. Madec) add the ice-shelf cavities boundary condition
+ !! 4.5 ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
@@ -34,7 +35,9 @@ MODULE traadv_mus
IMPLICIT NONE
PRIVATE
- PUBLIC tra_adv_mus ! routine called by traadv.F90
+ PUBLIC tra_adv_mus ! routine called by traadv.F90
+ PUBLIC tra_adv_mus_hls1 ! routine called by traadv.F90
+
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: upsmsk !: mixed upstream/centered scheme near some straits
! ! and in closed seas (orca 2 and 1 configurations)
@@ -55,6 +58,217 @@ MODULE traadv_mus
CONTAINS
SUBROUTINE tra_adv_mus( kt, kit000, cdtype, p2dt, pU, pV, pW, &
+ & Kbb, Kmm, pt, kjpt, Krhs, ld_msc_ups )
+ !!----------------------------------------------------------------------
+ !! *** ROUTINE tra_adv_mus ***
+ !!
+ !! ** Purpose : Compute the now trend due to total advection of tracers
+ !! using a MUSCL scheme (Monotone Upstream-centered Scheme for
+ !! Conservation Laws) and add it to the general tracer trend.
+ !!
+ !! ** Method : MUSCL scheme plus centered scheme at ocean boundaries
+ !! ld_msc_ups=T :
+ !!
+ !! ** Action : - update pt(:,:,:,:,Krhs) with the now advective tracer trends
+ !! - send trends to trdtra module for further diagnostcs (l_trdtra=T)
+ !! - poleward advective heat and salt transport (ln_diaptr=T)
+ !!
+ !! References : Estubier, A., and M. Levy, Notes Techn. Pole de Modelisation
+ !! IPSL, Sept. 2000 (http://www.lodyc.jussieu.fr/opa)
+ !!----------------------------------------------------------------------
+ INTEGER , INTENT(in ) :: kt ! ocean time-step index
+ INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs ! ocean time level indices
+ INTEGER , INTENT(in ) :: kit000 ! first time step index
+ CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
+ INTEGER , INTENT(in ) :: kjpt ! number of tracers
+ LOGICAL , INTENT(in ) :: ld_msc_ups ! use upstream scheme within muscl
+ REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step
+ ! TEMP: [tiling] This can be A2D(nn_hls) after all lbc_lnks removed in the nn_hls = 2 case in tra_adv_fct
+ REAL(wp), DIMENSION(jpi,jpj,jpk ), INTENT(in ) :: pU, pV, pW ! 3 ocean volume flux components
+ REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation
+ !
+ INTEGER :: ji, jj, jk, jn ! dummy loop indices
+ INTEGER :: ierr ! local integer
+ INTEGER :: ik
+ REAL(wp) :: zu, z0u, zzslpx, zzwx, zw , zalpha ! local scalars
+ REAL(wp) :: zv, z0v, zzslpy, zzwy, z0w ! - -
+ REAL(wp) :: zdzt_kp2, zslpz_kp1, zfW_kp1
+ REAL(wp), DIMENSION(A2D(2)) :: zdxt, zslpx, zwx ! 2D workspace
+ REAL(wp), DIMENSION(A2D(2)) :: zdyt, zslpy, zwy ! - -
+ !!----------------------------------------------------------------------
+ !
+ IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
+ IF( kt == kit000 ) THEN
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*) 'tra_adv : MUSCL advection scheme on ', cdtype
+ IF(lwp) WRITE(numout,*) ' : mixed up-stream ', ld_msc_ups
+ IF(lwp) WRITE(numout,*) '~~~~~~~'
+ IF(lwp) WRITE(numout,*)
+ !
+ ! Upstream / MUSCL scheme indicator
+ !
+ ALLOCATE( xind(jpi,jpj,jpk), STAT=ierr )
+ xind(:,:,:) = 1._wp ! set equal to 1 where up-stream is not needed
+ IF( ld_msc_ups ) THEN ! define the upstream indicator (if asked)
+ DO jk = 1, jpkm1
+ xind(:,:,jk) = 1._wp & ! =>1 where up-stream is not needed
+ & - rnfmsk(:,:) * rnfmsk_z(jk) * tmask(:,:,jk) ! =>0 near runoff mouths (& closed sea outflows)
+ END DO
+ ENDIF
+ !
+ ENDIF
+ !
+ l_trd = .FALSE.
+ l_hst = .FALSE.
+ l_ptr = .FALSE.
+ IF( ( cdtype == 'TRA' .AND. l_trdtra ) .OR. ( cdtype == 'TRC' .AND. l_trdtrc ) ) l_trd = .TRUE.
+ IF( l_diaptr .AND. cdtype == 'TRA' .AND. ( iom_use( 'sophtadv' ) .OR. iom_use( 'sophtadv' ) ) ) l_ptr = .TRUE.
+ IF( l_iom .AND. cdtype == 'TRA' .AND. ( iom_use("uadv_heattr") .OR. iom_use("vadv_heattr") .OR. &
+ & iom_use("uadv_salttr") .OR. iom_use("vadv_salttr") ) ) l_hst = .TRUE.
+ ENDIF
+ !
+ !
+ DO jn = 1, kjpt !== loop over the tracers ==!
+ !
+ DO jk = 1, jpkm1
+ ! !* Horizontal advective fluxes
+ !
+ ! !-- first guess of the slopes
+ DO_2D( 2, 1, 2, 1 )
+ zdxt(ji,jj) = ( pt(ji+1,jj ,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) ) * umask(ji,jj,jk)
+ zdyt(ji,jj) = ( pt(ji ,jj+1,jk,jn,Kbb) - pt(ji,jj,jk,jn,Kbb) ) * vmask(ji,jj,jk)
+ END_2D
+ ! !-- Slopes of tracer
+ DO_2D( 1, 1, 1, 1 )
+ ! ! 1/2 Slopes at T-point (set to 0 if adjectent slopes are of opposite sign)
+ zzslpx = ( zdxt(ji,jj) + zdxt(ji-1,jj ) ) &
+ & * ( 0.25_wp + SIGN( 0.25_wp, zdxt(ji,jj) * zdxt(ji-1,jj ) ) )
+ zzslpy = ( zdyt(ji,jj) + zdyt(ji ,jj-1) ) &
+ & * ( 0.25_wp + SIGN( 0.25_wp, zdyt(ji,jj) * zdyt(ji ,jj-1) ) )
+ ! ! Slopes limitation
+ zslpx(ji,jj) = SIGN( 1.0_wp, zzslpx ) * MIN( ABS( zzslpx ), &
+ & 2._wp*ABS( zdxt (ji-1,jj) ), &
+ & 2._wp*ABS( zdxt (ji ,jj) ) )
+ zslpy(ji,jj) = SIGN( 1.0_wp, zzslpy ) * MIN( ABS( zzslpy ), &
+ & 2._wp*ABS( zdyt (ji,jj-1) ), &
+ & 2._wp*ABS( zdyt (ji,jj ) ) )
+ END_2D
+!!gm + !!st ticket ? comparaison pommes et carrottes ABS(zzslpx) et 2._wp*ABS( zdxt (ji-1,jj) )
+ !
+ DO_2D( 1, 0, 1, 0 ) !-- MUSCL horizontal advective fluxes
+ z0u = SIGN( 0.5_wp, pU(ji,jj,jk) )
+ zalpha = 0.5_wp - z0u
+ zu = z0u - 0.5_wp * pU(ji,jj,jk) * p2dt * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm)
+ zzwx = pt(ji+1,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji+1,jj)
+ zzwy = pt(ji ,jj,jk,jn,Kbb) + xind(ji,jj,jk) * zu * zslpx(ji ,jj)
+ zwx(ji,jj) = pU(ji,jj,jk) * ( zalpha * zzwx + (1._wp-zalpha) * zzwy )
+ !
+ z0v = SIGN( 0.5_wp, pV(ji,jj,jk) )
+ zalpha = 0.5_wp - z0v
+ zv = z0v - 0.5_wp * pV(ji,jj,jk) * p2dt * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm)
+ zzwx = pt(ji,jj+1,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj+1)
+ zzwy = pt(ji,jj ,jk,jn,Kbb) + xind(ji,jj,jk) * zv * zslpy(ji,jj )
+ zwy(ji,jj) = pV(ji,jj,jk) * ( zalpha * zzwx + (1._wp-zalpha) * zzwy )
+ END_2D
+ !
+ DO_2D( 0, 0, 0, 0 ) !-- Tracer advective trend
+ pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj) - zwx(ji-1,jj ) &
+ & + zwy(ji,jj) - zwy(ji ,jj-1) ) &
+ & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+ END_2D
+ END DO
+!!gm + !!st to be done with the whole rewritting of trd
+!! trd routine arguments MUST be changed adding jk and zwx, zwy in 2D
+!!
+!! ! ! trend diagnostics
+!! IF( l_trd ) THEN
+!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jk, jptra_xad, zwx(:,:), pU, pt(:,:,:,jn,Kbb) )
+!! CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jk, jptra_yad, zwy(:,:), pV, pt(:,:,:,jn,Kbb) )
+!! END IF
+!! ! ! "Poleward" heat and salt transports
+!! IF( l_ptr ) CALL dia_ptr_hst( jn, jk, 'adv', zwy(:,:) )
+!! ! ! heat transport
+!! IF( l_hst ) CALL dia_ar5_hst( jn, jk, 'adv', zwx(:,:), zwy(:,:) )
+ !
+ ! !* Vertical advective fluxes
+ !
+#define zdzt_kp1 zdxt
+#define zslpz zslpx
+#define zfW zwx
+
+ zfW (A2D(0)) = 0._wp ! anciennement zwx at jk = 1
+ ! ! anciennement zwx at jk = 2
+ DO_2D( 0, 0, 0, 0 )
+ zdzt_kp1(ji,jj) = tmask(ji,jj,2) * ( pt(ji,jj,1,jn,Kbb) - pt(ji,jj,2,jn,Kbb) )
+ END_2D
+ zslpz (A2D(0)) = 0._wp ! anciennement zslpx at jk = 1
+ !
+ IF( ln_linssh ) THEN !-- linear ssh : non zero top values
+ DO_2D( 0, 0, 0, 0 ) ! at the ocean surface
+ zfW(ji,jj) = pW(ji,jj,1) * pt(ji,jj,1,jn,Kbb) ! surface flux
+ END_2D
+ IF( ln_isfcav ) THEN ! ice-shelf cavities (top of the ocean)
+ DO_2D( 0, 0, 0, 0 ) ! update pt(Krhs) under the ice-shelf
+ ik = mikt(ji,jj) ! the flux at ik-1 is zero ( inside ice-shelf )
+ IF( ik > 1 ) THEN
+ pt(ji,jj,ik,jn,Krhs) = pt(ji,jj,ik,jn,Krhs) - pW(ji,jj,ik) * pt(ji,jj,ik,jn,Kbb) &
+ & * r1_e1e2t(ji,jj) / e3t(ji,jj,ik,Kmm)
+ ENDIF
+ END_2D
+ ENDIF
+ ENDIF
+ !
+ ! wmask usage for computing zw and zwk is needed in isf case and linear ssh
+ !
+ !
+ DO jk = 1, jpkm1
+ IF( jk < jpkm1 ) THEN
+ DO_2D( 0, 0, 0, 0 )
+ ! !-- Slopes of tracer
+ ! ! masked vertical gradient at jk+2
+ zdzt_kp2 = ( pt(ji,jj,jk+1,jn,Kbb) - pt(ji,jj,jk+2,jn,Kbb) ) * tmask(ji,jj,jk+2) !!st wmask(ji,jj,jk+2)
+ ! ! vertical slope at jk+1
+ zslpz_kp1 = ( zdzt_kp1(ji,jj) + zdzt_kp2 ) &
+ & * ( 0.25_wp + SIGN( 0.25_wp, zdzt_kp1(ji,jj) * zdzt_kp2 ) )
+ ! ! slope limitation
+ zslpz_kp1 = SIGN( 1.0_wp, zslpz_kp1 ) * MIN( ABS( zslpz_kp1 ), &
+ & 2.*ABS( zdzt_kp2 ), &
+ & 2.*ABS( zdzt_kp1(ji,jj) ) )
+ ! !-- vertical advective flux at jk+1
+ ! ! caution: zfW_kp1 is masked for ice-shelf cavities
+ ! ! since top fluxes already added to pt(Krhs) before the vertical loop
+ z0w = SIGN( 0.5_wp, pW(ji,jj,jk+1) )
+ zalpha = 0.5_wp + z0w
+ zw = z0w - 0.5_wp * pW(ji,jj,jk+1) * p2dt * r1_e1e2t(ji,jj) / e3w(ji,jj,jk+1,Kmm)
+ zzwx = pt(ji,jj,jk+1,jn,Kbb) + xind(ji,jj,jk) * zw * zslpz_kp1
+ zzwy = pt(ji,jj,jk ,jn,Kbb) + xind(ji,jj,jk) * zw * zslpz(ji,jj)
+ zfW_kp1 = pW(ji,jj,jk+1) * ( zalpha * zzwx + (1.-zalpha) * zzwy ) * wmask(ji,jj,jk)!!st * wmask(ji,jj,jk+1)
+ ! !-- vertical advective trend at jk
+ pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zfW(ji,jj) - zfW_kp1 ) &
+ & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+ ! ! updates for next level
+ zdzt_kp1(ji,jj) = zdzt_kp2
+ zslpz (ji,jj) = zslpz_kp1
+ zfW (ji,jj) = zfW_kp1
+ END_2D
+ ELSE
+ DO_2D( 0, 0, 0, 0 ) !-- vertical advective trend at jpkm1
+ pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - zfW(ji,jj) &
+ & * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
+ END_2D
+ ENDIF
+ END DO ! end of jk loop
+ !
+!!gm + !!st idem see above
+!! ! ! send trends for diagnostic
+!! IF( l_trd ) CALL trd_tra( kt, Kmm, Krhs, cdtype, jn, jptra_zad, zwx, pW, pt(:,:,:,jn,Kbb) )
+ !
+ END DO ! end of tracer loop
+ !
+ END SUBROUTINE tra_adv_mus
+
+
+ SUBROUTINE tra_adv_mus_hls1( kt, kit000, cdtype, p2dt, pU, pV, pW, &
& Kbb, Kmm, pt, kjpt, Krhs, ld_msc_ups )
!!----------------------------------------------------------------------
!! *** ROUTINE tra_adv_mus ***
@@ -178,7 +392,7 @@ CONTAINS
!
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !-- Tracer advective trend
pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) - ( zwx(ji,jj,jk) - zwx(ji-1,jj ,jk ) &
- & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) &
+ & + zwy(ji,jj,jk) - zwy(ji ,jj-1,jk ) ) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
END_3D
! ! trend diagnostics
@@ -239,7 +453,7 @@ CONTAINS
!
END DO ! end of tracer loop
!
- END SUBROUTINE tra_adv_mus
+ END SUBROUTINE tra_adv_mus_hls1
!!======================================================================
END MODULE traadv_mus
diff --git a/src/OCE/TRA/trazdf.F90 b/src/OCE/TRA/trazdf.F90
index 6e7dd9a249bebc44102a67bf2ff140da9282c1de..7123aeba37c7fe1419e207dd89407a28d569e25e 100644
--- a/src/OCE/TRA/trazdf.F90
+++ b/src/OCE/TRA/trazdf.F90
@@ -6,6 +6,7 @@ MODULE trazdf
!! History : 1.0 ! 2005-11 (G. Madec) Original code
!! 3.0 ! 2008-01 (C. Ethe, G. Madec) merge TRC-TRA
!! 4.0 ! 2017-06 (G. Madec) remove explict time-stepping option
+ !! 4.5 ! 2022-06 (G. Madec) refactoring to reduce memory usage (j-k-i loops)
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
@@ -22,7 +23,7 @@ MODULE trazdf
USE ldfslp ! lateral diffusion: iso-neutral slope
USE trd_oce ! trends: ocean variables
USE trdtra ! trends: tracer trend manager
- USE eosbn2, ONLY: ln_SEOS, rn_b0
+ USE eosbn2 , ONLY: ln_SEOS, rn_b0
!
USE in_out_manager ! I/O manager
USE prtctl ! Print control
@@ -77,7 +78,7 @@ CONTAINS
ENDIF
!
! !* compute lateral mixing trend and add it to the general trend
- CALL tra_zdf_imp( kt, nit000, 'TRA', rDt, Kbb, Kmm, Krhs, pts, Kaa, jpts )
+ CALL tra_zdf_imp( 'TRA', Kbb, Kmm, Krhs, pts, Kaa, jpts )
!!gm WHY here ! and I don't like that !
! DRAKKAR SSS control {
@@ -116,7 +117,7 @@ CONTAINS
END SUBROUTINE tra_zdf
- SUBROUTINE tra_zdf_imp( kt, kit000, cdtype, p2dt, Kbb, Kmm, Krhs, pt, Kaa, kjpt )
+ SUBROUTINE tra_zdf_imp( cdtype, Kbb, Kmm, Krhs, pt, Kaa, kjpt )
!!----------------------------------------------------------------------
!! *** ROUTINE tra_zdf_imp ***
!!
@@ -136,128 +137,177 @@ CONTAINS
!!
!! ** Action : - pt(:,:,:,:,Kaa) becomes the after tracer
!!---------------------------------------------------------------------
- INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kbb, Kmm, Krhs, Kaa ! ocean time level indices
- INTEGER , INTENT(in ) :: kit000 ! first time step index
CHARACTER(len=3) , INTENT(in ) :: cdtype ! =TRA or TRC (tracer indicator)
INTEGER , INTENT(in ) :: kjpt ! number of tracers
- REAL(wp) , INTENT(in ) :: p2dt ! tracer time-step
REAL(wp), DIMENSION(jpi,jpj,jpk,kjpt,jpt), INTENT(inout) :: pt ! tracers and RHS of tracer equation
!
INTEGER :: ji, jj, jk, jn ! dummy loop indices
REAL(wp) :: zrhs, zzwi, zzws ! local scalars
- REAL(wp), DIMENSION(A2D(0),jpk) :: zwi, zwt, zwd, zws
+ REAL(wp), DIMENSION(A1Di(0),jpk) :: zwi, zwt, zwd, zws
!!---------------------------------------------------------------------
!
- ! ! ============= !
- DO jn = 1, kjpt ! tracer loop !
- ! ! ============= !
- ! Matrix construction
- ! --------------------
- ! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer
- !
- IF( ( cdtype == 'TRA' .AND. ( jn == jp_tem .OR. ( jn == jp_sal .AND. ln_zdfddm ) ) ) .OR. &
- & ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN
+ ! ! ================= !
+ DO_1Dj( 0, 0 ) ! i-k slices loop !
+ ! ! ================= !
+ DO jn = 1, kjpt ! tracer loop !
+ ! ! ================= !
!
- ! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers
- IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN
- DO_3D( 0, 0, 0, 0, 2, jpk )
- zwt(ji,jj,jk) = avt(ji,jj,jk)
- END_3D
- ELSE
- DO_3D( 0, 0, 0, 0, 2, jpk )
- zwt(ji,jj,jk) = avs(ji,jj,jk)
- END_3D
- ENDIF
- zwt(:,:,1) = 0._wp
+ ! Matrix construction
+ ! --------------------
+ ! Build matrix if temperature or salinity (only in double diffusion case) or first passive tracer
!
- IF( l_ldfslp ) THEN ! isoneutral diffusion: add the contribution
- IF( ln_traldf_msc ) THEN ! MSC iso-neutral operator
- DO_3D( 0, 0, 0, 0, 2, jpkm1 )
- zwt(ji,jj,jk) = zwt(ji,jj,jk) + akz(ji,jj,jk)
- END_3D
- ELSE ! standard or triad iso-neutral operator
- DO_3D( 0, 0, 0, 0, 2, jpkm1 )
- zwt(ji,jj,jk) = zwt(ji,jj,jk) + ah_wslp2(ji,jj,jk)
- END_3D
+ IF( ( cdtype == 'TRA' .AND. ( jn == jp_tem .OR. ( jn == jp_sal .AND. ln_zdfddm ) ) ) .OR. &
+ & ( cdtype == 'TRC' .AND. jn == 1 ) ) THEN
+ !
+ ! vertical mixing coef.: avt for temperature, avs for salinity and passive tracers
+ !
+ IF( cdtype == 'TRA' .AND. jn == jp_tem ) THEN ! use avt for temperature
+ !
+ IF( l_ldfslp ) THEN ! use avt + isoneutral diffusion contribution
+ IF( ln_traldf_msc ) THEN ! MSC iso-neutral operator
+ DO_2Dik( 0, 0, 2, jpk, 1 )
+ zwt(ji,jk) = avt(ji,jj,jk) + akz(ji,jj,jk)
+ END_2D
+ ELSE ! standard or triad iso-neutral operator
+ DO_2Dik( 0, 0, 2, jpk, 1 )
+ zwt(ji,jk) = avt(ji,jj,jk) + ah_wslp2(ji,jj,jk)
+ END_2D
+ ENDIF
+ ELSE ! use avt only
+ DO_2Dik( 0, 0, 2, jpk, 1 )
+ zwt(ji,jk) = avt(ji,jj,jk)
+ END_2D
+ ENDIF
+ !
+ ELSE ! use avs for salinty or passive tracers
+ !
+ IF( l_ldfslp ) THEN ! use avs + isoneutral diffusion contribution
+ IF( ln_traldf_msc ) THEN ! MSC iso-neutral operator
+ DO_2Dik( 0, 0, 2, jpk, 1 )
+ zwt(ji,jk) = avs(ji,jj,jk) + akz(ji,jj,jk)
+ END_2D
+ ELSE ! standard or triad iso-neutral operator
+ DO_2Dik( 0, 0, 2, jpk, 1 )
+ zwt(ji,jk) = avs(ji,jj,jk) + ah_wslp2(ji,jj,jk)
+ END_2D
+ ENDIF
+ ELSE !
+ DO_2Dik( 0, 0, 2, jpk, 1 )
+ zwt(ji,jk) = avs(ji,jj,jk)
+ END_2D
+ ENDIF
ENDIF
+ zwt(:,1) = 0._wp
+ !
+ ! Diagonal, lower (i), upper (s) (including the bottom boundary condition since avt is masked)
+ IF( ln_zad_Aimp ) THEN ! Adaptive implicit vertical advection
+ DO_2Dik( 0, 0, 1, jpkm1, 1 )
+ zzwi = - rDt * zwt(ji,jk ) / e3w(ji,jj,jk ,Kmm)
+ zzws = - rDt * zwt(ji,jk+1) / e3w(ji,jj,jk+1,Kmm)
+ zwd(ji,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws &
+ & + rDt * ( MAX( wi(ji,jj,jk ) , 0._wp ) &
+ & - MIN( wi(ji,jj,jk+1) , 0._wp ) )
+ zwi(ji,jk) = zzwi + rDt * MIN( wi(ji,jj,jk ) , 0._wp )
+ zws(ji,jk) = zzws - rDt * MAX( wi(ji,jj,jk+1) , 0._wp )
+ END_2D
+ ELSE
+ DO_2Dik( 0, 0, 1, jpkm1, 1 )
+ zwi(ji,jk) = - rDt * zwt(ji,jk ) / e3w(ji,jj,jk,Kmm)
+ zws(ji,jk) = - rDt * zwt(ji,jk+1) / e3w(ji,jj,jk+1,Kmm)
+ zwd(ji,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jk) - zws(ji,jk)
+ END_2D
+ ENDIF
+ !
+!!gm BUG?? : if edmfm is equivalent to a w ==>>> just add +/- rDt * edmfm(ji,jj,jk+1/jk )
+!! but edmfm is at t-point !!!! crazy??? why not keep it at w-point????
+ !
+ IF( ln_zdfmfc ) THEN ! add upward Mass Flux in the matrix
+ DO_2Dik( 0, 0, 1, jpkm1, 1 )
+ zws(ji,jk) = zws(ji,jk) + e3t(ji,jj,jk,Kaa) * rDt * edmfm(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
+ zwd(ji,jk) = zwd(ji,jk) - e3t(ji,jj,jk,Kaa) * rDt * edmfm(ji,jj,jk ) / e3w(ji,jj,jk+1,Kmm)
+ END_2D
+ ENDIF
+! DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+! edmfa(ji,jj,jk) = 0._wp
+! edmfb(ji,jj,jk) = -edmfm(ji,jj,jk ) / e3w(ji,jj,jk+1,Kmm)
+! edmfc(ji,jj,jk) = edmfm(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
+! END_3D
+!!gm BUG : level jpk never used in the inversion
+! DO_2D( 0, 0, 0, 0 )
+! edmfa(ji,jj,jpk) = -edmfm(ji,jj,jpk-1) / e3w(ji,jj,jpk,Kmm)
+! edmfb(ji,jj,jpk) = edmfm(ji,jj,jpk ) / e3w(ji,jj,jpk,Kmm)
+! edmfc(ji,jj,jpk) = 0._wp
+! END_2D
+!!
+!!gm BUG ??? below e3t_Kmm should be used ?
+!! or even no multiplication by e3t unless there is a bug in wi calculation
+!!
+! DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+!!gm edmfa = 0._wp except at jpk which is not used ==>> zdiagi update is useless !
+! zdiagi(ji,jj,jk) = zdiagi(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfa(ji,jj,jk)
+! zdiags(ji,jj,jk) = zdiags(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfc(ji,jj,jk)
+! zdiagd(ji,jj,jk) = zdiagd(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfb(ji,jj,jk)
+! END_3D
+!!gm CALL diag_mfc( zwi, zwd, zws, rDt, Kaa )
+!!gm SUBROUTINE diag_mfc( zdiagi, zdiagd, zdiags, p2dt, Kaa )
+ !
+ !! Matrix inversion from the first level
+ !!----------------------------------------------------------------------
+ ! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
+ !
+ ! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
+ ! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
+ ! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
+ ! ( ... )( ... ) ( ... )
+ ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
+ !
+ ! m is decomposed in the product of an upper and lower triangular matrix.
+ ! The 3 diagonal terms are in 3d arrays: zwd, zws, zwi.
+ ! Suffices i,s and d indicate "inferior" (below diagonal), diagonal
+ ! and "superior" (above diagonal) components of the tridiagonal system.
+ ! The solution will be in the 4d array pta.
+ ! The 3d array zwt is used as a work space array.
+ ! En route to the solution pt(:,:,:,:,Kaa) is used a to evaluate the rhs and then
+ ! used as a work space array: its value is modified.
+ !
+ DO_1Di( 0, 0 ) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k) ! done one for all passive tracers (so included in the IF instruction)
+ zwt(ji,1) = zwd(ji,1)
+ END_1D
+ DO_2Dik( 0, 0, 2, jpkm1, 1 )
+ zwt(ji,jk) = zwd(ji,jk) - zwi(ji,jk) * zws(ji,jk-1) / zwt(ji,jk-1)
+ END_2D
+ !
ENDIF
!
- ! Diagonal, lower (i), upper (s) (including the bottom boundary condition since avt is masked)
- IF( ln_zad_Aimp ) THEN ! Adaptive implicit vertical advection
- DO_3D( 0, 0, 0, 0, 1, jpkm1 )
- zzwi = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk ,Kmm)
- zzws = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
- zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zzwi - zzws &
- & + p2dt * ( MAX( wi(ji,jj,jk ) , 0._wp ) - MIN( wi(ji,jj,jk+1) , 0._wp ) )
- zwi(ji,jj,jk) = zzwi + p2dt * MIN( wi(ji,jj,jk ) , 0._wp )
- zws(ji,jj,jk) = zzws - p2dt * MAX( wi(ji,jj,jk+1) , 0._wp )
- END_3D
- ELSE
- DO_3D( 0, 0, 0, 0, 1, jpkm1 )
- zwi(ji,jj,jk) = - p2dt * zwt(ji,jj,jk ) / e3w(ji,jj,jk,Kmm)
- zws(ji,jj,jk) = - p2dt * zwt(ji,jj,jk+1) / e3w(ji,jj,jk+1,Kmm)
- zwd(ji,jj,jk) = e3t(ji,jj,jk,Kaa) - zwi(ji,jj,jk) - zws(ji,jj,jk)
- END_3D
+ IF( ln_zdfmfc ) THEN ! add Mass Flux to the RHS
+ DO_2Dik( 0, 0, 1, jpkm1, 1 )
+ pt(ji,jj,jk,jn,Krhs) = pt(ji,jj,jk,jn,Krhs) + edmftra(ji,jj,jk,jn)
+ END_2D
+!!gm CALL rhs_mfc( pt(:,:,:,jn,Krhs), jn )
ENDIF
!
- ! Modification of diagonal to add MF scheme
- IF ( ln_zdfmfc ) THEN
- CALL diag_mfc( zwi, zwd, zws, p2dt, Kaa )
- END IF
- !
- !! Matrix inversion from the first level
- !!----------------------------------------------------------------------
- ! solve m.x = y where m is a tri diagonal matrix ( jpk*jpk )
- !
- ! ( zwd1 zws1 0 0 0 )( zwx1 ) ( zwy1 )
- ! ( zwi2 zwd2 zws2 0 0 )( zwx2 ) ( zwy2 )
- ! ( 0 zwi3 zwd3 zws3 0 )( zwx3 )=( zwy3 )
- ! ( ... )( ... ) ( ... )
- ! ( 0 0 0 zwik zwdk )( zwxk ) ( zwyk )
- !
- ! m is decomposed in the product of an upper and lower triangular matrix.
- ! The 3 diagonal terms are in 3d arrays: zwd, zws, zwi.
- ! Suffices i,s and d indicate "inferior" (below diagonal), diagonal
- ! and "superior" (above diagonal) components of the tridiagonal system.
- ! The solution will be in the 4d array pta.
- ! The 3d array zwt is used as a work space array.
- ! En route to the solution pt(:,:,:,:,Kaa) is used a to evaluate the rhs and then
- ! used as a work space array: its value is modified.
- !
- DO_2D( 0, 0, 0, 0 ) !* 1st recurrence: Tk = Dk - Ik Sk-1 / Tk-1 (increasing k) ! done one for all passive tracers (so included in the IF instruction)
- zwt(ji,jj,1) = zwd(ji,jj,1)
+ DO_1Di( 0, 0 ) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1
+ pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb ) &
+ & + rDt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs)
+ END_1D
+ DO_2Dik( 0, 0, 2, jpkm1, 1 )
+ zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb ) &
+ & + rDt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs) ! zrhs=right hand side
+ pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jk) / zwt(ji,jk-1) * pt(ji,jj,jk-1,jn,Kaa)
END_2D
- DO_3D( 0, 0, 0, 0, 2, jpkm1 )
- zwt(ji,jj,jk) = zwd(ji,jj,jk) - zwi(ji,jj,jk) * zws(ji,jj,jk-1) / zwt(ji,jj,jk-1)
- END_3D
!
- ENDIF
- !
- ! Modification of rhs to add MF scheme
- IF ( ln_zdfmfc ) THEN
- CALL rhs_mfc( pt(:,:,:,jn,Krhs), jn )
- END IF
- !
- DO_2D( 0, 0, 0, 0 ) !* 2nd recurrence: Zk = Yk - Ik / Tk-1 Zk-1
- pt(ji,jj,1,jn,Kaa) = e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb) &
- & + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs)
- END_2D
- DO_3D( 0, 0, 0, 0, 2, jpkm1 )
- zrhs = e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) &
- & + p2dt * e3t(ji,jj,jk,Kmm) * pt(ji,jj,jk,jn,Krhs) ! zrhs=right hand side
- pt(ji,jj,jk,jn,Kaa) = zrhs - zwi(ji,jj,jk) / zwt(ji,jj,jk-1) * pt(ji,jj,jk-1,jn,Kaa)
- END_3D
- !
- DO_2D( 0, 0, 0, 0 ) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk (result is the after tracer)
- pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jj,jpkm1) * tmask(ji,jj,jpkm1)
- END_2D
- DO_3DS( 0, 0, 0, 0, jpk-2, 1, -1 )
- pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jj,jk) * pt(ji,jj,jk+1,jn,Kaa) ) &
- & / zwt(ji,jj,jk) * tmask(ji,jj,jk)
- END_3D
+ DO_1Di( 0, 0 ) !* 3d recurrence: Xk = (Zk - Sk Xk+1 ) / Tk (result is the after tracer)
+ pt(ji,jj,jpkm1,jn,Kaa) = pt(ji,jj,jpkm1,jn,Kaa) / zwt(ji,jpkm1) * tmask(ji,jj,jpkm1)
+ END_1D
+ DO_2Dik( 0, 0, jpk-2, 1, -1 )
+ pt(ji,jj,jk,jn,Kaa) = ( pt(ji,jj,jk,jn,Kaa) - zws(ji,jk) * pt(ji,jj,jk+1,jn,Kaa) ) &
+ & / zwt(ji,jk) * tmask(ji,jj,jk)
+ END_2D
+ ! ! ================= !
+ END DO ! tracer loop !
! ! ================= !
- END DO ! end tracer loop !
+ END_1D ! i-k slices loop !
! ! ================= !
END SUBROUTINE tra_zdf_imp
diff --git a/src/OCE/do_loop_substitute.h90 b/src/OCE/do_loop_substitute.h90
index f957d0741b54399e90fba13f4876ae152b2ea041..c6db00d51aef101b5ca421abf373c001deeffedc 100644
--- a/src/OCE/do_loop_substitute.h90
+++ b/src/OCE/do_loop_substitute.h90
@@ -58,7 +58,10 @@
!
#endif
-#define DO_2D(L, R, B, T) DO jj = ntsj-(B), ntej+(T) ; DO ji = ntsi-(L), ntei+(R)
+#define DO_1Di(L, R) DO ji = ntsi-(L), ntei+(R)
+#define DO_1Dj(B, T) DO jj = ntsj-(B), ntej+(T)
+#define DO_2Dik(L, R, ks, ke, ki) DO jk = ks, ke, ki ; DO_1Di(L, R)
+#define DO_2D(L, R, B, T) DO_1Dj(B, T) ; DO_1Di(L, R)
#define DO_2D_OVR(L, R, B, T) DO_2D(L-(L+R)*nthl, R-(R+L)*nthr, B-(B+T)*nthb, T-(T+B)*ntht)
#define A1Di(H) ntsi-(H):ntei+(H)
#define A1Dj(H) ntsj-(H):ntej+(H)
@@ -76,5 +79,6 @@
#define DO_3DS(L, R, B, T, ks, ke, ki) DO jk = ks, ke, ki ; DO_2D(L, R, B, T)
#define DO_3DS_OVR(L, R, B, T, ks, ke, ki) DO jk = ks, ke, ki ; DO_2D_OVR(L, R, B, T)
+#define END_1D END DO
#define END_2D END DO ; END DO
-#define END_3D END DO ; END DO ; END DO
+#define END_3D END DO ; END DO ; END DO
\ No newline at end of file
diff --git a/src/TOP/TRP/trcadv.F90 b/src/TOP/TRP/trcadv.F90
index e98590e1693612ef16395d623624c477273d89d9..fd1b0aeaac3911aee54cb4ec5cd942eaca372cc1 100644
--- a/src/TOP/TRP/trcadv.F90
+++ b/src/TOP/TRP/trcadv.F90
@@ -8,6 +8,7 @@ MODULE trcadv
!! 3.7 ! 2014-05 (G. Madec, C. Ethe) Add 2nd/4th order cases for CEN and FCT schemes
!! 4.0 ! 2017-09 (G. Madec) remove vertical time-splitting option
!! 4.5 ! 2021-08 (G. Madec, S. Techene) add advective velocities as optional arguments
+ !! - ! 2022-06 (S. Techene, G, Madec) refactorization to reduce local memory usage
!!----------------------------------------------------------------------
#if defined key_top
!!----------------------------------------------------------------------
@@ -156,11 +157,19 @@ CONTAINS
SELECT CASE ( nadv ) !== compute advection trend and add it to general trend ==!
!
CASE ( np_CEN ) ! Centered : 2nd / 4th order
- CALL tra_adv_cen ( kt, nittrc000,'TRC', zuu, zvv, zww, Kmm, ptr, jptra, Krhs, nn_cen_h, nn_cen_v )
+ IF( nn_hls == 1 ) THEN
+ CALL tra_adv_cen_hls1( kt, nittrc000,'TRC', zuu, zvv, zww, Kmm, ptr, jptra, Krhs, nn_cen_h, nn_cen_v )
+ ELSE
+ CALL tra_adv_cen ( kt, nittrc000,'TRC', zuu, zvv, zww, Kmm, ptr, jptra, Krhs, nn_cen_h, nn_cen_v )
+ ENDIF
CASE ( np_FCT ) ! FCT : 2nd / 4th order
CALL tra_adv_fct( kt, nittrc000,'TRC', rDt_trc, zuu, zvv, zww, Kbb, Kmm, ptr, jptra, Krhs, nn_fct_h, nn_fct_v )
CASE ( np_MUS ) ! MUSCL
+ IF( nn_hls == 1 ) THEN
+ CALL tra_adv_mus_hls1( kt, nittrc000,'TRC', rDt_trc, zuu, zvv, zww, Kbb, Kmm, ptr, jptra, Krhs, ln_mus_ups )
+ ELSE
CALL tra_adv_mus( kt, nittrc000,'TRC', rDt_trc, zuu, zvv, zww, Kbb, Kmm, ptr, jptra, Krhs, ln_mus_ups )
+ ENDIF
CASE ( np_UBS ) ! UBS
CALL tra_adv_ubs ( kt, nittrc000,'TRC', rDt_trc, zuu, zvv, zww, Kbb, Kmm, ptr, jptra, Krhs, nn_ubs_v )
CASE ( np_QCK ) ! QUICKEST
diff --git a/src/TOP/TRP/trczdf.F90 b/src/TOP/TRP/trczdf.F90
index 2f0b5a8f2dbebb0ba09b1467a6b1b7831777b741..b5a9dd4ccd34296c7f3ba759cc8830496b836d95 100644
--- a/src/TOP/TRP/trczdf.F90
+++ b/src/TOP/TRP/trczdf.F90
@@ -53,7 +53,7 @@ CONTAINS
!
IF( l_trdtrc ) ztrtrd(:,:,:,:) = ptr(:,:,:,:,Krhs)
!
- CALL tra_zdf_imp( kt, nittrc000, 'TRC', rDt_trc, Kbb, Kmm, Krhs, ptr, Kaa, jptra ) ! implicit scheme
+ CALL tra_zdf_imp( 'TRC', Kbb, Kmm, Krhs, ptr, Kaa, jptra ) ! implicit scheme
!
IF( l_trdtrc ) THEN ! save the vertical diffusive trends for further diagnostics
DO jn = 1, jptra