Name and subject of the action

Last edition: Wikinfo(changed_ts)? by Wikinfo(changed_by)?

The PI is responsible to closely follow the progress of the action, and especially to contact NEMO project manager if the delay on preview (or review) are longer than the 2 weeks expected.

1. Summary

Summary

Action	RK3 stage 1
PI(S)	Gurvan et Sibylle
Digest	Run a GYRE configuration with new RK3 scheme
Dependencies	If any
Branch	source:/NEMO/branches/2021/dev_r14318_RK3_stage1
Previewer(s)	Gurvan
Reviewer(s)	Names
Ticket	#2605 #2715

Description

Error: Failed to load processor box

No macro or processor named 'box' found

RK3 time stepping implementation for NEMO includes at this stage dynamic and active tracers implementation, time spitting single first with 2D mode integration.

...

Implementation

Error: Failed to load processor box

No macro or processor named 'box' found

RK3 implementation is splitted up into :

• code preparation
• dynamic and active tracers (barocline)
• vertical physics (TKE) ?
• barotropic mode (barotrope)
• mass forcing
• passive tracers

Code preparation In order to preserve constancy property velocity for momentum and active tracers must be the same. Advection routines in flux form are modified to take (u,v,w) as an input argument. In order to use advection routines for the barotropic mode we need the possibility to de-activate vertical advection computation. Advection routines in flux and vector form are modified to take an optional argument (no_zad) to do so.

Barocline part For sake of simplicity we started to implement RK3 regarding a GYRE configuration validation with no barotrope mode (ssh, uu_b, un_adv are set to zero at each time step). Forcing have been removed except winds and heat flux. key_qco is active and vertical physics is modeled as constant with high viscosity coefficients.

• Prepare routines
  • Change eos divhor and sshwzv interface.
• Add RK3 time stepping routines
  • rk3stg deals with time integration at N+1/3, N+1/2 and N+1
  • stprk3 orchestrates

Barotrope part In order to validate 2D mode implementation we remove above zero forcing for barotropic variables mass forcing remains to zero.

• Prepare routines
  • Change dynadv, dynvor, dynspg_ts
• Add RK3 2D mode time stepping routines
  • rk3ssh prepare 2D forcing, get dynamics 2D RHS from 3D trends, integrate 2D mode

...

Implementation details : Code preparation

r14418 Allow an advective velocity to be passed as an argument.
3D velocity can be a pointer.

OCE
 |-- oce.F90
     REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:), TARGET   ::   uu   ,  vv     !: horizontal velocities        [m/s] 
     REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  , TARGET   ::   ww             !: vertical velocity            [m/s]

3D velocity added as an input argument of advective routines passed through dyn_adv

OCE
 |--DYN
     |-- dynadv.F90
         SUBROUTINE dyn_adv( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )
         ...
         CALL dyn_adv_cen2( kt                , Kmm, puu, pvv, Krhs, pau, pav, paw )   ! 2nd order centered scheme 
         CALL dyn_adv_ubs ( kt           , Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw )   ! 3rd order UBS      scheme (UP3) 
     |-- dynadv_cen2.F90
         SUBROUTINE dyn_adv_cen2( kt, Kmm, puu, pvv, Krhs, pau, pav, paw )  
         ...       
         IF( PRESENT( pau ) ) THEN     ! RK3: advective velocity (pau,pav,paw) /= advected velocity (puu,pvv,ww) 
            zptu => pau(:,:,:) 
            ...
         zfu(:,:,jk) = 0.25_wp * e2u(:,:) * e3u(:,:,jk,Kmm) * zptu(:,:,jk) 
     |-- dynadv_ubs.F90
         SUBROUTINE dyn_adv_ubs( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw ) 
         ...
         IF( PRESENT( pau ) ) THEN     ! RK3: advective velocity (pau,pav,paw) /= advected velocity (puu,pvv,ww) 
            zptu => pau(:,:,:) 
            ...
         zfu(:,:,jk) = e2u(:,:) * e3u(:,:,jk,Kmm) * zptu(:,:,jk) 
 |--TRA  
     |-- traadv.F90
         SUBROUTINE tra_adv( kt, Kbb, Kmm, pts, Krhs, pau, pav, paw ) 
         ...
         IF( PRESENT( pau ) ) THEN     ! RK3: advective velocity (pau,pav,paw) /= advected velocity (puu,pvv,ww) 
            zptu => pau(:,:,:) 
            ...
         zuu(ji,jj,jk) = e2u  (ji,jj) * e3u(ji,jj,jk,Kmm) * ( zptu(ji,jj,jk) + usd(ji,jj,jk) ) 

Finally this new structure is used in step and tested with usual velocities

OCE
 |-- stpmlf.F90
     REAL(wp), TARGET     , DIMENSION(jpi,jpj,jpk) ::   zau, zav, zaw   ! advective velocity 
     ...
     zau(:,:,:) = uu(:,:,:,Nnn)        !!st   patch for MLF will be computed in RK3
     ...
     CALL dyn_adv( kstp, Nbb, Nnn      , uu, vv, Nrhs, zau, zav, zaw )  ! advection (VF or FF)   ==> RHS
     ...
     CALL tra_adv    ( kstp, Nbb, Nnn, ts, Nrhs, zau, zav, zaw )  ! hor. + vert. advection ==> RHS 

Results should be exactly the same as the ones from from the trunk. It was not the case for an OVERFLOW. The use of ln_wAimp=T changes ww at the truncature in diawri.F90, and that produces a small error. This has been corrected.

r14428 Allow vertical advection to be de-activated with an optionnal input argument : no_zad.

3D velocity added as an input argument of advective routines passed through dyn_adv

OCE
 |--DYN
     |-- dynadv.F90
         SUBROUTINE dyn_adv( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )
         ...
         CALL dyn_adv_cen2( kt                , Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )   ! 2nd order centered scheme 
         CALL dyn_adv_ubs ( kt           , Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )   ! 3rd order UBS      scheme (UP3) 
     |-- dynadv_cen2.F90
         SUBROUTINE dyn_adv_cen2( kt, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad )  
         ...       
         IF( PRESENT( no_zad ) ) THEN  !==  No vertical advection  ==!   (except if linear free surface) 
            IF( ln_linssh ) THEN                ! linear free surface: advection through the surface z=0 
               DO_2D( 0, 0, 0, 0 ) 
                  zzu = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm) 
                  zzv = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm) 
                  puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) - zzu * r1_e1e2u(ji,jj)   & 
                     &                                        / e3u(ji,jj,1,Kmm) 
                  pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) - zzv * r1_e1e2v(ji,jj)   & 
                     &                                        / e3v(ji,jj,1,Kmm) 
               END_2D 
            ENDIF 
         ! 
         ELSE                          !==  Vertical advection  ==! 
            ...
     |-- dynadv_ubs.F90
         SUBROUTINE dyn_adv_ubs( kt, Kbb, Kmm, puu, pvv, Krhs, pau, pav, paw, no_zad ) 
         ...
         IF( PRESENT( no_zad ) ) THEN  !==  No vertical advection  ==!   (except if linear free surface) 
            IF( ln_linssh ) THEN                ! linear free surface: advection through the surface z=0 
               DO_2D( 0, 0, 0, 0 ) 
                  zzu = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji+1,jj) * zpt_w(ji+1,jj,1) ) * puu(ji,jj,1,Kmm)
                  zzv = 0.5_wp * ( e1e2t(ji,jj) * zpt_w(ji,jj,1) + e1e2t(ji,jj+1) * zpt_w(ji,jj+1,1) ) * pvv(ji,jj,1,Kmm) 
                  puu(ji,jj,1,Krhs) = puu(ji,jj,1,Krhs) - zzu * r1_e1e2u(ji,jj)   & 
                     &                                        / e3u(ji,jj,1,Kmm) 
                  pvv(ji,jj,1,Krhs) = pvv(ji,jj,1,Krhs) - zzv * r1_e1e2v(ji,jj)   & 
                     &                                        / e3v(ji,jj,1,Kmm) 
               END_2D 
            ENDIF 
         ! 
         ELSE                          !==  Vertical advection  ==! 

Gurvan added a loop optimisation for dynzad.F90

OCE
 |--DYN
     |-- dynzad.F90
All the loops are now gather in a single one.  

Implementation details : barocline processing

r14547 Allow RK3 time-stepping with 2D mode damped.
div_hor interface and sshwzv interface have been changed accordingly for RK3. eos also changed in order to avoid gdep to be used as an input argument in the key_qco framework.

OCE
 |--DYN
     |-- divhor.F90
         SUBROUTINE div_hor_RK3( kt, Kbb, Kmm, puu, pvv, pe3divUh ) 
     |-- sshwzv.F90
         SUBROUTINE wzv_RK3( kt, Kbb, Kmm, Kaa, puu, pvv, pww )
 |--TRA
     |-- eosbn2.F90
         SUBROUTINE eos_insitu_New( pts, Knn, prd ) 

Time step no longer need to be doubled. rk3 routines are added to the code and stprk3 is called through nemogcm when key_RK3 is active.

OCE
 |--DOM
     |-- domain.F90
         #if defined key_RK3 
              rDt   =         rn_Dt 
 	      r1_Dt = 1._wp / rDt
          ...
 |-- nemogcm.F90
     # if defined key_RK3 
             USE stprk3 
     ...
 |-- stprk3.F90
 |-- stprk3-stg.F90

Has been tested and validated against an modified leap frog GYRE in the same configuration with the same namelist.

r14549 Allow RK3 time-stepping with 2D mode.
Prepare forcings and barotropic 2D fields. dynspg_ts remains for 2D mode integration. dyn_vor_RK3 only computes 2D relative vorticity and metric term from 3D to 2D. stp_2D is called by stprk3 in single first.

OCE
 |--DYN
     |-- dynspg_ts.F90
         #remove k_only_ADV
         PUBLIC dyn_drg_init      ! called by rk3ssh
         !  Phase 1 : Coupling between general trend and barotropic estimates (1st step)
         IF( kt == nit000 ) THEN 
 	    IF( .NOT.ln_bt_fw  .OR. ln_bt_av )   CALL ctl_stop( 'dyn_spg_ts: RK3 requires ln_bt_fw=T AND ln_bt_av=F') 
 	 ENDIF
         !                          ! set values computed in RK3_ssh 
 	 zssh_frc(:,:) = sshe_rhs(:,:) 
 	 zu_frc(:,:) =   Ue_rhs(:,:) 
 	 zv_frc(:,:) =   Ve_rhs(:,:) 
 	 zCdU_u  (:,:) = CdU_u   (:,:) 
 	 zCdU_v  (:,:) = CdU_v   (:,:)
         !
 	 IF( kt == nit000 .OR. .NOT. ln_linssh )   CALL dyn_cor_2D_init( Kmm )   ! Set zwz, the barotropic Coriolis force coefficient
         ! Phase 3. update the general trend with the barotropic trend 
         IF(.NOT.ln_bt_av ) THEN                          !* Update Kaa barotropic external mode  
            uu_b(:,:,Kaa) = ua_e  (:,:) 
  	    pvv_b(:,:,Kaa) = va_e  (:,:) 
  	    pssh (:,:,Kaa) = ssha_e(:,:) 
  	 ENDIF
         un_adv...
         ubar...
     |-- dynspg.F90
         #remove k_only_ADV
     |-- dynvor.F90
         SUBROUTINE dyn_vor_3D( kt, Kmm, puu, pvv, Krhs ) 
     |-- dynzdf.F90
         zDt_2 = rDt * 0.5_wp # small cosmetic optim
 |-- stprk3_stg.F90
 |-- stp2d.F90
 |-- stprk3.F90
179	      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
180	      !  RK3 : single first external mode computation
181	      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
182	
183	      CALL stp_2D( kstp, Nbb, Nbb, Naa, Nrhs )   ! out: ssh, (uu_b,vv_b) and (un_adv,vn_adv) at Naa
184	
185	
186	      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
187	      !  RK3 time integration
188	      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
189	
190	      ! Stage 1 :
191	      CALL stp_RK3_stg( 1, kstp, Nbb, Nbb, Nrhs, Naa )
192	      !
193	      Nrhs = Nnn   ;   Nnn  = Naa   ;   Naa  = Nrhs    ! Swap: Nbb unchanged, Nnn <==> Naa
194	      !
195	      ! Stage 2 :
196	      CALL stp_RK3_stg( 2, kstp, Nbb, Nnn, Nrhs, Naa )
197	      !
198	      Nrhs = Nnn   ;   Nnn  = Naa   ;   Naa  = Nrhs    ! Swap: Nbb unchanged, Nnn <==> Naa
199	      !
200	      ! Stage 3 :
201	      CALL stp_RK3_stg( 3, kstp, Nbb, Nnn, Nrhs, Naa )
202	      !
203	      Nrhs = Nbb   ;   Nbb  = Naa   ;   Naa  = Nrhs    ! Swap: Nnn unchanged, Nbb <==> Naa

Implementation of AGRIF with RK3

r14800 introduces the major changes to activate AGRIF with RK3

• Add provision of ssh data for setting bcs during barotropic loop in agrif_dta_ts subroutine. Indeed, MLF makes use of sshwzv where a first guess of ssh was done and corrected at AGRIF bdys. This routine is not called anymore with RK3.

NST
 |-- agrif_oce_interp.F90
     SUBROUTINE Agrif_dta_ts( kt )
     #if defined key_RK3
        Agrif_SpecialValue    = 0._wp
        Agrif_UseSpecialValue = .TRUE.
        CALL Agrif_Bc_variable(sshn_id, procname=interpsshn )
        Agrif_UseSpecialValue = .FALSE.
     #endif

• Important change of time indexes on Parent grid in stprk3.F90: "Now" arrays on parent refers to "Nbb" index for RK3 instead of "Nnn" with MLF:

  OCE
   |-- stprk3.F90
       Kbb_a = Nbb; Kmm_a = Nbb; Krhs_a = Nrhs      ! agrif_oce module copies of time level indices
  

• Added (optional) time interpolation at intermediate stages in Agrif_dyn and Agrif_tra:

  NST
   |-- agrif_oce_interp.F90
        SUBROUTINE Agrif_tra( kt, kstg )
        !!----------------------------------------------------------------------
        !!                  ***  ROUTINE Agrif_tra  ***
        !!----------------------------------------------------------------------
        INTEGER, INTENT(in) ::   kt
        INTEGER, OPTIONAL, INTENT(in) :: kstg
        REAL(wp) :: ztindex
        !
        IF( Agrif_Root() )   RETURN
        !
        ! Set time index depending on stage in case of RK3 time stepping:
        IF ( PRESENT( kstg ) ) THEN
           ztindex = REAL(Agrif_Nbstepint(), wp)
           IF     ( kstg == 1 ) THEN
              ztindex = ztindex + 1._wp / 3._wp
           ELSEIF ( kstg == 2 ) THEN
              ztindex = ztindex + 1._wp / 2._wp
           ELSEIF ( kstg == 3 ) THEN
              ztindex = ztindex + 1._wp
           ENDIF
           ztindex = ztindex / Agrif_Rhot()
        ELSE
           ztindex = REAL(Agrif_Nbstepint()+1, wp) / Agrif_Rhot()
        ENDIF
        !
        Agrif_SpecialValue    = 0._wp
        Agrif_UseSpecialValue = .TRUE.
        l_vremap              = ln_vert_remap
        !
        CALL Agrif_Bc_variable( ts_interp_id, calledweight=ztindex, procname=interptsn )
        !
        Agrif_UseSpecialValue = .FALSE.
        l_vremap              = .FALSE.
        !
     END SUBROUTINE Agrif_tra
  

Functionnality activated in stprk3_stg.F90

OCE
 |-- stprk3_stg.F90
         !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
         ! Set boundary conditions
         !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
         !
     # if defined key_agrif
            CALL Agrif_tra( kstp, kstg )             !* AGRIF zoom boundaries
            CALL Agrif_dyn( kstp, kstg )
     # endif

• Removed updates of Parent "before" fields coming from asselin filtering. cpp directives "#if ! defined key_RK3" have been used in agrif_oce_update.F90 module.

• Added Agrif sponge trend in barotropic mode (stp2d.F90). It is taken as constant over the barotropic loop, but that's the only source of damping near grid boundaries... Should add time dependent nudging during barotropic loop which could help damping instabilities as reported below without time averaging of barotropic variables.

• The changes above have also been made for agrif routines related to TOP. At that stage, TOP was nevertheless not functional with RK3.

• Changes validated (qualitavely) in VORTEX test case. Successful doubling of time step. However results are substantiallly different in term of VORTEX amplitude around grid interface => impact of unchanged adimensionnalized sponge/nudging coefficients ?

