From b60262ce33cdbfc1b87f47bb6b5e48b08433fefe Mon Sep 17 00:00:00 2001
From: Sibylle TECHENE <techenes@irene191.c-irene.tgcc.ccc.cea.fr>
Date: Thu, 24 Mar 2022 14:36:05 +0100
Subject: [PATCH] Squashed commit of the following:

commit 87d91638e9b25f7b45192068b43db371c26a8ed6
Author: Sibylle TECHENE <techenes@irene191.c-irene.tgcc.ccc.cea.fr>
Date:   Thu Mar 24 14:31:02 2022 +0100

    comments are still not clear

commit 3e306fb89db460675f1a18e15f1d8910c8b9dded
Author: Sibylle TECHENE <techenes@irene191.c-irene.tgcc.ccc.cea.fr>
Date:   Thu Mar 24 14:24:56 2022 +0100

    RK3: fix linssh

commit 92d543ecce2c9b401b31c8dfb17c75d1a33fc0ac
Author: Sibylle TECHENE <techenes@irene191.c-irene.tgcc.ccc.cea.fr>
Date:   Thu Mar 24 14:24:07 2022 +0100

    RK3: management of u/v_bf, u/vb2_b and u/vb2_i_b allocation

commit 7381d7650a1c356406078a86e12c5d462afdc5cb
Author: Sibylle TECHENE <techenes@irene191.c-irene.tgcc.ccc.cea.fr>
Date:   Thu Mar 24 14:21:02 2022 +0100

    RK3: no use of before forcing fields except for SWE (u/vtau_b)

commit 1d53bab023d5ebfcf2f3e4e5b7d3a835a0795890
Author: Sibylle TECHENE <techenes@irene191.c-irene.tgcc.ccc.cea.fr>
Date:   Thu Mar 24 14:18:32 2022 +0100

    RK3: cleanning

commit 683b0c199ee34e621d27de16e688c4c81fb6b0b0
Author: Sibylle TECHENE <techenes@irene191.c-irene.tgcc.ccc.cea.fr>
Date:   Thu Mar 24 14:17:36 2022 +0100

    RK3: no use of ts_rst, ts_rst turned into ts_rst_MLF
---
 src/OCE/DYN/dynspg_ts.F90 | 55 ++++++++++++++++++++++++++-------------
 src/OCE/IOM/restart.F90   |  9 -------
 src/OCE/SBC/sbcmod.F90    | 16 ++++++++++--
 src/OCE/nemogcm.F90       | 12 ++++-----
 src/OCE/oce.F90           | 32 +++++++++++++++++------
 src/OCE/stp2d.F90         |  8 +++++-
 6 files changed, 88 insertions(+), 44 deletions(-)

diff --git a/src/OCE/DYN/dynspg_ts.F90 b/src/OCE/DYN/dynspg_ts.F90
index b47a85a7..efef39e9 100644
--- a/src/OCE/DYN/dynspg_ts.F90
+++ b/src/OCE/DYN/dynspg_ts.F90
@@ -349,9 +349,9 @@ CONTAINS
          END_2D
       ENDIF  
       !
-      !              !----------------!
-      !              !==  sssh_frc  ==!   Right-Hand-Side of the barotropic ssh equation   (over the FULL domain)
-      !              !----------------!
+      !              !---------------!
+      !              !==  ssh_frc  ==!   Right-Hand-Side of the barotropic ssh equation   (over the FULL domain)
+      !              !---------------!
       !                                   !=  Net water flux forcing applied to a water column  =!
       !                                   ! ---------------------------------------------------  !
       IF (ln_bt_fw) THEN                          ! FORWARD integration: use kt+1/2 fluxes (NOW+1/2)
@@ -424,7 +424,8 @@ CONTAINS
          zhtp2_e(:,:) = ht_0(:,:)
       ENDIF
       !
-      IF( ln_bt_fw ) THEN                 ! FORWARD integration: start from NOW fields                    
+      IF( ln_bt_fw ) THEN                 ! FORWARD integration: start from NOW fields
+         !                                   ! RK3: Kmm = Kbb when calling dynspg_ts
          sshn_e(:,:) =    pssh (:,:,Kmm)            
          un_e  (:,:) =    puu_b(:,:,Kmm)            
          vn_e  (:,:) =    pvv_b(:,:,Kmm)
@@ -821,10 +822,6 @@ CONTAINS
          !
       ENDIF
       !
