Skip to content
Snippets Groups Projects
Commit b2ec917b authored by Sebastien Masson's avatar Sebastien Masson
Browse files

Merge branch '50-correct-restartability-issues-for-rk3-isf' into 'main'

Resolve "Fix restartability issues for RK3+ISF"

Closes #50

See merge request nemo/nemo!85
parents 86aacda3 77465cf6
No related branches found
No related tags found
No related merge requests found
......@@ -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
......
......@@ -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)
!==========================================
......
......@@ -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 ***
......
......@@ -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) )
!
......
......@@ -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)
......
......@@ -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)
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment