Index: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/step.F90
===================================================================
--- NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/step.F90	(revision 14219)
+++ NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/step.F90	(revision 14644)
@@ -33,4 +33,5 @@
    !!            4.1  !  2019-08  (A. Coward, D. Storkey) rewrite in preparation for new timestepping scheme
    !!----------------------------------------------------------------------
+
 #if defined key_qco   ||   defined key_linssh
    !!----------------------------------------------------------------------
@@ -41,5 +42,5 @@
 #else
    !!----------------------------------------------------------------------
-   !!   stp             : OPA system time-stepping
+   !!   stp             : OCE system time-stepping
    !!----------------------------------------------------------------------
    USE step_oce         ! time stepping definition modules
@@ -50,13 +51,10 @@
    PUBLIC   stp   ! called by nemogcm.F90
 
-   !!----------------------------------------------------------------------
-   !! time level indices
-   !!----------------------------------------------------------------------
-   INTEGER, PUBLIC :: Nbb, Nnn, Naa, Nrhs          !! used by nemo_init
+   !                                          !**  time level indices  **!
+   INTEGER, PUBLIC ::   Nbb, Nnn, Naa, Nrhs   !: used by nemo_init
 
    !! * Substitutions
 #  include "do_loop_substitute.h90"
 #  include "single_precision_substitute.h90"
-
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
@@ -76,5 +74,5 @@
       !!                     ***  ROUTINE stp  ***
       !!
-      !! ** Purpose : - Time stepping of OPA  (momentum and active tracer eqs.)
+      !! ** Purpose : - Time stepping of OCE  (momentum and active tracer eqs.)
       !!              - Time stepping of SI3 (dynamic and thermodynamic eqs.)
       !!              - Time stepping of TRC  (passive tracer eqs.)
@@ -90,6 +88,4 @@
       !!----------------------------------------------------------------------
       INTEGER ::   ji, jj, jk, jtile   ! dummy loop indice
-!!gm kcall can be removed, I guess
-      INTEGER ::   kcall        ! optional integer argument (dom_vvl_sf_nxt)
       !! ---------------------------------------------------------------------
 #if defined key_agrif
@@ -102,5 +98,5 @@
       ENDIF
       IF( kstp == nit000 + 1 )   lk_agrif_fstep = .FALSE.
-# if defined key_iomput
+# if defined key_xios
       IF( Agrif_Nbstepint() == 0 )   CALL iom_swap( cxios_context )
 # endif
@@ -113,7 +109,6 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       !
-      IF( l_1st_euler ) THEN
-         ! start or restart with Euler 1st time-step
-         rDt =  rn_Dt
+      IF( l_1st_euler ) THEN     ! start or restart with Euler 1st time-step
+         rDt   = rn_Dt   
          r1_Dt = 1._wp / rDt
       ENDIF
@@ -122,4 +117,5 @@
       ! update I/O and calendar
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+      !
       IF( kstp == nit000 ) THEN                       ! initialize IOM context (must be done after nemo_init for AGRIF+XIOS+OASIS)
                              CALL iom_init( cxios_context, ld_closedef=.FALSE. )   ! for model grid (including possible AGRIF zoom)