-      ! advective transport from N to N+1 (i.e. Kbb to Kaa) 
-      ub2_b(:,:) = un_adv(:,:)      ! Save integrated transport for next computation (NOT USED)
-      vb2_b(:,:) = vn_adv(:,:)
-      ! 
       IF( iom_use("ubar") ) THEN    ! RK3 single first: hu[N+1/2] = 1/2 ( hu[N] + hu[N+1] ) 
          ALLOCATE( z2d(jpi,jpj) )
          z2d(:,:) = 2._wp / ( hu_e(:,:) + hu(:,:,Kbb) + 1._wp - ssumask(:,:) ) 
@@ -933,6 +930,12 @@ CONTAINS
 
       !
 #if defined key_agrif
+# if defined key_RK3
+      ! advective transport from N to N+1 (i.e. Kbb to Kaa) 
+      ub2_b(:,:) = un_adv(:,:)      ! Save integrated transport for AGRIF
+      vb2_b(:,:) = vn_adv(:,:)
+# endif
+      !
       ! Save time integrated fluxes during child grid integration
       ! (used to update coarse grid transports at next time step)
       !
@@ -946,9 +949,11 @@ CONTAINS
          ub2_i_b(:,:) = ub2_i_b(:,:) + za1 * ub2_b(:,:)
          vb2_i_b(:,:) = vb2_i_b(:,:) + za1 * vb2_b(:,:)
       ENDIF
-#endif      
-      !                                   !* write time-spliting arrays in the restart
-      IF( lrst_oce .AND.ln_bt_fw )   CALL ts_rst( kt, 'WRITE' )
+#endif
+#if ! defined key_RK3
+      !                                   !* MLF: write time-spliting arrays in the restart
+      IF( lrst_oce .AND.ln_bt_fw )   CALL ts_rst_MLF( kt, 'WRITE' )
+#endif
       !
       IF( ln_wd_il )   DEALLOCATE( zcpx, zcpy )
       IF( ln_wd_dl )   DEALLOCATE( ztwdmask, zuwdmask, zvwdmask, zuwdav2, zvwdav2 )
@@ -1037,9 +1042,9 @@ CONTAINS
    END SUBROUTINE ts_wgt
 
 
-   SUBROUTINE ts_rst( kt, cdrw )
+   SUBROUTINE ts_rst_MLF( kt, cdrw )
       !!---------------------------------------------------------------------
-      !!                   ***  ROUTINE ts_rst  ***
+      !!                   ***  ROUTINE ts_rst_MLF  ***
       !!
       !! ** Purpose : Read or write time-splitting arrays in restart file
       !!----------------------------------------------------------------------
@@ -1071,12 +1076,14 @@ CONTAINS
                ub2_i_b(:,:) = 0._wp   ;   vb2_i_b(:,:) = 0._wp   ! used in the 1st update of agrif
             ENDIF
 #endif
-         ELSE                                   !* Start from rest
+         ELSE
+            !                      !* Start from rest or use RK3 time-step
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set barotropic values to 0'
             ub2_b  (:,:) = 0._wp   ;   vb2_b  (:,:) = 0._wp   ! used in the 1st interpol of agrif
             un_adv (:,:) = 0._wp   ;   vn_adv (:,:) = 0._wp   ! used in the 1st interpol of agrif
             un_bf  (:,:) = 0._wp   ;   vn_bf  (:,:) = 0._wp   ! used in the 1st update   of agrif
+
 #if defined key_agrif
             ub2_i_b(:,:) = 0._wp   ;   vb2_i_b(:,:) = 0._wp   ! used in the 1st update of agrif
 #endif
@@ -1084,7 +1091,7 @@ CONTAINS
          !
       ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
          !                                   ! -------------------
-         IF(lwp) WRITE(numout,*) '---- ts_rst ----'
+         IF(lwp) WRITE(numout,*) '---- ts_rst_MLF ----'
          CALL iom_rstput( kt, nitrst, numrow, 'ub2_b'   , ub2_b  (:,:) )
          CALL iom_rstput( kt, nitrst, numrow, 'vb2_b'   , vb2_b  (:,:) )
          CALL iom_rstput( kt, nitrst, numrow, 'un_bf'   , un_bf  (:,:) )
@@ -1107,7 +1114,7 @@ CONTAINS
 #endif
       ENDIF
       !
-   END SUBROUTINE ts_rst
+   END SUBROUTINE ts_rst_MLF
 
 
    SUBROUTINE dyn_spg_ts_init
@@ -1193,8 +1200,20 @@ CONTAINS
       !                             ! Allocate time-splitting arrays
       IF( dyn_spg_ts_alloc() /= 0    )   CALL ctl_stop('STOP', 'dyn_spg_init: failed to allocate dynspg_ts  arrays' )
       !