• NB1: I noticed instabilities preventing runs to finish when switching off barotropic time averaging (and using Demange's dissipative time stepping). This occurs also with MLF. Is the periodic data exchange between grids (at baroclinic time step) suitable with such an approach ?

• NB2: Restartability not checked yet but it is likely that some changes on Agrif time indexes (such as the one reported above) will be needed at initialization.

Solar flux penetration optimization (r14990 start r15083 fully debugged rxxxxx optimized)

Tracer surface boundary condition optimization (r14993)

Ice Shelve now compatible with qco (r15083)

Adapt TOP for RK3

Merge with trunk r15557

Trunk has new management of the halos, RK3 additionnal routines have to be updated accordingly.

• lbc_lnk_multi routines no longer exist : they have been integrated into lbc_lnk generic routines, then they have to be replaced in :
  

=> src/OCE/DYN/dynspg_ts.F90
=> src/OCE/stprk3_stg.F90

• loop have changed and optimized for limiting communications : when relevant RK3 routines will mimic MLF ones
  

=> src/OCE/DYN/divhor.F90 : div_hor_RK3
=> src/OCE/DYN/sshwzv.F90 : wzv_RK3
=> src/OCE/TRA/trasbc.F90 : tra_sbc_RK3

• we have optimized tra_qsr routine by limiting intermidiate storing array, computations are now done only in the interior it works for all variables except fraqsr_1lev. Indeed it needs to be defined over the whole domain for the ICE. When not np_BIO fraqsr_1lev depends on qsr_hc, to minimize the code changed I decided to keep qsr_hc defined in the interior and to make a communication on fraqsr_1lev but it may not be the best solution. see with GM !
  

=> src/OCE/TRA/traqsr.F90 (fraqsr_1lev)

• src/OCE/stprk3_stg.F90 also has to be updated.
  

=> A communication is required on r3. Indeed r3. at Kbb and Kmm can be used over the whole domain. In MLF case r3. at Kbb and Kmm are lbc-ised when starting a time step, while in RK3 case because of the stages within a step and r3. at Kmm is never lbc-ised. That is why this communication is required !
=> Advective speeds at Kmm are computed on the whole domain and also required r3u to be defined over the whole domain.
=> As far as I understand ts, tr @ Kaa are not only computed in the interior in order to limit the number of loops, and sometime (ex pt update in src/OCE/TRA/traadv_fct.F90) values are computed but not used. In traadv case pt value outside the interior in erased with a lbc_lnk at the end of the time step. That is why ts/tr have to be initialised over the whole domain. Maybe it should be done at the init only as we discussed with GM !
=> As what has been done for MLF "IF( nn_hls == 2 .AND. l_zdfsh2 ) CALL lbc_lnk( 'stp_RK3_stg', avm_k, 'W', 1.0_wp )" needs to be externalised from the tke routine.

• rhop removal needs to be adressed properly. For now it has been quickly removed from RK3 but it is required for some DIA routines and over all for tra_mle (activated in ORCA2_ICE_PISCES configuration). The initial target was to compute rhop locally when needed on order to avoid to store a useless rhop array. But is means about 6 times. see with GM !!!
  

=> src/OCE/TRD/trdglo.F90 : glo_dyn_wri changes maye useless since glo_dyn_wri never called !
=> src/OCE/stprk3.F90 rhop@Nbb now computed with rhd at each time step
=> src/OCE/IOM/restart.F90 : remove rhop from restart

• initialisation with restart sometimes requires Kmm field with are not intialised in RK3, as a matter of simplicity we just copy Kbb into Kmm for ssh, ts, uu, vv fields.
  

=> src/OCE/IOM/restart.F90 : rst_read
=> src/OCE/IOM/restart.F90 : rst_read_ssh

• finally we wrote additionnal comments and did some cleaning !

DEC_dbg_RK3/src/OCE/stprk3.F90

@@ -156 +156 @@
 !!st                         CALL zdf_phy( kstp, Nbb, Nnn, Nrhs )   ! vertical physics update (top/bot drag, avt, avs, avm + MLD)
                          CALL zdf_phy( kstp, Nbb, Nbb, Nrhs )   ! vertical physics update (top/bot drag, avt, avs, avm + MLD)
 !!gm
+!!st check with gm !!!!!
+                         CALL eos( ts(:,:,:,:,Nbb), rhd, rhop, gdept_0(:,:,:) ) ! now in situ density for tramle slope computation
       !  LATERAL  PHYSICS
       !
       IF( l_ldfslp ) THEN                             ! slope of lateral mixing
 !!gm gdep
-                         CALL eos( ts(:,:,:,:,Nbb), rhd, gdept_0(:,:,:) )               ! before in situ density
+!!st         CALL eos( ts(:,:,:,:,Nbb), rhd, gdept_0(:,:,:) )               ! before in situ density
 
          IF( ln_zps .AND. .NOT. ln_isfcav)                                    &
             &            CALL zps_hde    ( kstp, Nbb, jpts, ts(:,:,:,:,Nbb), gtsu, gtsv,  &  ! Partial steps: before horizontal gradient

DEC_dbg_RK3/src/OCE/IOM/restart.F90

@@ -176 +176 @@
          CALL iom_rstput( kt, nitrst, numrow, 'vn'  , vv(:,:,:       ,Kmm) )
          CALL iom_rstput( kt, nitrst, numrow, 'tn'  , ts(:,:,:,jp_tem,Kmm) )
          CALL iom_rstput( kt, nitrst, numrow, 'sn'  , ts(:,:,:,jp_sal,Kmm) )
-         IF( .NOT.lk_SWE )   CALL iom_rstput( kt, nitrst, numrow, 'rhop', rhop )
+!!$         IF( .NOT.lk_SWE )   CALL iom_rstput( kt, nitrst, numrow, 'rhop', rhop )
 #endif
       ENDIF
 
@@ -290 +290 @@
       CALL iom_get( numror, jpdom_auto, 'sb'   , ts(:,:,:,jp_sal,Kbb) )
       CALL iom_get( numror, jpdom_auto, 'uu_b' , uu_b(:,:       ,Kbb), cd_type = 'U', psgn = -1._wp )
       CALL iom_get( numror, jpdom_auto, 'vv_b' , vv_b(:,:       ,Kbb), cd_type = 'V', psgn = -1._wp )
+      uu(:,:,:  ,Kmm) = uu(:,:,:  ,Kbb)         ! all now value set to before for initialisation (sbc_ssm_init)
+      vv(:,:,:  ,Kmm) = vv(:,:,:  ,Kbb)
+      ts(:,:,:,:,Kmm) = ts(:,:,:,:,Kbb)
 #else
       !                             !*  Read Kmm fields   (MLF only)
       IF(lwp) WRITE(numout,*)    '           Kmm u, v and T-S fields read in the restart file'
@@ -313 +316 @@
       ENDIF
 #endif
       !
-      IF( .NOT.lk_SWE ) THEN
-         IF( iom_varid( numror, 'rhop', ldstop = .FALSE. ) > 0 ) THEN
-            CALL iom_get( numror, jpdom_auto, 'rhop'   , rhop )   ! now    potential density
-         ELSE
-#if defined key_RK3
-            CALL eos( ts, Kbb, rhop )
-#else
-            CALL eos( ts, Kmm, rhop )
-#endif
-!!#if defined key_qco
-!!            ALLOCATE( zgdept(jpi,jpj,jpk) )
-!!            DO jk = 1, jpk
-!!               zgdept(:,:,jk) = gdept(:,:,jk,Kmm)
-!!            END DO
-!!            CALL eos( ts(:,:,:,:,Kmm), rhd, rhop, zgdept )
-!!            DEALLOCATE( zgdept )
-!!#else
-!!            CALL eos( ts(:,:,:,:,Kmm), rhd, rhop, gdept(:,:,:,Kmm) )
-!!#endif
-         ENDIF
-      ENDIF
+!!$      IF( .NOT.lk_SWE ) THEN
+!!$         IF( iom_varid( numror, 'rhop', ldstop = .FALSE. ) > 0 ) THEN
+!!$            CALL iom_get( numror, jpdom_auto, 'rhop'   , rhop )   ! now    potential density
+!!$         ELSE
+!!$#if defined key_RK3
+!!$            CALL eos( ts, Kbb, rhop ) ! not rhop but rhd 
+!!$#else
+!!$            CALL eos( ts, Kmm, rhop )
+!!$#endif
+!!$         ENDIF
+!!$      ENDIF
       !
    END SUBROUTINE rst_read
 
@@ -377 +370 @@
          CALL iom_get( numror, jpdom_auto, 'sshb'   , ssh(:,:,Kbb) )
          !
          !                                     !*  RK3: Set ssh at Kmm for AGRIF
-         ssh(:,:,Kmm) = 0._wp
+         ssh(:,:,Kmm) = ssh(:,:,Kbb)
 #else
          !                                     !*  MLF: Read ssh at Kmm
          IF(lwp) WRITE(numout,*)

DEC_dbg_RK3/src/OCE/DYN/divhor.F90

@@ -83 +83 @@
       ! 
       pe3divUh(:,:,:) = 0._wp    !!gm to be applied to the halos only
       !
-      DO_3D( 0, 0, 0, 0, 1, jpkm1 )                                                  
+      DO_3D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls, 1, jpkm1 )
          hdiv(ji,jj,jk) = (   e2u(ji  ,jj) * e3u(ji  ,jj,jk,Kmm) * puu(ji  ,jj,jk)      &
             &               - e2u(ji-1,jj) * e3u(ji-1,jj,jk,Kmm) * puu(ji-1,jj,jk)      &
             &               + e1v(ji,jj  ) * e3v(ji,jj  ,jk,Kmm) * pvv(ji,jj  ,jk)      &
@@ -100 +100 @@
       !
       IF( ln_isf )   CALL isf_hdiv( kt, Kmm, hdiv )            !==  + ice-shelf mass exchange ==!
       !
-      CALL lbc_lnk( 'divhor', hdiv, 'T', 1._wp )   !   (no sign change)
+      IF( nn_hls==1 )   CALL lbc_lnk( 'divhor', hdiv, 'T', 1._wp )   !   (no sign change)
       !
 !!gm Patch before suppression of hdiv from all modules that use it
 !      DO_3D( 0, 0, 0, 0, 1, jpkm1 )                            !==  e3t * Horizontal divergence  ==!
@@ -107 +107 @@
 !         pe3divUh(ji,jj,jk) = hdiv(ji,jj,jk) * e3t(ji,jj,jk,Kmm)
 !      END_3D
 !JC: over whole domain, and after lbclnk on hdiv to prevent from reproducibility issues
-      DO jk=1, jpkm1
-         pe3divUh(:,:,jk) = hdiv(:,:,jk) * e3t(:,:,jk,Kmm)
-      END DO
+      DO_3D_OVR( nn_hls-1, nn_hls, nn_hls-1, nn_hls, 1, jpkm1 )
+         pe3divUh(ji,jj,jk) = hdiv(ji,jj,jk) * e3t(ji,jj,jk,Kmm)
+      END_3D
 !!gm end
       !
       !

DEC_dbg_RK3/src/OCE/DYN/sshwzv.F90

@@ -328 +328 @@
          DO jk = 1, jpkm1
             ! horizontal divergence of thickness diffusion transport ( velocity multiplied by e3t)
             ! - ML - note: computation already done in dom_vvl_sf_nxt. Could be optimized (not critical and clearer this way)
-            DO_2D( 0, 0, 0, 0 )
+            DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls ) 
                zhdiv(ji,jj,jk) = r1_e1e2t(ji,jj) * ( un_td(ji,jj,jk) - un_td(ji-1,jj,jk) + vn_td(ji,jj,jk) - vn_td(ji,jj-1,jk) )
             END_2D
          END DO
-         CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
+         IF( nn_hls == 1)   CALL lbc_lnk('sshwzv', zhdiv, 'T', 1.0_wp)  ! - ML - Perhaps not necessary: not used for horizontal "connexions"
          !                             ! Is it problematic to have a wrong vertical velocity in boundary cells?
          !                             ! Same question holds for hdiv. Perhaps just for security
-         DO jk = jpkm1, 1, -1                       ! integrate from the bottom the hor. divergence
-            ! computation of w
-            pww(:,:,jk) = pww(:,:,jk+1) - (   ze3div(:,:,jk) + zhdiv(:,:,jk)   &
-               &                            + r1_Dt * (  e3t(:,:,jk,Kaa)       &
-               &                                       - e3t(:,:,jk,Kbb) )   ) * tmask(:,:,jk)
-         END DO
-         !          IF( ln_vvl_layer ) pww(:,:,:) = 0.e0
+         DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 )     ! integrate from the bottom the hor. divergence
+            pww(ji,jj,jk) = pww(ji,jj,jk+1) - (   ze3div(ji,jj,jk) + zhdiv(ji,jj,jk)   &
+                 &                            + r1_Dt * (  e3t(ji,jj,jk,Kaa)       &
+                 &                                       - e3t(ji,jj,jk,Kbb) )   ) * tmask(ji,jj,jk)
+         END_3D
+         !
          DEALLOCATE( zhdiv ) 
          !                                            !=================================!
       ELSEIF( ln_linssh )   THEN                      !==  linear free surface cases  ==!
          !                                            !=================================!