@@ -128,21 +124,26 @@
          IF( ln_crs      )   CALL iom_init( TRIM(cxios_context)//"_crs" )  ! for coarse grid
       ENDIF
-      IF((kstp == nitrst) .AND. lwxios) THEN
-         CALL iom_swap(      cw_ocerst_cxt          )
-         CALL iom_init_closedef(cw_ocerst_cxt)
-         CALL iom_setkt( kstp - nit000 + 1,      cw_ocerst_cxt          )
+      IF( kstp == nitrst .AND. lwxios ) THEN
+                             CALL iom_swap(                     cw_ocerst_cxt )
+                             CALL iom_init_closedef(            cw_ocerst_cxt )
+                             CALL iom_setkt( kstp - nit000 + 1, cw_ocerst_cxt )
 #if defined key_top
-         CALL iom_swap(      cw_toprst_cxt          )
-         CALL iom_init_closedef(cw_toprst_cxt)
-         CALL iom_setkt( kstp - nit000 + 1,      cw_toprst_cxt          )
-#endif
-      ENDIF
+                             CALL iom_swap(                     cw_toprst_cxt )
+                             CALL iom_init_closedef(            cw_toprst_cxt )
+                             CALL iom_setkt( kstp - nit000 + 1, cw_toprst_cxt )
+#endif
+      ENDIF
+      IF( kstp + nn_fsbc - 1 == nitrst .AND. lwxios ) THEN
 #if defined key_si3
-      IF(((kstp + nn_fsbc - 1) == nitrst) .AND. lwxios) THEN
-         CALL iom_swap(      cw_icerst_cxt          )
-         CALL iom_init_closedef(cw_icerst_cxt)
-         CALL iom_setkt( kstp - nit000 + 1,      cw_icerst_cxt          )
-      ENDIF
-#endif
+                             CALL iom_swap(                     cw_icerst_cxt )
+                             CALL iom_init_closedef(            cw_icerst_cxt )
+                             CALL iom_setkt( kstp - nit000 + 1, cw_icerst_cxt )
+#endif
+         IF( ln_abl      ) THEN
+                             CALL iom_swap(                     cw_ablrst_cxt )
+                             CALL iom_init_closedef(            cw_ablrst_cxt )
+                             CALL iom_setkt( kstp - nit000 + 1, cw_ablrst_cxt )
+         ENDIF
+      ENDIF
       IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init)
                              CALL iom_setkt( kstp - nit000 + 1,      cxios_context          )   ! tell IOM we are at time step kstp
@@ -161,6 +162,6 @@
       ! Update stochastic parameters and random T/S fluctuations
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-      IF( ln_sto_eos ) CALL sto_par( kstp )          ! Stochastic parameters
-      IF( ln_sto_eos ) CALL sto_pts( ts(:,:,:,:,Nnn)  )          ! Random T/S fluctuations
+      IF( ln_sto_eos )   CALL sto_par( kstp )                         ! Stochastic parameters
+      IF( ln_sto_eos )   CALL sto_pts( ts(:,:,:,:,Nnn)  )             ! Random T/S fluctuations
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -202,9 +203,10 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
-                            CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )    ! after ssh (includes call to div_hor)
-      IF( .NOT.ln_linssh )  CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn,      Naa )    ! after vertical scale factors
-                            CALL wzv           ( kstp, Nbb, Nnn, Naa, ww  )    ! now cross-level velocity
-      IF( ln_zad_Aimp )     CALL wAimp         ( kstp,      Nnn           )  ! Adaptive-implicit vertical advection partitioning
-                            CALL eos    ( CASTWP(ts(:,:,:,:,Nnn)), rhd, rhop, CASTWP(gdept(:,:,:,Nnn)) )  ! now in situ density for hpg computation
+                         CALL ssh_nxt       ( kstp, Nbb, Nnn, ssh, Naa )   ! after ssh (includes call to div_hor)
+      IF( .NOT.ln_linssh )   &
+                       & CALL dom_vvl_sf_nxt( kstp, Nbb, Nnn,      Naa )   ! after vertical scale factors
+                         CALL wzv           ( kstp, Nbb, Nnn, Naa, ww  )   ! now cross-level velocity
+      IF( ln_zad_Aimp )  CALL wAimp         ( kstp,      Nnn           )   ! Adaptive-implicit vertical advection partitioning
+                         CALL eos    ( ts(:,:,:,:,Nnn), rhd, rhop, gdept(:,:,:,Nnn) )  ! now in situ density for hpg computation
 
 
@@ -233,5 +235,5 @@
                             CALL dyn_zdf    ( kstp, Nbb, Nnn, Nrhs, uu, vv, Naa  )  ! vertical diffusion
       IF( ln_dynspg_ts ) THEN                                                       ! vertical scale factors and vertical velocity need to be updated
-                            CALL wzv        ( kstp, Nbb, Nnn, Naa, ww )             ! now cross-level velocity
+                            CALL wzv        ( kstp, Nbb, Nnn, Naa, ww )             ! Nnn cross-level velocity
          IF( ln_zad_Aimp )  CALL wAimp      ( kstp,      Nnn )                      ! Adaptive-implicit vertical advection partitioning
       ENDIF
@@ -261,5 +263,5 @@
       ! Passive Tracer Model
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-                         CALL trc_stp       ( kstp, Nbb, Nnn, Nrhs, Naa )  ! time-stepping
+                         CALL trc_stp    ( kstp, Nbb, Nnn, Nrhs, Naa )           ! time-stepping
 #endif
 
@@ -269,8 +271,8 @@
       ! Loop over tile domains
       DO jtile = 1, nijtile