-      !                             ! read restart when needed
-      CALL ts_rst( nit000, 'READ' )
+#if defined key_RK3
+      !                      !* RK3: no restart
+# if defined key_agrif
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) '   ==>>>   RK3: set barotropic values to 0'
+      ub2_b  (:,:) = 0._wp   ;   vb2_b  (:,:) = 0._wp   ! used in the 1st interpol of agrif
+      un_adv (:,:) = 0._wp   ;   vn_adv (:,:) = 0._wp   ! used in the 1st interpol of agrif
+      un_bf  (:,:) = 0._wp   ;   vn_bf  (:,:) = 0._wp   ! used in the 1st update   of agrif
+# endif
+#else
+      !                      !* MLF: restart/initialise
+      !
+      CALL ts_rst_MLF( nit000, 'READ' )
+#endif
       !
    END SUBROUTINE dyn_spg_ts_init
 
diff --git a/src/OCE/IOM/restart.F90 b/src/OCE/IOM/restart.F90
index 0d93b821..01fcd453 100644
--- a/src/OCE/IOM/restart.F90
+++ b/src/OCE/IOM/restart.F90
@@ -174,7 +174,6 @@ CONTAINS
 #if defined key_RK3
          CALL iom_rstput( kt, nitrst, numrow, 'uu_b'   , uu_b(:,:       ,Kbb) )     ! before fields
          CALL iom_rstput( kt, nitrst, numrow, 'vv_b'   , vv_b(:,:       ,Kbb) )     ! before fields
-!!st-10/03         CALL iom_rstput( kt, nitrst, numrow, 'ssha '  , ssh(:,:        ,Kaa) )     ! after field post swap (n-1)
 #else
          CALL iom_rstput( kt, nitrst, numrow, 'sshn', ssh(:,:        ,Kmm) )     ! now fields : n
          CALL iom_rstput( kt, nitrst, numrow, 'un'  , uu(:,:,:       ,Kmm) )
@@ -380,8 +379,6 @@ CONTAINS
          !                                     !*  RK3: Set ssh at Kmm for AGRIF
          ssh(:,:,Kmm) = ssh(:,:,Kbb)
          !
-!!st-10/03         !                                     !*  RK3: Set ssh at Kaa (n-1) for ww computation
-!!st         CALL iom_get( numror, jpdom_auto, 'ssha'   , ssh(:,:,Kaa) )
 #else
          !                                     !*  MLF: Read ssh at Kmm
          IF(lwp) WRITE(numout,*)
@@ -438,15 +435,9 @@ CONTAINS
          ssh(:,:,Kmm) = ssh(:,:,Kbb)              !* set now values from to before ones
       ENDIF
       !
-#if defined key_RK3
-      IF(.NOT. ln_rstart ) THEN
-#endif
       !                            !==========================!
       ssh(:,:,Kaa) = 0._wp         !==  Set to 0 for AGRIF  ==!
       !                            !==========================!
-#if defined key_RK3
-      ENDIF
-#endif
       !
    END SUBROUTINE rst_read_ssh
 
diff --git a/src/OCE/SBC/sbcmod.F90 b/src/OCE/SBC/sbcmod.F90
index 12fcf66d..22069336 100644
--- a/src/OCE/SBC/sbcmod.F90
+++ b/src/OCE/SBC/sbcmod.F90
@@ -530,7 +530,11 @@ CONTAINS
       !
       IF( kt == nit000 ) THEN                          !   set the forcing field at nit000 - 1    !
          !                                             ! ---------------------------------------- !
-         IF( ln_rstart .AND. .NOT.l_1st_euler ) THEN            !* Restart: read in restart file
+#if defined key_RK3
+         IF( ln_rstart .AND. lk_SWE ) THEN                      !* RK3 + SWE: Restart: read in restart file
+#else
+         IF( ln_rstart .AND. .NOT.l_1st_euler ) THEN            !* MLF: Restart: read in restart file
+#endif
             IF(lwp) WRITE(numout,*) '          nit000-1 surface forcing fields read in the restart file'
             CALL iom_get( numror, jpdom_auto, 'utau_b', utau_b )   ! i-stress
             CALL iom_get( numror, jpdom_auto, 'vtau_b', vtau_b )   ! j-stress
@@ -552,9 +556,17 @@ CONTAINS
             sfx_b (:,:) = sfx (:,:)
          ENDIF
       ENDIF