-         DO jk = jpkm1, 1, -1                               ! integrate from the bottom the hor. divergence
-            pww(:,:,jk) = pww(:,:,jk+1) - ze3div(:,:,jk)
-         END DO
+         DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 )     ! integrate from the bottom the hor. divergence
+            pww(ji,jj,jk) = pww(ji,jj,jk+1) - ze3div(ji,jj,jk) 
+         END_3D
          !                                            !==========================================!
       ELSE                                            !==  Quasi-Eulerian vertical coordinate  ==!   ('key_qco')
          !                                            !==========================================!
-         DO jk = jpkm1, 1, -1                               ! integrate from the bottom the hor. divergence
-            !                                               ! NB: [e3t[a] -e3t[b] ]=e3t_0*[r3t[a]-r3t[b]]
-            pww(:,:,jk) = pww(:,:,jk+1) - (  ze3div(:,:,jk)                            &
-               &                            + r1_Dt * e3t_0(:,:,jk) * ( r3t(:,:,Kaa) - r3t(:,:,Kbb) )  ) * tmask(:,:,jk)
-         END DO
+         DO_3DS( nn_hls-1, nn_hls, nn_hls-1, nn_hls, jpkm1, 1, -1 )     ! integrate from the bottom the hor. divergence
+         !                                                              ! NB: [e3t[a] -e3t[b] ]=e3t_0*[r3t[a]-r3t[b]]
+            pww(ji,jj,jk) = pww(ji,jj,jk+1) - (  ze3div(ji,jj,jk)                            &
+               &                               + r1_Dt * e3t_0(ji,jj,jk) * ( r3t(ji,jj,Kaa) - r3t(ji,jj,Kbb) )  ) * tmask(ji,jj,jk)
+         END_3D
       ENDIF
 
       IF( ln_bdy ) THEN
          DO jk = 1, jpkm1
