From f1728754d2173725b7cc6ce0c8b2efcde516b7e5 Mon Sep 17 00:00:00 2001
From: Sibylle TECHENE <techenes@irene190.c-irene.tgcc.ccc.cea.fr>
Date: Wed, 5 Jan 2022 23:32:18 +0100
Subject: [PATCH] merged debug_rk3
---
src/OCE/DOM/domqco.F90 | 91 ++++++++++++++++------
src/OCE/DOM/istate.F90 | 33 ++++----
src/OCE/DYN/divhor.F90 | 10 +--
src/OCE/DYN/dynspg_ts.F90 | 2 +-
src/OCE/DYN/sshwzv.F90 | 37 ++++-----
src/OCE/IOM/restart.F90 | 21 ++---
src/OCE/SBC/sbcssm.F90 | 7 +-
src/OCE/TRA/traatf.F90 | 4 +
src/OCE/TRA/traqsr.F90 | 11 +--
src/OCE/TRA/trasbc.F90 | 2 +-
src/OCE/ZDF/zdfphy.F90 | 4 +-
src/OCE/ZDF/zdftke.F90 | 3 +-
src/OCE/stprk3_stg.F90 | 82 ++++++++++++++------
src/TOP/trcini.F90 | 7 --
src/TOP/trcrst.F90 | 3 -
tests/ISOMIP+/EXPREF/namelist_cfg | 3 +-
tests/ISOMIP+/MY_SRC/eosbn2.F90 | 4 +-
tests/ISOMIP+/MY_SRC/isf_oce.F90 | 1 -
tests/ISOMIP+/MY_SRC/isfcavgam.F90 | 4 +-
tests/ISOMIP+/MY_SRC/isfstp.F90 | 11 ++-
tests/ISOMIP+/MY_SRC/istate.F90 | 46 +++++++----
tests/ISOMIP+/MY_SRC/sbcfwb.F90 | 118 +++++++++++++++++------------
tests/ISOMIP+/MY_SRC/tradmp.F90 | 42 ++++++----
23 files changed, 335 insertions(+), 211 deletions(-)
diff --git a/src/OCE/DOM/domqco.F90 b/src/OCE/DOM/domqco.F90
index c0e8fb93..2e981d71 100644
--- a/src/OCE/DOM/domqco.F90
+++ b/src/OCE/DOM/domqco.F90
@@ -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
@@ -117,7 +118,7 @@ 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) )
+ CALL dom_qco_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm) ) !!st needed for Agrif_Grid call in nemo_gcm
#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(:,:) )
@@ -130,6 +131,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 ***
!!
@@ -158,7 +214,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 ) &
@@ -166,52 +222,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
diff --git a/src/OCE/DOM/istate.F90 b/src/OCE/DOM/istate.F90
index 5c9b2b78..ad97be1b 100644
--- a/src/OCE/DOM/istate.F90
+++ b/src/OCE/DOM/istate.F90
@@ -141,20 +141,6 @@ CONTAINS
ENDIF
#endif
!
- ! Initialize "now" barotropic velocities:
- ! Do it whatever the free surface method, these arrays being used eventually
- !
-!!gm the use of umask & vmask is not necessary below as uu(:,:,:,Kmm), vv(:,:,:,Kmm), uu(:,:,:,Kbb), vv(:,:,:,Kbb) are always masked
-#if ! defined key_RK3
- 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)
- END_3D
- uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm)
- vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv(:,:,Kmm)
-#endif
- !
#if defined key_RK3
IF( .NOT. ln_rstart ) THEN
#endif
@@ -171,6 +157,25 @@ CONTAINS
!
#if defined key_RK3
ENDIF
+#endif
+ !
+ ! Initialize "now" barotropic velocities:
+ ! Do it whatever the free surface method, these arrays being used eventually
+ !
+#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)
+ END_3D
+ uu_b(:,:,Kmm) = uu_b(:,:,Kmm) * r1_hu(:,:,Kmm)
+ vv_b(:,:,Kmm) = vv_b(:,:,Kmm) * r1_hv(:,:,Kmm)
#endif
!
END SUBROUTINE istate_init
diff --git a/src/OCE/DYN/divhor.F90 b/src/OCE/DYN/divhor.F90
index edfccb67..2cdf53d1 100644
--- a/src/OCE/DYN/divhor.F90
+++ b/src/OCE/DYN/divhor.F90
@@ -83,7 +83,7 @@ CONTAINS
!
pe3divUh(:,:,:) = 0._wp !!gm to be applied to the halos only
!
- DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+ 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) &
@@ -100,16 +100,16 @@ CONTAINS
!
IF( ln_isf ) CALL isf_hdiv( kt, Kmm, hdiv ) !== + ice-shelf mass exchange ==!
!
- CALL lbc_lnk( 'divhor', hdiv, 'T', 1._wp ) ! (no sign change)
+ 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 jk=1, jpkm1
- pe3divUh(:,:,jk) = hdiv(:,:,jk) * e3t(:,:,jk,Kmm)
- END DO
+ 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
!
!
diff --git a/src/OCE/DYN/dynspg_ts.F90 b/src/OCE/DYN/dynspg_ts.F90
index 911cf67f..69d04f57 100644
--- a/src/OCE/DYN/dynspg_ts.F90
+++ b/src/OCE/DYN/dynspg_ts.F90
@@ -802,7 +802,7 @@ CONTAINS
END_2D
# endif
!
- CALL lbc_lnk_multi( 'dynspg_ts', puu_b, 'U', -1._wp, pvv_b, 'V', -1._wp ) ! Boundary conditions
+ CALL lbc_lnk( 'dynspg_ts', puu_b, 'U', -1._wp, pvv_b, 'V', -1._wp ) ! Boundary conditions
!
ENDIF
!
diff --git a/src/OCE/DYN/sshwzv.F90 b/src/OCE/DYN/sshwzv.F90
index 4c0b9778..26a0e789 100644
--- a/src/OCE/DYN/sshwzv.F90
+++ b/src/OCE/DYN/sshwzv.F90
@@ -328,40 +328,41 @@ CONTAINS
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( 0, 0, 0, 0 )
+ 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
- CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp) ! - ML - Perhaps not necessary: not used for horizontal "connexions"
+ 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 jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
- ! computation of w
- pww(:,:,jk) = pww(:,:,jk+1) - ( ze3div(:,:,jk) + zhdiv(:,:,jk) &
- & + r1_Dt * ( e3t(:,:,jk,Kaa) &
- & - e3t(:,:,jk,Kbb) ) ) * tmask(:,:,jk)
- END DO
- ! IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
+ 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 jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
- pww(:,:,jk) = pww(:,:,jk+1) - ze3div(:,:,jk)
- END DO
+ 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 jk = jpkm1, 1, -1 ! integrate from the bottom the hor. divergence
- ! ! NB: [e3t[a] -e3t[b] ]=e3t_0*[r3t[a]-r3t[b]]
- pww(:,:,jk) = pww(:,:,jk+1) - ( ze3div(:,:,jk) &
- & + r1_Dt * e3t_0(:,:,jk) * ( r3t(:,:,Kaa) - r3t(:,:,Kbb) ) ) * tmask(:,:,jk)
- END DO
+ 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
IF( ln_bdy ) THEN
DO jk = 1, jpkm1
- pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:)
+ 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
!
diff --git a/src/OCE/IOM/restart.F90 b/src/OCE/IOM/restart.F90
index 0aec277c..60e2710a 100644
--- a/src/OCE/IOM/restart.F90
+++ b/src/OCE/IOM/restart.F90
@@ -176,7 +176,6 @@ CONTAINS
CALL iom_rstput( kt, nitrst, numrow, 'vn' , vv(:,:,: ,Kmm) )
CALL iom_rstput( kt, nitrst, numrow, 'tn' , ts(:,:,:,jp_tem,Kmm) )
CALL iom_rstput( kt, nitrst, numrow, 'sn' , ts(:,:,:,jp_sal,Kmm) )
- IF( .NOT.lk_SWE ) CALL iom_rstput( kt, nitrst, numrow, 'rhop', rhop )
#endif
ENDIF
@@ -290,6 +289,12 @@ CONTAINS
CALL iom_get( numror, jpdom_auto, 'sb' , ts(:,:,:,jp_sal,Kbb) )
CALL iom_get( numror, jpdom_auto, 'uu_b' , uu_b(:,: ,Kbb), cd_type = 'U', psgn = -1._wp )
CALL iom_get( numror, jpdom_auto, 'vv_b' , vv_b(:,: ,Kbb), cd_type = 'V', psgn = -1._wp )
+ uu(:,:,: ,Kmm) = uu(:,:,: ,Kbb) ! Kmm values set to Kbb for initialisation (sbc_ssm_init)
+ vv(:,:,: ,Kmm) = vv(:,:,: ,Kbb)
+ ts(:,:,:,:,Kmm) = ts(:,:,:,:,Kbb)
+ !
+ uu_b(:,:,Kmm) = uu_b(:,:,Kbb) ! Kmm value set to Kbb for initialisation in Agrif_Regrid
+ vv_b(:,:,Kmm) = vv_b(:,:,Kbb)
#else
! !* Read Kmm fields (MLF only)
IF(lwp) WRITE(numout,*) ' Kmm u, v and T-S fields read in the restart file'
@@ -313,18 +318,6 @@ CONTAINS
ENDIF
#endif
!
- IF( .NOT.lk_SWE ) THEN
- IF( iom_varid( numror, 'rhop', ldstop = .FALSE. ) > 0 ) THEN
- CALL iom_get( numror, jpdom_auto, 'rhop' , rhop ) ! now potential density
- ELSE
-#if defined key_RK3
- CALL eos( ts, Kbb, rhop )
-#else
- CALL eos( ts, Kmm, rhop )
-#endif
- ENDIF
- ENDIF
- !
END SUBROUTINE rst_read
@@ -367,7 +360,7 @@ CONTAINS
CALL iom_get( numror, jpdom_auto, 'sshb' , ssh(:,:,Kbb) )
!
! !* RK3: Set ssh at Kmm for AGRIF
- ssh(:,:,Kmm) = 0._wp
+ ssh(:,:,Kmm) = ssh(:,:,Kbb)
#else
! !* MLF: Read ssh at Kmm
IF(lwp) WRITE(numout,*)
diff --git a/src/OCE/SBC/sbcssm.F90 b/src/OCE/SBC/sbcssm.F90
index 1c108df0..5f3daebf 100644
--- a/src/OCE/SBC/sbcssm.F90
+++ b/src/OCE/SBC/sbcssm.F90
@@ -64,7 +64,12 @@ CONTAINS
zts(:,:,jp_tem) = ts(:,:,1,jp_tem,Kmm)
zts(:,:,jp_sal) = ts(:,:,1,jp_sal,Kmm)
!
- ! ! ---------------------------------------- !
+ ! !===>>> CAUTION: lbc_lnk is required on fraqsr_lev since sea ice computes on the full domain
+ ! ! otherwise restartability and reproducibility are broken
+ ! ! computed in traqsr only on the inner domain
+ CALL lbc_lnk( 'sbc_ssm', fraqsr_1lev(:,:), 'T', 1._wp )
+ !
+ ! ! ---------------------------------------- !
IF( nn_fsbc == 1 ) THEN ! Instantaneous surface fields !
! ! ---------------------------------------- !
ssu_m(:,:) = uu(:,:,1,Kbb)
diff --git a/src/OCE/TRA/traatf.F90 b/src/OCE/TRA/traatf.F90
index 855e2cbc..57a9ca0a 100644
--- a/src/OCE/TRA/traatf.F90
+++ b/src/OCE/TRA/traatf.F90
@@ -268,6 +268,10 @@ CONTAINS
ALLOCATE( ztrd_atf(jpi,jpj,jpk,kjpt) )
ztrd_atf(:,:,:,:) = 0.0_wp
ENDIF
+ !
+!!st variables only computed in the interior by traqsr
+ IF( ll_traqsr ) CALL lbc_lnk( 'traatf', qsr_hc_b(:,:,:) , 'T', 1.0_wp, qsr_hc(:,:,:) , 'T', 1.0_wp )
+ !
zfact = 1._wp / p2dt
zfact1 = rn_atfp * p2dt
zfact2 = zfact1 * r1_rho0
diff --git a/src/OCE/TRA/traqsr.F90 b/src/OCE/TRA/traqsr.F90
index 7e76160a..717fd6b1 100644
--- a/src/OCE/TRA/traqsr.F90
+++ b/src/OCE/TRA/traqsr.F90
@@ -145,13 +145,13 @@ CONTAINS
ENDIF
ELSE ! No restart or Euler forward at 1st time step
z1_2 = 1._wp
- DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
+ DO_3D_OVR( 0, 0, 0, 0, 1, jpk )
qsr_hc_b(ji,jj,jk) = 0._wp
END_3D
ENDIF
ELSE !== Swap of qsr heat content ==!
z1_2 = 0.5_wp
- DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
+ DO_3D_OVR( 0, 0, 0, 0, 1, jpk )
qsr_hc_b(ji,jj,jk) = qsr_hc(ji,jj,jk)
END_3D
ENDIF
@@ -207,17 +207,12 @@ CONTAINS
& / e3t(ji,jj,jk,Kmm)
END_3D
!
-!!st7-2
! sea-ice: store the 1st ocean level attenuation coefficient
- DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
+ DO_2D( 0, 0, 0, 0 )
IF( qsr(ji,jj) /= 0._wp ) THEN ; fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) )
ELSE ; fraqsr_1lev(ji,jj) = 1._wp
ENDIF
END_2D
- ! !===>>> CAUTION: lbc_lnk is required on fraqsr_lev since sea ice computes on the full domain
- ! ! otherwise restartability and reproducibility are broken
- CALL lbc_lnk( 'tra_qsr', fraqsr_1lev(:,:), 'T', 1._wp )
-!!st CALL lbc_lnk( 'tra_qsr', qsr_hc(:,:,:), 'T', 1._wp )
!
IF( iom_use('qsr3d') ) THEN ! output the shortwave Radiation distribution
ALLOCATE( zetot(A2D(nn_hls),jpk) )
diff --git a/src/OCE/TRA/trasbc.F90 b/src/OCE/TRA/trasbc.F90
index b23c6062..0f8472a3 100644
--- a/src/OCE/TRA/trasbc.F90
+++ b/src/OCE/TRA/trasbc.F90
@@ -275,7 +275,7 @@ CONTAINS
!!gm This should be moved into sbcmod.F90 module ? (especially now that ln_traqsr is read in namsbc namelist)
IF( .NOT.ln_traqsr .AND. kstg == 1) THEN ! no solar radiation penetration
- DO_2D( 0, 0, 0, 0 )
+ DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
qns(ji,jj) = qns(ji,jj) + qsr(ji,jj) ! total heat flux in qns
qsr(ji,jj) = 0._wp ! qsr set to zero
END_2D
diff --git a/src/OCE/ZDF/zdfphy.F90 b/src/OCE/ZDF/zdfphy.F90
index 66beff36..5fbf2322 100644
--- a/src/OCE/ZDF/zdfphy.F90
+++ b/src/OCE/ZDF/zdfphy.F90
@@ -209,7 +209,7 @@ CONTAINS
ioptio = 0
IF( ln_zdfcst ) THEN ; ioptio = ioptio + 1 ; nzdf_phy = np_CST ; ENDIF
IF( ln_zdfric ) THEN ; ioptio = ioptio + 1 ; nzdf_phy = np_RIC ; CALL zdf_ric_init ; ENDIF
- IF( ln_zdftke ) THEN ; ioptio = ioptio + 1 ; nzdf_phy = np_TKE ; CALL zdf_tke_init( Kmm ) ; ENDIF
+ IF( ln_zdftke ) THEN ; ioptio = ioptio + 1 ; nzdf_phy = np_TKE ; CALL zdf_tke_init ; ENDIF
IF( ln_zdfgls ) THEN ; ioptio = ioptio + 1 ; nzdf_phy = np_GLS ; CALL zdf_gls_init ; ENDIF
IF( ln_zdfosm ) THEN ; ioptio = ioptio + 1 ; nzdf_phy = np_OSM ; CALL zdf_osm_init( Kmm ) ; ENDIF
!
@@ -350,7 +350,7 @@ CONTAINS
IF( ln_zdfiwm ) CALL zdf_iwm( kt, Kmm, avm, avt, avs ) ! internal wave (de Lavergne et al 2017)
! !* Lateral boundary conditions (sign unchanged)
- IF(nn_hls==1) THEN
+ IF(nn_hls==1) THEN ! if nn_hls==2 lbc_lnk done in stp routines
IF( l_zdfsh2 ) THEN
CALL lbc_lnk( 'zdfphy', avm_k, 'W', 1.0_wp , avt_k, 'W', 1.0_wp, &
& avm , 'W', 1.0_wp , avt , 'W', 1.0_wp , avs , 'W', 1.0_wp )
diff --git a/src/OCE/ZDF/zdftke.F90 b/src/OCE/ZDF/zdftke.F90
index 653dfd0d..8bbdc4d8 100644
--- a/src/OCE/ZDF/zdftke.F90
+++ b/src/OCE/ZDF/zdftke.F90
@@ -698,7 +698,7 @@ CONTAINS
END SUBROUTINE tke_avn
- SUBROUTINE zdf_tke_init( Kmm )
+ SUBROUTINE zdf_tke_init
!!----------------------------------------------------------------------
!! *** ROUTINE zdf_tke_init ***
!!
@@ -714,7 +714,6 @@ CONTAINS
!!----------------------------------------------------------------------
USE zdf_oce , ONLY : ln_zdfiwm ! Internal Wave Mixing flag
!!
- INTEGER, INTENT(in) :: Kmm ! time level index
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: ios
!!
diff --git a/src/OCE/stprk3_stg.F90 b/src/OCE/stprk3_stg.F90
index 4ac8e54b..6853de23 100644
--- a/src/OCE/stprk3_stg.F90
+++ b/src/OCE/stprk3_stg.F90
@@ -21,6 +21,7 @@ MODULE stprk3_stg
USE domqco ! quasi-eulerian coordinate (dom_qco_r3c routine)
USE dynspg_ts, ONLY: un_adv , vn_adv ! advective transport from N to N+1
USE bdydyn ! ocean open boundary conditions (define bdy_dyn)
+ USE lbclnk ! ocean lateral boundary conditions (or mpp link)
# if defined key_top
USE trc ! ocean passive tracers variables
USE trcadv ! passive tracers advection (trc_adv routine)
@@ -48,6 +49,8 @@ MODULE stprk3_stg
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ssha ! sea-surface height at N+1
REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: ua_b, va_b ! barotropic velocity at N+1
+ REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: r3ta, r3ua, r3va ! ssh/h_0 ratio at t,u,v-column at N+1
+ REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: r3fb, r3fa ! ssh/h_0 ratio at f-column at N and N+1
!! * Substitutions
# include "do_loop_substitute.h90"
@@ -117,6 +120,24 @@ CONTAINS
uu_b(:,:,Kaa) = r2_3 * uu_b(:,:,Kbb) + r1_3 * ua_b(:,:)
vv_b(:,:,Kaa) = r2_3 * vv_b(:,:,Kbb) + r1_3 * va_b(:,:)
!
+ !
+ ! !== ssh/h0 ratio at Kaa ==!
+ !
+ IF( .NOT.lk_linssh ) THEN ! "after" ssh/h_0 ratio at t,u,v-column
+ !
+ ALLOCATE( r3ta(jpi,jpj) , r3ua(jpi,jpj) , r3va(jpi,jpj) , r3fa(jpi,jpj) , r3fb(jpi,jpj) )
+ !
+ r3fb(:,:) = r3f(:,:)
+ CALL dom_qco_r3c_RK3( ssha, r3ta, r3ua, r3va, r3fa )
+ !
+ CALL lbc_lnk( 'stprk3_stg', r3ua, 'U', 1._wp, r3va, 'V', 1._wp, r3fa, 'F', 1._wp )
+ !
+ r3t(:,:,Kaa) = r2_3 * r3t(:,:,Kbb) + r1_3 * r3ta(:,:) ! at N+1/3 (Kaa)
+ r3u(:,:,Kaa) = r2_3 * r3u(:,:,Kbb) + r1_3 * r3ua(:,:)
+ r3v(:,:,Kaa) = r2_3 * r3v(:,:,Kbb) + r1_3 * r3va(:,:)
+ ! r3f already properly set up ! at N (Kmm)
+ ENDIF
+ !
! !---------------!
CASE ( 2 ) !== Stage 2 ==! Kbb = N ; Kmm = N+1/3 ; Kaa = N+1/2
! !---------------!
@@ -127,7 +148,14 @@ CONTAINS
! ! set ssh and (uu_b,vv_b) at N+1/2 (Kaa)
ssh (:,:,Kaa) = r1_2 * ( ssh (:,:,Kbb) + ssha(:,:) )
uu_b(:,:,Kaa) = r1_2 * ( uu_b(:,:,Kbb) + ua_b(:,:) )
- vv_b(:,:,Kaa) = r1_2 * ( vv_b(:,:,Kbb) + va_b(:,:) )
+ vv_b(:,:,Kaa) = r1_2 * ( vv_b(:,:,Kbb) + va_b(:,:) )
+ !
+ IF( .NOT.lk_linssh ) THEN
+ r3t(:,:,Kaa) = r1_2 * ( r3t(:,:,Kbb) + r3ta(:,:) ) ! at N+1/2 (Kaa)
+ r3u(:,:,Kaa) = r1_2 * ( r3u(:,:,Kbb) + r3ua(:,:) )
+ r3v(:,:,Kaa) = r1_2 * ( r3v(:,:,Kbb) + r3va(:,:) )
+ r3f(:,:) = r2_3 * r3fb(:,:) + r1_3 * r3fa(:,:) ! at N+1/3 (Kmm)
+ ENDIF
!
! !---------------!
CASE ( 3 ) !== Stage 3 ==! Kbb = N ; Kmm = N+1/2 ; Kaa = N+1
@@ -142,22 +170,29 @@ CONTAINS
!
DEALLOCATE( ssha , ua_b , va_b )
!
+ IF( .NOT.lk_linssh ) THEN
+ r3t(:,:,Kaa) = r3ta(:,:) ! at N+1 (Kaa)
+ r3u(:,:,Kaa) = r3ua(:,:)
+ r3v(:,:,Kaa) = r3va(:,:)
+ r3f(:,: ) = r1_2 * ( r3fb(:,:) + r3fa(:,:) ) ! at N+1/2 (Kmm)
+ DEALLOCATE( r3ta, r3ua, r3va, r3fb ) ! deallocate all r3. except r3fa which will be
+ ! ! saved in r3f at the end of the time integration and then deallocated
+ !
+ ENDIF
+ !
END SELECT
!
- ! !== ssh/h0 ratio at Kaa ==!
- !
- IF( .NOT.lk_linssh ) CALL dom_qco_r3c( ssh(:,:,Kaa), r3t(:,:,Kaa), r3u(:,:,Kaa), r3v(:,:,Kaa), r3f(:,:) ) ! "after" ssh/h_0 ratio at t,u,v-column
- !
- !
! !== advective velocity at Kmm ==!
!
! !- horizontal components -! (zaU,zaV)
- zub(:,:) = un_adv(:,:)*r1_hu(:,:,Kmm) - uu_b(:,:,Kmm) ! barotropic velocity correction
- zvb(:,:) = vn_adv(:,:)*r1_hv(:,:,Kmm) - vv_b(:,:,Kmm)
- DO jk = 1, jpkm1 ! horizontal advective velocity
- zaU(:,:,jk) = uu(:,:,jk,Kmm) + zub(:,:)*umask(:,:,jk)
- zaV(:,:,jk) = vv(:,:,jk,Kmm) + zvb(:,:)*vmask(:,:,jk)
- END DO
+ DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
+ zub(ji,jj) = un_adv(ji,jj)*r1_hu(ji,jj,Kmm) - uu_b(ji,jj,Kmm) ! barotropic velocity correction
+ zvb(ji,jj) = vn_adv(ji,jj)*r1_hv(ji,jj,Kmm) - vv_b(ji,jj,Kmm)
+ END_2D
+ DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 ) ! horizontal advective velocity
+ zaU(ji,jj,jk) = uu(ji,jj,jk,Kmm) + zub(ji,jj)*umask(ji,jj,jk)
+ zaV(ji,jj,jk) = vv(ji,jj,jk,Kmm) + zvb(ji,jj)*vmask(ji,jj,jk)
+ END_3D
! !- vertical components -! ww
CALL wzv ( kstp, Kbb, Kmm, Kaa, zaU, zaV, ww ) ! ww cross-level velocity
@@ -223,9 +258,7 @@ CONTAINS
IF(.NOT. ln_trcadv_mus ) THEN
!
DO jn = 1, jptra
- DO_3D( 0, 0, 0, 0, 1, jpkm1 )
- tr(ji,jj,jk,jn,Krhs) = 0._wp ! set tracer trends to zero
- END_3D
+ tr(:,:,:,jn,Krhs) = 0._wp ! set tracer trends to zero !!st ::: required because of tra_adv new loops
END DO
! !== advection of passive tracers ==!
rDt_trc = rDt
@@ -253,9 +286,7 @@ CONTAINS
! !---------------!
!
DO jn = 1, jptra
- DO_3D( 0, 0, 0, 0, 1, jpkm1 )
- tr(ji,jj,jk,jn,Krhs) = 0._wp ! set tracer trends to zero
- END_3D
+ tr(:,:,:,jn,Krhs) = 0._wp
END DO
! !== advection of passive tracers ==!
rDt_trc = rDt
@@ -284,6 +315,7 @@ CONTAINS
ts(:,:,:,jn,Krhs) = 0._wp
END DO
+ CALL eos( ts(:,:,:,:,Kmm), rhd, rhop, gdept_0 ) ! now in potential density for tra_mle computation
!===>>> CAUTION here may be without GM velocity but stokes drift required ! 0 barotropic divergence for GM != 0 barotropic divergence for SD
!!st consistence 2D / 3D - flux de masse
CALL tra_adv( kstp, Kbb, Kmm, ts, Krhs, zaU, zaV, ww ) ! hor. + vert. advection ==> RHS
@@ -386,6 +418,9 @@ CONTAINS
CALL tra_zdf( kstp, Kbb, Kmm, Krhs, ts , Kaa ) ! vertical mixing and after tracer fields
IF( ln_zdfnpc ) CALL tra_npc( kstp, Kmm, Krhs, ts , Kaa ) ! update after fields by non-penetrative convection
!
+ r3f(:,:) = r3fa(:,:) ! save r3fa in r3f before deallocation
+ DEALLOCATE( r3fa ) ! (r3f = r3f(Kbb) of the next time step)
+ !
END SELECT
! !== correction of the barotropic (all stages) ==! at Kaa = N+1/3, N+1/2 or N+1
! ! barotropic velocity correction
@@ -393,10 +428,9 @@ CONTAINS
zvb(A2D(0)) = vv_b(A2D(0),Kaa) - SUM( e3v_0(A2D(0),:)*vv(A2D(0),:,Kaa), 3 ) * r1_hv_0(A2D(0))
!
DO jk = 1, jpkm1 ! corrected horizontal velocity
- uu(:,:,jk,Kaa) = uu(:,:,jk,Kaa) + zub(:,:)*umask(:,:,jk)
- vv(:,:,jk,Kaa) = vv(:,:,jk,Kaa) + zvb(:,:)*vmask(:,:,jk)
+ uu(A2D(0),jk,Kaa) = uu(A2D(0),jk,Kaa) + zub(A2D(0))*umask(A2D(0),jk)
+ vv(A2D(0),jk,Kaa) = vv(A2D(0),jk,Kaa) + zvb(A2D(0))*vmask(A2D(0),jk)
END DO
-!!st pourquoi ne pas mettre A2D(0) ici ?
!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -408,8 +442,10 @@ CONTAINS
CALL Agrif_dyn( kstp, kstg )
# endif
! !* local domain boundaries (T-point, unchanged sign)
- CALL lbc_lnk_multi( 'stp_RK3_stg', uu(:,:,:, Kaa), 'U', -1., vv(:,:,: ,Kaa), 'V', -1. &
- & , ts(:,:,:,jp_tem,Kaa), 'T', 1., ts(:,:,:,jp_sal,Kaa), 'T', 1. )
+ CALL lbc_lnk( 'stp_RK3_stg', uu(:,:,:, Kaa), 'U', -1., vv(:,:,: ,Kaa), 'V', -1. &
+ & , ts(:,:,:,jp_tem,Kaa), 'T', 1., ts(:,:,:,jp_sal,Kaa), 'T', 1. )
+ !
+ IF( nn_hls == 2 .AND. l_zdfsh2 ) CALL lbc_lnk( 'stp_RK3_stg', avm_k, 'W', 1.0_wp ) ! lbc_lnk needed for zdf_sh2 when using nn_hls = 2, moved here to allow tiling in zdf_phy
!
! !* BDY open boundaries
IF( ln_bdy ) THEN
diff --git a/src/TOP/trcini.F90 b/src/TOP/trcini.F90
index 5953dd54..486f5208 100644
--- a/src/TOP/trcini.F90
+++ b/src/TOP/trcini.F90
@@ -246,13 +246,6 @@ CONTAINS
!
tr(:,:,:,:,Kaa) = 0._wp
!
- IF( ln_trcbc .AND. lltrcbc ) THEN
- CALL trc_bc_ini ( jptra, Kmm ) ! set tracers Boundary Conditions
- CALL trc_bc ( nit000, Kmm, tr, Kaa ) ! tracers: surface and lateral Boundary Conditions
- ENDIF
- !
- IF( ln_trcais ) CALL trc_ais_ini ! set tracers from Antarctic Ice Sheet
- !
IF( ln_rsttr ) THEN ! restart from a file
!
CALL trc_rst_read( Kbb, Kmm )
diff --git a/src/TOP/trcrst.F90 b/src/TOP/trcrst.F90
index 05179e6e..78ff2bbc 100644
--- a/src/TOP/trcrst.F90
+++ b/src/TOP/trcrst.F90
@@ -159,9 +159,6 @@ CONTAINS
! prognostic variables
! --------------------
#if defined key_RK3
- DO jn = 1, jptra ! MLF : After time step (before the swap) put in TRN
- CALL iom_rstput( kt, nitrst, numrtw, 'TRN'//ctrcnm(jn), tr(:,:,:,jn,Kaa) )
- END DO
DO jn = 1, jptra ! RK3 : After time step (before the swap) put in TRB
CALL iom_rstput( kt, nitrst, numrtw, 'TRB'//ctrcnm(jn), tr(:,:,:,jn,Kaa) )
END DO
diff --git a/tests/ISOMIP+/EXPREF/namelist_cfg b/tests/ISOMIP+/EXPREF/namelist_cfg
index d2b10f0c..0b8ae4c9 100644
--- a/tests/ISOMIP+/EXPREF/namelist_cfg
+++ b/tests/ISOMIP+/EXPREF/namelist_cfg
@@ -193,7 +193,6 @@ rn_Dt = 720.
! ! vel_stab = velocity and stability dependent transfert coeficient (Holland et al. 1999 for a complete description)
rn_gammat0 = 0.0215 ! gammat coefficient used in blk formula
rn_gammas0 = 0.614e-3 ! gammas coefficient used in blk formula
- rn_vtide = 0.01 ! tidal velocity [m/s]
!
rn_htbl = 20. ! thickness of the top boundary layer (Losh et al. 2008)
! ! 0 => thickness of the tbl = thickness of the first wet cell
@@ -279,7 +278,7 @@ rn_Dt = 720.
&namdrg_top ! TOP friction (ln_drg_OFF =F & ln_isfcav=T)
!-----------------------------------------------------------------------
rn_Cd0 = 2.5e-3 ! drag coefficient [-]
- rn_ke0 = 0.0e-3 ! background kinetic energy [m2/s2] (non-linear cases)
+ rn_ke0 = 1.0e-4 ! background kinetic energy [m2/s2] (non-linear cases)
/
!-----------------------------------------------------------------------
&namdrg_bot ! BOTTOM friction (ln_drg_OFF =F)
diff --git a/tests/ISOMIP+/MY_SRC/eosbn2.F90 b/tests/ISOMIP+/MY_SRC/eosbn2.F90
index b83d87e1..d2da8cfd 100644
--- a/tests/ISOMIP+/MY_SRC/eosbn2.F90
+++ b/tests/ISOMIP+/MY_SRC/eosbn2.F90
@@ -16,7 +16,7 @@ MODULE eosbn2
!! 3.0 ! 2006-08 (G. Madec) add tfreez function (now eos_fzp function)
!! 3.3 ! 2010-05 (C. Ethe, G. Madec) merge TRC-TRA
!! - ! 2010-10 (G. Nurser, G. Madec) add alpha/beta used in ldfslp
- !! 3.7 ! 2012-0 3 (F. Roquet, G. Madec) add primitive of alpha and beta used in PE computation
+ !! 3.7 ! 2012-03 (F. Roquet, G. Madec) add primitive of alpha and beta used in PE computation
!! - ! 2012-05 (F. Roquet) add Vallis and original JM95 equation of state
!! - ! 2013-04 (F. Roquet, G. Madec) add eos_rab, change bn2 computation and reorganize the module
!! - ! 2014-09 (F. Roquet) add TEOS-10, S-EOS, and modify EOS-80
@@ -294,7 +294,7 @@ CONTAINS
!
CASE( np_leos ) !== linear ISOMIP EOS ==!
!
- DO_3D( 1, 1, 1, 1, 1, jpkm1 )
+ DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )
zt = pts (ji,jj,jk,jp_tem,Knn) - (-1._wp)
zs = pts (ji,jj,jk,jp_sal,Knn) - 34.2_wp
zh = gdept(ji,jj,jk, Knn)
diff --git a/tests/ISOMIP+/MY_SRC/isf_oce.F90 b/tests/ISOMIP+/MY_SRC/isf_oce.F90
index 9f5b8ebd..3bf1b30d 100644
--- a/tests/ISOMIP+/MY_SRC/isf_oce.F90
+++ b/tests/ISOMIP+/MY_SRC/isf_oce.F90
@@ -37,7 +37,6 @@ MODULE isf_oce
LOGICAL , PUBLIC :: ln_isfcav_mlt !: logical for the use of ice shelf parametrisation
REAL(wp) , PUBLIC :: rn_gammat0 !: temperature exchange coeficient []
REAL(wp) , PUBLIC :: rn_gammas0 !: salinity exchange coeficient []
- REAL(wp) , PUBLIC :: rn_vtide !: tidal background velocity (can be different to what is used in the
REAL(wp) , PUBLIC :: rn_htbl !: Losch top boundary layer thickness [m]
REAL(wp) , PUBLIC :: rn_isfload_T !:
REAL(wp) , PUBLIC :: rn_isfload_S !:
diff --git a/tests/ISOMIP+/MY_SRC/isfcavgam.F90 b/tests/ISOMIP+/MY_SRC/isfcavgam.F90
index ac0697a8..6c0ac2a4 100644
--- a/tests/ISOMIP+/MY_SRC/isfcavgam.F90
+++ b/tests/ISOMIP+/MY_SRC/isfcavgam.F90
@@ -94,9 +94,9 @@ CONTAINS
pgt(:,:) = rn_gammat0
pgs(:,:) = rn_gammas0
CASE ( 'vel' ) ! gamma is proportional to u*
- CALL gammats_vel ( zutbl, zvtbl, rCd0_top, rn_vtide**2, pgt, pgs )
+ CALL gammats_vel ( zutbl, zvtbl, rCd0_top, r_ke0_top, pgt, pgs )
CASE ( 'vel_stab' ) ! gamma depends of stability of boundary layer and u*
- CALL gammats_vel_stab (Kmm, pttbl, pstbl, zutbl, zvtbl, rCd0_top, rn_vtide**2, pqoce, pqfwf, pRc, pgt, pgs )
+ CALL gammats_vel_stab (Kmm, pttbl, pstbl, zutbl, zvtbl, rCd0_top, r_ke0_top, pqoce, pqfwf, pRc, pgt, pgs )
CASE DEFAULT
CALL ctl_stop('STOP','method to compute gamma (cn_gammablk) is unknown (should not see this)')
END SELECT
diff --git a/tests/ISOMIP+/MY_SRC/isfstp.F90 b/tests/ISOMIP+/MY_SRC/isfstp.F90
index 58b388ac..2ab92084 100644
--- a/tests/ISOMIP+/MY_SRC/isfstp.F90
+++ b/tests/ISOMIP+/MY_SRC/isfstp.F90
@@ -38,7 +38,7 @@ MODULE isfstp
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
- !! $Id: isfstp.F90 11876 2019-11-08 11:26:42Z mathiot $
+ !! $Id: isfstp.F90 15574 2021-12-03 19:32:50Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
@@ -194,6 +194,11 @@ CONTAINS
WRITE(numout,*)
!
IF ( ln_isf ) THEN
+#if key_qco
+# if ! defined key_isf
+ CALL ctl_stop( 'STOP', 'isf_ctl: ice shelf requires both ln_isf=T AND key_isf activated' )
+# endif
+#endif
WRITE(numout,*) ' Add debug print in isf module ln_isfdebug = ', ln_isfdebug
WRITE(numout,*)
WRITE(numout,*) ' melt inside the cavity ln_isfcav_mlt = ', ln_isfcav_mlt
@@ -204,7 +209,7 @@ CONTAINS
IF ( TRIM(cn_gammablk) .NE. 'spe' ) THEN
WRITE(numout,*) ' gammat coefficient rn_gammat0 = ', rn_gammat0
WRITE(numout,*) ' gammas coefficient rn_gammas0 = ', rn_gammas0
- WRITE(numout,*) ' top background ke used (from namdrg_top) rn_vtide**2 = ', rn_vtide**2
+ WRITE(numout,*) ' top background ke used (from namdrg_top) rn_ke0 = ', r_ke0_top
WRITE(numout,*) ' top drag coef. used (from namdrg_top) rn_Cd0 = ', r_Cdmin_top
END IF
END IF
@@ -299,7 +304,7 @@ CONTAINS
& ln_isfcav_mlt , cn_isfcav_mlt , sn_isfcav_fwf , &
& ln_isfpar_mlt , cn_isfpar_mlt , sn_isfpar_fwf , &
& sn_isfpar_zmin, sn_isfpar_zmax, sn_isfpar_Leff, &
- & ln_isfcpl , nn_drown , ln_isfcpl_cons, ln_isfdebug, rn_vtide, &
+ & ln_isfcpl , nn_drown , ln_isfcpl_cons, ln_isfdebug, &
& cn_isfload , rn_isfload_T , rn_isfload_S , cn_isfdir , &
& rn_isfpar_bg03_gt0
!!----------------------------------------------------------------------
diff --git a/tests/ISOMIP+/MY_SRC/istate.F90 b/tests/ISOMIP+/MY_SRC/istate.F90
index 35af0fd2..0f32b25a 100644
--- a/tests/ISOMIP+/MY_SRC/istate.F90
+++ b/tests/ISOMIP+/MY_SRC/istate.F90
@@ -32,6 +32,7 @@ MODULE istate
USE in_out_manager ! I/O manager
USE iom ! I/O library
USE lib_mpp ! MPP library
+ USE lbclnk ! lateal boundary condition / mpp exchanges
USE restart ! restart
#if defined key_agrif
@@ -49,7 +50,7 @@ MODULE istate
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
- !! $Id: istate.F90 11423 2019-08-08 14:02:49Z mathiot $
+ !! $Id: istate.F90 15581 2021-12-07 13:08:22Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
@@ -59,6 +60,8 @@ CONTAINS
!! *** ROUTINE istate_init ***
!!
!! ** Purpose : Initialization of the dynamics and tracer fields.
+ !!
+ !! ** Method :
!!----------------------------------------------------------------------
INTEGER, INTENT( in ) :: Kbb, Kmm, Kaa ! ocean time level indices
!
@@ -86,7 +89,7 @@ CONTAINS
#endif
#if defined key_agrif
- IF ( (.NOT.Agrif_root()).AND.ln_init_chfrpar ) THEN
+ IF ( .NOT.Agrif_root() .AND. ln_init_chfrpar ) THEN
numror = 0 ! define numror = 0 -> no restart file to read
ln_1st_euler = .true. ! Set time-step indicator at nit000 (euler forward)
CALL day_init
@@ -125,37 +128,50 @@ CONTAINS
DO jk = 1, jpk
zgdept(:,:,jk) = gdept(:,:,jk,Kbb)
END DO
- CALL usr_def_istate( zgdept, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb) )
+ CALL usr_def_istate( zgdept, tmask, ts(:,:,:,:,Kbb), uu(:,:,:,Kbb), vv(:,:,:,Kbb) )
+ ! make sure that periodicities are properly applied
+ CALL lbc_lnk( 'istate', ts(:,:,:,jp_tem,Kbb), 'T', 1._wp, ts(:,:,:,jp_sal,Kbb), 'T', 1._wp, &
+ & uu(:,:,:, Kbb), 'U', -1._wp, vv(:,:,:, Kbb), 'V', -1._wp )
ENDIF
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
+ ! Initialize "now" barotropic velocities:
+ ! Do it whatever the free surface method, these arrays being used eventually
!
+!!gm the use of umask & vmask is not necessary below as uu(:,:,:,Kmm), vv(:,:,:,Kmm), uu(:,:,:,Kbb), vv(:,:,:,Kbb) are always masked
+#if ! defined key_RK3
uu_b(:,:,Kmm) = 0._wp ; vv_b(:,:,Kmm) = 0._wp
- uu_b(:,:,Kbb) = 0._wp ; vv_b(:,:,Kbb) = 0._wp
- !
-!!gm the use of umsak & vmask is not necessary below as uu(:,:,:,Kmm), vv(:,:,:,Kmm), uu(:,:,:,Kbb), vv(:,:,:,Kbb) are always masked
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)
+#endif
!
- uu_b(:,:,Kbb) = uu_b(:,:,Kbb) * r1_hu(:,:,Kbb)
- vv_b(:,:,Kbb) = vv_b(:,:,Kbb) * r1_hv(:,:,Kbb)
+#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
!
END SUBROUTINE istate_init
diff --git a/tests/ISOMIP+/MY_SRC/sbcfwb.F90 b/tests/ISOMIP+/MY_SRC/sbcfwb.F90
index 5aee8bb6..0102aa12 100644
--- a/tests/ISOMIP+/MY_SRC/sbcfwb.F90
+++ b/tests/ISOMIP+/MY_SRC/sbcfwb.F90
@@ -35,8 +35,9 @@ MODULE sbcfwb
PUBLIC sbc_fwb ! routine called by step
REAL(wp) :: rn_fwb0 ! initial freshwater adjustment flux [kg/m2/s] (nn_fwb = 2 only)
- REAL(wp) :: a_fwb ! annual domain averaged freshwater budget from the
- ! previous year
+ REAL(wp) :: a_fwb ! annual domain averaged freshwater budget from the previous year
+ REAL(wp) :: a_fwb_b ! annual domain averaged freshwater budget from the year before or at initial state
+ REAL(wp) :: a_fwb_ini ! initial domain averaged freshwater budget
REAL(wp) :: area ! global mean ocean surface (interior domain)
!!----------------------------------------------------------------------
@@ -128,68 +129,65 @@ CONTAINS
IF( iom_use('vflx_fwb_cea') ) CALL iom_put( 'vflx_fwb_cea', z_fwfprv(1) * tmask(:,:,1) )
ENDIF
!
- CASE ( 4 ) !== global mean fwf set to zero (ISOMIP case) ==!
- !
- IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
- z_fwf = glob_sum( 'sbcfwb', e1e2t(:,:) * ( emp(:,:) - rnf(:,:) - fwfisf_cav(:,:) - fwfisf_par(:,:) - snwice_fmass(:,:) ) )
- !
- ! correction for ice sheet coupling testing (ie remove the excess through the surface)
- ! test impact on the melt as conservation correction made in depth
- ! test conservation level as sbcfwb is conserving
- ! avoid the model to blow up for large ssh drop (isomip OCEAN3 with melt switch off and uniform T/S)
- IF (ln_isfcpl .AND. ln_isfcpl_cons) THEN
- z_fwf = z_fwf + glob_sum( 'sbcfwb', e1e2t(:,:) * risfcpl_cons_ssh(:,:) * rho0 )
- END IF
- !
- z_fwf = z_fwf / area
- zcoef = z_fwf * rcp
- emp(:,:) = emp(:,:) - z_fwf * tmask(:,:,1) ! (Eq. 34 AD2015)
- qns(:,:) = qns(:,:) + zcoef * sst_m(:,:) * tmask(:,:,1) ! (Eq. 35 AD2015) ! use sst_m to avoid generation of any bouyancy fluxes
- sfx(:,:) = sfx(:,:) + z_fwf * sss_m(:,:) * tmask(:,:,1) ! (Eq. 36 AD2015) ! use sss_m to avoid generation of any bouyancy fluxes
- ENDIF
- !
CASE ( 2 ) !== fw adjustment based on fw budget at the end of the previous year ==!
- !
- IF( kt == nit000 ) THEN ! initialisation
- ! ! set the fw adjustment (a_fwb)
- IF ( ln_rstart .AND. iom_varid( numror, 'a_fwb', ldstop = .FALSE. ) > 0 ) THEN ! as read from restart file
- IF(lwp) WRITE(numout,*) 'sbc_fwb : reading FW-budget adjustment from restart file'
- CALL iom_get( numror, 'a_fwb', a_fwb )
- ELSE ! as specified in namelist
- a_fwb = rn_fwb0
+ ! simulation is supposed to start 1st of January
+ IF( kt == nit000 ) THEN ! initialisation
+ ! ! set the fw adjustment (a_fwb)
+ IF ( ln_rstart .AND. iom_varid( numror, 'a_fwb_b', ldstop = .FALSE. ) > 0 & ! as read from restart file
+ & .AND. iom_varid( numror, 'a_fwb', ldstop = .FALSE. ) > 0 ) THEN
+ IF(lwp) WRITE(numout,*) 'sbc_fwb : reading freshwater-budget from restart file'
+ CALL iom_get( numror, 'a_fwb_b', a_fwb_b )
+ CALL iom_get( numror, 'a_fwb' , a_fwb )
+ !
+ a_fwb_ini = a_fwb_b
+ ELSE ! as specified in namelist
+ IF(lwp) WRITE(numout,*) 'sbc_fwb : setting freshwater-budget from namelist rn_fwb0'
+ a_fwb = rn_fwb0
+ a_fwb_b = 0._wp ! used only the first year then it is replaced by a_fwb_ini
+ !
+ a_fwb_ini = glob_sum( 'sbcfwb', e1e2t(:,:) * ( ssh(:,:,Kmm) + snwice_mass(:,:) * r1_rho0 ) ) &
+ & * rho0 / ( area * rday * REAL(nyear_len(1), wp) )
END IF
!
- IF(lwp)WRITE(numout,*)
- IF(lwp)WRITE(numout,*)'sbc_fwb : initial freshwater-budget adjustment = ', a_fwb, 'kg/m2/s'
+ IF(lwp) WRITE(numout,*)
+ IF(lwp) WRITE(numout,*)'sbc_fwb : freshwater-budget at the end of previous year = ', a_fwb , 'kg/m2/s'
+ IF(lwp) WRITE(numout,*)' freshwater-budget at initial state = ', a_fwb_ini, 'kg/m2/s'
+ !
+ ELSE
+ ! at the end of year n:
+ ikty = nyear_len(1) * 86400 / NINT(rn_Dt)
+ IF( MOD( kt, ikty ) == 0 ) THEN ! Update a_fwb at the last time step of a year
+ ! It should be the first time step of a year MOD(kt-1,ikty) but then the restart would be wrong
+ ! Hence, we make a small error here but the code is restartable
+ a_fwb_b = a_fwb_ini
+ ! mean sea level taking into account ice+snow
+ a_fwb = glob_sum( 'sbcfwb', e1e2t(:,:) * ( ssh(:,:,Kmm) + snwice_mass(:,:) * r1_rho0 ) )
+ a_fwb = a_fwb * rho0 / ( area * rday * REAL(nyear_len(1), wp) ) ! convert in kg/m2/s
+ ENDIF
!
- ENDIF
- ! ! Update a_fwb if new year start
- ikty = 365 * 86400 / rn_Dt !!bug use of 365 days leap year or 360d year !!!!!!!
- IF( MOD( kt, ikty ) == 0 ) THEN
- ! mean sea level taking into account the ice+snow
- ! sum over the global domain
- a_fwb = glob_sum( 'sbcfwb', e1e2t(:,:) * ( ssh(:,:,Kmm) + snwice_mass(:,:) * r1_rho0 ) )
- a_fwb = a_fwb * 1.e+3 / ( area * rday * 365. ) ! convert in Kg/m3/s = mm/s
-!!gm ! !!bug 365d year
ENDIF
- !
- IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN ! correct the freshwater fluxes
- zcoef = a_fwb * rcp
- emp(:,:) = emp(:,:) + a_fwb * tmask(:,:,1)
- qns(:,:) = qns(:,:) - zcoef * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction
+ !
+ IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN ! correct the freshwater fluxes using previous year budget minus initial state
+ zcoef = ( a_fwb - a_fwb_b )
+ emp(:,:) = emp(:,:) + zcoef * tmask(:,:,1)
+ qns(:,:) = qns(:,:) - zcoef * rcp * sst_m(:,:) * tmask(:,:,1) ! account for change to the heat budget due to fw correction
! outputs
- IF( iom_use('hflx_fwb_cea') ) CALL iom_put( 'hflx_fwb_cea', -zcoef * sst_m(:,:) * tmask(:,:,1) )
- IF( iom_use('vflx_fwb_cea') ) CALL iom_put( 'vflx_fwb_cea', -a_fwb * tmask(:,:,1) )
+ IF( iom_use('hflx_fwb_cea') ) CALL iom_put( 'hflx_fwb_cea', -zcoef * rcp * sst_m(:,:) * tmask(:,:,1) )
+ IF( iom_use('vflx_fwb_cea') ) CALL iom_put( 'vflx_fwb_cea', -zcoef * tmask(:,:,1) )
ENDIF
! Output restart information
IF( lrst_oce ) THEN
IF(lwp) WRITE(numout,*)
IF(lwp) WRITE(numout,*) 'sbc_fwb : writing FW-budget adjustment to ocean restart file at it = ', kt
IF(lwp) WRITE(numout,*) '~~~~'
- CALL iom_rstput( kt, nitrst, numrow, 'a_fwb', a_fwb )
+ CALL iom_rstput( kt, nitrst, numrow, 'a_fwb_b', a_fwb_b )
+ CALL iom_rstput( kt, nitrst, numrow, 'a_fwb', a_fwb )
END IF
!
- IF( kt == nitend .AND. lwp ) WRITE(numout,*) 'sbc_fwb : final freshwater-budget adjustment = ', a_fwb, 'kg/m2/s'
+ IF( kt == nitend .AND. lwp ) THEN
+ WRITE(numout,*) 'sbc_fwb : freshwater-budget at the end of simulation (year now) = ', a_fwb , 'kg/m2/s'
+ WRITE(numout,*) ' freshwater-budget at initial state = ', a_fwb_b, 'kg/m2/s'
+ ENDIF
!
CASE ( 3 ) !== global fwf set to zero and spread out over erp area ==!
!
@@ -248,6 +246,26 @@ CONTAINS
ENDIF
DEALLOCATE( ztmsk_neg , ztmsk_pos , ztmsk_tospread , z_wgt , zerp_cor )
!
+ CASE ( 4 ) !== global mean fwf set to zero (ISOMIP case) ==!
+ !
+ IF( MOD( kt-1, kn_fsbc ) == 0 ) THEN
+ z_fwf = glob_sum( 'sbcfwb', e1e2t(:,:) * ( emp(:,:) - rnf(:,:) - fwfisf_cav(:,:) - fwfisf_par(:,:) - snwice_fmass(:,:) ) )
+ !
+ ! correction for ice sheet coupling testing (ie remove the excess through the surface)
+ ! test impact on the melt as conservation correction made in depth
+ ! test conservation level as sbcfwb is conserving
+ ! avoid the model to blow up for large ssh drop (isomip OCEAN3 with melt switch off and uniform T/S)
+ IF (ln_isfcpl .AND. ln_isfcpl_cons) THEN
+ z_fwf = z_fwf + glob_sum( 'sbcfwb', e1e2t(:,:) * risfcpl_cons_ssh(:,:) * rho0 )
+ END IF
+ !
+ z_fwf = z_fwf / area
+ zcoef = z_fwf * rcp
+ emp(:,:) = emp(:,:) - z_fwf * tmask(:,:,1) ! (Eq. 34 AD2015)
+ qns(:,:) = qns(:,:) + zcoef * sst_m(:,:) * tmask(:,:,1) ! (Eq. 35 AD2015) ! use sst_m to avoid generation of any bouyancy fluxes
+ sfx(:,:) = sfx(:,:) + z_fwf * sss_m(:,:) * tmask(:,:,1) ! (Eq. 36 AD2015) ! use sss_m to avoid generation of any bouyancy fluxes
+ ENDIF
+ !
CASE DEFAULT !== you should never be there ==!
CALL ctl_stop( 'sbc_fwb : wrong nn_fwb value for the FreshWater Budget correction, choose either 1, 2 or 3' )
!
diff --git a/tests/ISOMIP+/MY_SRC/tradmp.F90 b/tests/ISOMIP+/MY_SRC/tradmp.F90
index 734d1767..05ed5b09 100644
--- a/tests/ISOMIP+/MY_SRC/tradmp.F90
+++ b/tests/ISOMIP+/MY_SRC/tradmp.F90
@@ -23,7 +23,6 @@ MODULE tradmp
!!----------------------------------------------------------------------
USE oce ! ocean: variables
USE dom_oce ! ocean: domain variables
- USE c1d ! 1D vertical configuration
USE trd_oce ! trends: ocean variables
USE trdtra ! trends manager: tracers
USE zdf_oce ! ocean: vertical physics
@@ -55,7 +54,7 @@ MODULE tradmp
# include "domzgr_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2018)
- !! $Id: tradmp.F90 10425 2018-12-19 21:54:16Z smasson $
+ !! $Id: tradmp.F90 15574 2021-12-03 19:32:50Z techene $
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
@@ -96,15 +95,19 @@ CONTAINS
!
INTEGER :: ji, jj, jk, jn ! dummy loop indices
REAL(wp), DIMENSION(A2D(nn_hls),jpk,jpts) :: zts_dta
- REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t
+ REAL(wp), DIMENSION(:,:,:) , ALLOCATABLE :: zwrk
REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE :: ztrdts
!!----------------------------------------------------------------------
!
IF( ln_timing ) CALL timing_start('tra_dmp')
!
IF( l_trdtra .OR. iom_use('hflx_dmp_cea') .OR. iom_use('sflx_dmp_cea') ) THEN !* Save ta and sa trends
- ALLOCATE( ztrdts(jpi,jpj,jpk,jpts) )
- ztrdts(:,:,:,:) = pts(:,:,:,:,Krhs)
+ ALLOCATE( ztrdts(A2D(nn_hls),jpk,jpts) )
+ DO jn = 1, jpts
+ DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
+ ztrdts(ji,jj,jk,jn) = pts(ji,jj,jk,jn,Krhs)
+ END_3D
+ END DO
ENDIF
! !== input T-S data at kt ==!
CALL dta_tsd( kt, 'dmp', zts_dta ) ! read and interpolates T-S data at kt
@@ -142,16 +145,25 @@ CONTAINS
END SELECT
!
! outputs (clem trunk)
- DO jk = 1, jpk
- ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
- END DO
- !
- IF( iom_use('hflx_dmp_cea') ) &
- & CALL iom_put('hflx_dmp_cea', &
- & SUM( ( pts(:,:,:,jp_tem,Krhs) - ztrdts(:,:,:,jp_tem) ) * ze3t(:,:,:), dim=3 ) * rcp * rho0 ) ! W/m2
- IF( iom_use('sflx_dmp_cea') ) &
- & CALL iom_put('sflx_dmp_cea', &
- & SUM( ( pts(:,:,:,jp_sal,Krhs) - ztrdts(:,:,:,jp_sal) ) * ze3t(:,:,:), dim=3 ) * rho0 ) ! g/m2/s
+ IF( iom_use('hflx_dmp_cea') .OR. iom_use('sflx_dmp_cea') ) THEN
+ ALLOCATE( zwrk(A2D(nn_hls),jpk) ) ! Needed to handle expressions containing e3t when using key_qco or key_linssh
+ zwrk(:,:,:) = 0._wp
+
+ IF( iom_use('hflx_dmp_cea') ) THEN
+ DO_3D( 0, 0, 0, 0, 1, jpk )
+ zwrk(ji,jj,jk) = ( pts(ji,jj,jk,jp_tem,Krhs) - ztrdts(ji,jj,jk,jp_tem) ) * e3t(ji,jj,jk,Kmm)
+ END_3D
+ CALL iom_put('hflx_dmp_cea', SUM( zwrk(:,:,:), dim=3 ) * rcp * rho0 ) ! W/m2
+ ENDIF
+ IF( iom_use('sflx_dmp_cea') ) THEN
+ DO_3D( 0, 0, 0, 0, 1, jpk )
+ zwrk(ji,jj,jk) = ( pts(ji,jj,jk,jp_sal,Krhs) - ztrdts(ji,jj,jk,jp_sal) ) * e3t(ji,jj,jk,Kmm)
+ END_3D
+ CALL iom_put('sflx_dmp_cea', SUM( zwrk(:,:,:), dim=3 ) * rho0 ) ! g/m2/s
+ ENDIF
+
+ DEALLOCATE( zwrk )
+ ENDIF
!
IF( l_trdtra ) THEN ! trend diagnostic
ztrdts(:,:,:,:) = pts(:,:,:,:,Krhs) - ztrdts(:,:,:,:)
--
GitLab