From a044e9c28afaef407419df3ed6c4c55ab3bcbb6e Mon Sep 17 00:00:00 2001
From: Sibylle TECHENE <techenes@irene190.c-irene.tgcc.ccc.cea.fr>
Date: Mon, 11 Apr 2022 16:34:53 +0200
Subject: [PATCH] #2 RK3 time-stepping currently uses forcing read at n+1/2
during a time-step
---
src/OCE/ISF/isfhdiv.F90 | 5 +++++
src/OCE/SBC/sbcrnf.F90 | 18 +++++++++++++++---
src/OCE/SBC/sbcssm.F90 | 14 +++++++++++++-
3 files changed, 33 insertions(+), 4 deletions(-)
diff --git a/src/OCE/ISF/isfhdiv.F90 b/src/OCE/ISF/isfhdiv.F90
index d8b8264d8..772ab9a74 100644
--- a/src/OCE/ISF/isfhdiv.F90
+++ b/src/OCE/ISF/isfhdiv.F90
@@ -71,6 +71,7 @@ CONTAINS
!
END SUBROUTINE isf_hdiv
+
SUBROUTINE isf_hdiv_mlt(ktop, kbot, phtbl, pfrac, pfwf, pfwf_b, phdiv)
!!----------------------------------------------------------------------
!! *** SUBROUTINE sbc_isf_div ***
@@ -97,7 +98,11 @@ CONTAINS
!
! compute integrated divergence correction
DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
+#if defined key_RK3
+ zhdiv(ji,jj) = pfwf(ji,jj) * r1_rho0 / phtbl(ji,jj)
+#else
zhdiv(ji,jj) = 0.5_wp * ( pfwf(ji,jj) + pfwf_b(ji,jj) ) * r1_rho0 / phtbl(ji,jj)
+#endif
END_2D
!
! update divergence at each level affected by ice shelf top boundary layer
diff --git a/src/OCE/SBC/sbcrnf.F90 b/src/OCE/SBC/sbcrnf.F90
index 785ba0d96..be3f221b8 100644
--- a/src/OCE/SBC/sbcrnf.F90
+++ b/src/OCE/SBC/sbcrnf.F90
@@ -212,7 +212,11 @@ CONTAINS
IF( ln_linssh ) THEN !* constant volume case : just apply the runoff input flow
DO_2D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
DO jk = 1, nk_rnf(ji,jj)
- phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact / h_rnf(ji,jj)
+#if defined key_RK3
+ phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - rnf(ji,jj) * r1_rho0 / h_rnf(ji,jj) ! RK3: rnf forcing at n+1/2
+#else
+ phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact / h_rnf(ji,jj) ! MLF: rnf forcing at Kmm (n)
+#endif
END DO
END_2D
ELSE !* variable volume case
@@ -224,7 +228,11 @@ CONTAINS
END_2D
DO_2D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls ) ! apply the runoff input flow
DO jk = 1, nk_rnf(ji,jj)
- phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact / h_rnf(ji,jj)
+#if defined key_RK3
+ phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - rnf(ji,jj) * r1_rho0 / h_rnf(ji,jj) ! RK3: rnf forcing at n+1/2
+#else
+ phdivn(ji,jj,jk) = phdivn(ji,jj,jk) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact / h_rnf(ji,jj) ! MLF: rnf forcing at Kmm (n)
+#endif
END DO
END_2D
ENDIF
@@ -233,7 +241,11 @@ CONTAINS
h_rnf (ji,jj) = e3t(ji,jj,1,Kmm) ! update h_rnf to be depth of top box
END_2D
DO_2D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
- phdivn(ji,jj,1) = phdivn(ji,jj,1) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact / e3t(ji,jj,1,Kmm)
+#if defined key_RK3
+ phdivn(ji,jj,1) = phdivn(ji,jj,1) - rnf(ji,jj) * r1_rho0 / e3t(ji,jj,1,Kmm) ! RK3: rnf forcing at n+1/2
+#else
+ phdivn(ji,jj,1) = phdivn(ji,jj,1) - ( rnf(ji,jj) + rnf_b(ji,jj) ) * zfact / e3t(ji,jj,1,Kmm) ! MLF: rnf forcing at Kmm (n)
+#endif
END_2D
ENDIF
!
diff --git a/src/OCE/SBC/sbcssm.F90 b/src/OCE/SBC/sbcssm.F90
index 5f3daebf6..7ff88f827 100644
--- a/src/OCE/SBC/sbcssm.F90
+++ b/src/OCE/SBC/sbcssm.F90
@@ -79,7 +79,11 @@ CONTAINS
ENDIF
sss_m(:,:) = zts(:,:,jp_sal)
! ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics)
- IF( ln_apr_dyn ) THEN ; ssh_m(:,:) = ssh(:,:,Kmm) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
+#if defined key_RK3
+ IF( ln_apr_dyn ) THEN ; ssh_m(:,:) = ssh(:,:,Kmm) - ssh_ib(:,:) ! RK3: forcing at n+1/2
+#else
+ IF( ln_apr_dyn ) THEN ; ssh_m(:,:) = ssh(:,:,Kmm) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) ! MLF: forcing at n (Kmm)
+#endif
ELSE ; ssh_m(:,:) = ssh(:,:,Kmm)
ENDIF
!
@@ -102,7 +106,11 @@ CONTAINS
ENDIF
sss_m(:,:) = zcoef * zts(:,:,jp_sal)
! ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics)
+#if defined key_RK3
+ IF( ln_apr_dyn ) THEN ; ssh_m(:,:) = zcoef * ( ssh(:,:,Kmm) - ssh_ib(:,:) )
+#else
IF( ln_apr_dyn ) THEN ; ssh_m(:,:) = zcoef * ( ssh(:,:,Kmm) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) ) )
+#endif
ELSE ; ssh_m(:,:) = zcoef * ssh(:,:,Kmm)
ENDIF
!
@@ -130,7 +138,11 @@ CONTAINS
ENDIF
sss_m(:,:) = sss_m(:,:) + zts(:,:,jp_sal)
! ! removed inverse barometer ssh when Patm forcing is used (for sea-ice dynamics)
+#if defined key_RK3
+ IF( ln_apr_dyn ) THEN ; ssh_m(:,:) = ssh_m(:,:) + ssh(:,:,Kmm) - ssh_ib(:,:)
+#else
IF( ln_apr_dyn ) THEN ; ssh_m(:,:) = ssh_m(:,:) + ssh(:,:,Kmm) - 0.5 * ( ssh_ib(:,:) + ssh_ibb(:,:) )
+#endif
ELSE ; ssh_m(:,:) = ssh_m(:,:) + ssh(:,:,Kmm)
ENDIF
!
--
GitLab