-            pww(:,:,jk) = pww(:,:,jk) * bdytmask(:,:)
+            DO_2D( nn_hls-1, nn_hls, nn_hls-1, nn_hls )
+               pww(ji,jj,jk) = pww(ji,jj,jk) * bdytmask(ji,jj)
+            END_2D
          END DO
       ENDIF
       !

DEC_dbg_RK3/src/OCE/DYN/dynspg_ts.F90

@@ -802 +802 @@
          END_2D
 # endif   
          !
-         CALL lbc_lnk_multi( 'dynspg_ts', puu_b, 'U', -1._wp, pvv_b, 'V', -1._wp ) ! Boundary conditions
+         CALL lbc_lnk( 'dynspg_ts', puu_b, 'U', -1._wp, pvv_b, 'V', -1._wp ) ! Boundary conditions
          !
       ENDIF
       !

DEC_dbg_RK3/src/OCE/stprk3_stg.F90

@@ -33 +33 @@
 # if defined key_agrif
    USE agrif_oce_interp
 # endif
+!!st -dbg stp2d (r3u)
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    
    !
    USE prtctl         ! print control
@@ -147 +149 @@
       !                     !==  ssh/h0 ratio at Kaa  ==! 
       !
       IF( .NOT.lk_linssh )   CALL dom_qco_r3c( ssh(:,:,Kaa), r3t(:,:,Kaa), r3u(:,:,Kaa), r3v(:,:,Kaa), r3f(:,:) )   ! "after" ssh/h_0 ratio at t,u,v-column
+!!      SELECT CASE( kstg )
+         !
+!!      CASE ( 3 )
+         !!st required at each stage for div hor loops
+         CALL lbc_lnk( 'stprk3_stg', r3u(:,:,Kaa), 'U', 1._wp, r3v(:,:,Kaa), 'V', 1._wp, r3f(:,:), 'F', 1._wp )
+         !
+!!      END SELECT
       !
       !
       !                     !==  advective velocity at Kmm  ==!
       !
       !                                            !- horizontal components -!   (zaU,zaV) 
-      zub(:,:) = un_adv(:,:)*r1_hu(:,:,Kmm) - uu_b(:,:,Kmm)    ! barotropic velocity correction
-      zvb(:,:) = vn_adv(:,:)*r1_hv(:,:,Kmm) - vv_b(:,:,Kmm)
-      DO jk = 1, jpkm1                                         ! horizontal advective velocity
-         zaU(:,:,jk) = uu(:,:,jk,Kmm) + zub(:,:)*umask(:,:,jk)
-         zaV(:,:,jk) = vv(:,:,jk,Kmm) + zvb(:,:)*vmask(:,:,jk)
-      END DO
+      DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
+         zub(ji,jj) = un_adv(ji,jj)*r1_hu(ji,jj,Kmm) - uu_b(ji,jj,Kmm)    ! barotropic velocity correction
+         zvb(ji,jj) = vn_adv(ji,jj)*r1_hv(ji,jj,Kmm) - vv_b(ji,jj,Kmm)
+      END_2D
+      DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpkm1 )           ! horizontal advective velocity
+         zaU(ji,jj,jk) = uu(ji,jj,jk,Kmm) + zub(ji,jj)*umask(ji,jj,jk)
+         zaV(ji,jj,jk) = vv(ji,jj,jk,Kmm) + zvb(ji,jj)*vmask(ji,jj,jk)
+      END_3D
+
       !                                            !- vertical components -!   ww
                          CALL wzv  ( kstp, Kbb, Kmm, Kaa, zaU, zaV, ww )     ! ww cross-level velocity
 
@@ -223 +235 @@
          IF(.NOT. ln_trcadv_mus ) THEN
             !
             DO jn = 1, jptra
-               DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-               tr(ji,jj,jk,jn,Krhs) = 0._wp                              ! set tracer trends to zero
-               END_3D
+               tr(:,:,:,jn,Krhs) = 0._wp                              ! set tracer trends to zero !!st ::: required because of tra_adv new loops
             END DO
             !                                      !==  advection of passive tracers  ==!
             rDt_trc = rDt
@@ -253 +263 @@
          !                 !---------------!
          !
          DO jn = 1, jptra
-            DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-            tr(ji,jj,jk,jn,Krhs) = 0._wp                              ! set tracer trends to zero
-            END_3D
+            tr(:,:,:,jn,Krhs) = 0._wp                              ! set tracer trends to zero !!st ::: required because of tra_adv new loops
          END DO
          !                                         !==  advection of passive tracers  ==!
          rDt_trc = rDt
@@ -393 +401 @@
       zvb(A2D(0)) = vv_b(A2D(0),Kaa) - SUM( e3v_0(A2D(0),:)*vv(A2D(0),:,Kaa), 3 ) * r1_hv_0(A2D(0))
       !
       DO jk = 1, jpkm1                                            ! corrected horizontal velocity