+      !
+#if defined key_RK3
+      !                                                ! ---------------------------------------- !
+      IF( lrst_oce .AND. lk_SWE ) THEN                 !   RK3: Write in the ocean restart file   !
+         !                                             ! ---------------------------------------- !
+#else
       !                                                ! ---------------------------------------- !
-      IF( lrst_oce ) THEN                              !      Write in the ocean restart file     !
+      IF( lrst_oce ) THEN                              !   MLF: Write in the ocean restart file   !
          !                                             ! ---------------------------------------- !
+#endif
+         !
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'sbc : ocean surface forcing fields written in ocean restart file ',   &
             &                    'at it= ', kt,' date= ', ndastp
diff --git a/src/OCE/nemogcm.F90 b/src/OCE/nemogcm.F90
index deab2b74..9e5ebf36 100644
--- a/src/OCE/nemogcm.F90
+++ b/src/OCE/nemogcm.F90
@@ -126,9 +126,9 @@ CONTAINS
       !                            !-----------------------!
 #if defined key_agrif
 # if defined key_RK3
-      Kbb_a = Nbb; Kmm_a = Nbb; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices
+      Kbb_a = Nbb; Kmm_a = Nbb; Krhs_a = Nrhs   ! RK3: agrif_oce module copies of time level indices
 # else
-      Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices
+      Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! MLF: agrif_oce module copies of time level indices
 # endif
       CALL Agrif_Declare_Var       !  "      "   "   "      "  DYN/TRA
 # if defined key_top
@@ -155,9 +155,9 @@ CONTAINS
       !
       ! Recursive update from highest nested level to lowest:
 # if defined key_RK3
-      Kbb_a = Nbb; Kmm_a = Nbb; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices
+      Kbb_a = Nbb; Kmm_a = Nbb; Krhs_a = Nrhs   ! RK3: agrif_oce module copies of time level indices
 # else
-      Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! agrif_oce module copies of time level indices
+      Kbb_a = Nbb; Kmm_a = Nnn; Krhs_a = Nrhs   ! MLF: agrif_oce module copies of time level indices
 # endif
       CALL Agrif_step_child_adj(Agrif_Update_All)
       CALL Agrif_step_child_adj(Agrif_Check_parent_bat)
@@ -412,9 +412,9 @@ CONTAINS
       Nbb = 1   ;   Nnn = 2   ;   Naa = 3   ;   Nrhs = Naa
 #if defined key_agrif
 # if defined key_RK3
-      Kbb_a = Nbb   ;   Kmm_a = Nbb   ;   Krhs_a = Nrhs   ! agrif_oce module copies of time level indices
+      Kbb_a = Nbb   ;   Kmm_a = Nbb   ;   Krhs_a = Nrhs   ! RK3: agrif_oce module copies of time level indices
 # else
-      Kbb_a = Nbb   ;   Kmm_a = Nnn   ;   Krhs_a = Nrhs   ! agrif_oce module copies of time level indices
+      Kbb_a = Nbb   ;   Kmm_a = Nnn   ;   Krhs_a = Nrhs   ! MLF: agrif_oce module copies of time level indices
 # endif
 #endif
       !                             !-------------------------------!
diff --git a/src/OCE/oce.F90 b/src/OCE/oce.F90
index 8feef4e0..8f5e2d45 100644
--- a/src/OCE/oce.F90
+++ b/src/OCE/oce.F90
@@ -85,35 +85,51 @@ CONTAINS
       !!----------------------------------------------------------------------
       !!                   ***  FUNCTION oce_alloc  ***
       !!----------------------------------------------------------------------
-      INTEGER :: ierr(6)
+      INTEGER               ::   ii
+      INTEGER, DIMENSION(7) :: ierr
       !!----------------------------------------------------------------------
       !
-      ierr(:) = 0 
+      ii = 0   ;   ierr(:) = 0 
+      !
+      ii = ii+1
       ALLOCATE( uu   (jpi,jpj,jpk,jpt)  , vv   (jpi,jpj,jpk,jpt)           ,     &          
          &      ww   (jpi,jpj,jpk)      , hdiv(jpi,jpj,jpk)                ,     &
          &      ts   (jpi,jpj,jpk,jpts,jpt)                                ,     &
          &      rab_b(jpi,jpj,jpk,jpts) , rab_n(jpi,jpj,jpk,jpts)          ,     &
          &      rn2b (jpi,jpj,jpk)      , rn2  (jpi,jpj,jpk)               ,     &
-         &      rhd  (jpi,jpj,jpk)      , rhop (jpi,jpj,jpk)               , STAT=ierr(1) )
+         &      rhd  (jpi,jpj,jpk)      , rhop (jpi,jpj,jpk)               , STAT=ierr(ii) )
          !
