Index: /NEMO/trunk/src/OCE/DIA/diawri.F90
===================================================================
--- /NEMO/trunk/src/OCE/DIA/diawri.F90	(revision 14178)
+++ /NEMO/trunk/src/OCE/DIA/diawri.F90	(revision 14179)
@@ -215,5 +215,5 @@
       ENDIF
 
-      CALL iom_put( "rhop", rhop(:,:,:) )          ! 3D potential density (sigma0)
+      IF( .NOT.lk_SWE )   CALL iom_put( "rhop", rhop(:,:,:) )          ! 3D potential density (sigma0)
 
       IF ( iom_use("taubot") ) THEN                ! bottom stress
Index: /NEMO/trunk/src/OCE/DOM/domqco.F90
===================================================================
--- /NEMO/trunk/src/OCE/DOM/domqco.F90	(revision 14178)
+++ /NEMO/trunk/src/OCE/DOM/domqco.F90	(revision 14179)
@@ -116,6 +116,11 @@
       !                    !== Set of all other vertical scale factors  ==!  (now and before)
       !                                ! Horizontal interpolation of e3t
+#if defined key_RK3
+      CALL dom_qco_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb), r3f(:,:) )
+      CALL dom_qco_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm)           )
+#else
       CALL dom_qco_r3c( ssh(:,:,Kbb), r3t(:,:,Kbb), r3u(:,:,Kbb), r3v(:,:,Kbb)           )
       CALL dom_qco_r3c( ssh(:,:,Kmm), r3t(:,:,Kmm), r3u(:,:,Kmm), r3v(:,:,Kmm), r3f(:,:) )
+#endif
       !
    END SUBROUTINE dom_qco_zgr
Index: /NEMO/trunk/src/OCE/IOM/restart.F90
===================================================================
--- /NEMO/trunk/src/OCE/IOM/restart.F90	(revision 14178)
+++ /NEMO/trunk/src/OCE/IOM/restart.F90	(revision 14179)
@@ -47,4 +47,5 @@
    !! * Substitutions
 #  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -169,4 +170,5 @@
          ENDIF
          !
+#if ! defined key_RK3
          CALL iom_rstput( kt, nitrst, numrow, 'sshn'   , ssh(:,:        ,Kmm) )     ! now fields
          CALL iom_rstput( kt, nitrst, numrow, 'un'     , uu(:,:,:       ,Kmm) )
@@ -177,4 +179,5 @@
             CALL iom_rstput( kt, nitrst, numrow, 'rhop', rhop                 )
          ENDIF
+#endif
       ENDIF
 
@@ -264,4 +267,5 @@
       INTEGER  ::   jk
       REAL(wp), DIMENSION(jpi, jpj, jpk) :: w3d
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   :: zgdept       ! 3D workspace for QCO
       !!----------------------------------------------------------------------
       !
@@ -279,5 +283,15 @@
       ENDIF
       !
-      !                             !*  Read Kmm fields
+#if defined key_RK3
+      !                             !*  Read Kbb fields   (NB: in RK3 Kmm = Kbb = Nbb)
+      IF(lwp) WRITE(numout,*) '           Kbb u, v and T-S fields read in the restart file'
+      CALL iom_get( numror, jpdom_auto, 'ub'   , uu(:,:,:       ,Kbb), cd_type = 'U', psgn = -1._wp )
+      CALL iom_get( numror, jpdom_auto, 'vb'   , vv(:,:,:       ,Kbb), cd_type = 'V', psgn = -1._wp )
+      IF( .NOT.lk_SWE ) THEN
+         CALL iom_get( numror, jpdom_auto, 'tb', ts(:,:,:,jp_tem,Kbb) )
+         CALL iom_get( numror, jpdom_auto, 'sb', ts(:,:,:,jp_sal,Kbb) )
+      ENDIF
+#else
+      !                             !*  Read Kmm fields   (MLF only)
       IF(lwp) WRITE(numout,*)    '           Kmm u, v and T-S fields read in the restart file'
       CALL iom_get( numror, jpdom_auto, 'un'   , uu(:,:,:       ,Kmm), cd_type = 'U', psgn = -1._wp )
@@ -288,5 +302,5 @@
       ENDIF
       !
