Index: branches/UKMO/dev_r5107_hadgem3_mct/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/UKMO/dev_r5107_hadgem3_mct/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 5284)
+++ branches/UKMO/dev_r5107_hadgem3_mct/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 5568)
@@ -83,16 +83,17 @@
       IF ( kstp == (nit000 + 1) ) lk_agrif_fstep = .FALSE.
 # if defined key_iomput
-      IF( Agrif_Nbstepint() == 0 )   CALL iom_swap( "nemo" )
+      IF( Agrif_Nbstepint() == 0 )   CALL iom_swap( cxios_context )
 # endif
 #endif
                              indic = 0           ! reset to no error condition
       IF( kstp == nit000 ) THEN
-                      CALL iom_init( "nemo" )      ! iom_put initialization (must be done after nemo_init for AGRIF+XIOS+OASIS)
-         IF( ln_crs ) CALL iom_init( "nemo_crs" )  ! initialize context for coarse grid
+         ! must be done after nemo_init for AGRIF+XIOS+OASIS
+                      CALL iom_init(      cxios_context          )  ! iom_put initialization
+         IF( ln_crs ) CALL iom_init( TRIM(cxios_context)//"_crs" )  ! initialize context for coarse grid
       ENDIF
 
       IF( kstp /= nit000 )   CALL day( kstp )         ! Calendar (day was already called at nit000 in day_init)
-                             CALL iom_setkt( kstp - nit000 + 1, "nemo"     )   ! say to iom that we are at time step kstp
-      IF( ln_crs     )       CALL iom_setkt( kstp - nit000 + 1, "nemo_crs" )   ! say to iom that we are at time step kstp
+                             CALL iom_setkt( kstp - nit000 + 1,      cxios_context          )   ! tell iom we are at time step kstp
+      IF( ln_crs     )       CALL iom_setkt( kstp - nit000 + 1, TRIM(cxios_context)//"_crs" )   ! tell iom we are at time step kstp
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -100,8 +101,14 @@
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
       IF( lk_tide    )   CALL sbc_tide( kstp )
-      IF( lk_bdy     )   CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
-
+      IF( lk_bdy     )  THEN
+         IF( ln_apr_dyn) CALL sbc_apr( kstp )   ! bdy_dta needs ssh_ib 
+                         CALL bdy_dta ( kstp, time_offset=+1 )   ! update dynamic & tracer data at open boundaries
+      ENDIF
                          CALL sbc    ( kstp )         ! Sea Boundary Condition (including sea-ice)
                                                       ! clem: moved here for bdy ice purpose
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Update stochastic parameters and random T/S fluctuations
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+                        CALL sto_par( kstp )          ! Stochastic parameters
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -122,7 +129,7 @@
       IF( lk_zdfkpp  )   CALL zdf_kpp( kstp )            ! KPP closure scheme for Kz
       IF( lk_zdfcst  ) THEN                              ! Constant Kz (reset avt, avm[uv] to the background value)
-         avt (:,:,:) = rn_avt0 * tmask(:,:,:)
-         avmu(:,:,:) = rn_avm0 * umask(:,:,:)
-         avmv(:,:,:) = rn_avm0 * vmask(:,:,:)
+         avt (:,:,:) = rn_avt0 * wmask (:,:,:)
+         avmu(:,:,:) = rn_avm0 * wumask(:,:,:)
+         avmv(:,:,:) = rn_avm0 * wvmask(:,:,:)
       ENDIF
       IF( ln_rnf_mouth ) THEN                         ! increase diffusivity at rivers mouths
@@ -145,8 +152,13 @@
       !
       IF( lk_ldfslp ) THEN                            ! slope of lateral mixing
-                         CALL eos( tsb, rhd, gdept_0(:,:,:) )             ! before in situ density
-         IF( ln_zps )    CALL zps_hde( kstp, jpts, tsb, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient
-            &                                      rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &             !
-            &                               gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )  ! of t, s, rd at the last ocean level
+         IF(ln_sto_eos ) CALL sto_pts( tsn )          ! Random T/S fluctuations
+                         CALL eos( tsb, rhd, gdept_0(:,:,:) )               ! before in situ density
+         IF( ln_zps .AND. .NOT. ln_isfcav)                               &
+            &            CALL zps_hde    ( kstp, jpts, tsb, gtsu, gtsv,  &  ! Partial steps: before horizontal gradient
+            &                                          rhd, gru , grv    )  ! of t, s, rd at the last ocean level
+         IF( ln_zps .AND.       ln_isfcav)                               &
+            &            CALL zps_hde_isf( kstp, jpts, tsb, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
+            &                                          rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
+            &                                   gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the first ocean level
          IF( ln_traldf_grif ) THEN                           ! before slope for Griffies operator
                          CALL ldf_slp_grif( kstp )
@@ -158,8 +170,8 @@
       IF( lk_traldf_eiv )   CALL ldf_eiv( kstp )      ! eddy induced velocity coefficient
 #endif
-#if defined key_traldf_c3d && key_traldf_smag
+#if defined key_traldf_c3d && defined key_traldf_smag
                           CALL ldf_tra_smag( kstp )      ! eddy induced velocity coefficient
 #  endif
-#if defined key_dynldf_c3d && key_dynldf_smag
+#if defined key_dynldf_c3d && defined key_dynldf_smag
                           CALL ldf_dyn_smag( kstp )      ! eddy induced velocity coefficient
 #  endif
@@ -176,8 +188,13 @@
           ! Note that the computation of vertical velocity above, hence "after" sea level
           ! is necessary to compute momentum advection for the rhs of barotropic loop:
+            IF(ln_sto_eos ) CALL sto_pts( tsn )                             ! Random T/S fluctuations
                             CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) ) ! now in situ density for hpg computation
-            IF( ln_zps )    CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient
-               &                                      rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &             !
-               &                               gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )  ! of t, s, rd at the last ocean level
+            IF( ln_zps .AND. .NOT. ln_isfcav)                               &
+               &            CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
+               &                                          rhd, gru , grv    )  ! of t, s, rd at the last ocean level
+            IF( ln_zps .AND.       ln_isfcav)                               &
+               &            CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
+               &                                          rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
+               &                                   gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
 
                                   ua(:,:,:) = 0.e0             ! set dynamics trends to zero