+      ii = ii+1
       ALLOCATE( ssh (jpi,jpj,jpt)  , uu_b(jpi,jpj,jpt) , vv_b(jpi,jpj,jpt) ,     &
          &      ssh_frc(jpi,jpj)                                           ,     &
          &      gtsu(jpi,jpj,jpts) , gtsv(jpi,jpj,jpts)                    ,     &
          &      gru (jpi,jpj)      , grv (jpi,jpj)                         ,     &
          &      gtui(jpi,jpj,jpts) , gtvi(jpi,jpj,jpts)                    ,     &
          &      grui(jpi,jpj)      , grvi(jpi,jpj)                         ,     &
-         &      riceload(jpi,jpj)                                          , STAT=ierr(2) )
+         &      riceload(jpi,jpj)                                          , STAT=ierr(ii) )
          !
-      ALLOCATE( fraqsr_1lev(jpi,jpj) , STAT=ierr(3) )
+      ii = ii+1
+      ALLOCATE( fraqsr_1lev(jpi,jpj)                                       , STAT=ierr(ii) )
          !
+      ii = ii+1
       ALLOCATE( ssha_e(jpi,jpj),  sshn_e(jpi,jpj), sshb_e(jpi,jpj), sshbb_e(jpi,jpj), &
          &        ua_e(jpi,jpj),    un_e(jpi,jpj),   ub_e(jpi,jpj),   ubb_e(jpi,jpj), &
          &        va_e(jpi,jpj),    vn_e(jpi,jpj),   vb_e(jpi,jpj),   vbb_e(jpi,jpj), &
-         &        hu_e(jpi,jpj),   hur_e(jpi,jpj),   hv_e(jpi,jpj),   hvr_e(jpi,jpj), STAT=ierr(4) )
+         &        hu_e(jpi,jpj),   hur_e(jpi,jpj),   hv_e(jpi,jpj),   hvr_e(jpi,jpj), STAT=ierr(ii) )
          !
-      ALLOCATE( ub2_b(jpi,jpj), vb2_b(jpi,jpj), un_bf(jpi,jpj), vn_bf(jpi,jpj)      , STAT=ierr(5) )
+#if defined key_RK3
+# if defined key_agrif
+      ii = ii+1                             ! RK3+AGRIF: ???
+      ALLOCATE( ub2_b(jpi,jpj), vb2_b(jpi,jpj)                                      , STAT=ierr(ii) )
+# endif
+#else
+      ii = ii+1                             ! MLF: arrays related to Asselin filter + ???
+      ALLOCATE( un_bf(jpi,jpj), vn_bf(jpi,jpj), ub2_b(jpi,jpj), vb2_b(jpi,jpj)      , STAT=ierr(ii) )
+#endif
+
 #if defined key_agrif
-      ALLOCATE( ub2_i_b(jpi,jpj), vb2_i_b(jpi,jpj)                                  , STAT=ierr(6) )
+      ii = ii+1                             ! AGRIF: ???
+      ALLOCATE( ub2_i_b(jpi,jpj), vb2_i_b(jpi,jpj)                                  , STAT=ierr(ii) )
 #endif
          !
       oce_alloc = MAXVAL( ierr )
diff --git a/src/OCE/stp2d.F90 b/src/OCE/stp2d.F90
index 6ca69903..96170dc2 100644
--- a/src/OCE/stp2d.F90
+++ b/src/OCE/stp2d.F90
@@ -109,7 +109,13 @@ CONTAINS
       !
       !                             !*  compute advection + coriolis *!
       !
-      r3t(:,:,Kaa) =  ssh(:,:,Kaa) * r1_ht_0(:,:)                           ! ratio at t-point at Kaa (n-1)
+      CALL ssh_nxt( kt, Kbb, Kbb, ssh, Kaa )
+      !
+      IF( .NOT.lk_linssh ) THEN
+         DO_2D_OVR( 1, nn_hls, 1, nn_hls )      ! loop bounds limited by ssh definition in ssh_nxt
+            r3t(ji,jj,Kaa) =  ssh(ji,jj,Kaa) * r1_ht_0(ji,jj)               ! "after" ssh/h_0 ratio guess at t-column at Kaa (n+1)
+         END_2D
+      ENDIF
       !
       CALL wzv    ( kt, Kbb, Kbb, Kaa , uu(:,:,:,Kbb), vv(:,:,:,Kbb), ww )  ! ww guess at Kbb (n)
       !
-- 
GitLab