-      IF( l_1st_euler ) THEN        !*  Euler restart
+      IF( l_1st_euler ) THEN        !*  Euler restart   (MLF only)
          IF(lwp) WRITE(numout,*) '           Kbb u, v and T-S fields set to Kmm values'
          uu(:,:,:  ,Kbb) = uu(:,:,:  ,Kmm)         ! all before fields set to now values
@@ -294,5 +308,5 @@
          IF( .NOT.lk_SWE ) ts(:,:,:,:,Kbb) = ts(:,:,:,:,Kmm)
          !
-      ELSE                          !* Leap frog restart
+      ELSE                          !* Leap frog restart   (MLF only)
          IF(lwp) WRITE(numout,*) '           Kbb u, v and T-S fields read in the restart file'
          CALL iom_get( numror, jpdom_auto, 'ub'   , uu(:,:,:       ,Kbb), cd_type = 'U', psgn = -1._wp )
@@ -303,4 +317,5 @@
          ENDIF
       ENDIF
+#endif
       !
       IF( .NOT.lk_SWE ) THEN
@@ -308,5 +323,14 @@
             CALL iom_get( numror, jpdom_auto, 'rhop'   , rhop )   ! now    potential density
          ELSE
+#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
@@ -347,19 +371,29 @@
          !                         !=============================!
          !
-         !                                     !*  Read ssh at Kmm
+#if defined key_RK3
+         !                                     !*  RK3: Read ssh at Kbb
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*)    '      Kbb sea surface height read in the restart file'
+         CALL iom_get( numror, jpdom_auto, 'sshb'   , ssh(:,:,Kbb) )
+         !
+         !                                     !*  RK3: Set ssh at Kmm for AGRIF
+         ssh(:,:,Kmm) = 0._wp
+#else
+         !                                     !*  MLF: Read ssh at Kmm
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*)    '      Kmm sea surface height read in the restart file'
          CALL iom_get( numror, jpdom_auto, 'sshn'   , ssh(:,:,Kmm) )
          !
-         IF( l_1st_euler ) THEN                !* Euler at first time-step
+         IF( l_1st_euler ) THEN                !*  MLF: Euler at first time-step
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '      Euler first time step : ssh(Kbb) = ssh(Kmm)'
             ssh(:,:,Kbb) = ssh(:,:,Kmm)
             !
-         ELSE                                  !* read ssh at Kbb
+         ELSE                                  !*  MLF: read ssh at Kbb
             IF(lwp) WRITE(numout,*)
             IF(lwp) WRITE(numout,*) '      Kbb sea surface height read in the restart file'
             CALL iom_get( numror, jpdom_auto, 'sshb', ssh(:,:,Kbb) )
          ENDIF
+#endif
          !                         !============================!
       ELSE                         !==  Initialize at "rest"  ==!
@@ -390,9 +424,13 @@
          ENDIF
          !
-         ssh(:,:,Kmm) = ssh(:,:,Kbb)           !* set now values from to before ones 
+#if defined key_RK3
+         ssh(:,:,Kmm) = 0._wp                  !* RK3: set Kmm to 0 for AGRIF
+#else
+         ssh(:,:,Kmm) = ssh(:,:,Kbb)           !* MLF: set now values from to before ones 
+#endif
       ENDIF
       !
       !                            !==========================!
-      ssh(:,:,Kaa) = 0._wp         !==  Set to 0 for agrif  ==!
+      ssh(:,:,Kaa) = 0._wp         !==  Set to 0 for AGRIF  ==!
       !                            !==========================!
       !
Index: /NEMO/trunk/src/SWE/stprk3.F90
===================================================================
--- /NEMO/trunk/src/SWE/stprk3.F90	(revision 14178)
+++ /NEMO/trunk/src/SWE/stprk3.F90	(revision 14179)
@@ -129,22 +129,4 @@
       CALL dyn_ldf( kstp, Nbb, Nbb, uu, vv, Nrhs )  ! lateral mixing
 #endif