-         uu(:,:,jk,Kaa) = uu(:,:,jk,Kaa) + zub(:,:)*umask(:,:,jk)
-         vv(:,:,jk,Kaa) = vv(:,:,jk,Kaa) + zvb(:,:)*vmask(:,:,jk)
+         uu(A2D(0),jk,Kaa) = uu(A2D(0),jk,Kaa) + zub(A2D(0))*umask(A2D(0),jk)
+         vv(A2D(0),jk,Kaa) = vv(A2D(0),jk,Kaa) + zvb(A2D(0))*vmask(A2D(0),jk)
       END DO
-!!st pourquoi ne pas mettre A2D(0) ici ? 
-         
 
+
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Set boundary conditions
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -408 +415 @@
             CALL Agrif_dyn( kstp, kstg )
 # endif
       !                                              !* local domain boundaries  (T-point, unchanged sign)
-      CALL lbc_lnk_multi( 'stp_RK3_stg', uu(:,:,:,       Kaa), 'U', -1., vv(:,:,:       ,Kaa), 'V', -1.   &
-         &                             , ts(:,:,:,jp_tem,Kaa), 'T',  1., ts(:,:,:,jp_sal,Kaa), 'T',  1. )
+      CALL lbc_lnk( 'stp_RK3_stg', uu(:,:,:,       Kaa), 'U', -1., vv(:,:,:       ,Kaa), 'V', -1.   &
+         &                       , ts(:,:,:,jp_tem,Kaa), 'T',  1., ts(:,:,:,jp_sal,Kaa), 'T',  1. )
+      ! lbc_lnk needed for zdf_sh2 when using nn_hls = 2, moved here to allow tiling in zdf_phy
+      IF( nn_hls == 2 .AND. l_zdfsh2 ) CALL lbc_lnk( 'stp_RK3_stg', avm_k, 'W', 1.0_wp )
       !
+      !
       !                                              !* BDY open boundaries
       IF( ln_bdy )   THEN
                                CALL bdy_tra( kstp, Kbb, ts,     Kaa )

DEC_dbg_RK3/src/OCE/ZDF/zdfphy.F90