@@ -208,15 +225,13 @@
       ! diagnostics and outputs             (ua, va, tsa used as workspace)
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( lk_floats  )   CALL flo_stp( kstp )         ! drifting Floats
-      IF( lk_diahth  )   CALL dia_hth( kstp )         ! Thermocline depth (20 degres isotherm depth)
-      IF( .NOT. lk_cpl ) CALL dia_fwb( kstp )         ! Fresh water budget diagnostics
-      IF( ln_diaptr  )   CALL dia_ptr( kstp )         ! Poleward TRansports diagnostics
-      IF( lk_diadct  )   CALL dia_dct( kstp )         ! Transports
-      IF( lk_diaar5  )   CALL dia_ar5( kstp )         ! ar5 diag
-      IF( lk_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis
-                         CALL dia_wri( kstp )         ! ocean model: outputs
-      !
-      IF( ln_crs     )   CALL crs_fld( kstp )         ! ocean model: online field coarsening & output
-
+      IF( lk_floats  )      CALL flo_stp( kstp )         ! drifting Floats
+      IF( lk_diahth  )      CALL dia_hth( kstp )         ! Thermocline depth (20 degres isotherm depth)
+      IF( .NOT. ln_cpl )    CALL dia_fwb( kstp )         ! Fresh water budget diagnostics
+      IF( lk_diadct  )      CALL dia_dct( kstp )         ! Transports
+      IF( lk_diaar5  )      CALL dia_ar5( kstp )         ! ar5 diag
+      IF( lk_diaharm )      CALL dia_harm( kstp )        ! Tidal harmonic analysis
+                            CALL dia_wri( kstp )         ! ocean model: outputs
+      !
+      IF( ln_crs     )      CALL crs_fld( kstp )         ! ocean model: online field coarsening & output
 
 #if defined key_top
@@ -244,4 +259,7 @@
       IF( lk_zdfkpp      )   CALL tra_kpp    ( kstp )       ! KPP non-local tracer fluxes
                              CALL tra_ldf    ( kstp )       ! lateral mixing
+
+      IF( ln_diaptr      )   CALL dia_ptr                   ! Poleward adv/ldf TRansports diagnostics
+
 #if defined key_agrif
       IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_tra          ! tracers sponge
@@ -252,14 +270,24 @@
          IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
                              CALL tra_nxt( kstp )                ! tracer fields at next time step
+            IF( ln_sto_eos ) CALL sto_pts( tsn )                 ! Random T/S fluctuations
                              CALL eos    ( tsa, rhd, rhop, fsdept_n(:,:,:) )  ! Time-filtered in situ density for hpg computation
-         IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsa, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient
-            &                                      rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &             !
-            &                               gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )  ! of t, s, rd at the last ocean level
+            IF( ln_zps .AND. .NOT. ln_isfcav)                                &
+               &             CALL zps_hde    ( kstp, jpts, tsa, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
+               &                                           rhd, gru , grv    )  ! of t, s, rd at the last ocean level
+            IF( ln_zps .AND.       ln_isfcav)                                &
+               &             CALL zps_hde_isf( kstp, jpts, tsa, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
+               &                                           rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
+               &                                    gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
       ELSE                                                  ! centered hpg  (eos then time stepping)
          IF ( .NOT. lk_dynspg_ts ) THEN                     ! eos already called in time-split case
+            IF( ln_sto_eos ) CALL sto_pts( tsn )    ! Random T/S fluctuations
                              CALL eos    ( tsn, rhd, rhop, fsdept_n(:,:,:) )  ! now in situ density for hpg computation
-         IF( ln_zps      )   CALL zps_hde( kstp, jpts, tsn, gtsu, gtsv,  &        ! Partial steps: before horizontal gradient
-         &                                      rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &             !
-         &                               gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi  )  ! of t, s, rd at the last ocean level
+         IF( ln_zps .AND. .NOT. ln_isfcav)                                   &
+               &             CALL zps_hde    ( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
+               &                                           rhd, gru , grv    )  ! of t, s, rd at the last ocean level
+         IF( ln_zps .AND.       ln_isfcav)                                   & 
+               &             CALL zps_hde_isf( kstp, jpts, tsn, gtsu, gtsv,  &    ! Partial steps for top cell (ISF)
+               &                                           rhd, gru , grv , aru , arv , gzu , gzv , ge3ru , ge3rv ,   &
+               &                                    gtui, gtvi, grui, grvi, arui, arvi, gzui, gzvi, ge3rui, ge3rvi    ) ! of t, s, rd at the last ocean level
          ENDIF
          IF( ln_zdfnpc   )   CALL tra_npc( kstp )                ! update after fields by non-penetrative convection
@@ -323,5 +351,5 @@
                  CALL iom_close( numror )     ! close input  ocean restart file
          IF(lwm) CALL FLUSH    ( numond )     ! flush output namelist oce
-         IF( lwm.AND.numoni /= -1 ) CALL FLUSH    ( numoni )     ! flush output namelist ice    
+         IF( lwm.AND.numoni /= -1 ) CALL FLUSH    ( numoni )     ! flush output namelist ice
       ENDIF
       IF( lrst_oce         )   CALL rst_write    ( kstp )   ! write output ocean restart file
@@ -330,10 +358,10 @@
       ! Coupled mode
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      IF( lk_cpl           )   CALL sbc_cpl_snd( kstp )     ! coupled mode : field exchanges
+      IF( lk_oasis         )   CALL sbc_cpl_snd( kstp )     ! coupled mode : field exchanges
       !
 #if defined key_iomput
       IF( kstp == nitend .OR. indic < 0 ) THEN 
-                      CALL iom_context_finalize( "nemo"     ) ! needed for XIOS+AGRIF
-         IF( ln_crs ) CALL iom_context_finalize( "nemo_crs" ) ! 
+                      CALL iom_context_finalize(      cxios_context          ) ! needed for XIOS+AGRIF
+         IF( ln_crs ) CALL iom_context_finalize( trim(cxios_context)//"_crs" ) ! 
       ENDIF
 #endif