-!!st       !
-!!st       DO_3D( 0,0, 0,0, 1,jpkm1 )
-!!st          !                                          ! horizontal pressure gradient
-!!st          uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nbb) - ssh(ji,jj,Nbb) ) * r1_e1u(ji,jj)
-!!st          vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nbb) - ssh(ji,jj,Nbb) ) * r1_e2v(ji,jj)
-!!st       END_3D
-!!st       !
-!!st #if defined key_RK3all 
-!!st       !                                             ! wind stress and layer friction
-!!st       z5_6 = 5._wp/6._wp
-!!st       DO_3D( 0, 0, 0, 0,1,jpkm1)
-!!st          uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + r1_rho0 * ( z5_6*utau_b(ji,jj) + (1._wp - z5_6)*utau(ji,jj) ) / e3u(ji,jj,jk,Nbb)   &
-!!st             &                                  - rn_rfr * uu(ji,jj,jk,Nbb)
-!!st          vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + r1_rho0 * ( z5_6*vtau_b(ji,jj) + (1._wp - z5_6)*vtau(ji,jj) ) / e3v(ji,jj,jk,Nbb)   &
-!!st             &                                  - rn_rfr * vv(ji,jj,jk,Nbb)
-!!st       END_3D
-!!st #endif
-!!st why not ?
       z5_6 = 5._wp/6._wp
       DO_3D( 0, 0, 0, 0, 1, jpkm1 )
@@ -163,5 +145,4 @@
          vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + zrhs_v
       END_3D
-!!st end
       !
       !                                 !==  Time stepping of ssh Eq.  ==!   (and update r3_Naa)
@@ -229,21 +210,4 @@
          vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + zrhs_v
       END_3D
-!!st       !
-!!st       DO_3D( 0, 0, 0, 0, 1, jpkm1 )
-!!st          !                                          ! horizontal pressure gradient
-!!st          uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) - grav * ( ssh(ji+1,jj,Nnn) - ssh(ji,jj,Nnn) ) * r1_e1u(ji,jj)
-!!st          vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) - grav * ( ssh(ji,jj+1,Nnn) - ssh(ji,jj,Nnn) ) * r1_e2v(ji,jj)
-!!st       END_3D
-!!st       !
-!!st #if defined key_RK3all
-!!st       !                                             ! wind stress and layer friction
-!!st       z3_4 = 3._wp/4._wp
-!!st       DO_3D( 0, 0, 0, 0,1,jpkm1)
-!!st          uu(ji,jj,jk,Nrhs) = uu(ji,jj,jk,Nrhs) + r1_rho0 * ( z3_4*utau_b(ji,jj) + (1._wp - z3_4)*utau(ji,jj) ) / e3u(ji,jj,jk,Nbb)   &
-!!st             &                                  - rn_rfr * uu(ji,jj,jk,Nbb)
-!!st          vv(ji,jj,jk,Nrhs) = vv(ji,jj,jk,Nrhs) + r1_rho0 * ( z3_4*vtau_b(ji,jj) + (1._wp - z3_4)*vtau(ji,jj) ) / e3v(ji,jj,jk,Nbb)   &
-!!st             &                                  - rn_rfr * vv(ji,jj,jk,Nbb)
-!!st       END_3D
-!!st #endif
       !
       !                                 !==  Time stepping of ssh Eq.  ==!   (and update r3_Naa)
@@ -344,16 +308,16 @@
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      ! diagnostics and outputs
+      ! diagnostics and outputs at Nbb (i.e. the just computed time step)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       
-      IF( ln_diacfl  )   CALL dia_cfl      ( kstp,      Nnn )      ! Courant number diagnostics
-                         CALL dia_wri      ( kstp,      Nnn )      ! ocean model: outputs
-      !
-      IF( lrst_oce   )   CALL rst_write    ( kstp, Nbb, Nnn )   ! write output ocean restart file
+      IF( ln_diacfl  )   CALL dia_cfl      ( kstp,      Nbb )      ! Courant number diagnostics
+                         CALL dia_wri      ( kstp,      Nbb )      ! ocean model: outputs
+      !
+      IF( lrst_oce   )   CALL rst_write    ( kstp, Nbb, Nbb )   ! write output ocean restart file
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Control
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL stp_ctl_SWE  ( kstp     , Nnn )
+                         CALL stp_ctl_SWE  ( kstp     , Nbb )
 
       IF( kstp == nit000 ) THEN                          ! 1st time step only