@@ -209 +209 @@
       ioptio = 0
       IF( ln_zdfcst ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_CST   ;   ENDIF
       IF( ln_zdfric ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_RIC   ;   CALL zdf_ric_init          ;   ENDIF
-      IF( ln_zdftke ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_TKE   ;   CALL zdf_tke_init( Kmm )   ;   ENDIF
+      IF( ln_zdftke ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_TKE   ;   CALL zdf_tke_init          ;   ENDIF
       IF( ln_zdfgls ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_GLS   ;   CALL zdf_gls_init          ;   ENDIF
       IF( ln_zdfosm ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_OSM   ;   CALL zdf_osm_init( Kmm )   ;   ENDIF
       !
@@ -350 +350 @@
       IF( ln_zdfiwm )   CALL zdf_iwm( kt, Kmm, avm, avt, avs )   ! internal wave (de Lavergne et al 2017)
 
       !                                         !* Lateral boundary conditions (sign unchanged)
-      IF(nn_hls==1) THEN
+      IF(nn_hls==1) THEN                                         ! if nn_hls==2 lbc_lnk done in stpxxx routines
          IF( l_zdfsh2 ) THEN
             CALL lbc_lnk( 'zdfphy', avm_k, 'W', 1.0_wp , avt_k, 'W', 1.0_wp,   &
                   &                 avm  , 'W', 1.0_wp , avt  , 'W', 1.0_wp , avs , 'W', 1.0_wp )

DEC_dbg_RK3/src/OCE/ZDF/zdftke.F90

@@ -698 +698 @@
    END SUBROUTINE tke_avn
 
 
-   SUBROUTINE zdf_tke_init( Kmm )
+   SUBROUTINE zdf_tke_init
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE zdf_tke_init  ***
       !!
@@ -714 +714 @@
       !!----------------------------------------------------------------------
       USE zdf_oce , ONLY : ln_zdfiwm   ! Internal Wave Mixing flag
       !!
-      INTEGER, INTENT(in) ::   Kmm          ! time level index
       INTEGER             ::   ji, jj, jk   ! dummy loop indices
       INTEGER             ::   ios
       !!

DEC_dbg_RK3/src/OCE/TRA/traqsr.F90

@@ -179 +179 @@
 #endif
          END_3D
          !                                                     !- sea-ice : store the 1st level attenuation coefficient
-         WHERE( etot3(A2D(0),1) /= 0._wp )   ;   fraqsr_1lev(A2D(0)) = 1._wp - etot3(A2D(0),2) / etot3(A2D(0),1)
-         ELSEWHERE                           ;   fraqsr_1lev(A2D(0)) = 1._wp
+         WHERE( etot3(A2D(nn_hls),1) /= 0._wp )   ;   fraqsr_1lev(A2D(nn_hls)) = 1._wp - etot3(A2D(nn_hls),2) / etot3(A2D(nn_hls),1)
+         ELSEWHERE                                ;   fraqsr_1lev(A2D(nn_hls)) = 1._wp
          END WHERE
          !
       END SELECT
@@ -207 +207 @@
             &                             / e3t(ji,jj,jk,Kmm)
       END_3D
       !
-!!st7-2
       ! sea-ice: store the 1st ocean level attenuation coefficient
-      DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
+      DO_2D( 0, 0, 0, 0 )
          IF( qsr(ji,jj) /= 0._wp ) THEN   ;   fraqsr_1lev(ji,jj) = qsr_hc(ji,jj,1) / ( r1_rho0_rcp * qsr(ji,jj) )
          ELSE                             ;   fraqsr_1lev(ji,jj) = 1._wp
          ENDIF
@@ -217 +216 @@
       !                                       !===>>> CAUTION: lbc_lnk is required on fraqsr_lev since sea ice computes on the full domain
       !                                       !                otherwise restartability and reproducibility are broken 
       CALL lbc_lnk( 'tra_qsr', fraqsr_1lev(:,:), 'T', 1._wp )
-!!st      CALL lbc_lnk( 'tra_qsr', qsr_hc(:,:,:), 'T', 1._wp )
       !
       IF( iom_use('qsr3d') ) THEN      ! output the shortwave Radiation distribution
          ALLOCATE( zetot(A2D(nn_hls),jpk) )

DEC_dbg_RK3/src/OCE/TRA/traadv_fct.F90

@@ -186 +186 @@
                &      + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk  )   &
                &      + zwz(ji,jj,jk) - zwz(ji  ,jj  ,jk+1) ) * r1_e1e2t(ji,jj)
             !                               ! update and guess with monotonic sheme
+!!st pt should be updated with a DO_3D( 0,0,0,0, 1, jpkm1 ) values outside the interior are wrong but unused.
             pt(ji,jj,jk,jn,Krhs) =                   pt(ji,jj,jk,jn,Krhs) +       ztra   &
                &                                  / e3t(ji,jj,jk,Kmm ) * tmask(ji,jj,jk)
             zwi(ji,jj,jk)    = ( e3t(ji,jj,jk,Kbb) * pt(ji,jj,jk,jn,Kbb) + p2dt * ztra ) &

DEC_dbg_RK3/src/OCE/TRA/trasbc.F90

@@ -275 +275 @@
             
 !!gm  This should be moved into sbcmod.F90 module ? (especially now that ln_traqsr is read in namsbc namelist)
       IF( .NOT.ln_traqsr  .AND. kstg == 1) THEN     ! no solar radiation penetration
-         DO_2D( 0, 0, 0, 0 )
+         DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
             qns(ji,jj) = qns(ji,jj) + qsr(ji,jj)         ! total heat flux in qns
             qsr(ji,jj) = 0._wp                           ! qsr set to zero
          END_2D

DEC_dbg_RK3/src/OCE/TRD/trdglo.F90

@@ -201 +201 @@
          zkz  (:,:,:) = 0._wp
          zkepe(:,:,:) = 0._wp
    
-         CALL eos( ts(:,:,:,:,Kmm), rhd, rhop )       ! now potential density
+!!clem st dbg         CALL eos( ts(:,:,:,:,Kmm), rhd, rhop )       ! now potential density
 
          zcof = 0.5_wp / rho0             ! Density flux at w-point
          zkz(:,:,1) = 0._wp

DEC_dbg_RK3/src/TOP/trcini.F90

@@ -246 +246 @@
       !
       tr(:,:,:,:,Kaa) = 0._wp
       !
-      IF( ln_trcbc .AND. lltrcbc )  THEN 
-        CALL trc_bc_ini ( jptra, Kmm  )            ! set tracers Boundary Conditions
-        CALL trc_bc     ( nit000, Kmm, tr, Kaa )   ! tracers: surface and lateral Boundary Conditions
-      ENDIF
-      !
-      IF( ln_trcais ) CALL trc_ais_ini   ! set tracers from Antarctic Ice Sheet
-      !
       IF( ln_rsttr ) THEN              ! restart from a file
         !
         CALL trc_rst_read( Kbb, Kmm )

DEC_dbg_RK3/src/TOP/trcrst.F90

@@ -159 +159 @@
       ! prognostic variables 
       ! --------------------
 #if defined key_RK3
-      DO jn = 1, jptra      ! MLF : After time step (before the swap) put in TRN
-         CALL iom_rstput( kt, nitrst, numrtw, 'TRN'//ctrcnm(jn), tr(:,:,:,jn,Kaa) )
-      END DO
       DO jn = 1, jptra      ! RK3 : After time step (before the swap) put in TRB
          CALL iom_rstput( kt, nitrst, numrtw, 'TRB'//ctrcnm(jn), tr(:,:,:,jn,Kaa) )
       END DO

sette.sh

@@ -4 +4 @@
 MAIN_DIR=$(dirname $SETTE_DIR)
 export SETTE_TIMING='no'
 export NOT_USING_QCO='no'
+export USING_RK3='no'
 
 # Parse command-line arguments
 if [ $# -gt 0 ]; then
-  while getopts t:x:cshTq option; do
+  while getopts t:x:cshTqQ option; do
      case $option in
         c) export SETTE_CLEAN_CONFIGS='yes'
            export SETTE_SYNC_CONFIGS='yes'
@@ -36 +37 @@
            echo ""
            echo "key_qco and key_linssh will NOT be activated"
            echo "";;
+       Q) export USING_RK3='yes'
+           echo ""
+           echo "key_qco and key_RK3 will be activated"
+           echo "";;
         h | *) echo 'sette.sh with no arguments (in this case all configuration will be tested)'
                echo '-t "CFG1_to_test CFG2_to_test ..." to test some specific configurations'
                echo '-x "TEST_type TEST_type ..." to specify particular types of test (RESTART is mandatory)'

sette_reference-configurations.sh

@@ -143 +143 @@
    export DEL_KEYS="${DEL_KEYS} key_qco key_linssh"
 fi
 #
+if [ ${USING_RK3} == "yes" ]
+ then
+   export ADD_KEYS="${ADD_KEYS} key_qco key_RK3"
+fi
+#
 # Settings which control the use of stand alone servers (only relevant if using xios)
 #
 export NUM_XIOSERVERS=4

sette_test-cases.sh

@@ -140 +140 @@
    export DEL_KEYS="${DEL_KEYS} key_qco key_linssh"
 fi
 #
+if [ ${USING_RK3} == "yes" ]
+ then
+   export ADD_KEYS="${ADD_KEYS} key_qco key_RK3"
+fi
+#
 # Settings which control the use of stand alone servers (only relevant if using xios)
 #
 export NUM_XIOSERVERS=4

OCE/stp2d.F90

@@ -26 +26 @@
    USE sbc_ice , ONLY : snwice_mass, snwice_mass_b
    USE sbcapr         ! surface boundary condition: atmospheric pressure
    USE sbcwave,  ONLY : bhd_wave
+#if defined key_agrif
+   USE agrif_oce_interp
+#endif
 
    PRIVATE
 

OCE/stprk3_stg.F90

@@ -354 +354 @@
       !
 # if defined key_agrif
             CALL Agrif_tra                     !* AGRIF zoom boundaries
-            CALL Agrif_dyn( kt )
+            CALL Agrif_dyn( kstp )
 # endif
       !                                        ! local domain boundaries  (T-point, unchanged sign)
       CALL lbc_lnk_multi( 'stp_RK3_stg', uu(:,:,:,       Kaa), 'U', -1., vv(:,:,:       ,Kaa), 'V', -1.   &