Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
No results found
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 1182 additions and 595 deletions
src/NST/agrif_user.F90
View file @ d1c6ae27
......@@ -27,7 +27,11 @@
!!----------------------------------------------------------------------
!
CALL nemo_init !* Initializations of each fine grid
# if defined key_RK3
Kbb_a = Nbb; Kmm_a = Nbb; Krhs_a = Nrhs
# else
Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs ! agrif_oce module copies of time level indices
# endif
!
! !* Agrif initialization
CALL Agrif_InitValues_cont
......@@ -410,29 +414,17 @@
hbdy(:,:) = 0._wp
ssh(:,:,Krhs_a) = 0._wp
IF ( ln_dynspg_ts ) THEN
Agrif_UseSpecialValue = ln_spc_dyn
use_sign_north = .TRUE.
sign_north = -1.
CALL Agrif_Bc_variable(ub2b_interp_id,calledweight=1.,procname=interpub2b) ! must be called before unb_id to define ubdy
CALL Agrif_Bc_variable(vb2b_interp_id,calledweight=1.,procname=interpvb2b) ! must be called before vnb_id to define vbdy
CALL Agrif_Bc_variable( unb_interp_id,calledweight=1.,procname=interpunb )
CALL Agrif_Bc_variable( vnb_interp_id,calledweight=1.,procname=interpvnb )
use_sign_north = .FALSE.
ubdy(:,:) = 0._wp
vbdy(:,:) = 0._wp
ELSEIF ( ln_dynspg_EXP ) THEN
Agrif_UseSpecialValue = ln_spc_dyn
use_sign_north = .TRUE.
sign_north = -1.
ubdy(:,:) = 0._wp
vbdy(:,:) = 0._wp
CALL Agrif_Bc_variable( unb_interp_id,calledweight=1.,procname=interpunb )
CALL Agrif_Bc_variable( vnb_interp_id,calledweight=1.,procname=interpvnb )
use_sign_north = .FALSE.
ubdy(:,:) = 0._wp
vbdy(:,:) = 0._wp
ENDIF
Agrif_UseSpecialValue = ln_spc_dyn
use_sign_north = .TRUE.
sign_north = -1.
ubdy(:,:) = 0._wp
vbdy(:,:) = 0._wp
CALL Agrif_Bc_variable( unb_interp_id,calledweight=1.,procname=interpunb )
CALL Agrif_Bc_variable( vnb_interp_id,calledweight=1.,procname=interpvnb )
use_sign_north = .FALSE.
ubdy(:,:) = 0._wp
vbdy(:,:) = 0._wp
Agrif_UseSpecialValue = .FALSE.
l_vremap = .FALSE.
......
src/OCE/DIA/diaar5.F90
View file @ d1c6ae27
......@@ -79,7 +79,6 @@ CONTAINS
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z2d, zpe ! 2D workspace
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: z3d, zrhd, ztpot, zgdept ! 3D workspace (zgdept: needed to use the substitute)
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: ztsn ! 4D workspace
!!--------------------------------------------------------------------
IF( ln_timing ) CALL timing_start('dia_ar5')
......@@ -318,6 +317,7 @@ CONTAINS
!
INTEGER :: ji, jj, jk
REAL(wp), DIMENSION(A2D(nn_hls)) :: z2d
!!----------------------------------------------------------------------
z2d(:,:) = puflx(:,:,1)
DO_3D( 0, 0, 0, 0, 1, jpkm1 )
......@@ -355,7 +355,7 @@ CONTAINS
!! ** Purpose : initialization for AR5 diagnostic computation
!!----------------------------------------------------------------------
INTEGER :: inum
INTEGER :: ik, idep
INTEGER :: ik
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zztmp
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) :: zsaldta ! Jan/Dec levitus salinity
......@@ -384,9 +384,9 @@ CONTAINS
zvol0 (:,:) = 0._wp
thick0(:,:) = 0._wp
DO_3D( 1, 1, 1, 1, 1, jpkm1 ) ! interpolation of salinity at the last ocean level (i.e. the partial step)
idep = tmask(ji,jj,jk) * e3t_0(ji,jj,jk)
zvol0 (ji,jj) = zvol0 (ji,jj) + idep * e1e2t(ji,jj)
thick0(ji,jj) = thick0(ji,jj) + idep
zztmp = tmask(ji,jj,jk) * e3t_0(ji,jj,jk)
zvol0 (ji,jj) = zvol0 (ji,jj) + zztmp * e1e2t(ji,jj)
thick0(ji,jj) = thick0(ji,jj) + zztmp
END_3D
vol0 = glob_sum( 'diaar5', zvol0 )
DEALLOCATE( zvol0 )
......
src/OCE/DIA/diaptr.F90
View file @ d1c6ae27
......@@ -41,6 +41,7 @@ MODULE diaptr
END INTERFACE
PUBLIC dia_ptr ! call in step module
PUBLIC dia_ptr_init ! call in stprk3 module
PUBLIC dia_ptr_hst ! called from tra_ldf/tra_adv routines
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: hstr_adv, hstr_ldf, hstr_eiv !: Heat/Salt TRansports(adv, diff, Bolus.)
......@@ -65,7 +66,7 @@ MODULE diaptr
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: diaptr.F90 14834 2021-05-11 09:24:44Z hadcv $
!! $Id: diaptr.F90 15513 2021-11-15 17:31:29Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -82,7 +83,9 @@ CONTAINS
!
IF( ln_timing ) CALL timing_start('dia_ptr')
#if ! defined key_RK3
IF( kt == nit000 .AND. ll_init ) CALL dia_ptr_init ! -> will define l_diaptr and nbasin
#endif
!
IF( l_diaptr ) THEN
! Calculate zonal integrals
......
src/OCE/DOM/dom_oce.F90
View file @ d1c6ae27
......@@ -152,13 +152,12 @@ MODULE dom_oce
! ! reference depths of cells
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdept_0 !: t- depth [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gdepw_0 !: w- depth [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gde3w_0 !: w- depth (sum of e3w) [m]
! ! time-dependent depths of cells (domvvl)
#if defined key_qco || defined key_linssh
#else
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: gdept, gdepw
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gde3w
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: gde3w_0, gde3w !: w- depth (sum of e3w) [m]
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: gdept, gdepw
#endif
! ! reference heights of ocean water column and its inverse
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ht_0, r1_ht_0 !: t-depth [m] and [1/m]
......@@ -286,7 +285,7 @@ CONTAINS
& ff_f (jpi,jpj) , ff_t (jpi,jpj) , STAT=ierr(ii) )
!
ii = ii+1
ALLOCATE( gdept_0 (jpi,jpj,jpk) , gdepw_0 (jpi,jpj,jpk) , gde3w_0(jpi,jpj,jpk) , &
ALLOCATE( gdept_0 (jpi,jpj,jpk) , gdepw_0 (jpi,jpj,jpk) , &
& gdept_1d( jpk) , gdepw_1d( jpk) , STAT=ierr(ii) )
!
ii = ii+1
......@@ -311,7 +310,8 @@ CONTAINS
#else
! vvl : time varation for all vertical coordinate variables
ii = ii+1
ALLOCATE( gdept (jpi,jpj,jpk,jpt) , gdepw (jpi,jpj,jpk,jpt) , gde3w (jpi,jpj,jpk) , STAT=ierr(ii) )
ALLOCATE( gdept (jpi,jpj,jpk,jpt) , gdepw (jpi,jpj,jpk,jpt) , &
& gde3w_0(jpi,jpj,jpk) , gde3w (jpi,jpj,jpk) , STAT=ierr(ii) )
!
ii = ii+1
ALLOCATE( e3t(jpi,jpj,jpk,jpt) , e3u (jpi,jpj,jpk,jpt) , e3v (jpi,jpj,jpk,jpt) , e3f(jpi,jpj,jpk) , &
......
src/OCE/DOM/domain.F90
View file @ d1c6ae27
......@@ -64,7 +64,7 @@ MODULE domain
# include "do_loop_substitute.h90"
!!-------------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: domain.F90 15270 2021-09-17 14:27:55Z smasson $
!! $Id: domain.F90 14547 2021-02-25 17:07:15Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!-------------------------------------------------------------------------
CONTAINS
......@@ -88,7 +88,6 @@ CONTAINS
!
INTEGER :: ji, jj, jk, jt ! dummy loop indices
INTEGER :: iconf = 0 ! local integers
REAL(wp):: zrdt
CHARACTER (len=64) :: cform = "(A12, 3(A13, I7))"
INTEGER , DIMENSION(jpi,jpj) :: ik_top , ik_bot ! top and bottom ocean level
REAL(wp), DIMENSION(jpi,jpj) :: z1_hu_0, z1_hv_0
......@@ -310,8 +309,27 @@ CONTAINS
ENDIF
!
! set current model timestep rDt = 2*rn_Dt if MLF or rDt = rn_Dt if RK3
#if defined key_RK3
rDt = rn_Dt
r1_Dt = 1._wp / rDt
!
IF(lwp) THEN
WRITE(numout,*)
WRITE(numout,*) ' ===>>> Runge Kutta 3rd order (RK3) : rDt = ', rDt
WRITE(numout,*)
ENDIF
!
#else
rDt = 2._wp * rn_Dt
r1_Dt = 1._wp / rDt
!
IF(lwp) THEN
WRITE(numout,*)
WRITE(numout,*) ' ===>>> Modified Leap-Frog (MLF) : rDt = ', rDt
WRITE(numout,*)
ENDIF
!
#endif
!
IF( l_SAS .AND. .NOT.ln_linssh ) THEN
CALL ctl_warn( 'SAS requires linear ssh : force ln_linssh = T' )
......@@ -336,6 +354,9 @@ CONTAINS
IF( .NOT. Agrif_Root() ) THEN
nn_it000 = (Agrif_Parent(nn_it000)-1)*Agrif_IRhot() + 1
nn_itend = Agrif_Parent(nn_itend) *Agrif_IRhot()
nn_date0 = Agrif_Parent(nn_date0)
nn_time0 = Agrif_Parent(nn_time0)
nn_leapy = Agrif_Parent(nn_leapy)
ENDIF
#endif
!
......@@ -392,7 +413,16 @@ CONTAINS
IF( nn_wxios > 0 ) lwxios = .TRUE. !* set output file type for XIOS based on NEMO namelist
nxioso = nn_wxios
ENDIF
! !== Check consistency between ln_rstart and ln_1st_euler ==! (i.e. set l_1st_euler)
!
#if defined key_RK3
! !== RK3: Open the restart file ==!
IF( ln_rstart ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' open the restart file'
CALL rst_read_open
ENDIF
#else
! !== MLF: Check consistency between ln_rstart and ln_1st_euler ==! (i.e. set l_1st_euler)
l_1st_euler = ln_1st_euler
!
IF( ln_rstart ) THEN !* Restart case
......@@ -427,6 +457,7 @@ CONTAINS
IF(lwp) WRITE(numout,*)' an Euler initial time step is used : l_1st_euler is forced to .true. '
l_1st_euler = .TRUE.
ENDIF
#endif
!
! !== control of output frequency ==!
!
......
src/OCE/DOM/domqco.F90
View file @ d1c6ae27
......@@ -39,6 +39,7 @@ MODULE domqco
PUBLIC dom_qco_init ! called by domain.F90
PUBLIC dom_qco_zgr ! called by isfcpl.F90
PUBLIC dom_qco_r3c ! called by steplf.F90
PUBLIC dom_qco_r3c_RK3 ! called by stprk3_stg.F90
! !!* Namelist nam_vvl
LOGICAL , PUBLIC :: ln_vvl_zstar = .FALSE. ! zstar vertical coordinate
......@@ -123,7 +124,9 @@ CONTAINS
! ! Horizontal interpolation of e3t
#if defined key_RK3
CALL dom_qco_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb), r3f(:,:) )
CALL dom_qco_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm) )
r3t(:,:,Kmm) = r3t(:,:,Kbb) !!st r3 at Kmm needed to be initialised for Agrid_Grid call in nemo_gcm
r3u(:,:,Kmm) = r3u(:,:,Kbb) !! maybe we only need zeros ???
r3v(:,:,Kmm) = r3v(:,:,Kbb)
#else
CALL dom_qco_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb) )
CALL dom_qco_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm), r3f(:,:) )
......@@ -136,6 +139,61 @@ CONTAINS
SUBROUTINE dom_qco_r3c( pssh, pr3t, pr3u, pr3v, pr3f )
!!---------------------------------------------------------------------
!! *** ROUTINE r3c ***
!!
!! ** Purpose : compute the filtered ratio ssh/h_0 at t-,u-,v-,f-points
!!
!! ** Method : - compute the ssh at u- and v-points (f-point optional)
!! Vector Form : surface weighted averaging
!! Flux Form : simple averaging
!! - compute the ratio ssh/h_0 at t-,u-,v-pts, (f-pt optional)
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(:,:) , INTENT(in ) :: pssh ! sea surface height [m]
REAL(wp), DIMENSION(:,:) , INTENT( out) :: pr3t, pr3u, pr3v ! ssh/h0 ratio at t-, u-, v-,points [-]
REAL(wp), DIMENSION(:,:), OPTIONAL, INTENT( out) :: pr3f ! ssh/h0 ratio at f-point [-]
!
INTEGER :: ji, jj ! dummy loop indices
!!----------------------------------------------------------------------
!
!
DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
pr3t(ji,jj) = pssh(ji,jj) * r1_ht_0(ji,jj) !== ratio at t-point ==!
END_2D
!
! !== ratio at u-,v-point ==!
!
DO_2D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
pr3u(ji,jj) = 0.5_wp * ( e1e2t(ji ,jj) * pssh(ji ,jj) &
& + e1e2t(ji+1,jj) * pssh(ji+1,jj) ) * r1_hu_0(ji,jj) * r1_e1e2u(ji,jj)
pr3v(ji,jj) = 0.5_wp * ( e1e2t(ji,jj ) * pssh(ji,jj ) &
& + e1e2t(ji,jj+1) * pssh(ji,jj+1) ) * r1_hv_0(ji,jj) * r1_e1e2v(ji,jj)
END_2D
!
IF( .NOT.PRESENT( pr3f ) ) THEN !- lbc on ratio at u-, v-points only
IF (nn_hls==1) CALL lbc_lnk( 'dom_qco_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp )
!
!
ELSE !== ratio at f-point ==!
!
DO_2D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility
pr3f(ji,jj) = 0.25_wp * ( ( e1e2t(ji ,jj ) * pssh(ji ,jj ) &
& + e1e2t(ji+1,jj ) * pssh(ji+1,jj ) ) & ! bracket for halo 1 - halo 2 compatibility
& + ( e1e2t(ji ,jj+1) * pssh(ji ,jj+1) &
& + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1) ) & ! bracket for halo 1 - halo 2 compatibility
& ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj)
END_2D
! ! lbc on ratio at u-,v-,f-points
IF (nn_hls==1) CALL lbc_lnk( 'dom_qco_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp, pr3f, 'F', 1._wp )
!
ENDIF
!
END SUBROUTINE dom_qco_r3c
SUBROUTINE dom_qco_r3c_RK3( pssh, pr3t, pr3u, pr3v, pr3f )
!!---------------------------------------------------------------------
!! *** ROUTINE r3c ***
!!
......@@ -164,7 +222,7 @@ CONTAINS
!!st IF( ln_dynadv_vec ) THEN !- Vector Form (thickness weighted averaging)
#if ! defined key_qcoTest_FluxForm
! ! no 'key_qcoTest_FluxForm' : surface weighted ssh average
DO_2D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
DO_2D( 0, 0, 0, 0 )
pr3u(ji,jj) = 0.5_wp * ( e1e2t(ji ,jj) * pssh(ji ,jj) &
& + e1e2t(ji+1,jj) * pssh(ji+1,jj) ) * r1_hu_0(ji,jj) * r1_e1e2u(ji,jj)
pr3v(ji,jj) = 0.5_wp * ( e1e2t(ji,jj ) * pssh(ji,jj ) &
......@@ -172,52 +230,43 @@ CONTAINS
END_2D
!!st ELSE !- Flux Form (simple averaging)
#else
DO_2D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
DO_2D( 0, 0, 0, 0 )
pr3u(ji,jj) = 0.5_wp * ( pssh(ji,jj) + pssh(ji+1,jj ) ) * r1_hu_0(ji,jj)
pr3v(ji,jj) = 0.5_wp * ( pssh(ji,jj) + pssh(ji ,jj+1) ) * r1_hv_0(ji,jj)
END_2D
!!st ENDIF
#endif
!
IF( .NOT.PRESENT( pr3f ) ) THEN !- lbc on ratio at u-, v-points only
IF (nn_hls==1) CALL lbc_lnk( 'dom_qco_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp )
!
!
ELSE !== ratio at f-point ==!
IF( PRESENT( pr3f ) ) THEN !== ratio at f-point ==!
!
!!st IF( ln_dynadv_vec ) THEN !- Vector Form (thickness weighted averaging)
#if ! defined key_qcoTest_FluxForm
! ! no 'key_qcoTest_FluxForm' : surface weighted ssh average
DO_2D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
DO_2D( 0, 0, 0, 0 )
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility
pr3f(ji,jj) = 0.25_wp * ( ( e1e2t(ji ,jj ) * pssh(ji ,jj ) &
& + e1e2t(ji+1,jj ) * pssh(ji+1,jj ) &
& ) & ! bracket for halo 1 - halo 2 compatibility
& + ( e1e2t(ji ,jj+1) * pssh(ji ,jj+1) &
& + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1) &
& ) & ! bracket for halo 1 - halo 2 compatibility
pr3f(ji,jj) = 0.25_wp * ( ( e1e2t(ji ,jj ) * pssh(ji ,jj ) &
& + e1e2t(ji+1,jj ) * pssh(ji+1,jj ) ) & ! bracket for halo 1 - halo 2 compatibility
& + ( e1e2t(ji ,jj+1) * pssh(ji ,jj+1) &
& + e1e2t(ji+1,jj+1) * pssh(ji+1,jj+1) ) & ! bracket for halo 1 - halo 2 compatibility
& ) * r1_hf_0(ji,jj) * r1_e1e2f(ji,jj)
END_2D
!!st ELSE !- Flux Form (simple averaging)
#else
DO_2D_OVR( nn_hls, nn_hls-1, nn_hls, nn_hls-1 )
DO_2D( 0, 0, 0, 0 )
! round brackets added to fix the order of floating point operations
! needed to ensure halo 1 - halo 2 compatibility
pr3f(ji,jj) = 0.25_wp * ( ( pssh(ji,jj ) + pssh(ji+1,jj ) ) &
& + ( pssh(ji,jj+1) + pssh(ji+1,jj+1) &
& ) & ! bracket for halo 1 - halo 2 compatibility
pr3f(ji,jj) = 0.25_wp * ( ( pssh(ji,jj ) + pssh(ji+1,jj ) ) &
& + ( pssh(ji,jj+1) + pssh(ji+1,jj+1) ) & ! bracket for halo 1 - halo 2 compatibility
& ) * r1_hf_0(ji,jj)
END_2D
!!st ENDIF
#endif
! ! lbc on ratio at u-,v-,f-points
IF (nn_hls==1) CALL lbc_lnk( 'dom_qco_r3c', pr3u, 'U', 1._wp, pr3v, 'V', 1._wp, pr3f, 'F', 1._wp )
!
ENDIF
!
END SUBROUTINE dom_qco_r3c
END SUBROUTINE dom_qco_r3c_RK3
SUBROUTINE qco_ctl
......
src/OCE/DOM/domzgr.F90
View file @ d1c6ae27
......@@ -46,7 +46,7 @@ MODULE domzgr
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: domzgr.F90 15556 2021-11-29 15:23:06Z jchanut $
!! $Id: domzgr.F90 15157 2021-07-29 08:28:32Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -142,12 +142,17 @@ CONTAINS
CALL lbc_lnk( 'usrdef_zgr', zmsk, 'T', 1. ) ! set halos
k_top(:,:) = k_top(:,:) * NINT( zmsk(:,:) )
!
!!gm to be remove when removing the OLD definition of e3 scale factors so that gde3w disappears
#if ! defined key_qco && ! defined key_linssh
! OLD implementation of coordinate (not with 'key_qco' or 'key_linssh')
! gde3w_0 has to be defined
!!gm to be remove when removing the OLD definition of e3 scale factors so that gde3w_0=gdept_0
!!gm therefore gde3w_0 disappears
! Compute gde3w_0 (vertical sum of e3w)
gde3w_0(:,:,1) = 0.5_wp * e3w_0(:,:,1)
DO jk = 2, jpk
gde3w_0(:,:,jk) = gde3w_0(:,:,jk-1) + e3w_0(:,:,jk)
END DO
#endif
!
! Any closed seas (defined by closea_mask > 0 in domain_cfg file) to be filled
! in at runtime if ln_closea=.false.
......@@ -200,14 +205,20 @@ CONTAINS
WRITE(numout,*) ' MIN val k_top ', MINVAL( k_top(:,:) ), ' MAX ', MAXVAL( k_top(:,:) )
WRITE(numout,*) ' MIN val k_bot ', MINVAL( k_bot(:,:) ), ' MAX ', MAXVAL( k_bot(:,:) )
WRITE(numout,*) ' MIN val depth t ', MINVAL( gdept_0(:,:,:) ), &
& ' w ', MINVAL( gdepw_0(:,:,:) ), '3w ', MINVAL( gde3w_0(:,:,:) )
#if ! defined key_qco && ! defined key_linssh
& '3w ', MINVAL( gde3w_0(:,:,:) ), &
#endif
& ' w ', MINVAL( gdepw_0(:,:,:) )
WRITE(numout,*) ' MIN val e3 t ', MINVAL( e3t_0(:,:,:) ), ' f ', MINVAL( e3f_0(:,:,:) ), &
& ' u ', MINVAL( e3u_0(:,:,:) ), ' u ', MINVAL( e3v_0(:,:,:) ), &
& ' uw', MINVAL( e3uw_0(:,:,:) ), ' vw', MINVAL( e3vw_0(:,:,:)), &
& ' w ', MINVAL( e3w_0(:,:,:) )
WRITE(numout,*) ' MAX val depth t ', MAXVAL( gdept_0(:,:,:) ), &
& ' w ', MAXVAL( gdepw_0(:,:,:) ), '3w ', MAXVAL( gde3w_0(:,:,:) )
#if ! defined key_qco && ! defined key_linssh
& '3w ', MINVAL( gde3w_0(:,:,:) ), &
#endif
& ' w ', MINVAL( gdepw_0(:,:,:) )
WRITE(numout,*) ' MAX val e3 t ', MAXVAL( e3t_0(:,:,:) ), ' f ', MAXVAL( e3f_0(:,:,:) ), &
& ' u ', MAXVAL( e3u_0(:,:,:) ), ' u ', MAXVAL( e3v_0(:,:,:) ), &
& ' uw', MAXVAL( e3uw_0(:,:,:) ), ' vw', MAXVAL( e3vw_0(:,:,:) ), &
......
src/OCE/DOM/domzgr_substitute.h90
View file @ d1c6ae27
......@@ -51,4 +51,3 @@
# define gde3w(i,j,k) gdept_0(i,j,k)
#endif
!!----------------------------------------------------------------------
src/OCE/DOM/istate.F90
View file @ d1c6ae27
......@@ -50,7 +50,7 @@ MODULE istate
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: istate.F90 15052 2021-06-24 14:39:14Z smasson $
!! $Id: istate.F90 14991 2021-06-14 19:52:31Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -138,32 +138,47 @@ CONTAINS
ts (:,:,:,:,Kmm) = ts (:,:,:,:,Kbb) ! set now values from to before ones
uu (:,:,:,Kmm) = uu (:,:,:,Kbb)
vv (:,:,:,Kmm) = vv (:,:,:,Kbb)
ENDIF
#if defined key_agrif
ENDIF
#endif
!
! Initialize "now" and "before" barotropic velocities:
! Do it whatever the free surface method, these arrays being eventually used
#if defined key_RK3
IF( .NOT. ln_rstart ) THEN
#endif
! Initialize "before" barotropic velocities. "now" values are always set but
! "before" values may have been read from a restart to ensure restartability.
! In the non-restart or non-RK3 cases they need to be initialised here:
uu_b(:,:,Kbb) = 0._wp ; vv_b(:,:,Kbb) = 0._wp
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
uu_b(ji,jj,Kbb) = uu_b(ji,jj,Kbb) + e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb) * umask(ji,jj,jk)
vv_b(ji,jj,Kbb) = vv_b(ji,jj,Kbb) + e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb) * vmask(ji,jj,jk)
END_3D
uu_b(:,:,Kbb) = uu_b(:,:,Kbb) * r1_hu(:,:,Kbb)
vv_b(:,:,Kbb) = vv_b(:,:,Kbb) * r1_hv(:,:,Kbb)
!
#if defined key_RK3
ENDIF
#endif
!
uu_b(:,:,Kmm) = 0._wp ; vv_b(:,:,Kmm) = 0._wp
uu_b(:,:,Kbb) = 0._wp ; vv_b(:,:,Kbb) = 0._wp
! Initialize "now" barotropic velocities:
! Do it whatever the free surface method, these arrays being used eventually
!
!!gm the use of umsak & vmask is not necessary below as uu(:,:,:,Kmm), vv(:,:,:,Kmm), uu(:,:,:,Kbb), vv(:,:,:,Kbb) are always masked
#if defined key_RK3
IF( .NOT. ln_rstart ) THEN
uu_b(:,:,Kmm) = uu_b(:,:,Kbb) ! Kmm value set to Kbb for initialisation in Agrif_Regrid in namo_gcm
vv_b(:,:,Kmm) = vv_b(:,:,Kbb)
ENDIF
#else
!!gm the use of umask & vmask is not necessary below as uu(:,:,:,Kmm), vv(:,:,:,Kmm), uu(:,:,:,Kbb), vv(:,:,:,Kbb) are always masked
uu_b(:,:,Kmm) = 0._wp ; vv_b(:,:,Kmm) = 0._wp
DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
uu_b(ji,jj,Kmm) = uu_b(ji,jj,Kmm) + e3u(ji,jj,jk,Kmm) * uu(ji,jj,jk,Kmm) * umask(ji,jj,jk)
vv_b(ji,jj,Kmm) = vv_b(ji,jj,Kmm) + e3v(ji,jj,jk,Kmm) * vv(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
!
uu_b(ji,jj,Kbb) = uu_b(ji,jj,Kbb) + e3u(ji,jj,jk,Kbb) * uu(ji,jj,jk,Kbb) * umask(ji,jj,jk)
vv_b(ji,jj,Kbb) = vv_b(ji,jj,Kbb) + e3v(ji,jj,jk,Kbb) * vv(ji,jj,jk,Kbb) * vmask(ji,jj,jk)
END_3D
!
uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm)
vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv(:,:,Kmm)
!
uu_b(:,:,Kbb) = uu_b(:,:,Kbb) * r1_hu(:,:,Kbb)
vv_b(:,:,Kbb) = vv_b(:,:,Kbb) * r1_hv(:,:,Kbb)
#endif
!
END SUBROUTINE istate_init
......
src/OCE/DYN/divhor.F90
View file @ d1c6ae27
......@@ -12,6 +12,7 @@ MODULE divhor
!! 3.7 ! 2014-01 (G. Madec) suppression of velocity curl from in-core memory
!! - ! 2014-12 (G. Madec) suppression of cross land advection option
!! - ! 2015-10 (G. Madec) add velocity and rnf flag in argument of div_hor
!! 4.5 ! 2015-10 (S. Techene, G. Madec) hdiv replaced by e3divUh
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -35,19 +36,89 @@ MODULE divhor
IMPLICIT NONE
PRIVATE
PUBLIC div_hor ! routine called by step.F90 and istate.F90
! !! * Interface
INTERFACE div_hor
MODULE PROCEDURE div_hor_RK3, div_hor_old
END INTERFACE
PUBLIC div_hor ! routine called by ssh_nxt.F90 and istate.F90
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: divhor.F90 15150 2021-07-27 10:38:24Z smasson $
!! $Id: divhor.F90 14808 2021-05-07 12:00:45Z jchanut $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE div_hor( kt, Kbb, Kmm )
SUBROUTINE div_hor_RK3( kt, Kbb, Kmm, puu, pvv, pe3divUh )
!!----------------------------------------------------------------------
!! *** ROUTINE div_hor_RK3 ***
!!
!! ** Purpose : compute the horizontal divergence at now time-step
!!
!! ** Method : the now divergence is computed as :
!! hdiv = 1/(e1e2t*e3t) ( di[e2u*e3u un] + dj[e1v*e3v vn] )
!! and correct with runoff inflow (div_rnf) and cross land flow (div_cla)
!!
!! ** Action : - thickness weighted horizontal divergence of in input velocity (puu,pvv)
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt, Kbb, Kmm ! ocean time-step & time-level indices
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: puu, pvv ! horizontal velocity
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT( out) :: pe3divUh ! e3t*div[Uh]
!
INTEGER :: ji, jj, jk ! dummy loop indices
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('div_hor_RK3')
!
IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'div_hor_RK3 : thickness weighted horizontal divergence '
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~'
hdiv (:,:,:) = 0._wp ! initialize hdiv & pe3divUh for the halos and jpk level at the first time step
ENDIF
!
pe3divUh(:,:,:) = 0._wp !!gm to be applied to the halos only
!
DO_3D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls, 1, jpkm1 )
hdiv(ji,jj,jk) = ( e2u(ji ,jj) * e3u(ji ,jj,jk,Kmm) * puu(ji ,jj,jk) &
& - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * puu(ji-1,jj,jk) &
& + e1v(ji,jj ) * e3v(ji,jj ,jk,Kmm) * pvv(ji,jj ,jk) &
& - e1v(ji,jj-1) * e3v(ji,jj-1,jk,Kmm) * pvv(ji,jj-1,jk) ) &
& * r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm)
END_3D
!
IF( ln_rnf ) CALL sbc_rnf_div( hdiv, Kmm ) !== + runoffs divergence ==!
!
#if defined key_asminc
IF( ln_sshinc .AND. ln_asmiau ) & !== + SSH assimilation increment ==!
& CALL ssh_asm_div( kt, Kbb, Kmm, hdiv )
#endif
!
IF( ln_isf ) CALL isf_hdiv( kt, Kmm, hdiv ) !== + ice-shelf mass exchange ==!
!
IF( nn_hls==1 ) CALL lbc_lnk( 'divhor', hdiv, 'T', 1._wp ) ! (no sign change)
!
!!gm Patch before suppression of hdiv from all modules that use it
! DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !== e3t * Horizontal divergence ==!
! pe3divUh(ji,jj,jk) = hdiv(ji,jj,jk) * e3t(ji,jj,jk,Kmm)
! END_3D
!JC: over whole domain, and after lbclnk on hdiv to prevent from reproducibility issues
DO_3D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls, 1, jpkm1 )
pe3divUh(ji,jj,jk) = hdiv(ji,jj,jk) * e3t(ji,jj,jk,Kmm)
END_3D
!!gm end
!
!
IF( ln_timing ) CALL timing_stop('div_hor_RK3')
!
END SUBROUTINE div_hor_RK3
SUBROUTINE div_hor_old( kt, Kbb, Kmm )
!!----------------------------------------------------------------------
!! *** ROUTINE div_hor ***
!!
......@@ -102,7 +173,7 @@ CONTAINS
! ! needed for ww in sshwzv
IF( ln_timing ) CALL timing_stop('div_hor')
!
END SUBROUTINE div_hor
END SUBROUTINE div_hor_old
!!======================================================================
END MODULE divhor
src/OCE/DYN/dynadv.F90
View file @ d1c6ae27
......@@ -45,12 +45,12 @@ MODULE dynadv
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: dynadv.F90 14053 2020-12-03 13:48:38Z techene $
!! $Id: dynadv.F90 14419 2021-02-09 12:22:16Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dyn_adv( kt, Kbb, Kmm, puu, pvv, Krhs )
SUBROUTINE dyn_adv( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )
!!---------------------------------------------------------------------
!! *** ROUTINE dyn_adv ***
!!
......@@ -63,21 +63,22 @@ CONTAINS
!! it is the relative vorticity which is added to coriolis term
!! (see dynvor module).
!!----------------------------------------------------------------------
INTEGER , INTENT( in ) :: kt ! ocean time-step index
INTEGER , INTENT( in ) :: Kbb, 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 , 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.
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start( 'dyn_adv' )
!
SELECT CASE( n_dynadv ) !== compute advection trend and add it to general trend ==!
CASE( np_VEC_c2 )
CALL dyn_keg ( kt, nn_dynkeg, Kmm, puu, pvv, Krhs ) ! vector form : horizontal gradient of kinetic energy
CALL dyn_zad ( kt, Kmm, puu, pvv, Krhs ) ! vector form : vertical advection
CASE( np_FLX_c2 )
CALL dyn_adv_cen2( kt, Kmm, puu, pvv, Krhs ) ! 2nd order centered scheme
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
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 ) ! 3rd order UBS scheme (UP3)
CALL dyn_adv_ubs ( kt , Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad ) ! 3rd order UBS scheme (UP3)
END SELECT
!
IF( ln_timing ) CALL timing_stop( 'dyn_adv' )
......
src/OCE/DYN/dynadv_cen2.F90
View file @ d1c6ae27
......@@ -30,29 +30,36 @@ MODULE dynadv_cen2
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: dynadv_cen2.F90 14834 2021-05-11 09:24:44Z hadcv $
!! $Id: dynadv_cen2.F90 14419 2021-02-09 12:22:16Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dyn_adv_cen2( kt, Kmm, puu, pvv, Krhs )
SUBROUTINE dyn_adv_cen2( kt, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_adv_cen2 ***
!!
!! ** Purpose : Compute the now momentum advection trend in flux form
!! ** Purpose : Compute the momentum advection trend in flux form
!! and the general trend of the momentum equation.
!!
!! ** Method : Trend evaluated using now fields (centered in time)
!! ** Method : Trend evaluated with a 2nd order centered scheme
!! using fields at Kmm time-level.
!! 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(:,:,:,Krhs),pvv(:,:,:,Krhs)) updated with the now vorticity term trend
!! ** Action : (puu,pvv)(:,:,:,Krhs) updated with the advective trend
!!----------------------------------------------------------------------
INTEGER , INTENT( in ) :: kt ! ocean time-step index
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 , 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), 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, 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), DIMENSION(:,:,:) , POINTER :: zpt_u, zpt_v, zpt_w
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
......@@ -68,12 +75,22 @@ CONTAINS
zfv_vw(:,:,:) = pvv(:,:,:,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
!
! !== Horizontal advection ==!
!
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) * puu(ji,jj,jk,Kmm)
zfv(ji,jj,jk) = 0.25_wp * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Kmm)
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)
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) )
......@@ -83,11 +100,11 @@ CONTAINS
END_2D
DO_2D( 0, 0, 0, 0 ) ! divergence of horizontal momentum fluxes
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zfu_t(ji+1,jj,jk) - zfu_t(ji,jj ,jk) &
& + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * r1_e1e2u(ji,jj) &
& / e3u(ji,jj,jk,Kmm)
& + zfv_f(ji ,jj,jk) - zfv_f(ji,jj-1,jk) ) * r1_e1e2u(ji,jj) &
& / e3u(ji,jj,jk,Kmm)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zfu_f(ji,jj ,jk) - zfu_f(ji-1,jj,jk) &
& + zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,jk,Kmm)
& + zfv_t(ji,jj+1,jk) - zfv_t(ji ,jj,jk) ) * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,jk,Kmm)
END_2D
END DO
!
......@@ -99,42 +116,57 @@ CONTAINS
zfv_t(:,:,:) = pvv(:,:,:,Krhs)
ENDIF
!
! !== 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
END_2D
IF( ln_linssh ) THEN ! linear free surface: advection through the surface
DO_2D( 0, 0, 0, 0 )
zfu_uw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji+1,jj) * ww(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji,jj+1) * ww(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) * ww(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) )
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
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 ! 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 )
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) &
& / 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' )
!
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
......
src/OCE/DYN/dynadv_ubs.F90
View file @ d1c6ae27
......@@ -36,12 +36,12 @@ MODULE dynadv_ubs
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: dynadv_ubs.F90 14834 2021-05-11 09:24:44Z hadcv $
!! $Id: dynadv_ubs.F90 14419 2021-02-09 12:22:16Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dyn_adv_ubs( kt, Kbb, Kmm, puu, pvv, Krhs )
SUBROUTINE dyn_adv_ubs( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_adv_ubs ***
!!
......@@ -64,20 +64,26 @@ CONTAINS
!! Default value (hard coded in the begining of the module) are
!! gamma1=1/3 and gamma2=1/32.
!!
!! ** Action : - (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) updated with the 3D advective momentum trends
!! 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 ! ocean time-step index
INTEGER , INTENT( in ) :: Kbb, 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 , 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
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: zui, zvj, zfuj, zfvi, zl_u, zl_v ! local scalars
REAL(wp) :: zui, zvj, zfuj, zfvi, zl_u, zl_v, 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), DIMENSION(A2D(nn_hls),jpk,2) :: zlu_uu, zlu_uv
REAL(wp), DIMENSION(A2D(nn_hls),jpk,2) :: zlv_vv, zlv_vu
REAL(wp), DIMENSION(:,:,:), POINTER :: zpt_u, zpt_v, zpt_w
!!----------------------------------------------------------------------
!
IF( .NOT. l_istiled .OR. ntile == 1 ) THEN ! Do only on the first tile
......@@ -102,13 +108,24 @@ CONTAINS
zfu_uw(:,:,:) = puu(:,:,:,Krhs)
zfv_vw(:,:,:) = pvv(:,:,:,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( nn_hls, nn_hls, nn_hls, nn_hls )
zfu(ji,jj,jk) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Kmm)
zfv(ji,jj,jk) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Kmm)
zfu(ji,jj,jk) = e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * zpt_u(ji,jj,jk)
zfv(ji,jj,jk) = e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * zpt_v(ji,jj,jk)
END_2D
!
DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 ) ! laplacian
......@@ -157,8 +174,8 @@ CONTAINS
DO jk = 1, jpkm1 ! ====================== !
! ! horizontal volume fluxes
DO_2D( 1, 1, 1, 1 )
zfu(ji,jj,jk) = 0.25_wp * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * puu(ji,jj,jk,Kmm)
zfv(ji,jj,jk) = 0.25_wp * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * pvv(ji,jj,jk,Kmm)
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)
END_2D
!
DO_2D( 1, 0, 1, 0 ) ! horizontal momentum fluxes at T- and F-point
......@@ -212,42 +229,62 @@ CONTAINS
! ! ==================== !
! ! 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
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) * ww(ji,jj,1) + e1e2t(ji+1,jj) * ww(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
zfv_vw(ji,jj,1) = 0.5_wp * ( e1e2t(ji,jj) * ww(ji,jj,1) + e1e2t(ji,jj+1) * ww(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) * ww(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 )
! ! ======================== !
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
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
! ! 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
......
src/OCE/DYN/dynhpg.F90
View file @ d1c6ae27
......@@ -83,7 +83,7 @@ MODULE dynhpg
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: dynhpg.F90 15529 2021-11-23 15:00:19Z techene $
!! $Id: dynhpg.F90 15157 2021-07-29 08:28:32Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......
src/OCE/DYN/dynspg.F90
View file @ d1c6ae27
......@@ -51,12 +51,12 @@ MODULE dynspg
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: dynspg.F90 14225 2020-12-19 14:58:39Z smasson $
!! $Id: dynspg.F90 14547 2021-02-25 17:07:15Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dyn_spg( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa, k_only_ADV )
SUBROUTINE dyn_spg( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa )
!!----------------------------------------------------------------------
!! *** ROUTINE dyn_spg ***
!!
......@@ -78,10 +78,9 @@ CONTAINS
INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs, Kaa ! ocean time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
REAL(wp), DIMENSION(jpi,jpj,jpt) , INTENT(inout) :: pssh, puu_b, pvv_b ! SSH and barotropic velocities at main time levels
INTEGER , OPTIONAL , INTENT( in ) :: k_only_ADV ! only Advection in the RHS
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp) :: z2dt, zg_2, zintp, zgrho0r, zld ! local scalars
REAL(wp) :: zg_2, zintp, zgrho0r, zld ! local scalars
REAL(wp) , DIMENSION(jpi,jpj) :: zpgu, zpgv ! 2D workspace
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zpice
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv
......@@ -150,8 +149,8 @@ CONTAINS
!
IF( ln_wave .and. ln_bern_srfc ) THEN !== Add J terms: depth-independent Bernoulli head
DO_2D( 0, 0, 0, 0 )
zpgu(ji,jj) = zpgu(ji,jj) + ( bhd_wave(ji+1,jj) - bhd_wave(ji,jj) ) / e1u(ji,jj) !++ bhd_wave from wave model in m2/s2 [BHD parameters in WW3]
zpgv(ji,jj) = zpgv(ji,jj) + ( bhd_wave(ji,jj+1) - bhd_wave(ji,jj) ) / e2v(ji,jj)
zpgu(ji,jj) = zpgu(ji,jj) + ( bhd_wave(ji+1,jj) - bhd_wave(ji,jj) ) * r1_e1u(ji,jj) !++ bhd_wave from wave model in m2/s2 [BHD parameters in WW3]
zpgv(ji,jj) = zpgv(ji,jj) + ( bhd_wave(ji,jj+1) - bhd_wave(ji,jj) ) * r1_e2v(ji,jj)
END_2D
ENDIF
!
......@@ -166,7 +165,7 @@ CONTAINS
!
SELECT CASE ( nspg ) !== surface pressure gradient computed and add to the general trend ==!
CASE ( np_EXP ) ; CALL dyn_spg_exp( kt, Kmm, puu, pvv, Krhs ) ! explicit
CASE ( np_TS ) ; CALL dyn_spg_ts ( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa, k_only_ADV ) ! time-splitting
CASE ( np_TS ) ; CALL dyn_spg_ts ( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa ) ! time-splitting
END SELECT
!
IF( l_trddyn ) THEN ! save the surface pressure gradient trends for further diagnostics
......
src/OCE/DYN/dynspg_ts.F90
View file @ d1c6ae27
......@@ -67,9 +67,8 @@ MODULE dynspg_ts
PUBLIC dyn_spg_ts ! called by dyn_spg
PUBLIC dyn_spg_ts_init ! - - dyn_spg_init
PUBLIC dyn_drg_init ! called by stp2d
!! Time filtered arrays at baroclinic time step:
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: un_adv , vn_adv !: Advection vel. at "now" barocl. step
!
INTEGER, SAVE :: icycle ! Number of barotropic sub-steps for each internal step nn_e <= 2.5 nn_e
REAL(wp),SAVE :: rDt_e ! Barotropic time step
......@@ -79,6 +78,10 @@ MODULE dynspg_ts
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftnw, ftne ! triad of coriolis parameter
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ftsw, ftse ! (only used with een vorticity scheme)
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sshe_rhs ! RHS of ssh Eq.
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: Ue_rhs, Ve_rhs ! RHS of barotropic velocity Eq.
REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: CdU_u, CdU_v ! top/bottom stress at u- & v-points
REAL(wp) :: r1_12 = 1._wp / 12._wp ! local ratios
REAL(wp) :: r1_8 = 0.125_wp !
REAL(wp) :: r1_4 = 0.25_wp !
......@@ -89,7 +92,7 @@ MODULE dynspg_ts
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: dynspg_ts.F90 15489 2021-11-10 09:18:39Z jchanut $
!! $Id: dynspg_ts.F90 14747 2021-04-26 08:47:14Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -98,7 +101,7 @@ CONTAINS
!!----------------------------------------------------------------------
!! *** routine dyn_spg_ts_alloc ***
!!----------------------------------------------------------------------
INTEGER :: ierr(3)
INTEGER :: ierr(2)
!!----------------------------------------------------------------------
ierr(:) = 0
!
......@@ -106,7 +109,6 @@ CONTAINS
IF( ln_dynvor_een .OR. ln_dynvor_eeT ) &
& ALLOCATE( ftnw(jpi,jpj) , ftne(jpi,jpj) , ftsw(jpi,jpj) , ftse(jpi,jpj), STAT=ierr(2) )
!
ALLOCATE( un_adv(jpi,jpj), vn_adv(jpi,jpj) , STAT=ierr(3) )
!
dyn_spg_ts_alloc = MAXVAL( ierr(:) )
!
......@@ -116,7 +118,7 @@ CONTAINS
END FUNCTION dyn_spg_ts_alloc
SUBROUTINE dyn_spg_ts( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa, k_only_ADV )
SUBROUTINE dyn_spg_ts( kt, Kbb, Kmm, Krhs, puu, pvv, pssh, puu_b, pvv_b, Kaa )
!!----------------------------------------------------------------------
!!
!! ** Purpose : - Compute the now trend due to the explicit time stepping
......@@ -144,27 +146,24 @@ CONTAINS
INTEGER , INTENT( in ) :: kt ! ocean time-step index
INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs, Kaa ! ocean time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocities and RHS of momentum equation
REAL(wp), DIMENSION(jpi,jpj,jpt) , INTENT(inout) :: pssh, puu_b, pvv_b ! SSH and barotropic velocities at main time levels
INTEGER , OPTIONAL , INTENT( in ) :: k_only_ADV ! only Advection in the RHS
REAL(wp), DIMENSION(jpi,jpj ,jpt), INTENT(inout) :: pssh, puu_b, pvv_b ! SSH and barotropic velocities at main time levels
!
INTEGER :: ji, jj, jk, jn ! dummy loop indices
LOGICAL :: ll_fw_start ! =T : forward integration
LOGICAL :: ll_init ! =T : special startup of 2d equations
INTEGER :: noffset ! local integers : time offset for bdy update
REAL(wp) :: r1_Dt_b, z1_hu, z1_hv ! local scalars
REAL(wp) :: za0, za1, za2, za3 ! - -
REAL(wp) :: z1_hu , z1_hv ! local scalars
REAL(wp) :: zzsshu, zzsshv ! - -
REAL(wp) :: za0, za1, za2, za3 ! - -
REAL(wp) :: zztmp, zldg ! - -
REAL(wp) :: zhu_bck, zhv_bck, zhdiv ! - -
REAL(wp) :: zun_save, zvn_save ! - -
REAL(wp), DIMENSION(jpi,jpj) :: zu_trd, zu_frc, zu_spg
REAL(wp) :: zhu_bck, zhv_bck, zhdiv ! - -
REAL(wp) :: zun_save, zvn_save ! - -
REAL(wp), DIMENSION(jpi,jpj) :: zu_trd, zu_frc, zu_spg !!st tests , zssh_frc
REAL(wp), DIMENSION(jpi,jpj) :: zv_trd, zv_frc, zv_spg
REAL(wp), DIMENSION(jpi,jpj) :: zsshu_a, zhup2_e, zhtp2_e
REAL(wp), DIMENSION(jpi,jpj) :: zsshv_a, zhvp2_e, zsshp2_e
REAL(wp), DIMENSION(jpi,jpj) :: zCdU_u, zCdU_v ! top/bottom stress at u- & v-points
REAL(wp), DIMENSION(jpi,jpj) :: zhU, zhV ! fluxes
!!st#if defined key_qco
!!st REAL(wp), DIMENSION(jpi, jpj, jpk) :: ze3u, ze3v
!!st#endif
!
REAL(wp) :: zwdramp ! local scalar - only used if ln_wd_dl = .True.
......@@ -174,6 +173,7 @@ CONTAINS
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zcpx, zcpy ! Wetting/Dying gravity filter coef.
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: ztwdmask, zuwdmask, zvwdmask ! ROMS wetting and drying masks at t,u,v points
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zuwdav2, zvwdav2 ! averages over the sub-steps of zuwdmask and zvwdmask
REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: z2d ! 2D workspace
REAL(wp) :: zt0substep ! Time of day at the beginning of the time substep
!!----------------------------------------------------------------------
!
......@@ -184,7 +184,6 @@ CONTAINS
zwdramp = r_rn_wdmin1 ! simplest ramp
! zwdramp = 1._wp / (rn_wdmin2 - rn_wdmin1) ! more general ramp
! ! inverse of baroclinic time step
r1_Dt_b = 1._wp / rDt
!
ll_init = ln_bt_av ! if no time averaging, then no specific restart
ll_fw_start = .FALSE.
......@@ -227,17 +226,58 @@ CONTAINS
! -----------------------------------------------------------------------------
! Phase 1 : Coupling between general trend and barotropic estimates (1st step)
! -----------------------------------------------------------------------------
#if defined key_RK3
! !========================================!
! !== Phase 1 for RK3 time integration ==!
! !========================================!
!
! ! Currently, RK3 requires the forward mode
IF( kt == nit000 ) THEN
IF( .NOT.ln_bt_fw ) CALL ctl_stop( 'dyn_spg_ts: RK3 requires forward mode (ln_bt_fw=T)' )
ENDIF
!
! ! set values computed in RK3_ssh
ssh_frc(:,:) = sshe_rhs(:,:)
zu_frc(:,:) = Ue_rhs(:,:)
zv_frc(:,:) = Ve_rhs(:,:)
zCdU_u (:,:) = CdU_u (:,:)
zCdU_v (:,:) = CdU_v (:,:)
!!gm ==>>> !!ts ISSUe her on en discute changer les cas ENS ENE et triad ?
!
! != remove 2D Coriolis trend =!
! ! -------------------------- !
!
IF( kt == nit000 .OR. .NOT. ln_linssh ) CALL dyn_cor_2D_init( Kmm ) ! Set zwz, the barotropic Coriolis force coefficient
! ! recompute zwz = f/depth at every time step for (.NOT.ln_linssh) as the water colomn height changes
!
zhU(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:) ! now fluxes
zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:) ! NB: FULL domain : put a value in last row and column
!
CALL dyn_cor_2D( ht(:,:), hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV, & ! <<== in
& zu_trd, zv_trd ) ! ==>> out
!
DO_2D( 0, 0, 0, 0 ) ! Remove coriolis term (and possibly spg) from barotropic trend
zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
END_2D
#else
! !========================================!
! !== Phase 1 for MLF time integration ==!
! !========================================!
!
!
!
! != zu_frc = 1/H e3*d/dt(Ua) =! (Vertical mean of Ua, the 3D trends)
! ! --------------------------- !
#if defined key_qco
# if defined key_qco
zu_frc(:,:) = SUM( e3u_0(:,:,: ) * puu(:,:,:,Krhs) * umask(:,:,:), DIM=3 ) * r1_hu_0(:,:)
zv_frc(:,:) = SUM( e3v_0(:,:,: ) * pvv(:,:,:,Krhs) * vmask(:,:,:), DIM=3 ) * r1_hv_0(:,:)
#else
# else
zu_frc(:,:) = SUM( e3u(:,:,:,Kmm) * puu(:,:,:,Krhs) * umask(:,:,:), DIM=3 ) * r1_hu(:,:,Kmm)
zv_frc(:,:) = SUM( e3v(:,:,:,Kmm) * pvv(:,:,:,Krhs) * vmask(:,:,:), DIM=3 ) * r1_hv(:,:,Kmm)
#endif
# endif
!
!
! != U(Krhs) => baroclinic trend =! (remove its vertical mean)
......@@ -255,29 +295,21 @@ CONTAINS
IF( kt == nit000 .OR. .NOT. ln_linssh ) CALL dyn_cor_2D_init( Kmm ) ! Set zwz, the barotropic Coriolis force coefficient
! ! recompute zwz = f/depth at every time step for (.NOT.ln_linssh) as the water colomn height changes
!
IF( .NOT. PRESENT(k_only_ADV) ) THEN !* remove the 2D Coriolis trend
zhU(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:) ! now fluxes
zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:) ! NB: FULL domain : put a value in last row and column
!
CALL dyn_cor_2d( ht(:,:), hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV, & ! <<== in
& zu_trd, zv_trd ) ! ==>> out
!
DO_2D( 0, 0, 0, 0 ) ! Remove coriolis term (and possibly spg) from barotropic trend
zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
END_2D
ENDIF
zhU(:,:) = puu_b(:,:,Kmm) * hu(:,:,Kmm) * e2u(:,:) ! now fluxes
zhV(:,:) = pvv_b(:,:,Kmm) * hv(:,:,Kmm) * e1v(:,:) ! NB: FULL domain : put a value in last row and column
!
CALL dyn_cor_2D( ht(:,:), hu(:,:,Kmm), hv(:,:,Kmm), puu_b(:,:,Kmm), pvv_b(:,:,Kmm), zhU, zhV, & ! <<== in
& zu_trd, zv_trd ) ! ==>> out
!
DO_2D( 0, 0, 0, 0 ) ! Remove coriolis term (and possibly spg) from barotropic trend
zu_frc(ji,jj) = zu_frc(ji,jj) - zu_trd(ji,jj) * ssumask(ji,jj)
zv_frc(ji,jj) = zv_frc(ji,jj) - zv_trd(ji,jj) * ssvmask(ji,jj)
END_2D
!
! != Add bottom stress contribution from baroclinic velocities =!
! ! ----------------------------------------------------------- !
IF( PRESENT(k_only_ADV) ) THEN !* only Advection in the RHS : provide the barotropic bottom drag coefficients
DO_2D( 0, 0, 0, 0 )
zCdU_u(ji,jj) = r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) )
zCdU_v(ji,jj) = r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) )
END_2D
ELSE !* remove baroclinic drag AND provide the barotropic drag coefficients
CALL dyn_drg_init( Kbb, Kmm, puu, pvv, puu_b, pvv_b, zu_frc, zv_frc, zCdU_u, zCdU_v )
ENDIF
CALL dyn_drg_init( Kbb, Kmm, puu , pvv , puu_b , pvv_b , & ! <<= IN
& zu_frc, zv_frc, zCdU_u, zCdU_v ) ! =>> OUT
!
! != Add atmospheric pressure forcing =!
! ! ---------------------------------- !
......@@ -313,9 +345,9 @@ CONTAINS
END_2D
ENDIF
!
! !----------------!
! !== sssh_frc ==! Right-Hand-Side of the barotropic ssh equation (over the FULL domain)
! !----------------!
! !---------------!
! !== ssh_frc ==! Right-Hand-Side of the barotropic ssh equation (over the FULL domain)
! !---------------!
! != Net water flux forcing applied to a water column =!
! ! --------------------------------------------------- !
IF (ln_bt_fw) THEN ! FORWARD integration: use kt+1/2 fluxes (NOW+1/2)
......@@ -347,13 +379,18 @@ CONTAINS
!
END IF
!
#if defined key_asminc
# if defined key_asminc
! != Add the IAU weighted SSH increment =!
! ! ------------------------------------ !
IF( lk_asminc .AND. ln_sshinc .AND. ln_asmiau ) THEN
ssh_frc(:,:) = ssh_frc(:,:) - ssh_iau(:,:)
ENDIF
# endif
! !== END of Phase 1 for MLF time integration ==!
#endif
! != Fill boundary data arrays for AGRIF
! ! ------------------------------------
#if defined key_agrif
......@@ -383,7 +420,8 @@ CONTAINS
zhtp2_e(:,:) = ht_0(:,:)
ENDIF
!
IF( ln_bt_fw ) THEN ! FORWARD integration: start from NOW fields
IF( ln_bt_fw ) THEN ! FORWARD integration: start from NOW fields
! ! RK3: Kmm = Kbb when calling dynspg_ts
sshn_e(:,:) = pssh (:,:,Kmm)
un_e (:,:) = puu_b(:,:,Kmm)
vn_e (:,:) = pvv_b(:,:,Kmm)
......@@ -596,7 +634,7 @@ CONTAINS
! zwz array below or triads normally depend on sea level with ln_linssh=F and should be updated
! at each time step. We however keep them constant here for optimization.
! Recall that zhU and zhV hold fluxes at jn+0.5 (extrapolated not backward interpolated)
CALL dyn_cor_2d( zhtp2_e, zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV, zu_trd, zv_trd )
CALL dyn_cor_2D( zhtp2_e, zhup2_e, zhvp2_e, ua_e, va_e, zhU, zhV, zu_trd, zv_trd )
!
! Add tidal astronomical forcing if defined
IF ( ln_tide .AND. ln_tide_pot ) THEN
......@@ -712,34 +750,90 @@ CONTAINS
sshbb_e(:,:) = sshb_e(:,:)
sshb_e (:,:) = sshn_e(:,:)
sshn_e (:,:) = ssha_e(:,:)
! !* Sum over whole bt loop
!
! !* Sum over whole bt loop (except in weight average)
! ! ----------------------
za1 = wgtbtp1(jn)
IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities
puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:)
pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:)
ELSE ! Sum transports
IF ( .NOT.ln_wd_dl ) THEN
puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:)
pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:)
ELSE
puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:)
pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:)
END IF
IF( ln_bt_av ) THEN
za1 = wgtbtp1(jn)
IF( ln_dynadv_vec .OR. ln_linssh ) THEN ! Sum velocities
puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:)
pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:)
ELSE ! Sum transports
IF ( .NOT.ln_wd_dl ) THEN
puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:)
pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:)
ELSE
puu_b (:,:,Kaa) = puu_b (:,:,Kaa) + za1 * ua_e (:,:) * hu_e (:,:) * zuwdmask(:,:)
pvv_b (:,:,Kaa) = pvv_b (:,:,Kaa) + za1 * va_e (:,:) * hv_e (:,:) * zvwdmask(:,:)
ENDIF
ENDIF
! ! Sum sea level
pssh(:,:,Kaa) = pssh(:,:,Kaa) + za1 * ssha_e(:,:)
ENDIF
! ! Sum sea level
pssh(:,:,Kaa) = pssh(:,:,Kaa) + za1 * ssha_e(:,:)
!
! ! ==================== !
END DO ! end loop !
! ! ==================== !
! -----------------------------------------------------------------------------
! Phase 3. update the general trend with the barotropic trend
! -----------------------------------------------------------------------------
!
! Set advection velocity correction:
IF (ln_bt_fw) THEN
IF(.NOT.ln_bt_av ) THEN !* Update Kaa barotropic external mode
puu_b(:,:,Kaa) = ua_e (:,:)
pvv_b(:,:,Kaa) = va_e (:,:)
pssh (:,:,Kaa) = ssha_e(:,:)
ENDIF
#if defined key_RK3
! !* RK3 case
!
IF(.NOT.ln_dynadv_vec .AND. ln_bt_av ) THEN ! at this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points
!
# if defined key_qcoTest_FluxForm
! ! 'key_qcoTest_FluxForm' : simple ssh average
DO_2D( 0, 0, 0, 0 )
zzsshu = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji+1,jj ,Kaa) ) * ssumask(ji,jj)
zzsshv = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji ,jj+1,Kaa) ) * ssvmask(ji,jj)
!
! ! Save barotropic velocities (not transport)
puu_b(ji,jj,Kaa) = puu_b(ji,jj,Kaa) / ( hu_0(ji,jj) + zzsshu + 1._wp - ssumask(ji,jj) )
pvv_b(ji,jj,Kaa) = pvv_b(ji,jj,Kaa) / ( hv_0(ji,jj) + zzsshv + 1._wp - ssvmask(ji,jj) )
END_2D
# else
DO_2D( 0, 0, 0, 0 )
zzsshu = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj) * pssh(ji ,jj,Kaa) &
& + e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) ) * ssumask(ji,jj)
zzsshv = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj ) * pssh(ji,jj ,Kaa) &
& + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) * ssvmask(ji,jj)
!
! ! Save barotropic velocities (not transport)
puu_b(ji,jj,Kaa) = puu_b(ji,jj,Kaa) / ( hu_0(ji,jj) + zzsshu + 1._wp - ssumask(ji,jj) )
pvv_b(ji,jj,Kaa) = pvv_b(ji,jj,Kaa) / ( hv_0(ji,jj) + zzsshv + 1._wp - ssvmask(ji,jj) )
END_2D
# endif
!
CALL lbc_lnk( 'dynspg_ts', puu_b, 'U', -1._wp, pvv_b, 'V', -1._wp ) ! Boundary conditions
!
ENDIF
!
IF( iom_use("ubar") ) THEN ! RK3 single first: hu[N+1/2] = 1/2 ( hu[N] + hu[N+1] )
ALLOCATE( z2d(jpi,jpj) )
z2d(:,:) = 2._wp / ( hu_e(:,:) + hu(:,:,Kbb) + 1._wp - ssumask(:,:) )
CALL iom_put( "ubar", un_adv(:,:)*z2d(:,:) ) ! barotropic i-current
z2d(:,:) = 2._wp / ( hv_e(:,:) + hv(:,:,Kbb) + 1._wp - ssvmask(:,:) )
CALL iom_put( "vbar", vn_adv(:,:)*z2d(:,:) ) ! barotropic i-current
DEALLOCATE( z2d )
ENDIF
!
! !== END Phase 3 for RK3 (forward mode) ==!
#else
! !* MLF case
!
! Set advective velocity correction:
IF( ln_bt_fw ) THEN
IF( .NOT.( kt == nit000 .AND. l_1st_euler ) ) THEN
DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
zun_save = un_adv(ji,jj)
......@@ -759,44 +853,53 @@ CONTAINS
vn_bf(:,:) = 0._wp
ub2_b(:,:) = un_adv(:,:) ! Save integrated transport for next computation
vb2_b(:,:) = vn_adv(:,:)
END IF
ENDIF
ENDIF
!
! Update barotropic trend:
IF( ln_dynadv_vec .OR. ln_linssh ) THEN
DO jk=1,jpkm1
puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) ) * r1_Dt_b
pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) ) * r1_Dt_b
puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) ) * r1_Dt
pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) ) * r1_Dt
END DO
ELSE
! At this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points
#if defined key_qcoTest_FluxForm
! ! 'key_qcoTest_FluxForm' : simple ssh average
DO_2D( 1, 0, 1, 0 )
zsshu_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji+1,jj ,Kaa) ) * ssumask(ji,jj)
zsshv_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji ,jj+1,Kaa) ) * ssvmask(ji,jj)
END_2D
#else
DO_2D( 1, 0, 1, 0 )
zsshu_a(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj) * pssh(ji ,jj,Kaa) &
& + e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) ) * ssumask(ji,jj)
zsshv_a(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj ) * pssh(ji,jj ,Kaa) &
& + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) * ssvmask(ji,jj)
END_2D
#endif
CALL lbc_lnk( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
!
DO jk=1,jpkm1
puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) &
& * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt_b
pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) &
& * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt_b
END DO
! Save barotropic velocities not transport:
puu_b(:,:,Kaa) = puu_b(:,:,Kaa) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
pvv_b(:,:,Kaa) = pvv_b(:,:,Kaa) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
IF(.NOT.ln_bt_av ) THEN ! (puu_b,pvv_b)_Kaa is a velocity (hu,hv)_Kaa = (hu_e,hv_e)
!
DO jk=1,jpkm1
puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) &
& * ( puu_b(:,:,Kaa)*hu_e(:,:) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt
pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) &
& * ( pvv_b(:,:,Kaa)*hv_e(:,:) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt
END DO
!
ELSE ! at this stage, pssh(:,:,:,Krhs) has been corrected: compute new depths at velocity points
!
# if defined key_qcoTest_FluxForm
! ! 'key_qcoTest_FluxForm' : simple ssh average
DO_2D( 1, 0, 1, 0 )
zsshu_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji+1,jj ,Kaa) ) * ssumask(ji,jj)
zsshv_a(ji,jj) = r1_2 * ( pssh(ji,jj,Kaa) + pssh(ji ,jj+1,Kaa) ) * ssvmask(ji,jj)
END_2D
# else
DO_2D( 1, 0, 1, 0 )
zsshu_a(ji,jj) = r1_2 * r1_e1e2u(ji,jj) * ( e1e2t(ji ,jj) * pssh(ji ,jj,Kaa) &
& + e1e2t(ji+1,jj) * pssh(ji+1,jj,Kaa) ) * ssumask(ji,jj)
zsshv_a(ji,jj) = r1_2 * r1_e1e2v(ji,jj) * ( e1e2t(ji,jj ) * pssh(ji,jj ,Kaa) &
& + e1e2t(ji,jj+1) * pssh(ji,jj+1,Kaa) ) * ssvmask(ji,jj)
END_2D
# endif
CALL lbc_lnk( 'dynspg_ts', zsshu_a, 'U', 1._wp, zsshv_a, 'V', 1._wp ) ! Boundary conditions
!
DO jk=1,jpkm1
puu(:,:,jk,Krhs) = puu(:,:,jk,Krhs) + r1_hu(:,:,Kmm) &
& * ( puu_b(:,:,Kaa) - puu_b(:,:,Kbb) * hu(:,:,Kbb) ) * r1_Dt
pvv(:,:,jk,Krhs) = pvv(:,:,jk,Krhs) + r1_hv(:,:,Kmm) &
& * ( pvv_b(:,:,Kaa) - pvv_b(:,:,Kbb) * hv(:,:,Kbb) ) * r1_Dt
END DO
! Save barotropic velocities not transport:
puu_b(:,:,Kaa) = puu_b(:,:,Kaa) / ( hu_0(:,:) + zsshu_a(:,:) + 1._wp - ssumask(:,:) )
pvv_b(:,:,Kaa) = pvv_b(:,:,Kaa) / ( hv_0(:,:) + zsshv_a(:,:) + 1._wp - ssvmask(:,:) )
ENDIF
ENDIF
......@@ -806,36 +909,46 @@ CONTAINS
pvv(:,:,jk,Kmm) = ( pvv(:,:,jk,Kmm) + vn_adv(:,:)*r1_hv(:,:,Kmm) - pvv_b(:,:,Kmm) ) * vmask(:,:,jk)
END DO
IF ( ln_wd_dl .and. ln_wd_dl_bc) THEN
IF( ln_wd_dl .AND. ln_wd_dl_bc) THEN
DO jk = 1, jpkm1
puu(:,:,jk,Kmm) = ( un_adv(:,:)*r1_hu(:,:,Kmm) &
& + zuwdav2(:,:)*(puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm)) ) * umask(:,:,jk)
& + zuwdav2(:,:)*(puu(:,:,jk,Kmm) - un_adv(:,:)*r1_hu(:,:,Kmm)) ) * umask(:,:,jk)
pvv(:,:,jk,Kmm) = ( vn_adv(:,:)*r1_hv(:,:,Kmm) &
& + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm)) ) * vmask(:,:,jk)
& + zvwdav2(:,:)*(pvv(:,:,jk,Kmm) - vn_adv(:,:)*r1_hv(:,:,Kmm)) ) * vmask(:,:,jk)
END DO
END IF
ENDIF
CALL iom_put( "ubar", un_adv(:,:)*r1_hu(:,:,Kmm) ) ! barotropic i-current
CALL iom_put( "vbar", vn_adv(:,:)*r1_hv(:,:,Kmm) ) ! barotropic i-current
CALL iom_put( "vbar", vn_adv(:,:)*r1_hv(:,:,Kmm) ) ! barotropic j-current
! !== END Phase 3 for MLF time integration ==!
#endif
!
#if defined key_agrif
!
! Save time integrated fluxes during child grid integration
! (used to update coarse grid transports at next time step)
!
IF( .NOT.Agrif_Root() .AND. ln_bt_fw .AND. ln_agrif_2way ) THEN
IF( .NOT.Agrif_Root() .AND. ln_agrif_2way ) THEN
IF( Agrif_NbStepint() == 0 ) THEN
ub2_i_b(:,:) = 0._wp
vb2_i_b(:,:) = 0._wp
END IF
!
za1 = 1._wp / REAL(Agrif_rhot(), wp)
# if defined key_RK3
ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * un_adv(:,:)
vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vn_adv(:,:)
# else
ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
# endif
ENDIF
#endif
! !* write time-spliting arrays in the restart
IF( lrst_oce .AND.ln_bt_fw ) CALL ts_rst( kt, 'WRITE' )
#endif
! !: write time-spliting arrays in the restart
IF( lrst_oce ) CALL ts_rst( kt, 'WRITE' )
!
IF( ln_wd_il ) DEALLOCATE( zcpx, zcpy )
IF( ln_wd_dl ) DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
......@@ -845,80 +958,80 @@ CONTAINS
!
END SUBROUTINE dyn_spg_ts
SUBROUTINE ts_wgt( ll_av, ll_fw, jpit, zwgt1, zwgt2)
SUBROUTINE ts_wgt( ll_av, ll_fw, Kpit, zwgt1, zwgt2)
!!---------------------------------------------------------------------
!! *** ROUTINE ts_wgt ***
!!
!! ** Purpose : Set time-splitting weights for temporal averaging (or not)
!!----------------------------------------------------------------------
LOGICAL, INTENT(in) :: ll_av ! temporal averaging=.true.
LOGICAL, INTENT(in) :: ll_fw ! forward time splitting =.true.
INTEGER, INTENT(inout) :: jpit ! cycle length
REAL(wp), DIMENSION(3*nn_e), INTENT(inout) :: zwgt1, & ! Primary weights
zwgt2 ! Secondary weights
INTEGER :: jic, jn, ji ! temporary integers
REAL(wp) :: za1, za2
LOGICAL, INTENT(in ) :: ll_av ! temporal averaging=.true.
LOGICAL, INTENT(in ) :: ll_fw ! forward time splitting =.true.
INTEGER, INTENT(inout) :: Kpit ! cycle length
!!
INTEGER :: jic, jn, ji ! local integers
REAL(wp) :: za1, za2 ! loca scalars
REAL(wp), DIMENSION(3*nn_e), INTENT(inout) :: zwgt1, zwgt2 ! Primary & Secondary weights
!!----------------------------------------------------------------------
!
zwgt1(:) = 0._wp
zwgt2(:) = 0._wp
! Set time index when averaged value is requested
IF (ll_fw) THEN
jic = nn_e
ELSE
jic = 2 * nn_e
!
! !== Set time index when averaged value is requested ==!
IF (ll_fw) THEN ; jic = nn_e
ELSE ; jic = 2 * nn_e
ENDIF
! Set primary weights:
IF (ll_av) THEN
! Define simple boxcar window for primary weights
! (width = nn_e, centered around jic)
SELECT CASE ( nn_bt_flt )
CASE( 0 ) ! No averaging
zwgt1(jic) = 1._wp
jpit = jic
CASE( 1 ) ! Boxcar, width = nn_e
DO jn = 1, 3*nn_e
za1 = ABS(float(jn-jic))/float(nn_e)
IF (za1 < 0.5_wp) THEN
zwgt1(jn) = 1._wp
jpit = jn
ENDIF
ENDDO
CASE( 2 ) ! Boxcar, width = 2 * nn_e
DO jn = 1, 3*nn_e
za1 = ABS(float(jn-jic))/float(nn_e)
IF (za1 < 1._wp) THEN
zwgt1(jn) = 1._wp
jpit = jn
ENDIF
ENDDO
CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
!
! !== Set primary weights ==!
!
IF (ll_av) THEN != Define simple boxcar window for primary weights
! ! (width = nn_e, centered around jic)
SELECT CASE( nn_bt_flt )
!
CASE( 0 ) ! No averaging
zwgt1(jic) = 1._wp
Kpit = jic
!
CASE( 1 ) ! Boxcar, width = nn_e
DO jn = 1, 3*nn_e
za1 = ABS( REAL( jn-jic, wp) ) / REAL( nn_e, wp )
IF( za1 < 0.5_wp ) THEN
zwgt1(jn) = 1._wp
Kpit = jn
ENDIF
END DO
!
CASE( 2 ) ! Boxcar, width = 2 * nn_e
DO jn = 1, 3*nn_e
za1 = ABS(REAL( jn-jic, wp) ) / REAL( nn_e, wp )
IF( za1 < 1._wp ) THEN
zwgt1(jn) = 1._wp
Kpit = jn
ENDIF
END DO
!
CASE DEFAULT ; CALL ctl_stop( 'unrecognised value for nn_bt_flt' )
!
END SELECT
ELSE ! No time averaging
ELSE != No time averaging
zwgt1(jic) = 1._wp
jpit = jic
Kpit = jic
ENDIF
! Set secondary weights
DO jn = 1, jpit
DO ji = jn, jpit
zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
END DO
!
! !== Set secondary weights ==!
DO jn = 1, Kpit
DO ji = jn, Kpit
zwgt2(jn) = zwgt2(jn) + zwgt1(ji)
END DO
END DO
! Normalize weigths:
za1 = 1._wp / SUM(zwgt1(1:jpit))
za2 = 1._wp / SUM(zwgt2(1:jpit))
DO jn = 1, jpit
zwgt1(jn) = zwgt1(jn) * za1
zwgt2(jn) = zwgt2(jn) * za2
!
! !== Normalize weigths ==!
za1 = 1._wp / SUM( zwgt1(1:Kpit) )
za2 = 1._wp / SUM( zwgt2(1:Kpit) )
DO jn = 1, Kpit
zwgt1(jn) = zwgt1(jn) * za1
zwgt2(jn) = zwgt2(jn) * za2
END DO
!
END SUBROUTINE ts_wgt
......@@ -936,11 +1049,15 @@ CONTAINS
!
IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise
! ! ---------------
IF( ln_rstart .AND. ln_bt_fw .AND. .NOT.l_1st_euler ) THEN !* Read the restart file
CALL iom_get( numror, jpdom_auto, 'ub2_b' , ub2_b (:,:), cd_type = 'U', psgn = -1._wp )
CALL iom_get( numror, jpdom_auto, 'vb2_b' , vb2_b (:,:), cd_type = 'V', psgn = -1._wp )
CALL iom_get( numror, jpdom_auto, 'un_bf' , un_bf (:,:), cd_type = 'U', psgn = -1._wp )
CALL iom_get( numror, jpdom_auto, 'vn_bf' , vn_bf (:,:), cd_type = 'V', psgn = -1._wp )
IF( ln_rstart .AND. .NOT.l_1st_euler ) THEN !* Read the restart file
# if ! defined key_RK3
IF ( ln_bt_fw ) THEN
CALL iom_get( numror, jpdom_auto, 'ub2_b' , ub2_b (:,:), cd_type = 'U', psgn = -1._wp )
CALL iom_get( numror, jpdom_auto, 'vb2_b' , vb2_b (:,:), cd_type = 'V', psgn = -1._wp )
CALL iom_get( numror, jpdom_auto, 'un_bf' , un_bf (:,:), cd_type = 'U', psgn = -1._wp )
CALL iom_get( numror, jpdom_auto, 'vn_bf' , vn_bf (:,:), cd_type = 'V', psgn = -1._wp )
ENDIF
# endif
IF( .NOT.ln_bt_av ) THEN
CALL iom_get( numror, jpdom_auto, 'sshbb_e' , sshbb_e(:,:), cd_type = 'T', psgn = 1._wp )
CALL iom_get( numror, jpdom_auto, 'ubb_e' , ubb_e(:,:), cd_type = 'U', psgn = -1._wp )
......@@ -957,13 +1074,21 @@ CONTAINS
ELSE
ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif
ENDIF
# if defined key_RK3
CALL iom_get( numror, jpdom_auto, 'un_adv' , un_adv(:,:), cd_type = 'U', psgn = -1._wp )
CALL iom_get( numror, jpdom_auto, 'vn_adv' , vn_adv(:,:), cd_type = 'V', psgn = -1._wp )
# endif
#endif
ELSE !* Start from rest
ELSE
! !* Start from rest or use RK3 time-step
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) ' ==>>> start from rest: set barotropic values to 0'
# if ! defined key_RK3
ub2_b (:,:) = 0._wp ; vb2_b (:,:) = 0._wp ! used in the 1st interpol of agrif
un_adv (:,:) = 0._wp ; vn_adv (:,:) = 0._wp ! used in the 1st interpol of agrif
un_bf (:,:) = 0._wp ; vn_bf (:,:) = 0._wp ! used in the 1st update of agrif
#else
un_adv (:,:) = 0._wp ; vn_adv (:,:) = 0._wp ! used in the 1st interpol of agrif
#endif
#if defined key_agrif
ub2_i_b(:,:) = 0._wp ; vb2_i_b(:,:) = 0._wp ! used in the 1st update of agrif
#endif
......@@ -972,10 +1097,14 @@ CONTAINS
ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN ! Create restart file
! ! -------------------
IF(lwp) WRITE(numout,*) '---- ts_rst ----'
CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:) )
CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:) )
CALL iom_rstput( kt, nitrst, numrow, 'un_bf' , un_bf (:,:) )
CALL iom_rstput( kt, nitrst, numrow, 'vn_bf' , vn_bf (:,:) )
# if ! defined key_RK3
IF ( ln_bt_fw ) THEN
CALL iom_rstput( kt, nitrst, numrow, 'ub2_b' , ub2_b (:,:) )
CALL iom_rstput( kt, nitrst, numrow, 'vb2_b' , vb2_b (:,:) )
CALL iom_rstput( kt, nitrst, numrow, 'un_bf' , un_bf (:,:) )
CALL iom_rstput( kt, nitrst, numrow, 'vn_bf' , vn_bf (:,:) )
ENDIF
# endif
!
IF (.NOT.ln_bt_av) THEN
CALL iom_rstput( kt, nitrst, numrow, 'sshbb_e' , sshbb_e(:,:) )
......@@ -991,6 +1120,10 @@ CONTAINS
CALL iom_rstput( kt, nitrst, numrow, 'ub2_i_b' , ub2_i_b(:,:) )
CALL iom_rstput( kt, nitrst, numrow, 'vb2_i_b' , vb2_i_b(:,:) )
ENDIF
# if defined key_RK3
CALL iom_rstput( kt, nitrst, numrow, 'un_adv' , un_adv(:,:) )
CALL iom_rstput( kt, nitrst, numrow, 'vn_adv' , vn_adv(:,:) )
# endif
#endif
ENDIF
!
......@@ -1034,19 +1167,23 @@ CONTAINS
ELSE
IF(lwp) WRITE(numout,*) ' ln_ts_auto=.false.: Use nn_e in namelist nn_e = ', nn_e
ENDIF
!
IF(ln_bt_av) THEN
IF(lwp) WRITE(numout,*) ' ln_bt_av =.true. ==> Time averaging over nn_e time steps is on '
ELSE
IF(lwp) WRITE(numout,*) ' ln_bt_av =.false. => No time averaging of barotropic variables '
ENDIF
!
!
# if ! defined key_RK3
IF(ln_bt_fw) THEN
IF(lwp) WRITE(numout,*) ' ln_bt_fw=.true. => Forward integration of barotropic variables '
ELSE
IF(lwp) WRITE(numout,*) ' ln_bt_fw =.false.=> Centred integration of barotropic variables '
ENDIF
# else
! Enforce ln_bt_fw = T with RK3
ln_bt_fw = .true.
# endif
!
#if defined key_agrif
! Restrict the use of Agrif to the forward case only
......@@ -1081,7 +1218,7 @@ CONTAINS
! ! Allocate time-splitting arrays
IF( dyn_spg_ts_alloc() /= 0 ) CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts arrays' )
!
! ! read restart when needed
! !: restart/initialise
CALL ts_rst( nit000, 'READ' )
!
END SUBROUTINE dyn_spg_ts_init
......@@ -1099,7 +1236,8 @@ CONTAINS
!! To remove this approximation, copy lines below inside barotropic loop
!! and update depths at T- points (ht) at each barotropic time step
!!
!! Compute zwz = f / ( height of the water colomn )
!! Compute zwz = f/h (potential planetary voricity)
!! Compute ftne, ftnw, ftse, ftsw (triad of potential planetary voricity)
!!----------------------------------------------------------------------
INTEGER, INTENT(in) :: Kmm ! Time index
INTEGER :: ji ,jj, jk ! dummy loop indices
......@@ -1149,13 +1287,13 @@ CONTAINS
END_2D
!
END SELECT
END SUBROUTINE dyn_cor_2d_init
!
END SUBROUTINE dyn_cor_2D_init
SUBROUTINE dyn_cor_2d( pht, phu, phv, punb, pvnb, zhU, zhV, zu_trd, zv_trd )
SUBROUTINE dyn_cor_2D( pht, phu, phv, punb, pvnb, zhU, zhV, zu_trd, zv_trd )
!!---------------------------------------------------------------------
!! *** ROUTINE dyn_cor_2d ***
!! *** ROUTINE dyn_cor_2D ***
!!
!! ** Purpose : Compute u and v coriolis trends
!!----------------------------------------------------------------------
......@@ -1298,11 +1436,13 @@ CONTAINS
!!
!! ** Purpose :
!!----------------------------------------------------------------------
INTEGER :: ji ,jj ! dummy loop indices
LOGICAL :: ll_tmp1, ll_tmp2
REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pshn
REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) :: zcpx, zcpy
!
INTEGER :: ji ,jj ! dummy loop indices
LOGICAL :: ll_tmp1, ll_tmp2
!!----------------------------------------------------------------------
!
DO_2D( 0, 0, 0, 0 )
ll_tmp1 = MIN( pshn(ji,jj) , pshn(ji+1,jj) ) > &
& MAX( -ht_0(ji,jj) , -ht_0(ji+1,jj) ) .AND. &
......@@ -1341,7 +1481,7 @@ CONTAINS
zcpy(ji,jj) = 0._wp
ENDIF
END_2D
!
END SUBROUTINE wad_spg
......@@ -1385,15 +1525,13 @@ CONTAINS
! !== BOTTOM stress contribution from baroclinic velocities ==!
!
IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW bottom baroclinic velocities
DO_2D( 0, 0, 0, 0 )
ikbu = mbku(ji,jj)
ikbv = mbkv(ji,jj)
zu_i(ji,jj) = puu(ji,jj,ikbu,Kmm) - puu_b(ji,jj,Kmm)
zv_i(ji,jj) = pvv(ji,jj,ikbv,Kmm) - pvv_b(ji,jj,Kmm)
END_2D
ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities
ELSE ! CENTRED integration: use BEFORE bottom baroclinic velocities
DO_2D( 0, 0, 0, 0 )
ikbu = mbku(ji,jj)
ikbv = mbkv(ji,jj)
......@@ -1411,7 +1549,6 @@ CONTAINS
& r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) , zztmp )
END_2D
ELSE ! use "unclipped" drag (even if explicit friction is used in 3D calculation)
DO_2D( 0, 0, 0, 0 )
pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * zu_i(ji,jj)
pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * zv_i(ji,jj)
......@@ -1423,15 +1560,13 @@ CONTAINS
IF( ln_isfcav.OR.ln_drgice_imp ) THEN
!
IF( ln_bt_fw ) THEN ! FORWARD integration: use NOW top baroclinic velocity
DO_2D( 0, 0, 0, 0 )
iktu = miku(ji,jj)
iktv = mikv(ji,jj)
zu_i(ji,jj) = puu(ji,jj,iktu,Kmm) - puu_b(ji,jj,Kmm)
zv_i(ji,jj) = pvv(ji,jj,iktv,Kmm) - pvv_b(ji,jj,Kmm)
END_2D
ELSE ! CENTRED integration: use BEFORE top baroclinic velocity
ELSE ! CENTRED integration: use BEFORE top baroclinic velocity
DO_2D( 0, 0, 0, 0 )
iktu = miku(ji,jj)
iktv = mikv(ji,jj)
......@@ -1439,9 +1574,7 @@ CONTAINS
zv_i(ji,jj) = pvv(ji,jj,iktv,Kbb) - pvv_b(ji,jj,Kbb)
END_2D
ENDIF
!
! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
! ! use "unclipped" top drag (even if explicit friction is used in 3D calculation)
DO_2D( 0, 0, 0, 0 )
pu_RHSi(ji,jj) = pu_RHSi(ji,jj) + r1_hu(ji,jj,Kmm) * r1_2*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * zu_i(ji,jj)
pv_RHSi(ji,jj) = pv_RHSi(ji,jj) + r1_hv(ji,jj,Kmm) * r1_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * zv_i(ji,jj)
......@@ -1451,6 +1584,7 @@ CONTAINS
!
END SUBROUTINE dyn_drg_init
SUBROUTINE ts_bck_interp( jn, ll_init, & ! <== in
& za0, za1, za2, za3 ) ! ==> out
!!----------------------------------------------------------------------
......@@ -1488,6 +1622,5 @@ CONTAINS
ENDIF
END SUBROUTINE ts_bck_interp
!!======================================================================
END MODULE dynspg_ts
src/OCE/DYN/dynvor.F90
View file @ d1c6ae27
......@@ -50,6 +50,10 @@ MODULE dynvor
IMPLICIT NONE
PRIVATE
INTERFACE dyn_vor
MODULE PROCEDURE dyn_vor_3D, dyn_vor_RK3
END INTERFACE
PUBLIC dyn_vor ! routine called by step.F90
PUBLIC dyn_vor_init ! routine called by nemogcm.F90
......@@ -73,8 +77,9 @@ MODULE dynvor
INTEGER, PUBLIC, PARAMETER :: np_EEN = 4 ! EEN scheme
INTEGER, PUBLIC, PARAMETER :: np_MIX = 5 ! MIX scheme
INTEGER :: ncor, nrvm, ntot ! choice of calculated vorticity
! ! associated indices:
! !: choice of calculated vorticity
INTEGER, PUBLIC :: ncor, nrvm, ntot ! Coriolis, relative vorticity, total vorticity
! ! associated indices:
INTEGER, PUBLIC, PARAMETER :: np_COR = 1 ! Coriolis (planetary)
INTEGER, PUBLIC, PARAMETER :: np_RVO = 2 ! relative vorticity
INTEGER, PUBLIC, PARAMETER :: np_MET = 3 ! metric term
......@@ -98,12 +103,12 @@ MODULE dynvor
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: dynvor.F90 14834 2021-05-11 09:24:44Z hadcv $
!! $Id: dynvor.F90 14547 2021-02-25 17:07:15Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dyn_vor( kt, Kmm, puu, pvv, Krhs )
SUBROUTINE dyn_vor_3D( kt, Kmm, puu, pvv, Krhs )
!!----------------------------------------------------------------------
!!
!! ** Purpose : compute the lateral ocean tracer physics.
......@@ -120,7 +125,7 @@ CONTAINS
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztrdu, ztrdv
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('dyn_vor')
IF( ln_timing ) CALL timing_start('dyn_vor_3D')
!
IF( l_trddyn ) THEN !== trend diagnostics case : split the added trend in two parts ==!
!
......@@ -208,9 +213,85 @@ CONTAINS
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' vor - Ua: ', mask1=umask, &
& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
IF( ln_timing ) CALL timing_stop('dyn_vor')
IF( ln_timing ) CALL timing_stop('dyn_vor_3D')
!
END SUBROUTINE dyn_vor_3D
SUBROUTINE dyn_vor_RK3( kt, Kmm, puu, pvv, Krhs, knoco )
!!----------------------------------------------------------------------
!!
!! ** Purpose : compute the lateral ocean tracer physics.
!!
!! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the now vorticity term trend
!! - save the trends in (ztrdu,ztrdv) in 2 parts (relative
!! and planetary vorticity trends) and send them to trd_dyn
!! for futher diagnostics (l_trddyn=T)
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! ocean time-step index
INTEGER , INTENT(in ) :: Kmm, Krhs ! ocean time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk,jpt), INTENT(inout) :: puu, pvv ! ocean velocity field and RHS of momentum equation
INTEGER , INTENT(in ) :: knoco ! specified vorticity trend (= np_MET or np_RVO)
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('dyn_vor_RK3')
!
! !== total vorticity trend added to the general trend ==!
!!st WARNING 22/02 !!!!!!!! stoke drift or not stoke drift ? => bar to do later !!!
!! stoke drift a garder pas vortex force a priori !!
!! ATTENTION déja appelé avec Kbb !!
!
SELECT CASE ( nvor_scheme ) !== vorticity trend added to the general trend ==!
CASE( np_ENT ) !* energy conserving scheme (T-pts)
CALL vor_enT( kt, Kmm, knoco, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN
CALL vor_enT( 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_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)
CALL vor_eeT( kt, Kmm, knoco, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
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
CALL vor_eeT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force
ENDIF
CASE( np_ENE ) !* energy conserving scheme
CALL vor_ene( kt, Kmm, knoco, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN
CALL vor_ene( 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_ene( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force
ENDIF
CASE( np_ENS ) !* enstrophy conserving scheme
CALL vor_ens( kt, Kmm, knoco, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN
CALL vor_ens( 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_ens( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force
ENDIF
CASE( np_MIX ) !* mixed ene-ens scheme
CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! relative vorticity or metric trend (ens)
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
CALL vor_een( kt, Kmm, knoco, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend
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
CALL vor_een( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force
ENDIF
END SELECT
!
! ! print sum trends (used for debugging)
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=puu(:,:,:,Krhs), clinfo1=' vor - Ua: ', mask1=umask, &
& tab3d_2=pvv(:,:,:,Krhs), clinfo2= ' Va: ', mask2=vmask, clinfo3='dyn' )
!
IF( ln_timing ) CALL timing_stop('dyn_vor_RK3')
!
END SUBROUTINE dyn_vor
END SUBROUTINE dyn_vor_RK3
SUBROUTINE vor_enT( kt, Kmm, kvor, pu, pv, pu_rhs, pv_rhs )
......
src/OCE/DYN/dynzad.F90
View file @ d1c6ae27
......@@ -3,8 +3,9 @@ MODULE dynzad
!! *** MODULE dynzad ***
!! Ocean dynamics : vertical advection trend
!!======================================================================
!! History : OPA ! 1991-01 (G. Madec) Original code
!! NEMO 0.5 ! 2002-07 (G. Madec) Free form, F90
!! History : OPA ! 1991-01 (G. Madec) Original code
!! NEMO 0.5 ! 2002-07 (G. Madec) Free form, F90
!! 4.5 ! 2021-01 (S. Techene, G. Madec) memory optimization
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
......@@ -32,7 +33,7 @@ MODULE dynzad
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: dynzad.F90 14834 2021-05-11 09:24:44Z hadcv $
!! $Id: dynzad.F90 14428 2021-02-10 18:12:36Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -53,15 +54,14 @@ CONTAINS
!! ** Action : - Update (puu(:,:,:,Krhs),pvv(:,:,:,Krhs)) with the vert. momentum adv. trends
!! - Send the trends to trddyn for diagnostics (l_trddyn=T)
!!----------------------------------------------------------------------
INTEGER , INTENT( in ) :: kt ! ocean time-step inedx
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 , INTENT(in ) :: kt, Kmm, Krhs ! ocean time-step & 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) :: zua, zva ! local scalars
REAL(wp), DIMENSION(A2D(nn_hls)) :: zww
REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zwuw, zwvw
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdu, ztrdv
REAL(wp) :: zWf, zWfi, zzWfu, zzWdzU ! local scalars
REAL(wp) :: zWfj, zzWfv, zzWdzV ! - -
REAL(wp), DIMENSION(A2D(0)) :: zWdzU, zWdzV ! 2D inner workspace
REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: ztrdu, ztrdv ! 3D workspace
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('dyn_zad')
......@@ -72,44 +72,53 @@ CONTAINS
IF(lwp) WRITE(numout,*) 'dyn_zad : 2nd order vertical advection scheme'
ENDIF
ENDIF
!
IF( l_trddyn ) THEN ! Save puu(:,:,:,Krhs) and pvv(:,:,:,Krhs) trends
ALLOCATE( ztrdu(jpi,jpj,jpk) , ztrdv(jpi,jpj,jpk) )
ztrdu(:,:,:) = puu(:,:,:,Krhs)
ztrdv(:,:,:) = pvv(:,:,:,Krhs)
ENDIF
DO jk = 2, jpkm1 ! Vertical momentum advection at level w and u- and v- vertical
IF( ln_vortex_force ) THEN ! vertical fluxes
DO_2D( 0, 1, 0, 1 )
zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * ( ww(ji,jj,jk) + wsd(ji,jj,jk) )
END_2D
ELSE
DO_2D( 0, 1, 0, 1 )
zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk)
END_2D
!
! !== vertical momentum advection ==! at u- and v-points
!
zWdzU(A2D(0)) = 0._wp ! set surface (jk=1) vertical advection to zero
zWdzV(A2D(0)) = 0._wp
!
DO_3D( 0, 0, 0, 0 , 1, jpk-2 ) != surface to jpk-2 vertical advection
! ! vertical transport at jk+1 uw/vw-level (x2): 2*mi/j[e1e2t*(We)]
IF( ln_vortex_force ) THEN ! We = ww+wsd
zWf = e1e2t(ji ,jj ) * ( ww(ji ,jj ,jk+1) + wsd(ji ,jj ,jk+1) )
zWfi = e1e2t(ji+1,jj ) * ( ww(ji+1,jj ,jk+1) + wsd(ji+1,jj ,jk+1) )
zWfj = e1e2t(ji ,jj+1) * ( ww(ji ,jj+1,jk+1) + wsd(ji ,jj+1,jk+1) )
ELSE ! We = ww
zWf = e1e2t(ji ,jj ) * ww(ji ,jj ,jk+1)
zWfi = e1e2t(ji+1,jj ) * ww(ji+1,jj ,jk+1)
zWfj = e1e2t(ji ,jj+1) * ww(ji ,jj+1,jk+1)
ENDIF
DO_2D( 0, 0, 0, 0 ) ! vertical momentum advection at w-point
zwuw(ji,jj,jk) = ( zww(ji+1,jj ) + zww(ji,jj) ) * ( puu(ji,jj,jk-1,Kmm) - puu(ji,jj,jk,Kmm) )
zwvw(ji,jj,jk) = ( zww(ji ,jj+1) + zww(ji,jj) ) * ( pvv(ji,jj,jk-1,Kmm) - pvv(ji,jj,jk,Kmm) )
END_2D
END DO
zzWfu = zWfi + zWf ! averaging at uw- and vw-points (x2)
zzWfv = zWfj + zWf
! ! vertical advection at jk+1 uw-level (x4): zzWfu/v*dk+1[u/v]
zzWdzU = zzWfu * ( puu(ji,jj,jk,Kmm) - puu(ji,jj,jk+1,Kmm) )
zzWdzV = zzWfv * ( pvv(ji,jj,jk,Kmm) - pvv(ji,jj,jk+1,Kmm) )
!
! ! vertical advection at jk u/v-level
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - 0.25_wp * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) &
& * ( zWdzU(ji,jj) + zzWdzU )
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - 0.25_wp * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) &
& * ( zWdzV(ji,jj) + zzWdzV )
!
zWdzU(ji,jj) = zzWdzU ! save for next level computation
zWdzV(ji,jj) = zzWdzV
END_3D
!
! Surface and bottom advective fluxes set to zero
DO_2D( 0, 0, 0, 0 )
zwuw(ji,jj, 1 ) = 0._wp
zwvw(ji,jj, 1 ) = 0._wp
zwuw(ji,jj,jpk) = 0._wp
zwvw(ji,jj,jpk) = 0._wp
jk = jpkm1
DO_2D( 0, 0, 0, 0 ) != bottom vertical advection at jpkm1
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - 0.25_wp * r1_e1e2u(ji,jj) / e3u(ji,jj,jk,Kmm) &
& * zWdzU(ji,jj)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - 0.25_wp * r1_e1e2v(ji,jj) / e3v(ji,jj,jk,Kmm) &
& * zWdzV(ji,jj)
END_2D
!
DO_3D( 0, 0, 0, 0, 1, jpkm1 ) ! Vertical momentum advection at u- and v-points
puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) - ( zwuw(ji,jj,jk) + zwuw(ji,jj,jk+1) ) * r1_e1e2u(ji,jj) &
& / e3u(ji,jj,jk,Kmm)
pvv(ji,jj,jk,Krhs) = pvv(ji,jj,jk,Krhs) - ( zwvw(ji,jj,jk) + zwvw(ji,jj,jk+1) ) * r1_e1e2v(ji,jj) &
& / e3v(ji,jj,jk,Kmm)
END_3D
IF( l_trddyn ) THEN ! save the vertical advection trends for diagnostic
ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:)
ztrdv(:,:,:) = pvv(:,:,:,Krhs) - ztrdv(:,:,:)
......
src/OCE/DYN/dynzdf.F90
View file @ d1c6ae27
......@@ -37,14 +37,12 @@ MODULE dynzdf
PUBLIC dyn_zdf ! routine called by step.F90
REAL(wp) :: r_vvl ! non-linear free surface indicator: =0 if ln_linssh=T, =1 otherwise
!! * Substitutions
# include "do_loop_substitute.h90"
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: dynzdf.F90 14834 2021-05-11 09:24:44Z hadcv $
!! $Id: dynzdf.F90 14547 2021-02-25 17:07:15Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -73,14 +71,14 @@ CONTAINS
INTEGER , INTENT( in ) :: Kbb, Kmm, Krhs, Kaa ! 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
INTEGER :: iku, ikv ! local integers
REAL(wp) :: zzwi, ze3ua, zdt ! local scalars
REAL(wp) :: zzws, ze3va ! - -
REAL(wp) :: z1_e3ua, z1_e3va ! - -
REAL(wp) :: zWu , zWv ! - -
REAL(wp) :: zWui, zWvi ! - -
REAL(wp) :: zWus, zWvs ! - -
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: iku, ikv ! local integers
REAL(wp) :: zzwi, ze3ua, zDt_2 ! local scalars
REAL(wp) :: zzws, ze3va ! - -
REAL(wp) :: z1_e3ua, z1_e3va ! - -
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(:,:,:), ALLOCATABLE :: ztrdu, ztrdv ! - -
!!---------------------------------------------------------------------
......@@ -92,12 +90,11 @@ CONTAINS
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'dyn_zdf_imp : vertical momentum diffusion implicit operator'
IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
!
If( ln_linssh ) THEN ; r_vvl = 0._wp ! non-linear free surface indicator
ELSE ; r_vvl = 1._wp
ENDIF
ENDIF
ENDIF
!
zDt_2 = rDt * 0.5_wp
!
! !* explicit top/bottom drag case
IF( .NOT.ln_drgimp ) CALL zdf_drg_exp( kt, Kmm, puu(:,:,:,Kbb), pvv(:,:,:,Kbb), puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add top/bottom friction trend to (puu(Kaa),pvv(Kaa))
!
......@@ -139,23 +136,19 @@ CONTAINS
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)
ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) &
& + r_vvl * e3u(ji,jj,iku,Kaa)
ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) &
& + r_vvl * e3v(ji,jj,ikv,Kaa)
puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + rDt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua
pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va
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)
ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) &
& + r_vvl * e3u(ji,jj,iku,Kaa)
ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) &
& + r_vvl * e3v(ji,jj,ikv,Kaa)
puu(ji,jj,iku,Kaa) = puu(ji,jj,iku,Kaa) + rDt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) * uu_b(ji,jj,Kaa) / ze3ua
pvv(ji,jj,ikv,Kaa) = pvv(ji,jj,ikv,Kaa) + rDt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) * vv_b(ji,jj,Kaa) / ze3va
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_2D
END IF
ENDIF
......@@ -163,70 +156,61 @@ CONTAINS
! !== Vertical diffusion on u ==!
!
! !* Matrix construction
zdt = rDt * 0.5
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 )
ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) &
& + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point
zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) &
& / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk )
zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) &
& / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) / ze3ua
zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua
zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp )
zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
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 )
ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) & ! after scale factor at U-point
& + r_vvl * e3u(ji,jj,jk,Kaa)
zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) &
& / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk )
zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) &
& / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+1)
zWui = ( wi(ji,jj,jk ) + wi(ji+1,jj,jk ) ) / ze3ua
zWus = ( wi(ji,jj,jk+1) + wi(ji+1,jj,jk+1) ) / ze3ua
zwi(ji,jj,jk) = zzwi + zdt * MIN( zWui, 0._wp )
zws(ji,jj,jk) = zzws - zdt * MAX( zWus, 0._wp )
zwd(ji,jj,jk) = 1._wp - zzwi - zzws + zdt * ( MAX( zWui, 0._wp ) - MIN( zWus, 0._wp ) )
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
ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) &
& + r_vvl * e3u(ji,jj,1,Kaa)
zzws = - zdt * ( avm(ji+1,jj,2) + avm(ji ,jj,2) ) &
& / ( ze3ua * e3uw(ji,jj,2,Kmm) ) * wumask(ji,jj,2)
zWus = ( wi(ji ,jj,2) + wi(ji+1,jj,2) ) / ze3ua
zws(ji,jj,1 ) = zzws - zdt * MAX( zWus, 0._wp )
zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWus, 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 )
ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) &
& + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point
zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) + akzu(ji,jj,jk ) ) &
& / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk )
zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) + akzu(ji,jj,jk+1) ) &
& / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+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,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 )
ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,jk,Kmm) &
& + r_vvl * e3u(ji,jj,jk,Kaa) ! after scale factor at U-point
zzwi = - zdt * ( avm(ji+1,jj,jk ) + avm(ji,jj,jk ) ) &
& / ( ze3ua * e3uw(ji,jj,jk ,Kmm) ) * wumask(ji,jj,jk )
zzws = - zdt * ( avm(ji+1,jj,jk+1) + avm(ji,jj,jk+1) ) &
& / ( ze3ua * e3uw(ji,jj,jk+1,Kmm) ) * wumask(ji,jj,jk+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,jj,jk) = zzwi
zws(ji,jj,jk) = zzws
zwd(ji,jj,jk) = 1._wp - zzwi - zzws
......@@ -248,17 +232,15 @@ CONTAINS
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
ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) &
& + r_vvl * e3u(ji,jj,iku,Kaa) ! after scale factor at T-point
zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_bot(ji+1,jj)+rCdU_bot(ji,jj) ) / ze3ua
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
ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,iku,Kmm) &
& + r_vvl * e3u(ji,jj,iku,Kaa) ! after scale factor at T-point
zwd(ji,jj,iku) = zwd(ji,jj,iku) - rDt * 0.5*( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) / ze3ua
zwd(ji,jj,iku) = zwd(ji,jj,iku) - zDt_2 *( rCdU_top(ji+1,jj)+rCdU_top(ji,jj) ) &
& / e3u(ji,jj,iku,Kaa)
END_2D
END IF
ENDIF
......@@ -283,10 +265,14 @@ CONTAINS
END_3D
!
DO_2D( 0, 0, 0, 0 ) !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==!
ze3ua = ( 1._wp - r_vvl ) * e3u(ji,jj,1,Kmm) &
& + r_vvl * e3u(ji,jj,1,Kaa)
puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * 0.5_wp * ( utau_b(ji,jj) + utau(ji,jj) ) &
& / ( ze3ua * rho0 ) * umask(ji,jj,1)
#if defined key_RK3
! ! RK3: use only utau (not utau_b)
puu(ji,jj,1,Kaa) = puu(ji,jj,1,Kaa) + rDt * utau(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) + utau(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)
......@@ -302,70 +288,61 @@ CONTAINS
! !== Vertical diffusion on v ==!
!
! !* Matrix construction
zdt = rDt * 0.5
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 )
ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) &
& + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point
zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) &
& / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk )
zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) &
& / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) / ze3va
zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va
zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp )
zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp )
zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
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 )
ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) &
& + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point
zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) &
& / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk )
zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) &
& / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+1)
zWvi = ( wi(ji,jj,jk ) + wi(ji,jj+1,jk ) ) / ze3va
zWvs = ( wi(ji,jj,jk+1) + wi(ji,jj+1,jk+1) ) / ze3va
zwi(ji,jj,jk) = zzwi + zdt * MIN( zWvi, 0._wp )
zws(ji,jj,jk) = zzws - zdt * MAX( zWvs, 0._wp )
zwd(ji,jj,jk) = 1._wp - zzwi - zzws - zdt * ( - MAX( zWvi, 0._wp ) + MIN( zWvs, 0._wp ) )
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
ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) &
& + r_vvl * e3v(ji,jj,1,Kaa)
zzws = - zdt * ( avm(ji,jj+1,2) + avm(ji,jj,2) ) &
& / ( ze3va * e3vw(ji,jj,2,Kmm) ) * wvmask(ji,jj,2)
zWvs = ( wi(ji,jj ,2) + wi(ji,jj+1,2) ) / ze3va
zws(ji,jj,1 ) = zzws - zdt * MAX( zWvs, 0._wp )
zwd(ji,jj,1 ) = 1._wp - zzws - zdt * ( MIN( zWvs, 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 )
ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) &
& + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point
zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) + akzv(ji,jj,jk ) ) &
& / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk )
zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) + akzv(ji,jj,jk+1) ) &
& / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+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,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 )
ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,jk,Kmm) &
& + r_vvl * e3v(ji,jj,jk,Kaa) ! after scale factor at V-point
zzwi = - zdt * ( avm(ji,jj+1,jk ) + avm(ji,jj,jk ) ) &
& / ( ze3va * e3vw(ji,jj,jk ,Kmm) ) * wvmask(ji,jj,jk )
zzws = - zdt * ( avm(ji,jj+1,jk+1) + avm(ji,jj,jk+1) ) &
& / ( ze3va * e3vw(ji,jj,jk+1,Kmm) ) * wvmask(ji,jj,jk+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,jj,jk) = zzwi
zws(ji,jj,jk) = zzws
zwd(ji,jj,jk) = 1._wp - zzwi - zzws
......@@ -386,16 +363,14 @@ CONTAINS
IF( ln_drgimp ) THEN
DO_2D( 0, 0, 0, 0 )
ikv = mbkv(ji,jj) ! (deepest ocean u- and v-points)
ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) &
& + r_vvl * e3v(ji,jj,ikv,Kaa) ! after scale factor at T-point
zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) / ze3va
zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - zDt_2*( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj) ) &
& / e3v(ji,jj,ikv,Kaa)
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)
ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,ikv,Kmm) &
& + r_vvl * e3v(ji,jj,ikv,Kaa) ! after scale factor at T-point
zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - rDt * 0.5*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) / ze3va
zwd(ji,jj,ikv) = zwd(ji,jj,ikv) - zDt_2*( rCdU_top(ji,jj+1)+rCdU_top(ji,jj) ) &
& / e3v(ji,jj,ikv,Kaa)
END_2D
ENDIF
ENDIF
......@@ -420,10 +395,14 @@ CONTAINS
END_3D
!
DO_2D( 0, 0, 0, 0 ) !== second recurrence: SOLk = RHSk - Lk / Dk-1 Lk-1 ==!
ze3va = ( 1._wp - r_vvl ) * e3v(ji,jj,1,Kmm) &
& + r_vvl * e3v(ji,jj,1,Kaa)
pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * 0.5_wp * ( vtau_b(ji,jj) + vtau(ji,jj) ) &
& / ( ze3va * rho0 ) * vmask(ji,jj,1)
#if defined key_RK3
! ! RK3: use only vtau (not vtau_b)
pvv(ji,jj,1,Kaa) = pvv(ji,jj,1,Kaa) + rDt * vtau(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) + vtau(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)
......@@ -437,8 +416,8 @@ CONTAINS
END_3D
!
IF( l_trddyn ) THEN ! save the vertical diffusive trends for further diagnostics
ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) ) / rDt - ztrdu(:,:,:)
ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) ) / rDt - ztrdv(:,:,:)
ztrdu(:,:,:) = ( puu(:,:,:,Kaa) - puu(:,:,:,Kbb) )*r1_Dt - ztrdu(:,:,:)
ztrdv(:,:,:) = ( pvv(:,:,:,Kaa) - pvv(:,:,:,Kbb) )*r1_Dt - ztrdv(:,:,:)
CALL trd_dyn( ztrdu, ztrdv, jpdyn_zdf, kt, Kmm )
DEALLOCATE( ztrdu, ztrdv )
ENDIF
......
src/OCE/DYN/sshwzv.F90
View file @ d1c6ae27
......@@ -16,7 +16,9 @@ MODULE sshwzv
!!----------------------------------------------------------------------
!! ssh_nxt : after ssh
!! ssh_atf : time filter the ssh arrays
!! wzv : compute now vertical velocity
!! wzv : generic interface of vertical velocity calculation
!! wzv_MLF : MLF: compute NOW vertical velocity
!! wzv_RK3 : RK3: Compute a vertical velocity
!!----------------------------------------------------------------------
USE oce ! ocean dynamics and tracers variables
USE isf_oce ! ice shelf
......@@ -43,6 +45,10 @@ MODULE sshwzv
IMPLICIT NONE
PRIVATE
! !! * Interface
INTERFACE wzv
MODULE PROCEDURE wzv_MLF, wzv_RK3
END INTERFACE
PUBLIC ssh_nxt ! called by step.F90
PUBLIC wzv ! called by step.F90
......@@ -54,7 +60,7 @@ MODULE sshwzv
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
!! $Id: sshwzv.F90 15150 2021-07-27 10:38:24Z smasson $
!! $Id: sshwzv.F90 14618 2021-03-19 14:42:32Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
......@@ -76,10 +82,11 @@ CONTAINS
INTEGER , INTENT(in ) :: kt ! time step
INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! time level index
REAL(wp), DIMENSION(jpi,jpj,jpt), INTENT(inout) :: pssh ! sea-surface height
!
!
INTEGER :: ji, jj, jk ! dummy loop index
REAL(wp) :: zcoef ! local scalar
REAL(wp), DIMENSION(jpi,jpj) :: zhdiv ! 2D workspace
REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d ! 3D workspace
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('ssh_nxt')
......@@ -91,12 +98,11 @@ CONTAINS
ENDIF
!
zcoef = 0.5_wp * r1_rho0
! !------------------------------!
! ! After Sea Surface Height !
! !------------------------------!
IF(ln_wd_il) THEN
CALL wad_lmt(pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), rDt, Kmm, uu, vv )
CALL wad_lmt( pssh(:,:,Kbb), zcoef * (emp_b(:,:) + emp(:,:)), rDt, Kmm, uu, vv )
ENDIF
CALL div_hor( kt, Kbb, Kmm ) ! Horizontal divergence
......@@ -110,7 +116,11 @@ CONTAINS
! compute the vertical velocity which can be used to compute the non-linear terms of the momentum equations.
!
DO_2D_OVR( 1, nn_hls, 1, nn_hls ) ! Loop bounds limited by hdiv definition in div_hor
pssh(ji,jj,Kaa) = ( pssh(ji,jj,Kbb) - rDt * ( zcoef * ( emp_b(ji,jj) + emp(ji,jj) ) + zhdiv(ji,jj) ) ) * ssmask(ji,jj)
#if defined key_RK3
pssh(ji,jj,Kaa) = ( pssh(ji,jj,Kbb) - rDt * ( r1_rho0 * emp(ji,jj) + zhdiv(ji,jj) ) ) * ssmask(ji,jj)
#else
pssh(ji,jj,Kaa) = ( pssh(ji,jj,Kbb) - rDt * ( zcoef * ( emp_b(ji,jj) + emp(ji,jj) ) + zhdiv(ji,jj) ) ) * ssmask(ji,jj)
#endif
END_2D
! pssh must be defined everywhere (true for dyn_spg_ts, not for dyn_spg_exp)
IF ( .NOT. ln_dynspg_ts .AND. nn_hls == 2 ) CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp )
......@@ -123,23 +133,21 @@ CONTAINS
IF ( .NOT.ln_dynspg_ts ) THEN
IF( ln_bdy ) THEN
IF (nn_hls==1) CALL lbc_lnk( 'sshwzv', pssh(:,:,Kaa), 'T', 1.0_wp ) ! Not sure that's necessary
CALL bdy_ssh( pssh(:,:,Kaa) ) ! Duplicate sea level across open boundaries
CALL bdy_ssh( pssh(:,:,Kaa) ) ! Duplicate sea level across open boundaries
ENDIF
ENDIF
! !------------------------------!
! ! outputs !
! !------------------------------!
!
! ! Control print
IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab2d_1=pssh(:,:,Kaa), clinfo1=' pssh(:,:,Kaa) - : ', mask1=tmask )
!
IF( ln_timing ) CALL timing_stop('ssh_nxt')
!
END SUBROUTINE ssh_nxt
SUBROUTINE wzv( kt, Kbb, Kmm, Kaa, pww )
SUBROUTINE wzv_MLF( kt, Kbb, Kmm, Kaa, pww )
!!----------------------------------------------------------------------
!! *** ROUTINE wzv ***
!! *** ROUTINE wzv_MLF ***
!!
!! ** Purpose : compute the now vertical velocity
!!
......@@ -160,12 +168,12 @@ CONTAINS
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zhdiv
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('wzv')
IF( ln_timing ) CALL timing_start('wzv_MLF')
!
IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'wzv : now vertical velocity '
IF(lwp) WRITE(numout,*) '~~~~~ '
IF(lwp) WRITE(numout,*) 'wzv_MLF : now vertical velocity '
IF(lwp) WRITE(numout,*) '~~~~~~~'
!
pww(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)
ENDIF
......@@ -198,7 +206,7 @@ CONTAINS
& - e3t(ji,jj,jk,Kbb) ) ) * tmask(ji,jj,jk)
END_3D
! IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
DEALLOCATE( zhdiv )
DEALLOCATE( zhdiv )
! !=================================!
ELSEIF( ln_linssh ) THEN !== linear free surface cases ==!
! !=================================!
......@@ -209,9 +217,146 @@ CONTAINS
ELSE !== Quasi-Eulerian vertical coordinate ==! ('key_qco')
! !==========================================!
DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 ) ! integrate from the bottom the hor. divergence
#if defined key_qco
!!gm slightly faster
pww(ji,jj,jk) = pww(ji,jj,jk+1) - ( e3t(ji,jj,jk,Kmm) * hdiv(ji,jj,jk) &
& + r1_Dt * e3t_0(ji,jj,jk) * ( r3t(ji,jj,Kaa) - r3t(ji,jj,Kbb) ) ) * tmask(ji,jj,jk)
#else
pww(ji,jj,jk) = pww(ji,jj,jk+1) - ( e3t(ji,jj,jk,Kmm) * hdiv(ji,jj,jk) &
& + r1_Dt * ( e3t(ji,jj,jk,Kaa) &
& - e3t(ji,jj,jk,Kbb) ) ) * tmask(ji,jj,jk)
#endif
END_3D
ENDIF
IF( ln_bdy ) THEN
DO jk = 1, jpkm1
DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
pww(ji,jj,jk) = pww(ji,jj,jk) * bdytmask(ji,jj)
END_2D
END DO
ENDIF
!
#if defined key_agrif
IF( .NOT. AGRIF_Root() ) THEN
!
! Mask vertical velocity at first/last columns/row
! inside computational domain (cosmetic)
DO jk = 1, jpkm1
IF( lk_west ) THEN ! --- West --- !
DO ji = mi0(2+nn_hls), mi1(2+nn_hls)
DO jj = 1, jpj
pww(ji,jj,jk) = 0._wp
END DO
END DO
ENDIF
IF( lk_east ) THEN ! --- East --- !
DO ji = mi0(jpiglo-1-nn_hls), mi1(jpiglo-1-nn_hls)
DO jj = 1, jpj
pww(ji,jj,jk) = 0._wp
END DO
END DO
ENDIF
IF( lk_south ) THEN ! --- South --- !
DO jj = mj0(2+nn_hls), mj1(2+nn_hls)
DO ji = 1, jpi
pww(ji,jj,jk) = 0._wp
END DO
END DO
ENDIF
IF( lk_north ) THEN ! --- North --- !
DO jj = mj0(jpjglo-1-nn_hls), mj1(jpjglo-1-nn_hls)
DO ji = 1, jpi
pww(ji,jj,jk) = 0._wp
END DO
END DO
ENDIF
!
END DO
!
ENDIF
#endif
!
IF( ln_timing ) CALL timing_stop('wzv_MLF')
!
END SUBROUTINE wzv_MLF
SUBROUTINE wzv_RK3( kt, Kbb, Kmm, Kaa, puu, pvv, pww )
!!----------------------------------------------------------------------
!! *** ROUTINE wzv_RK3 ***
!!
!! ** Purpose : compute the now vertical velocity
!!
!! ** Method : - Using the incompressibility hypothesis, the vertical
!! velocity is computed by integrating the horizontal divergence
!! from the bottom to the surface minus the scale factor evolution.
!! The boundary conditions are w=0 at the bottom (no flux) and.
!!
!! ** action : pww : now vertical velocity
!!
!! Reference : Leclair, M., and G. Madec, 2009, Ocean Modelling.
!!----------------------------------------------------------------------
INTEGER , INTENT(in ) :: kt ! time step
INTEGER , INTENT(in ) :: Kbb, Kmm, Kaa ! time level indices
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(in ) :: puu, pvv ! horizontal velocity at Kmm
REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pww ! vertical velocity at Kmm
!
INTEGER :: ji, jj, jk ! dummy loop indices
REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zhdiv
REAL(wp) , DIMENSION(jpi,jpj,jpk) :: ze3div
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('wzv_RK3')
!
IF( kt == nit000 ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'wzv_RK3 : now vertical velocity '
IF(lwp) WRITE(numout,*) '~~~~~ '
!
pww(:,:,jpk) = 0._wp ! bottom boundary condition: w=0 (set once for all)
ENDIF
!
CALL div_hor( kt, Kbb, Kmm, puu, pvv, ze3div )
! !------------------------------!
! ! Now Vertical Velocity !
! !------------------------------!
!
! !===============================!
IF( ln_vvl_ztilde .OR. ln_vvl_layer ) THEN !== z_tilde and layer cases ==!
! !===============================!
ALLOCATE( zhdiv(jpi,jpj,jpk) )
!
DO jk = 1, jpkm1
! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
END_2D
END DO
IF( nn_hls == 1) CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp) ! - ML - Perhaps not necessary: not used for horizontal "connexions"
! ! Is it problematic to have a wrong vertical velocity in boundary cells?
! ! Same question holds for hdiv. Perhaps just for security
DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 ) ! integrate from the bottom the hor. divergence
pww(ji,jj,jk) = pww(ji,jj,jk+1) - ( ze3div(ji,jj,jk) + zhdiv(ji,jj,jk) &
& + r1_Dt * ( e3t(ji,jj,jk,Kaa) &
& - e3t(ji,jj,jk,Kbb) ) ) * tmask(ji,jj,jk)
END_3D
!
DEALLOCATE( zhdiv )
! !=================================!
ELSEIF( ln_linssh ) THEN !== linear free surface cases ==!
! !=================================!
DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 ) ! integrate from the bottom the hor. divergence
pww(ji,jj,jk) = pww(ji,jj,jk+1) - ze3div(ji,jj,jk)
END_3D
! !==========================================!
ELSE !== Quasi-Eulerian vertical coordinate ==! ('key_qco')
! !==========================================!
DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 ) ! integrate from the bottom the hor. divergence
! ! NB: [e3t[a] -e3t[b] ]=e3t_0*[r3t[a]-r3t[b]]
pww(ji,jj,jk) = pww(ji,jj,jk+1) - ( ze3div(ji,jj,jk) &
& + r1_Dt * e3t_0(ji,jj,jk) * ( r3t(ji,jj,Kaa) - r3t(ji,jj,Kbb) ) ) * tmask(ji,jj,jk)
END_3D
ENDIF
......@@ -263,9 +408,9 @@ CONTAINS
ENDIF
#endif
!
IF( ln_timing ) CALL timing_stop('wzv')
IF( ln_timing ) CALL timing_stop('wzv_RK3')
!
END SUBROUTINE wzv
END SUBROUTINE wzv_RK3
SUBROUTINE ssh_atf( kt, Kbb, Kmm, Kaa, pssh )
......