From 77465cf6a52224417f58d9a7646afc33b7c7dea4 Mon Sep 17 00:00:00 2001
From: Sibylle Techene <sibylle.techene@locean-ipsl.upmc.fr>
Date: Mon, 13 Jun 2022 14:56:29 +0000
Subject: [PATCH] Resolve "Fix restartability issues for RK3+ISF"
---
src/OCE/ISF/isf_oce.F90 | 30 ++++++++++++++++++++++++------
src/OCE/ISF/isfcav.F90 | 9 +++++++--
src/OCE/ISF/isfhdiv.F90 | 29 ++++++++++++++++++++---------
src/OCE/ISF/isfpar.F90 | 8 ++++++--
src/OCE/ISF/isfstp.F90 | 6 ++++++
src/OCE/TRA/traisf.F90 | 36 +++++++++++++++++++++++++++---------
6 files changed, 90 insertions(+), 28 deletions(-)
diff --git a/src/OCE/ISF/isf_oce.F90 b/src/OCE/ISF/isf_oce.F90
index 16ff46027..eb62d83e9 100644
--- a/src/OCE/ISF/isf_oce.F90
+++ b/src/OCE/ISF/isf_oce.F90
@@ -80,8 +80,8 @@ MODULE isf_oce
!
! 2.2 -------- ice shelf cavity melt namelist parameter -------------
INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: mskisf_cav !:
- INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: misfkt_cav , misfkb_cav !:
- REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhisf_tbl_cav, rfrac_tbl_cav !:
+ INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: misfkt_cav , misfkb_cav !: shallowest and deepest level of the ice shelf
+ REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rhisf_tbl_cav, rfrac_tbl_cav !: thickness and fraction of deepest cell affected by the ice shelf
REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fwfisf_cav , fwfisf_cav_b !: before and now net fwf from the ice shelf [kg/m2/s]
REAL(wp) , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: risf_cav_tsc , risf_cav_tsc_b !: before and now T & S isf contents [K.m/s & PSU.m/s]
TYPE(FLD), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: sf_isfcav_fwf !:
@@ -237,17 +237,28 @@ CONTAINS
!
ierr = 0 ! set to zero if no array to be allocated
!
- ALLOCATE( fwfisf_par (jpi,jpj) , fwfisf_par_b(jpi,jpj) , &
- & fwfisf_cav (jpi,jpj) , fwfisf_cav_b(jpi,jpj) , &
+ ALLOCATE( fwfisf_par (jpi,jpj) , fwfisf_cav (jpi,jpj) , &
& fwfisf_oasis(jpi,jpj) , STAT=ialloc )
ierr = ierr + ialloc
!
- ALLOCATE( risf_par_tsc(jpi,jpj,jpts) , risf_par_tsc_b(jpi,jpj,jpts) , STAT=ialloc )
+ ALLOCATE( risf_par_tsc(jpi,jpj,jpts) , STAT=ialloc )
ierr = ierr + ialloc
!
- ALLOCATE( risf_cav_tsc(jpi,jpj,jpts) , risf_cav_tsc_b(jpi,jpj,jpts) , STAT=ialloc )
+ ALLOCATE( risf_cav_tsc(jpi,jpj,jpts) , STAT=ialloc )
ierr = ierr + ialloc
!
+#if ! defined key_RK3
+ ! MLF : need to allocate before arrays
+ ALLOCATE( fwfisf_par_b(jpi,jpj) , fwfisf_cav_b(jpi,jpj) , STAT=ialloc )
+ ierr = ierr + ialloc
+ !
+ ALLOCATE( risf_par_tsc_b(jpi,jpj,jpts) , STAT=ialloc )
+ ierr = ierr + ialloc
+ !
+ ALLOCATE( risf_cav_tsc_b(jpi,jpj,jpts) , STAT=ialloc )
+ ierr = ierr + ialloc
+#endif
+ !
ALLOCATE( risfload(jpi,jpj) , STAT=ialloc )
ierr = ierr + ialloc
!
@@ -260,10 +271,17 @@ CONTAINS
! initalisation of fwf and tsc array to 0
risfload (:,:) = 0._wp
fwfisf_oasis(:,:) = 0._wp
+#if defined key_RK3
+ fwfisf_par (:,:) = 0._wp
+ fwfisf_cav (:,:) = 0._wp
+ risf_cav_tsc(:,:,:) = 0._wp
+ risf_par_tsc(:,:,:) = 0._wp
+#else
fwfisf_par (:,:) = 0._wp ; fwfisf_par_b (:,:) = 0._wp
fwfisf_cav (:,:) = 0._wp ; fwfisf_cav_b (:,:) = 0._wp
risf_cav_tsc(:,:,:) = 0._wp ; risf_cav_tsc_b(:,:,:) = 0._wp
risf_par_tsc(:,:,:) = 0._wp ; risf_par_tsc_b(:,:,:) = 0._wp
+#endif
!
END SUBROUTINE isf_alloc
diff --git a/src/OCE/ISF/isfcav.F90 b/src/OCE/ISF/isfcav.F90
index 1a0597a83..27ca9912e 100644
--- a/src/OCE/ISF/isfcav.F90
+++ b/src/OCE/ISF/isfcav.F90
@@ -174,8 +174,10 @@ CONTAINS
! output fluxes
CALL isf_diags_flx( Kmm, misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, 'cav', pqfwf, zqoce, zqlat, zqhc)
!
- ! write restart variables (qoceisf, qhcisf, fwfisf for now and before)
+#if ! defined key_RK3
+ ! MLF: write restart variables (qoceisf, qhcisf, fwfisf for now and before)
IF (lrst_oce) CALL isfrst_write(kt, 'cav', ptsc, pqfwf)
+#endif
!
IF ( ln_isfdebug ) THEN
IF(lwp) WRITE(numout,*) ''
@@ -215,6 +217,8 @@ CONTAINS
!
! cavity mask
mskisf_cav(:,:) = (1._wp - tmask(:,:,1)) * ssmask(:,:)
+ !
+#if ! defined key_RK3
!================
! 2: activate restart
!================
@@ -223,9 +227,10 @@ CONTAINS
! 3: read restart
!================
!
- ! read cav variable from restart
+ ! MLF: read cav variable from restart
IF ( ln_rstart ) CALL isfrst_read('cav', risf_cav_tsc, fwfisf_cav, risf_cav_tsc_b, fwfisf_cav_b)
!
+#endif
!==========================================
! 3: specific allocation and initialisation (depending of scheme choice)
!==========================================
diff --git a/src/OCE/ISF/isfhdiv.F90 b/src/OCE/ISF/isfhdiv.F90
index 772ab9a74..58275063b 100644
--- a/src/OCE/ISF/isfhdiv.F90
+++ b/src/OCE/ISF/isfhdiv.F90
@@ -5,6 +5,7 @@ MODULE isfhdiv
!! with the ice shelf melt and coupling correction
!!======================================================================
!! History : 4.0 ! 2019-09 (P. Mathiot) Original code
+ !! 4.2 ! 2022-05 (S. Techene) Update wrt RK3 time-stepping
!!----------------------------------------------------------------------
!!----------------------------------------------------------------------
@@ -45,11 +46,19 @@ CONTAINS
!
IF ( ln_isf ) THEN
!
- ! ice shelf cavity contribution
- IF ( ln_isfcav_mlt ) CALL isf_hdiv_mlt(misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, fwfisf_cav, fwfisf_cav_b, phdiv)
+#if defined key_RK3
+ ! ice shelf cavity contribution (RK3)
+ IF ( ln_isfcav_mlt ) CALL isf_hdiv_mlt(misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, fwfisf_cav, phdiv)
+ !
+ ! ice shelf parametrisation contribution (RK3)
+ IF ( ln_isfpar_mlt ) CALL isf_hdiv_mlt(misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, fwfisf_par, phdiv)
+#else
+ ! ice shelf cavity contribution (MLF)
+ IF ( ln_isfcav_mlt ) CALL isf_hdiv_mlt(misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, fwfisf_cav, phdiv, fwfisf_cav_b)
!
- ! ice shelf parametrisation contribution
- IF ( ln_isfpar_mlt ) CALL isf_hdiv_mlt(misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, fwfisf_par, fwfisf_par_b, phdiv)
+ ! ice shelf parametrisation contribution (MLF)
+ IF ( ln_isfpar_mlt ) CALL isf_hdiv_mlt(misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, fwfisf_par, phdiv, fwfisf_par_b)
+#endif
!
! ice sheet coupling contribution
IF ( ln_isfcpl .AND. kt /= 0 ) THEN
@@ -72,7 +81,7 @@ CONTAINS
END SUBROUTINE isf_hdiv
- SUBROUTINE isf_hdiv_mlt(ktop, kbot, phtbl, pfrac, pfwf, pfwf_b, phdiv)
+ SUBROUTINE isf_hdiv_mlt(ktop, kbot, phtbl, pfrac, pfwf, phdiv, pfwf_b)
!!----------------------------------------------------------------------
!! *** SUBROUTINE sbc_isf_div ***
!!
@@ -83,11 +92,12 @@ CONTAINS
!!
!! ** Action : phdivn increased by the ice shelf outflow
!!----------------------------------------------------------------------
- REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: phdiv
+ REAL(wp), DIMENSION(jpi,jpj,jpk) , INTENT(inout) :: phdiv
!!----------------------------------------------------------------------
- INTEGER , DIMENSION(jpi,jpj), INTENT(in ) :: ktop , kbot
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pfrac, phtbl
- REAL(wp), DIMENSION(jpi,jpj), INTENT(in ) :: pfwf , pfwf_b
+ INTEGER , DIMENSION(jpi,jpj) , INTENT(in ) :: ktop , kbot
+ REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pfrac, phtbl
+ REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: pfwf
+ REAL(wp), DIMENSION(:,:) , OPTIONAL, INTENT(in ) :: pfwf_b
!!----------------------------------------------------------------------
INTEGER :: ji, jj, jk ! dummy loop indices
INTEGER :: ikt, ikb
@@ -119,6 +129,7 @@ CONTAINS
!
END SUBROUTINE isf_hdiv_mlt
+
SUBROUTINE isf_hdiv_cpl(Kmm, pqvol, phdiv)
!!----------------------------------------------------------------------
!! *** SUBROUTINE isf_hdiv_cpl ***
diff --git a/src/OCE/ISF/isfpar.F90 b/src/OCE/ISF/isfpar.F90
index 3dae222da..31f363687 100644
--- a/src/OCE/ISF/isfpar.F90
+++ b/src/OCE/ISF/isfpar.F90
@@ -93,9 +93,11 @@ CONTAINS
! output fluxes
CALL isf_diags_flx( Kmm, misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, 'par', pqfwf, zqoce, zqlat, zqhc)
!
- ! write restart variables (qoceisf, qhcisf, fwfisf for now and before)
+#if ! defined key_RK3
+ ! MLF: write restart variables (qoceisf, qhcisf, fwfisf for now and before)
IF (lrst_oce) CALL isfrst_write(kt, 'par', ptsc, pqfwf)
!
+#endif
IF ( ln_isfdebug ) THEN
IF(lwp) WRITE(numout,*)
CALL debug('isf_par: ptsc T',ptsc(:,:,1))
@@ -155,8 +157,10 @@ CONTAINS
mskisf_par(:,:) = 1._wp
END WHERE
!
- ! read par variable from restart
+#if ! defined key_RK3
+ ! MLF: read par variable from restart
IF ( ln_rstart ) CALL isfrst_read('par', risf_par_tsc, fwfisf_par, risf_par_tsc_b, fwfisf_par_b)
+#endif
!
SELECT CASE ( TRIM(cn_isfpar_mlt) )
!
diff --git a/src/OCE/ISF/isfstp.F90 b/src/OCE/ISF/isfstp.F90
index cef91bb94..edc5f92aa 100644
--- a/src/OCE/ISF/isfstp.F90
+++ b/src/OCE/ISF/isfstp.F90
@@ -74,11 +74,14 @@ CONTAINS
!
IF ( ln_isfcav_mlt ) THEN
!
+#if ! defined key_RK3
+ ! MLF : need risf_cav_tsc_b update
! 1.1: before time step
IF ( kt /= nit000 ) THEN
risf_cav_tsc_b (:,:,:) = risf_cav_tsc (:,:,:)
fwfisf_cav_b(:,:) = fwfisf_cav(:,:)
END IF
+#endif
!
! 1.2: compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl)
rhisf_tbl_cav(:,:) = rn_htbl * mskisf_cav(:,:)
@@ -102,11 +105,14 @@ CONTAINS
!
IF ( ln_isfpar_mlt ) THEN
!
+#if ! defined key_RK3
+ ! MLF : need risf_par_tsc_b update
! 2.1: before time step
IF ( kt /= nit000 ) THEN
risf_par_tsc_b(:,:,:) = risf_par_tsc(:,:,:)
fwfisf_par_b (:,:) = fwfisf_par (:,:)
END IF
+#endif
!
! 2.2: compute misfkb, rhisf_tbl, rfrac (deepest level, thickness, fraction of deepest cell affected by tbl)
! by simplicity, we assume the top level where param applied do not change with time (done in init part)
diff --git a/src/OCE/TRA/traisf.F90 b/src/OCE/TRA/traisf.F90
index 7f43285b9..610f9064f 100644
--- a/src/OCE/TRA/traisf.F90
+++ b/src/OCE/TRA/traisf.F90
@@ -17,6 +17,8 @@ MODULE traisf
USE isfutils, ONLY : debug ! debug option
USE timing , ONLY : timing_start, timing_stop ! Timing
USE in_out_manager ! I/O manager
+ !
+ USE prtctl ! print control
IMPLICIT NONE
PRIVATE
@@ -56,11 +58,19 @@ CONTAINS
ENDIF
ENDIF
!
- ! cavity case
- IF ( ln_isfcav_mlt ) CALL isf_mlt(misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, risf_cav_tsc, risf_cav_tsc_b, pts(:,:,:,:,Krhs))
+#if defined key_RK3
+ ! cavity case (RK3)
+ IF ( ln_isfcav_mlt ) CALL isf_mlt(misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, risf_cav_tsc, pts(:,:,:,:,Krhs))
!
- ! parametrisation case
- IF ( ln_isfpar_mlt ) CALL isf_mlt(misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, risf_par_tsc, risf_par_tsc_b, pts(:,:,:,:,Krhs))
+ ! parametrisation case (RK3)
+ IF ( ln_isfpar_mlt ) CALL isf_mlt(misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, risf_par_tsc, pts(:,:,:,:,Krhs))
+#else
+ ! cavity case (MLF)
+ IF ( ln_isfcav_mlt ) CALL isf_mlt(misfkt_cav, misfkb_cav, rhisf_tbl_cav, rfrac_tbl_cav, risf_cav_tsc, pts(:,:,:,:,Krhs), risf_cav_tsc_b)
+ !
+ ! parametrisation case (MLF)
+ IF ( ln_isfpar_mlt ) CALL isf_mlt(misfkt_par, misfkb_par, rhisf_tbl_par, rfrac_tbl_par, risf_par_tsc, pts(:,:,:,:,Krhs), risf_par_tsc_b)
+#endif
!
! ice sheet coupling case
IF ( ln_isfcpl ) THEN
@@ -86,12 +96,15 @@ CONTAINS
ENDIF
END IF
!
+ IF(sn_cfctl%l_prtctl) CALL prt_ctl( tab3d_1=pts(:,:,:,jp_tem,Krhs), clinfo1=' isf - Ta: ', mask1=tmask, &
+ & tab3d_2=pts(:,:,:,jp_sal,Krhs), clinfo2= ' Sa: ', mask2=tmask, clinfo3='tra' )
+ !
IF( ln_timing ) CALL timing_stop('tra_isf')
!
END SUBROUTINE tra_isf
- SUBROUTINE isf_mlt( ktop, kbot, phtbl, pfrac, ptsc, ptsc_b, pts )
+ SUBROUTINE isf_mlt( ktop, kbot, phtbl, pfrac, ptsc, pts, ptsc_b )
!!----------------------------------------------------------------------
!! *** ROUTINE isf_mlt ***
!!
@@ -100,10 +113,11 @@ CONTAINS
!! *** Action :: Update pts(:,:,:,:,Krhs) with the surface boundary condition trend
!!
!!----------------------------------------------------------------------
- INTEGER , DIMENSION(jpi,jpj) , INTENT(in ) :: ktop , kbot
- REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: phtbl, pfrac
- REAL(wp), DIMENSION(jpi,jpj,jpts) , INTENT(in ) :: ptsc , ptsc_b
- REAL(wp), DIMENSION(jpi,jpj,jpk,jpts), INTENT(inout) :: pts
+ REAL(wp), DIMENSION(jpi,jpj,jpk,jpts) , INTENT(inout) :: pts
+ INTEGER , DIMENSION(jpi,jpj) , INTENT(in ) :: ktop , kbot
+ REAL(wp), DIMENSION(jpi,jpj) , INTENT(in ) :: phtbl, pfrac
+ REAL(wp), DIMENSION(jpi,jpj,jpts) , INTENT(in ) :: ptsc
+ REAL(wp), DIMENSION(:,:,:) , OPTIONAL, INTENT(in ) :: ptsc_b
!!
INTEGER :: ji,jj,jk ! dummy loop index
INTEGER :: ikt, ikb ! top and bottom level of the tbl
@@ -112,7 +126,11 @@ CONTAINS
!
! compute 2d total trend due to isf
DO_2D( 0, 0, 0, 0 )
+#if defined key_RK3
+ ztc(ji,jj) = ptsc(ji,jj,jp_tem) / phtbl(ji,jj)
+#else
ztc(ji,jj) = 0.5_wp * ( ptsc(ji,jj,jp_tem) + ptsc_b(ji,jj,jp_tem) ) / phtbl(ji,jj)
+#endif
END_2D
!
! update pts(:,:,:,:,Krhs)
--
GitLab