From d1e52b888584fe31ea76addf7ca9edb0a47e47a8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=A9r=C3=B4me=20Chanut?=
<jerome.chanut@mercator-ocean.fr>
Date: Thu, 13 Apr 2023 09:25:40 +0000
Subject: [PATCH] Resolve "Bounds issue in AGRIF sponges preventing from
identical rotationally symmetric results"
---
src/NST/agrif_oce_sponge.F90 | 15 ++++++++-------
src/NST/agrif_oce_update.F90 | 3 +--
src/NST/agrif_top_sponge.F90 | 12 ++++++------
src/NST/agrif_user.F90 | 7 ++++---
4 files changed, 19 insertions(+), 18 deletions(-)
diff --git a/src/NST/agrif_oce_sponge.F90 b/src/NST/agrif_oce_sponge.F90
index fe2997c0..ce2d6575 100644
--- a/src/NST/agrif_oce_sponge.F90
+++ b/src/NST/agrif_oce_sponge.F90
@@ -395,7 +395,7 @@ CONTAINS
INTEGER :: iku, ikv
REAL(wp) :: ztsa, zabe1, zabe2, zbtr, zhtot
REAl(wp) :: zflag, zdmod, zdtot
- REAL(wp), DIMENSION(i1-1:i2,j1-1:j2,jpk) :: ztu, ztv
+ REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv
REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tsbdiff
! vertical interpolation:
REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child
@@ -551,17 +551,18 @@ CONTAINS
DO jn = 1, jpts
DO jk = 1, jpkm1
- ztu(i1-1:i2,j1-1:j2,jk) = 0._wp
+ ztu(i2,j1:j2,jk) = 0._wp
DO jj = j1,j2
DO ji = i1,i2-1
zabe1 = rn_sponge_tra * r1_Dt * umask(ji,jj,jk) * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
zdtot = tsbdiff(ji+1,jj,jk,jn) - tsbdiff(ji,jj,jk,jn)
zdmod = ts(ji+1,jj,jk,jn,Kbb_a) - ts(ji,jj,jk,jn,Kbb_a)
zflag = 0.5_wp + SIGN(0.5_wp, zdtot*zdmod)
+ IF (ABS(zdtot*zdmod)<=1.e-6) zflag = 1._wp ! for 1:1 repro
ztu(ji,jj,jk) = zabe1 * fspu(ji,jj) * ( zflag * zdtot + (1._wp - zflag) * zdmod )
END DO
END DO
- ztv(i1-1:i2,j1-1:j2,jk) = 0._wp
+ ztv(i1:i2,j2,jk) = 0._wp
DO ji = i1,i2
DO jj = j1,j2-1
zabe2 = rn_sponge_tra * r1_Dt * vmask(ji,jj,jk) * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
@@ -569,6 +570,7 @@ CONTAINS
zdtot = tsbdiff(ji,jj+1,jk,jn) - tsbdiff(ji,jj,jk,jn)
zdmod = ts(ji,jj+1,jk,jn,Kbb_a) - ts(ji,jj,jk,jn,Kbb_a)
zflag = 0.5_wp + SIGN(0.5_wp, zdtot*zdmod)
+ IF (ABS(zdtot*zdmod)<=1.e-6) zflag = 1._wp ! for 1:1 repro
ztv(ji,jj,jk) = zabe2 * fspv(ji,jj) * ( zflag * zdtot + (1._wp - zflag) * zdmod )
END DO
END DO
@@ -586,10 +588,9 @@ CONTAINS
ENDIF
END DO
!
-! JC: there is something wrong with the Laplacian in corners
DO jk = 1, jpkm1
- DO jj = j1,j2
- DO ji = i1,i2
+ DO jj = j1+1,j2-1
+ DO ji = i1+1,i2-1
IF (.NOT. tabspongedone_tsn(ji,jj)) THEN
zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a)
! horizontal diffusive trends
@@ -606,7 +607,7 @@ CONTAINS
!
END DO
!
- tabspongedone_tsn(i1:i2,j1:j2) = .TRUE.
+ tabspongedone_tsn(i1+1:i2-1,j1+1:j2-1) = .TRUE.
!
ENDIF
!
diff --git a/src/NST/agrif_oce_update.F90 b/src/NST/agrif_oce_update.F90
index e48bd2df..8e9142fd 100644
--- a/src/NST/agrif_oce_update.F90
+++ b/src/NST/agrif_oce_update.F90
@@ -202,15 +202,14 @@ CONTAINS
Agrif_UseSpecialValueInUpdate = .FALSE.
# if ! defined DECAL_FEEDBACK_2D
CALL Agrif_Update_Variable(r3t_id, locupdate=(/ nn_shift_bar,-2/), procname=update_r3t)
- CALL Agrif_Update_Variable(r3f_id, locupdate=(/ nn_shift_bar,-2/), procname=update_r3f)
CALL Agrif_Update_Variable(r3u_id, locupdate1=(/ nn_shift_bar,-2/), locupdate2=(/ nn_shift_bar,-2/), procname=update_r3u)
CALL Agrif_Update_Variable(r3v_id, locupdate1=(/ nn_shift_bar,-2/), locupdate2=(/ nn_shift_bar,-2/), procname=update_r3v)
# else
CALL Agrif_Update_Variable(r3t_id, locupdate=(/1+nn_shift_bar,-2/), procname=update_r3t)
- CALL Agrif_Update_Variable(r3f_id, locupdate=(/1+nn_shift_bar,-2/), procname=update_r3f)
CALL Agrif_Update_Variable(r3u_id, locupdate1=(/ nn_shift_bar,-2/), locupdate2=(/1+nn_shift_bar,-2/), procname=update_r3u)
CALL Agrif_Update_Variable(r3v_id, locupdate1=(/1+nn_shift_bar,-2/), locupdate2=(/ nn_shift_bar,-2/), procname=update_r3v)
# endif
+ CALL Agrif_Update_Variable(r3f_id, locupdate=(/1+nn_shift_bar,-2/), procname=update_r3f)
!
! Old way (update e3 at UVF-points everywhere on parent domain):
! CALL Agrif_ChildGrid_To_ParentGrid()
diff --git a/src/NST/agrif_top_sponge.F90 b/src/NST/agrif_top_sponge.F90
index a3706693..1ef72628 100644
--- a/src/NST/agrif_top_sponge.F90
+++ b/src/NST/agrif_top_sponge.F90
@@ -77,7 +77,7 @@ CONTAINS
INTEGER :: ji, jj, jk, jn ! dummy loop indices
INTEGER :: iku, ikv
REAL(wp) :: ztra, zabe1, zabe2, zbtr, zhtot
- REAL(wp), DIMENSION(i1-1:i2,j1-1:j2,jpk) :: ztu, ztv
+ REAL(wp), DIMENSION(i1:i2,j1:j2,jpk) :: ztu, ztv
REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::trbdiff
! vertical interpolation:
REAL(wp), DIMENSION(i1:i2,j1:j2,jpk,n1:n2) ::tabres_child
@@ -230,14 +230,14 @@ CONTAINS
DO jn = 1, jptra
DO jk = 1, jpkm1
- ztu(i1-1:i2,j1-1:j2,jk) = 0._wp
+ ztu(i2,j1:j2,jk) = 0._wp
DO jj = j1,j2
DO ji = i1,i2-1
zabe1 = rn_sponge_tra * r1_Dt * umask(ji,jj,jk) * e1e2u(ji,jj) * e3u(ji,jj,jk,Kmm_a)
ztu(ji,jj,jk) = zabe1 * fspu(ji,jj) * ( trbdiff(ji+1,jj ,jk,jn) - trbdiff(ji,jj,jk,jn) )
END DO
END DO
- ztv(i1-1:i2,j1-1:j2,jk) = 0._wp
+ ztv(i1:i2,j2,jk) = 0._wp
DO ji = i1,i2
DO jj = j1,j2-1
zabe2 = rn_sponge_tra * r1_Dt * vmask(ji,jj,jk) * e1e2v(ji,jj) * e3v(ji,jj,jk,Kmm_a)
@@ -260,8 +260,8 @@ CONTAINS
!
! JC: there is something wrong with the Laplacian in corners
DO jk = 1, jpkm1
- DO jj = j1,j2
- DO ji = i1,i2
+ DO jj = j1+1,j2-1
+ DO ji = i1+1,i2-1
IF (.NOT. tabspongedone_trn(ji,jj)) THEN
zbtr = r1_e1e2t(ji,jj) / e3t(ji,jj,jk,Kmm_a)
! horizontal diffusive trends
@@ -278,7 +278,7 @@ CONTAINS
!
END DO
!
- tabspongedone_trn(i1:i2,j1:j2) = .TRUE.
+ tabspongedone_trn(i1+1:i2-1,j1+1:j2-1) = .TRUE.
!
ENDIF
!
diff --git a/src/NST/agrif_user.F90 b/src/NST/agrif_user.F90
index 5f547049..5d4a83b8 100644
--- a/src/NST/agrif_user.F90
+++ b/src/NST/agrif_user.F90
@@ -592,9 +592,10 @@
CALL Agrif_Set_bc( un_interp_id, (/0,ind1-1/) )
CALL Agrif_Set_bc( vn_interp_id, (/0,ind1-1/) )
- CALL Agrif_Set_bc( ts_sponge_id, (/-nn_sponge_len*imaxrho-1,0/) ) ! if west, rhox=3, nn_sponge_len=2
- CALL Agrif_Set_bc( un_sponge_id, (/-nn_sponge_len*imaxrho-1,0/) ) ! and nbghost=3:
- CALL Agrif_Set_bc( vn_sponge_id, (/-nn_sponge_len*imaxrho-1,0/) ) ! columns 4 to 11
+ ! if west, rhox=3, nn_sponge_len=2 and and nbghost=3:
+ CALL Agrif_Set_bc( ts_sponge_id, (/-nn_sponge_len*imaxrho-1,0/) ) ! columns 4-11
+ CALL Agrif_Set_bc( un_sponge_id, (/-nn_sponge_len*imaxrho-2,0/) ) ! columns 4-12
+ CALL Agrif_Set_bc( vn_sponge_id, (/-nn_sponge_len*imaxrho-2,0/) ) ! columns 4-12
CALL Agrif_Set_bc( sshn_id, (/-imaxrho*nn_shift_bar,ind1-1/) )
CALL Agrif_Set_bc( sshn_frc_id, (/-imaxrho*nn_shift_bar,ind1-1/) )
--
GitLab