-         IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile )
+         IF( ln_tile    )   CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile )
 
          DO_3D( 0, 0, 0, 0, 1, jpk )
-            ts(ji,jj,jk,:,Nrhs) = 0._wp                                         ! set tracer trends to zero
+            ts(ji,jj,jk,:,Nrhs) = 0._wp                                   ! set tracer trends to zero
          END_3D
 
@@ -287,19 +289,20 @@
 
 #if defined key_agrif
-         IF(.NOT. Agrif_Root()) THEN
-            IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )
-            CALL Agrif_Sponge_tra        ! tracers sponge
-         ENDIF
+      IF(.NOT. Agrif_Root() )   THEN
+         IF( ln_tile    )   CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = 0 )
+                            CALL Agrif_Sponge_tra        ! tracers sponge
+      ENDIF
 #endif
 
       ! TEMP: [tiling] Separate loop over tile domains (due to tra_adv workarounds for tiling)
       DO jtile = 1, nijtile
-         IF( ln_tile ) CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile )
+         IF( ln_tile    )   CALL dom_tile( ntsi, ntsj, ntei, ntej, ktile = jtile )
 
                             CALL tra_adv    ( kstp, Nbb, Nnn, ts, Nrhs )  ! hor. + vert. advection	==> RHS
          IF( ln_zdfmfc  )   CALL tra_mfc    ( kstp, Nbb,      ts, Nrhs )  ! Mass Flux Convection
-         IF( ln_zdfosm  )   CALL tra_osm    ( kstp,      Nnn, ts, Nrhs )  ! OSMOSIS non-local tracer fluxes ==> RHS
-         IF( lrst_oce .AND. ln_zdfosm ) &
-              &             CALL osm_rst    ( kstp,      Nnn, 'WRITE'  )  ! write OSMOSIS outputs + ww (so must do here) to restarts
+         IF( ln_zdfosm  ) THEN
+                            CALL tra_osm    ( kstp,      Nnn, ts, Nrhs )  ! OSMOSIS non-local tracer fluxes ==> RHS
+            IF( lrst_oce )  CALL osm_rst    ( kstp,      Nnn, 'WRITE'  )  ! write OSMOSIS outputs + ww (so must do here) to restarts
+         ENDIF
                             CALL tra_ldf    ( kstp, Nbb, Nnn, ts, Nrhs )  ! lateral mixing
 
@@ -365,10 +368,10 @@
       ! AGRIF update
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( Agrif_NbStepint() == 0 .AND. nstop == 0 ) THEN
-                         CALL Agrif_update_all( )                  ! Update all components
-      ENDIF
-
-#endif
-      IF( ln_diaobs .AND. nstop == 0 )  CALL dia_obs( kstp, Nnn )  ! obs-minus-model (assimilation) diags (after dynamics update)
+      IF( Agrif_NbStepint() == 0 .AND. nstop == 0 )   &
+         &               CALL Agrif_update_all( )                  ! Update all components
+
+#endif
+      IF( ln_diaobs .AND. nstop == 0 )   &
+         &               CALL dia_obs( kstp, Nnn )  ! obs-minus-model (assimilation) diags (after dynamics update)
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -377,5 +380,5 @@
       IF( kstp == nit000 ) THEN                          ! 1st time step only
                                         CALL iom_close( numror )   ! close input  ocean restart file
-         IF( lrxios )                   CALL iom_context_finalize(      cr_ocerst_cxt         )
+         IF( lrxios )                   CALL iom_context_finalize( cr_ocerst_cxt )
          IF(lwm)                        CALL FLUSH    ( numond )   ! flush output namelist oce
          IF(lwm .AND. numoni /= -1 )    CALL FLUSH    ( numoni )   ! flush output namelist ice (if exist)
@@ -387,5 +390,5 @@
       IF( lk_oasis .AND. nstop == 0 )   CALL sbc_cpl_snd( kstp, Nbb, Nnn )     ! coupled mode : field exchanges
       !
-#if defined key_iomput
+#if defined key_xios
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Finalize contextes if end of simulation or error detected
@@ -398,5 +401,5 @@
       !
       IF( l_1st_euler ) THEN         ! recover Leap-frog timestep
-         rDt = 2._wp * rn_Dt
+         rDt   = 2._wp * rn_Dt
          r1_Dt = 1._wp / rDt
          l_1st_euler = .FALSE.
