Index: /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg	(revision 14912)
+++ /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg	(revision 14913)
@@ -42,7 +42,4 @@
 &namtile        !   parameters of the tiling
 !-----------------------------------------------------------------------
-   ln_tile = .false.     !  Use tiling (T) or not (F)
-   nn_ltile_i = 10       !  Length of tiles in i
-   nn_ltile_j = 10       !  Length of tiles in j
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/cfgs/SHARED/namelist_ref
===================================================================
--- /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/cfgs/SHARED/namelist_ref	(revision 14912)
+++ /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/cfgs/SHARED/namelist_ref	(revision 14913)
@@ -99,5 +99,5 @@
 !-----------------------------------------------------------------------
    ln_tile = .false.     !  Use tiling (T) or not (F)
-   nn_ltile_i = 10       !  Length of tiles in i
+   nn_ltile_i = 99999    !  Length of tiles in i
    nn_ltile_j = 10       !  Length of tiles in j
 /
Index: /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/ICE/icestp.F90
===================================================================
--- /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/ICE/icestp.F90	(revision 14912)
+++ /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/ICE/icestp.F90	(revision 14913)
@@ -460,4 +460,5 @@
             qcn_ice    (ji,jj,jl) = 0._wp   ! conductive flux (ln_cndflx=T & ln_cndemule=T)
             qtr_ice_bot(ji,jj,jl) = 0._wp   ! part of solar radiation transmitted through the ice needed at least for outputs
+            qml_ice    (ji,jj,jl) = 0._wp   ! surface melt heat flux
             ! Melt pond surface melt diagnostics (mv - more efficient: grouped into one water volume flux)
             dh_i_sum_2d(ji,jj,jl) = 0._wp
Index: /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/ICE/icethd.F90
===================================================================
--- /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/ICE/icethd.F90	(revision 14912)
+++ /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/ICE/icethd.F90	(revision 14913)
@@ -536,4 +536,5 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), qcn_ice_bot_1d(1:npti), qcn_ice_bot(:,:,kl) )
          CALL tab_1d_2d( npti, nptidx(1:npti), qcn_ice_top_1d(1:npti), qcn_ice_top(:,:,kl) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), qml_ice_1d    (1:npti), qml_ice    (:,:,kl) )
          ! extensive variables
          CALL tab_1d_2d( npti, nptidx(1:npti), v_i_1d (1:npti), v_i (:,:,kl) )
Index: /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfiwm.F90
===================================================================
--- /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfiwm.F90	(revision 14912)
+++ /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfiwm.F90	(revision 14913)
@@ -144,17 +144,11 @@
       ! Set to zero the 1st and last vertical levels of appropriate variables
       IF( iom_use("emix_iwm") ) THEN
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
-            zemx_iwm (ji,jj,1) = 0._wp   ;   zemx_iwm (ji,jj,jpk) = 0._wp
-         END_2D
+         zemx_iwm(:,:,:) = 0._wp
       ENDIF
       IF( iom_use("av_ratio") ) THEN
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
-            zav_ratio(ji,jj,1) = 0._wp   ;   zav_ratio(ji,jj,jpk) = 0._wp
-         END_2D
+         zav_ratio(:,:,:) = 0._wp
       ENDIF
       IF( iom_use("av_wave") .OR. sn_cfctl%l_prtctl ) THEN
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
-            zav_wave (ji,jj,1) = 0._wp   ;   zav_wave (ji,jj,jpk) = 0._wp
-         END_2D
+         zav_wave(:,:,:) = 0._wp
       ENDIF
       !
@@ -366,18 +360,9 @@
       IF( iom_use("bflx_iwm") .OR. iom_use("pcmap_iwm") ) THEN
          ALLOCATE( z2d(A2D(nn_hls)) , z3d(A2D(nn_hls),jpk) )
+         ! Initialisation for iom_put
          z2d(:,:) = 0._wp ; z3d(:,:,:) = 0._wp
-
-         ! Initialisation for iom_put
-         DO_2D( 0, 0, 0, 0 )
-            z3d(ji,jj,1) = 0._wp   ;   z3d(ji,jj,jpk) = 0._wp
-         END_2D
 
          DO_3D( 0, 0, 0, 0, 2, jpkm1 )
             z3d(ji,jj,jk) = MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk)
-         END_3D
-         DO_2D( 0, 0, 0, 0 )
-            z2d(ji,jj) = 0._wp
-         END_2D
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 ) 
             z2d(ji,jj) = z2d(ji,jj) + e3w(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * wmask(ji,jj,jk)
          END_3D
Index: /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90
===================================================================
--- /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90	(revision 14912)
+++ /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90	(revision 14913)
@@ -264,22 +264,22 @@
       zdf_osm_alloc = zdf_osm_alloc + ierr
       !
-      ALLOCATE( etmean(A2D((nn_hls-1)),jpk), dh(jpi,jpj), r1_ft(A2D((nn_hls-1))), STAT=ierr )
+      ALLOCATE( etmean(A2D(nn_hls-1),jpk), dh(jpi,jpj), r1_ft(A2D(nn_hls-1)), STAT=ierr )
       zdf_osm_alloc = zdf_osm_alloc + ierr
       !
-      ALLOCATE( nbld(jpi,jpj), nmld(A2D((nn_hls-1))), STAT=ierr )
+      ALLOCATE( nbld(jpi,jpj), nmld(A2D(nn_hls-1)), STAT=ierr )
       zdf_osm_alloc = zdf_osm_alloc + ierr
       !
-      ALLOCATE( n_ddh(A2D((nn_hls-1))), STAT=ierr )
+      ALLOCATE( n_ddh(A2D(nn_hls-1)), STAT=ierr )
       zdf_osm_alloc = zdf_osm_alloc + ierr
       !
-      ALLOCATE( l_conv(A2D((nn_hls-1))), l_shear(A2D((nn_hls-1))), l_coup(A2D((nn_hls-1))), l_pyc(A2D((nn_hls-1))),   &
-         &      l_flux(A2D((nn_hls-1))), l_mle(A2D((nn_hls-1))), STAT=ierr )
+      ALLOCATE( l_conv(A2D(nn_hls-1)), l_shear(A2D(nn_hls-1)), l_coup(A2D(nn_hls-1)), l_pyc(A2D(nn_hls-1)),   &
+         &      l_flux(A2D(nn_hls-1)), l_mle(A2D(nn_hls-1)),   STAT=ierr )
       zdf_osm_alloc = zdf_osm_alloc + ierr
       !
-      ALLOCATE( swth0(A2D((nn_hls-1))),  sws0(A2D((nn_hls-1))),      swb0(A2D((nn_hls-1))),      suw0(A2D((nn_hls-1))),      &
-         &      sustar(A2D((nn_hls-1))), scos_wind(A2D((nn_hls-1))), ssin_wind(A2D((nn_hls-1))), swthav(A2D((nn_hls-1))),    &
-         &      swsav(A2D((nn_hls-1))),  swbav(A2D((nn_hls-1))),     sustke(A2D((nn_hls-1))),    dstokes(A2D((nn_hls-1))),   &
-         &      swstrl(A2D((nn_hls-1))), swstrc(A2D((nn_hls-1))),    sla(A2D((nn_hls-1))),       svstr(A2D((nn_hls-1))),     &
-         &      shol(A2D((nn_hls-1))),   STAT=ierr )
+      ALLOCATE( swth0(A2D(nn_hls-1)),  sws0(A2D(nn_hls-1)),      swb0(A2D(nn_hls-1)),      suw0(A2D(nn_hls-1)),      &
+         &      sustar(A2D(nn_hls-1)), scos_wind(A2D(nn_hls-1)), ssin_wind(A2D(nn_hls-1)), swthav(A2D(nn_hls-1)),    &
+         &      swsav(A2D(nn_hls-1)),  swbav(A2D(nn_hls-1)),     sustke(A2D(nn_hls-1)),    dstokes(A2D(nn_hls-1)),   &
+         &      swstrl(A2D(nn_hls-1)), swstrc(A2D(nn_hls-1)),    sla(A2D(nn_hls-1)),       svstr(A2D(nn_hls-1)),     &
+         &      shol(A2D(nn_hls-1)),   STAT=ierr )
       zdf_osm_alloc = zdf_osm_alloc + ierr
       !
@@ -359,36 +359,36 @@
       REAL(wp) ::   zesh2, zri, zfri   ! Interior Richardson mixing
       !! Scales
-      REAL(wp), DIMENSION(A2D((nn_hls-1)))  ::   zrad0       ! Surface solar temperature flux (deg m/s)
-      REAL(wp), DIMENSION(A2D((nn_hls-1)))  ::   zradh       ! Radiative flux at bl base (Buoyancy units)
-      REAL(wp)                              ::   zradav      ! Radiative flux, bl average (Buoyancy Units)
-      REAL(wp)                              ::   zvw0        ! Surface v-momentum flux
-      REAL(wp), DIMENSION(A2D((nn_hls-1)))  ::   zwb0tot     ! Total surface buoyancy flux including insolation
-      REAL(wp), DIMENSION(A2D((nn_hls-1)))  ::   zwb_ent     ! Buoyancy entrainment flux
-      REAL(wp), DIMENSION(A2D((nn_hls-1)))  ::   zwb_min
-      REAL(wp), DIMENSION(A2D((nn_hls-1)))  ::   zwb_fk_b    ! MLE buoyancy flux averaged over OSBL
-      REAL(wp), DIMENSION(A2D((nn_hls-1)))  ::   zwb_fk      ! Max MLE buoyancy flux
-      REAL(wp), DIMENSION(A2D((nn_hls-1)))  ::   zdiff_mle   ! Extra MLE vertical diff
-      REAL(wp), DIMENSION(A2D((nn_hls-1)))  ::   zvel_mle    ! Velocity scale for dhdt with stable ML and FK
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zrad0       ! Surface solar temperature flux (deg m/s)
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zradh       ! Radiative flux at bl base (Buoyancy units)
+      REAL(wp)                           ::   zradav      ! Radiative flux, bl average (Buoyancy Units)
+      REAL(wp)                           ::   zvw0        ! Surface v-momentum flux
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zwb0tot     ! Total surface buoyancy flux including insolation
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zwb_ent     ! Buoyancy entrainment flux
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zwb_min
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zwb_fk_b    ! MLE buoyancy flux averaged over OSBL
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zwb_fk      ! Max MLE buoyancy flux
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zdiff_mle   ! Extra MLE vertical diff
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zvel_mle    ! Velocity scale for dhdt with stable ML and FK
       !! Mixed-layer variables
-      INTEGER,  DIMENSION(A2D((nn_hls-1))) ::   jk_nlev   ! Number of levels
-      INTEGER,  DIMENSION(A2D((nn_hls-1))) ::   jk_ext   ! Offset for external level
-      !!
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zhbl   ! BL depth - grid
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zhml   ! ML depth - grid
-      !!
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zhmle   ! MLE depth - grid
-      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zmld    ! ML depth on grid
-      !!
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zdh                          ! Pycnocline depth - grid
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zdhdt                        ! BL depth tendency
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zdtdz_bl_ext, zdsdz_bl_ext   ! External temperature/salinity gradients
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zdbdz_bl_ext                 ! External buoyancy gradients
-      REAL(wp), DIMENSION(A2D(nn_hls))     ::   zdtdx, zdtdy, zdsdx, zdsdy   ! Horizontal gradients for Fox-Kemper parametrization
-      !!
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zdbds_mle   ! Magnitude of horizontal buoyancy gradient
+      INTEGER,  DIMENSION(A2D(nn_hls-1)) ::   jk_nlev  ! Number of levels
+      INTEGER,  DIMENSION(A2D(nn_hls-1)) ::   jk_ext   ! Offset for external level
+      !!
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zhbl   ! BL depth - grid
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zhml   ! ML depth - grid
+      !!
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zhmle   ! MLE depth - grid
+      REAL(wp), DIMENSION(A2D(nn_hls))   ::   zmld    ! ML depth on grid
+      !!
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zdh                          ! Pycnocline depth - grid
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zdhdt                        ! BL depth tendency
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zdtdz_bl_ext, zdsdz_bl_ext   ! External temperature/salinity gradients
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zdbdz_bl_ext                 ! External buoyancy gradients
+      REAL(wp), DIMENSION(A2D(nn_hls))   ::   zdtdx, zdtdy, zdsdx, zdsdy   ! Horizontal gradients for Fox-Kemper parametrization
+      !!
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zdbds_mle   ! Magnitude of horizontal buoyancy gradient
       !! Flux-gradient relationship variables
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zshear   ! Shear production
-      !!
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zhbl_t   ! Holds boundary layer depth updated by full timestep
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zshear   ! Shear production
+      !!
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zhbl_t   ! Holds boundary layer depth updated by full timestep
       !! For calculating Ri#-dependent mixing
       REAL(wp), DIMENSION(A2D(nn_hls)) ::   z2du     ! u-shear^2
@@ -399,6 +399,6 @@
       REAL(wp)                                 ::   zz0, zz1, zfac
       REAL(wp)                                 ::   zus_x, zus_y     ! Temporary Stokes drift
-      REAL(wp), DIMENSION(A2D((nn_hls-1)),jpk) ::   zviscos          ! Viscosity
-      REAL(wp), DIMENSION(A2D((nn_hls-1)),jpk) ::   zdiffut          ! t-diffusivity
+      REAL(wp), DIMENSION(A2D(nn_hls-1),jpk)   ::   zviscos          ! Viscosity
+      REAL(wp), DIMENSION(A2D(nn_hls-1),jpk)   ::   zdiffut          ! t-diffusivity
       REAL(wp)                                 ::   zabsstke
       REAL(wp)                                 ::   zsqrtpi, z_two_thirds, zthickness
@@ -644,11 +644,11 @@
       !
       ! Averages over well-mixed and boundary layer, note BL averages use jk_ext=2 everywhere
-      jk_nlev(:,:) = nbld(A2D((nn_hls-1)))
+      jk_nlev(:,:) = nbld(A2D(nn_hls-1))
       jk_ext(:,:) = 1   ! ag 19/03
       CALL zdf_osm_vertical_average( Kbb,      Kmm,      jk_nlev,  av_t_bl,  av_s_bl,    &
          &                           av_b_bl,  av_u_bl,  av_v_bl,  jk_ext,   av_dt_bl,   &
          &                           av_ds_bl, av_db_bl, av_du_bl, av_dv_bl )
-      jk_nlev(:,:) = nmld(A2D((nn_hls-1))) - 1
-      jk_ext(:,:) = nbld(A2D((nn_hls-1))) - nmld(A2D((nn_hls-1))) + jk_ext(:,:) + 1   ! ag 19/03
+      jk_nlev(:,:) = nmld(A2D(nn_hls-1)) - 1
+      jk_ext(:,:) = nbld(A2D(nn_hls-1)) - nmld(A2D(nn_hls-1)) + jk_ext(:,:) + 1   ! ag 19/03
       CALL zdf_osm_vertical_average( Kbb,      Kmm,      jk_nlev,  av_t_ml,  av_s_ml,    &
          &                           av_b_ml,  av_u_ml,  av_v_ml,  jk_ext,   av_dt_ml,   &
@@ -672,5 +672,5 @@
             IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN( mbkt(ji,jj), jk)
          END_3D
-         jk_nlev(:,:) = mld_prof(A2D((nn_hls-1)))
+         jk_nlev(:,:) = mld_prof(A2D(nn_hls-1))
          CALL zdf_osm_vertical_average( Kbb,      Kmm,      jk_nlev,  av_t_mle, av_s_mle,   &
             &                           av_b_mle, av_u_mle, av_v_mle )
@@ -700,5 +700,5 @@
       !
       !! External gradient below BL needed both with and w/o FK
-      jk_ext(:,:) = nbld(A2D((nn_hls-1))) + 1
+      jk_ext(:,:) = nbld(A2D(nn_hls-1)) + 1
       CALL zdf_osm_external_gradients( Kmm, jk_ext, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext )   ! ag 19/03
       !
@@ -721,11 +721,11 @@
       !
       ! Recalculate bl averages using jk_ext & ml averages .... note no rotation of u & v here..
-      jk_nlev(:,:) = nbld(A2D((nn_hls-1)))
+      jk_nlev(:,:) = nbld(A2D(nn_hls-1))
       jk_ext(:,:) = 1   ! ag 19/03
       CALL zdf_osm_vertical_average( Kbb,      Kmm,      jk_nlev,  av_t_bl,  av_s_bl,    &
          &                           av_b_bl,  av_u_bl,  av_v_bl,  jk_ext,   av_dt_bl,   &
          &                           av_ds_bl, av_db_bl, av_du_bl, av_dv_bl )
-      jk_nlev(:,:) = nmld(A2D((nn_hls-1))) - 1
-      jk_ext(:,:) = nbld(A2D((nn_hls-1))) - nmld(A2D((nn_hls-1))) + jk_ext(:,:) + 1   ! ag 19/03
+      jk_nlev(:,:) = nmld(A2D(nn_hls-1)) - 1
+      jk_ext(:,:) = nbld(A2D(nn_hls-1)) - nmld(A2D(nn_hls-1)) + jk_ext(:,:) + 1   ! ag 19/03
       CALL zdf_osm_vertical_average( Kbb,      Kmm,      jk_nlev,  av_t_ml,  av_s_ml,    &
          &                           av_b_ml,  av_u_ml,  av_v_ml,  jk_ext,   av_dt_ml,   &
@@ -768,5 +768,5 @@
       !
       ! Recalculate BL averages and differences using new BL depth
-      jk_nlev(:,:) = nbld(A2D((nn_hls-1)))
+      jk_nlev(:,:) = nbld(A2D(nn_hls-1))
       jk_ext(:,:) = 1   ! ag 19/03
       CALL zdf_osm_vertical_average( Kbb,      Kmm,      jk_nlev,  av_t_bl,  av_s_bl,    &
@@ -800,11 +800,11 @@
       !      jk_ext = nbld - nmld + 1
       ! Recalculate ML averages and differences using new ML depth
-      jk_nlev(:,:) = nmld(A2D((nn_hls-1))) - 1
-      jk_ext(:,:) = nbld(A2D((nn_hls-1))) - nmld(A2D((nn_hls-1))) + jk_ext(:,:) + 1   ! ag 19/03
+      jk_nlev(:,:) = nmld(A2D(nn_hls-1)) - 1
+      jk_ext(:,:) = nbld(A2D(nn_hls-1)) - nmld(A2D(nn_hls-1)) + jk_ext(:,:) + 1   ! ag 19/03
       CALL zdf_osm_vertical_average( Kbb,      Kmm,      jk_nlev,  av_t_ml,  av_s_ml,    &
          &                           av_b_ml,  av_u_ml,  av_v_ml,  jk_ext,   av_dt_ml,   &
          &                           av_ds_ml, av_db_ml, av_du_ml, av_dv_ml )
       !
-      jk_ext(:,:) = nbld(A2D((nn_hls-1))) + 1
+      jk_ext(:,:) = nbld(A2D(nn_hls-1)) + 1
       CALL zdf_osm_external_gradients( Kmm, jk_ext, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext )
       ! Rotate mean currents and changes onto wind aligned co-ordinates
@@ -834,5 +834,5 @@
       !
       ! Rotate non-gradient velocity terms back to model reference frame
-      jk_nlev(:,:) = nbld(A2D((nn_hls-1)))
+      jk_nlev(:,:) = nbld(A2D(nn_hls-1))
       CALL zdf_osm_velocity_rotation( ghamu, ghamv, .FALSE.,  2, jk_nlev )
       !
@@ -1038,18 +1038,18 @@
       !!              knlev+kp_ext
       !!----------------------------------------------------------------------
-      INTEGER,                              INTENT(in   )           ::   Kbb, Kmm   ! Ocean time-level indices
-      INTEGER,  DIMENSION(A2D((nn_hls-1))), INTENT(in   )           ::   knlev      ! Number of levels to average over.
-      REAL(wp), DIMENSION(jpi,jpj),         INTENT(  out)           ::   pt, ps     ! Average temperature and salinity
-      REAL(wp), DIMENSION(jpi,jpj),         INTENT(  out)           ::   pb         ! Average buoyancy
-      REAL(wp), DIMENSION(jpi,jpj),         INTENT(  out)           ::   pu, pv     ! Average current components
-      INTEGER,  DIMENSION(A2D((nn_hls-1))), INTENT(in   ), OPTIONAL ::   kp_ext     ! External-level offsets
-      REAL(wp), DIMENSION(jpi,jpj),         INTENT(  out), OPTIONAL ::   pdt        ! Difference between average temperature,
-      REAL(wp), DIMENSION(jpi,jpj),         INTENT(  out), OPTIONAL ::   pds        !    salinity,
-      REAL(wp), DIMENSION(jpi,jpj),         INTENT(  out), OPTIONAL ::   pdb        !    buoyancy, and
-      REAL(wp), DIMENSION(jpi,jpj),         INTENT(  out), OPTIONAL ::   pdu, pdv   !    velocity components and the OSBL
+      INTEGER,                            INTENT(in   )           ::   Kbb, Kmm   ! Ocean time-level indices
+      INTEGER,  DIMENSION(A2D(nn_hls-1)), INTENT(in   )           ::   knlev      ! Number of levels to average over.
+      REAL(wp), DIMENSION(jpi,jpj),       INTENT(  out)           ::   pt, ps     ! Average temperature and salinity
+      REAL(wp), DIMENSION(jpi,jpj),       INTENT(  out)           ::   pb         ! Average buoyancy
+      REAL(wp), DIMENSION(jpi,jpj),       INTENT(  out)           ::   pu, pv     ! Average current components
+      INTEGER,  DIMENSION(A2D(nn_hls-1)), INTENT(in   ), OPTIONAL ::   kp_ext     ! External-level offsets
+      REAL(wp), DIMENSION(jpi,jpj),       INTENT(  out), OPTIONAL ::   pdt        ! Difference between average temperature,
+      REAL(wp), DIMENSION(jpi,jpj),       INTENT(  out), OPTIONAL ::   pds        !    salinity,
+      REAL(wp), DIMENSION(jpi,jpj),       INTENT(  out), OPTIONAL ::   pdb        !    buoyancy, and
+      REAL(wp), DIMENSION(jpi,jpj),       INTENT(  out), OPTIONAL ::   pdu, pdv   !    velocity components and the OSBL
       !!
       INTEGER                              ::   jk, jkflt, jkmax, ji, jj   ! Loop indices
       INTEGER                              ::   ibld_ext                   ! External-layer index
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zthick                     ! Layer thickness
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zthick                     ! Layer thickness
       REAL(wp)                             ::   zthermal                   ! Thermal expansion coefficient
       REAL(wp)                             ::   zbeta                      ! Haline contraction coefficient
@@ -1174,8 +1174,8 @@
       !!
       !!----------------------------------------------------------------------      
-      REAL(wp),           INTENT(inout), DIMENSION(jpi,jpj,jpk)     ::   pu, pv   ! Components of current
-      LOGICAL,  OPTIONAL, INTENT(in   )                             ::   fwd      ! Forward (default) or reverse rotation
-      INTEGER,  OPTIONAL, INTENT(in   )                             ::   ktop     ! Minimum depth index
-      INTEGER,  OPTIONAL, INTENT(in   ), DIMENSION(A2D((nn_hls-1))) ::   knlev    ! Array of maximum depth indices
+      REAL(wp),           INTENT(inout), DIMENSION(jpi,jpj,jpk)   ::   pu, pv   ! Components of current
+      LOGICAL,  OPTIONAL, INTENT(in   )                           ::   fwd      ! Forward (default) or reverse rotation
+      INTEGER,  OPTIONAL, INTENT(in   )                           ::   ktop     ! Minimum depth index
+      INTEGER,  OPTIONAL, INTENT(in   ), DIMENSION(A2D(nn_hls-1)) ::   knlev    ! Array of maximum depth indices
       !!
       INTEGER  ::   ji, jj, jk, jktop, jkmax   ! Loop indices
@@ -1221,18 +1221,18 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER,                              INTENT(in   ) ::   Kmm       ! Ocean time-level index
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(  out) ::   pwb_ent   ! Buoyancy fluxes at base
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(  out) ::   pwb_min   !    of well-mixed layer
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(  out) ::   pshear    ! Production of TKE due to shear across the pycnocline
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   phbl      ! BL depth
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   phml      ! ML depth
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pdh       ! Pycnocline depth
+      INTEGER,                            INTENT(in   ) ::   Kmm       ! Ocean time-level index
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(  out) ::   pwb_ent   ! Buoyancy fluxes at base
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(  out) ::   pwb_min   !    of well-mixed layer
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(  out) ::   pshear    ! Production of TKE due to shear across the pycnocline
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   phbl      ! BL depth
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   phml      ! ML depth
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   pdh       ! Pycnocline depth
       !!
       INTEGER :: jj, ji   ! Loop indices
       !!
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zekman
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zri_p, zri_b   ! Richardson numbers
-      REAL(wp)                             ::   zshear_u, zshear_v, zwb_shr
-      REAL(wp)                             ::   zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zekman
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zri_p, zri_b   ! Richardson numbers
+      REAL(wp)                           ::   zshear_u, zshear_v, zwb_shr
+      REAL(wp)                           ::   zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes
       !!
       REAL(wp), PARAMETER ::   pp_a_shr         = 0.4_wp,  pp_b_shr    = 6.5_wp,  pp_a_wb_s = 0.8_wp
@@ -1269,6 +1269,6 @@
          pshear(ji,jj) = 0.0_wp
       END_2D
-      zekman(:,:) = EXP( -1.0_wp * pp_ek * ABS( ff_t(A2D((nn_hls-1))) ) * phbl(A2D((nn_hls-1))) /   &
-         &               MAX( sustar(A2D((nn_hls-1))), 1.e-8 ) )
+      zekman(:,:) = EXP( -1.0_wp * pp_ek * ABS( ff_t(A2D(nn_hls-1)) ) * phbl(A2D(nn_hls-1)) /   &
+         &               MAX( sustar(A2D(nn_hls-1)), 1.e-8 ) )
       !
       DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
@@ -1384,8 +1384,8 @@
       !!
       !!----------------------------------------------------------------------   
-      INTEGER,                              INTENT(in   ) ::   Kmm            ! Ocean time-level index
-      INTEGER,  DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   kbase          ! OSBL base layer index
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(  out) ::   pdtdz, pdsdz   ! External gradients of temperature, salinity
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(  out) ::   pdbdz          !    and buoyancy
+      INTEGER,                            INTENT(in   ) ::   Kmm            ! Ocean time-level index
+      INTEGER,  DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   kbase          ! OSBL base layer index
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(  out) ::   pdtdz, pdsdz   ! External gradients of temperature, salinity
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(  out) ::   pdbdz          !    and buoyancy
       !!
       INTEGER  ::   ji, jj, jkb, jkb1
@@ -1429,13 +1429,13 @@
       !!
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(  out) ::   pdhdt          ! Rate of change of hbl
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   phbl           ! BL depth
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pdh            ! Pycnocline depth
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pwb_min
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(  out) ::   pwb_fk_b       ! MLE buoyancy flux averaged over OSBL
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pwb_fk         ! Max MLE buoyancy flux
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pvel_mle       ! Vvelocity scale for dhdt with stable ML and FK
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(  out) ::   pdhdt          ! Rate of change of hbl
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   phbl           ! BL depth
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   pdh            ! Pycnocline depth
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   pwb_min
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(  out) ::   pwb_fk_b       ! MLE buoyancy flux averaged over OSBL
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   pwb_fk         ! Max MLE buoyancy flux
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   pvel_mle       ! Vvelocity scale for dhdt with stable ML and FK
       !!
       INTEGER  ::   jj, ji
@@ -1598,10 +1598,10 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER,                              INTENT(in   ) ::   Kmm        ! Ocean time-level index
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) ::   pdhdt      ! Rates of change of hbl
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) ::   phbl       ! BL depth
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   phbl_t     ! BL depth
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pwb_ent    ! Buoyancy entrainment flux
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pwb_fk_b   ! MLE buoyancy flux averaged over OSBL
+      INTEGER,                            INTENT(in   ) ::   Kmm        ! Ocean time-level index
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(inout) ::   pdhdt      ! Rates of change of hbl
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(inout) ::   phbl       ! BL depth
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   phbl_t     ! BL depth
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   pwb_ent    ! Buoyancy entrainment flux
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   pwb_fk_b   ! MLE buoyancy flux averaged over OSBL
       !!
       INTEGER  ::   jk, jj, ji, jm
@@ -1698,12 +1698,12 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER,                              INTENT(in   ) ::   Kmm            ! Ocean time-level index
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) ::   pdh            ! Pycnocline thickness
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) ::   phml           ! ML depth
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pdhdt          ! BL depth tendency
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   phbl           ! BL depth
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pwb_fk_b       ! MLE buoyancy flux averaged over OSBL
+      INTEGER,                            INTENT(in   ) ::   Kmm            ! Ocean time-level index
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(inout) ::   pdh            ! Pycnocline thickness
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(inout) ::   phml           ! ML depth
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   pdhdt          ! BL depth tendency
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   phbl           ! BL depth
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   pwb_fk_b       ! MLE buoyancy flux averaged over OSBL
       !!
       INTEGER  ::   jj, ji
@@ -1858,13 +1858,13 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER,                                   INTENT(in   ) ::   Kmm            ! Ocean time-level index
-      INTEGER,  DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   kp_ext         ! External-level offsets
-      REAL(wp), DIMENSION(A2D((nn_hls-1)),jpk),  INTENT(  out) ::   pdbdz          ! Gradients in the pycnocline
-      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(  out) ::   palpha
-      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pdh            ! Pycnocline thickness
-      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   phbl           ! BL depth
-      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients
-      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   phml           ! ML depth
-      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pdhdt          ! Rates of change of hbl
+      INTEGER,                                 INTENT(in   ) ::   Kmm            ! Ocean time-level index
+      INTEGER,  DIMENSION(A2D(nn_hls-1)),      INTENT(in   ) ::   kp_ext         ! External-level offsets
+      REAL(wp), DIMENSION(A2D(nn_hls-1),jpk),  INTENT(  out) ::   pdbdz          ! Gradients in the pycnocline
+      REAL(wp), DIMENSION(A2D(nn_hls-1)),      INTENT(  out) ::   palpha
+      REAL(wp), DIMENSION(A2D(nn_hls-1)),      INTENT(in   ) ::   pdh            ! Pycnocline thickness
+      REAL(wp), DIMENSION(A2D(nn_hls-1)),      INTENT(in   ) ::   phbl           ! BL depth
+      REAL(wp), DIMENSION(A2D(nn_hls-1)),      INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients
+      REAL(wp), DIMENSION(A2D(nn_hls-1)),      INTENT(in   ) ::   phml           ! ML depth
+      REAL(wp), DIMENSION(A2D(nn_hls-1)),      INTENT(in   ) ::   pdhdt          ! Rates of change of hbl
       !!
       INTEGER  ::   jk, jj, ji
@@ -1969,22 +1969,22 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER,                                   INTENT(in   ) ::   Kbb, Kmm       ! Ocean time-level indices
-      REAL(wp), DIMENSION(A2D((nn_hls-1)),jpk),  INTENT(inout) ::   pdiffut        ! t-diffusivity
-      REAL(wp), DIMENSION(A2D((nn_hls-1)),jpk),  INTENT(inout) ::   pviscos        ! Viscosity
-      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   phbl           ! BL depth
-      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   phml           ! ML depth
-      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pdh            ! Pycnocline depth
-      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pdhdt          ! BL depth tendency
-      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pshear         ! Shear production
-      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux
-      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pwb_min
+      INTEGER,                                 INTENT(in   ) ::   Kbb, Kmm       ! Ocean time-level indices
+      REAL(wp), DIMENSION(A2D(nn_hls-1),jpk),  INTENT(inout) ::   pdiffut        ! t-diffusivity
+      REAL(wp), DIMENSION(A2D(nn_hls-1),jpk),  INTENT(inout) ::   pviscos        ! Viscosity
+      REAL(wp), DIMENSION(A2D(nn_hls-1)),      INTENT(in   ) ::   phbl           ! BL depth
+      REAL(wp), DIMENSION(A2D(nn_hls-1)),      INTENT(in   ) ::   phml           ! ML depth
+      REAL(wp), DIMENSION(A2D(nn_hls-1)),      INTENT(in   ) ::   pdh            ! Pycnocline depth
+      REAL(wp), DIMENSION(A2D(nn_hls-1)),      INTENT(in   ) ::   pdhdt          ! BL depth tendency
+      REAL(wp), DIMENSION(A2D(nn_hls-1)),      INTENT(in   ) ::   pshear         ! Shear production
+      REAL(wp), DIMENSION(A2D(nn_hls-1)),      INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux
+      REAL(wp), DIMENSION(A2D(nn_hls-1)),      INTENT(in   ) ::   pwb_min
       !!
       INTEGER ::   ji, jj, jk   ! Loop indices
       !! Scales used to calculate eddy diffusivity and viscosity profiles
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zdifml_sc,    zvisml_sc
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zdifpyc_n_sc, zdifpyc_s_sc
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zvispyc_n_sc, zvispyc_s_sc
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zbeta_d_sc,   zbeta_v_sc
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zb_coup,      zc_coup_vis,  zc_coup_dif
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zdifml_sc,    zvisml_sc
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zdifpyc_n_sc, zdifpyc_s_sc
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zvispyc_n_sc, zvispyc_s_sc
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zbeta_d_sc,   zbeta_v_sc
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zb_coup,      zc_coup_vis,  zc_coup_dif
       !!
       REAL(wp) ::   zvel_sc_pyc, zvel_sc_ml, zstab_fac, zz_b
@@ -2164,48 +2164,48 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER,                                   INTENT(in   ) ::   Kmm            ! Time-level index
-      INTEGER,  DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   kp_ext         ! Offset for external level
-      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   phbl           ! BL depth
-      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   phml           ! ML depth
-      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pdh            ! Pycnocline depth
-      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pdhdt          ! BL depth tendency
-      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pshear         ! Shear production
-      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pdtdz_bl_ext   ! External temperature gradients
-      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pdsdz_bl_ext   ! External salinity gradients
-      REAL(wp), DIMENSION(A2D((nn_hls-1))),      INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients
-      REAL(wp), DIMENSION(A2D((nn_hls-1)),jpk),  INTENT(in   ) ::   pdiffut        ! t-diffusivity
-      REAL(wp), DIMENSION(A2D((nn_hls-1)),jpk),  INTENT(in   ) ::   pviscos        ! Viscosity
-      !!
-      REAL(wp), DIMENSION(A2D((nn_hls-1)))     ::   zalpha_pyc   !
-      REAL(wp), DIMENSION(A2D((nn_hls-1)),jpk) ::   zdbdz_pyc    ! Parametrised gradient of buoyancy in the pycnocline
+      INTEGER,                                 INTENT(in   ) ::   Kmm            ! Time-level index
+      INTEGER,  DIMENSION(A2D(nn_hls-1)),      INTENT(in   ) ::   kp_ext         ! Offset for external level
+      REAL(wp), DIMENSION(A2D(nn_hls-1)),      INTENT(in   ) ::   phbl           ! BL depth
+      REAL(wp), DIMENSION(A2D(nn_hls-1)),      INTENT(in   ) ::   phml           ! ML depth
+      REAL(wp), DIMENSION(A2D(nn_hls-1)),      INTENT(in   ) ::   pdh            ! Pycnocline depth
+      REAL(wp), DIMENSION(A2D(nn_hls-1)),      INTENT(in   ) ::   pdhdt          ! BL depth tendency
+      REAL(wp), DIMENSION(A2D(nn_hls-1)),      INTENT(in   ) ::   pshear         ! Shear production
+      REAL(wp), DIMENSION(A2D(nn_hls-1)),      INTENT(in   ) ::   pdtdz_bl_ext   ! External temperature gradients
+      REAL(wp), DIMENSION(A2D(nn_hls-1)),      INTENT(in   ) ::   pdsdz_bl_ext   ! External salinity gradients
+      REAL(wp), DIMENSION(A2D(nn_hls-1)),      INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients
+      REAL(wp), DIMENSION(A2D(nn_hls-1),jpk),  INTENT(in   ) ::   pdiffut        ! t-diffusivity
+      REAL(wp), DIMENSION(A2D(nn_hls-1),jpk),  INTENT(in   ) ::   pviscos        ! Viscosity
+      !!
+      REAL(wp), DIMENSION(A2D(nn_hls-1))     ::   zalpha_pyc   !
+      REAL(wp), DIMENSION(A2D(nn_hls-1),jpk) ::   zdbdz_pyc    ! Parametrised gradient of buoyancy in the pycnocline
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   z3ddz_pyc_1, z3ddz_pyc_2   ! Pycnocline gradient/shear profiles
       !!
-      INTEGER                              ::   ji, jj, jk, jkm_bld, jkf_mld, jkm_mld   ! Loop indices
-      INTEGER                              ::   istat                                   ! Memory allocation status
-      REAL(wp)                             ::   zznd_d, zznd_ml, zznd_pyc, znd          ! Temporary non-dimensional depths
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zsc_wth_1,zsc_ws_1                      ! Temporary scales
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zsc_uw_1, zsc_uw_2                      ! Temporary scales
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zsc_vw_1, zsc_vw_2                      ! Temporary scales
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   ztau_sc_u                               ! Dissipation timescale at base of WML
-      REAL(wp)                             ::   zbuoy_pyc_sc, zdelta_pyc                !
-      REAL(wp)                             ::   zl_c,zl_l,zl_eps                        ! Used to calculate turbulence length scale
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   za_cubic, zb_cubic                      ! Coefficients in cubic polynomial specifying
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zc_cubic, zd_cubic                      !    diffusivity in pycnocline
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zwt_pyc_sc_1, zws_pyc_sc_1              !
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zzeta_pyc                               !
-      REAL(wp)                             ::   zomega, zvw_max                         !
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zuw_bse,zvw_bse                         ! Momentum, heat, and salinity fluxes
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zwth_ent,zws_ent                        !    at the top of the pycnocline
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   zsc_wth_pyc, zsc_ws_pyc                 ! Scales for pycnocline transport term
-      REAL(wp)                             ::   ztmp                                    !
-      REAL(wp)                             ::   ztgrad, zsgrad, zbgrad                  ! Variables used to calculate pycnocline
-      !!                                                                                !    gradients
-      REAL(wp)                             ::   zugrad, zvgrad                          ! Variables for calculating pycnocline shear
-      REAL(wp)                             ::   zdtdz_pyc                               ! Parametrized gradient of temperature in
-      !!                                                                                !    pycnocline
-      REAL(wp)                             ::   zdsdz_pyc                               ! Parametrised gradient of salinity in
-      !!                                                                                !    pycnocline
-      REAL(wp)                             ::   zdudz_pyc                               ! u-shear across the pycnocline
-      REAL(wp)                             ::   zdvdz_pyc                               ! v-shear across the pycnocline
+      INTEGER                            ::   ji, jj, jk, jkm_bld, jkf_mld, jkm_mld   ! Loop indices
+      INTEGER                            ::   istat                                   ! Memory allocation status
+      REAL(wp)                           ::   zznd_d, zznd_ml, zznd_pyc, znd          ! Temporary non-dimensional depths
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zsc_wth_1,zsc_ws_1                      ! Temporary scales
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zsc_uw_1, zsc_uw_2                      ! Temporary scales
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zsc_vw_1, zsc_vw_2                      ! Temporary scales
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   ztau_sc_u                               ! Dissipation timescale at base of WML
+      REAL(wp)                           ::   zbuoy_pyc_sc, zdelta_pyc                !
+      REAL(wp)                           ::   zl_c,zl_l,zl_eps                        ! Used to calculate turbulence length scale
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   za_cubic, zb_cubic                      ! Coefficients in cubic polynomial specifying
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zc_cubic, zd_cubic                      !    diffusivity in pycnocline
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zwt_pyc_sc_1, zws_pyc_sc_1              !
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zzeta_pyc                               !
+      REAL(wp)                           ::   zomega, zvw_max                         !
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zuw_bse,zvw_bse                         ! Momentum, heat, and salinity fluxes
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zwth_ent,zws_ent                        !    at the top of the pycnocline
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   zsc_wth_pyc, zsc_ws_pyc                 ! Scales for pycnocline transport term
+      REAL(wp)                           ::   ztmp                                    !
+      REAL(wp)                           ::   ztgrad, zsgrad, zbgrad                  ! Variables used to calculate pycnocline
+      !!                                                                              !    gradients
+      REAL(wp)                           ::   zugrad, zvgrad                          ! Variables for calculating pycnocline shear
+      REAL(wp)                           ::   zdtdz_pyc                               ! Parametrized gradient of temperature in
+      !!                                                                              !    pycnocline
+      REAL(wp)                           ::   zdsdz_pyc                               ! Parametrised gradient of salinity in
+      !!                                                                              !    pycnocline
+      REAL(wp)                           ::   zdudz_pyc                               ! u-shear across the pycnocline
+      REAL(wp)                           ::   zdvdz_pyc                               ! v-shear across the pycnocline
       !!----------------------------------------------------------------------
       !
@@ -2229,12 +2229,12 @@
       ! Stokes term in scalar flux, flux-gradient relationship
       ! ------------------------------------------------------
-      WHERE ( l_conv(A2D((nn_hls-1))) )
-         zsc_wth_1(:,:) = swstrl(A2D((nn_hls-1)))**3 * swth0(A2D((nn_hls-1))) /   &
-            &             ( svstr(A2D((nn_hls-1)))**3 + 0.5_wp * swstrc(A2D((nn_hls-1)))**3 + epsln )
-         zsc_ws_1(:,:)  = swstrl(A2D((nn_hls-1)))**3 * sws0(A2D((nn_hls-1)))  /   &
-            &             ( svstr(A2D((nn_hls-1)))**3 + 0.5_wp * swstrc(A2D((nn_hls-1)))**3 + epsln )
+      WHERE ( l_conv(A2D(nn_hls-1)) )
+         zsc_wth_1(:,:) = swstrl(A2D(nn_hls-1))**3 * swth0(A2D(nn_hls-1)) /   &
+            &             ( svstr(A2D(nn_hls-1))**3 + 0.5_wp * swstrc(A2D(nn_hls-1))**3 + epsln )
+         zsc_ws_1(:,:)  = swstrl(A2D(nn_hls-1))**3 * sws0(A2D(nn_hls-1))  /   &
+            &             ( svstr(A2D(nn_hls-1))**3 + 0.5_wp * swstrc(A2D(nn_hls-1))**3 + epsln )
       ELSEWHERE
-         zsc_wth_1(:,:) = 2.0_wp * swthav(A2D((nn_hls-1)))
-         zsc_ws_1(:,:)  = 2.0_wp * swsav(A2D((nn_hls-1)))
+         zsc_wth_1(:,:) = 2.0_wp * swthav(A2D(nn_hls-1))
+         zsc_ws_1(:,:)  = 2.0_wp * swsav(A2D(nn_hls-1))
       ENDWHERE
       DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) )
@@ -2266,18 +2266,18 @@
       ! svstr since term needs to go to zero as swstrl goes to zero)
       ! ---------------------------------------------------------------------
-      WHERE ( l_conv(A2D((nn_hls-1))) )
-         zsc_uw_1(:,:) = ( swstrl(A2D((nn_hls-1)))**3 +                                                &
-            &              0.5_wp * swstrc(A2D((nn_hls-1)))**3 )**pthird * sustke(A2D((nn_hls-1))) /   &
-            &              MAX( ( 1.0_wp - 1.0_wp * 6.5_wp * sla(A2D((nn_hls-1)))**( 8.0_wp / 3.0_wp ) ), 0.2_wp )
-         zsc_uw_2(:,:) = ( swstrl(A2D((nn_hls-1)))**3 +                                                &
-            &              0.5_wp * swstrc(A2D((nn_hls-1)))**3 )**pthird * sustke(A2D((nn_hls-1))) /   &
-            &              MIN( sla(A2D((nn_hls-1)))**( 8.0_wp / 3.0_wp ) + epsln, 0.12_wp )
-         zsc_vw_1(:,:) = ff_t(A2D((nn_hls-1))) * phml(A2D((nn_hls-1))) * sustke(A2D((nn_hls-1)))**3 *   &
-            &            MIN( sla(A2D((nn_hls-1)))**( 8.0_wp / 3.0_wp ), 0.12_wp ) /                    &
-            &            ( ( svstr(A2D((nn_hls-1)))**3 + 0.5_wp * swstrc(A2D((nn_hls-1)))**3 )**( 2.0_wp / 3.0_wp ) + epsln )
+      WHERE ( l_conv(A2D(nn_hls-1)) )
+         zsc_uw_1(:,:) = ( swstrl(A2D(nn_hls-1))**3 +                                                &
+            &              0.5_wp * swstrc(A2D(nn_hls-1))**3 )**pthird * sustke(A2D(nn_hls-1)) /   &
+            &              MAX( ( 1.0_wp - 1.0_wp * 6.5_wp * sla(A2D(nn_hls-1))**( 8.0_wp / 3.0_wp ) ), 0.2_wp )
+         zsc_uw_2(:,:) = ( swstrl(A2D(nn_hls-1))**3 +                                                &
+            &              0.5_wp * swstrc(A2D(nn_hls-1))**3 )**pthird * sustke(A2D(nn_hls-1)) /   &
+            &              MIN( sla(A2D(nn_hls-1))**( 8.0_wp / 3.0_wp ) + epsln, 0.12_wp )
+         zsc_vw_1(:,:) = ff_t(A2D(nn_hls-1)) * phml(A2D(nn_hls-1)) * sustke(A2D(nn_hls-1))**3 *   &
+            &            MIN( sla(A2D(nn_hls-1))**( 8.0_wp / 3.0_wp ), 0.12_wp ) /                    &
+            &            ( ( svstr(A2D(nn_hls-1))**3 + 0.5_wp * swstrc(A2D(nn_hls-1))**3 )**( 2.0_wp / 3.0_wp ) + epsln )
       ELSEWHERE
-         zsc_uw_1(:,:) = sustar(A2D((nn_hls-1)))**2
-         zsc_vw_1(:,:) = ff_t(A2D((nn_hls-1))) * phbl(A2D((nn_hls-1))) * sustke(A2D((nn_hls-1)))**3 *   &
-            &            MIN( sla(A2D((nn_hls-1)))**( 8.0_wp / 3.0_wp ), 0.12_wp ) / ( svstr(A2D((nn_hls-1)))**2 + epsln )
+         zsc_uw_1(:,:) = sustar(A2D(nn_hls-1))**2
+         zsc_vw_1(:,:) = ff_t(A2D(nn_hls-1)) * phbl(A2D(nn_hls-1)) * sustke(A2D(nn_hls-1))**3 *   &
+            &            MIN( sla(A2D(nn_hls-1))**( 8.0_wp / 3.0_wp ), 0.12_wp ) / ( svstr(A2D(nn_hls-1))**2 + epsln )
       ENDWHERE
       DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) )
@@ -2303,9 +2303,9 @@
       ! (X0.3) and pressure (X0.5)]
       ! ----------------------------------------------------------------------
-      WHERE ( l_conv(A2D((nn_hls-1))) )
-         zsc_wth_1(:,:) = swbav(A2D((nn_hls-1))) * swth0(A2D((nn_hls-1))) * ( 1.0_wp + EXP( 0.2_wp * shol(A2D((nn_hls-1))) ) ) *   &
-            &             phml(A2D((nn_hls-1))) / ( svstr(A2D((nn_hls-1)))**3 + 0.5_wp * swstrc(A2D((nn_hls-1)))**3 + epsln )
-         zsc_ws_1(:,:)  = swbav(A2D((nn_hls-1))) * sws0(A2D((nn_hls-1)))  * ( 1.0_wp + EXP( 0.2_wp * shol(A2D((nn_hls-1))) ) ) *   &
-            &             phml(A2D((nn_hls-1))) / ( svstr(A2D((nn_hls-1)))**3 + 0.5_wp * swstrc(A2D((nn_hls-1)))**3 + epsln )
+      WHERE ( l_conv(A2D(nn_hls-1)) )
+         zsc_wth_1(:,:) = swbav(A2D(nn_hls-1)) * swth0(A2D(nn_hls-1)) * ( 1.0_wp + EXP( 0.2_wp * shol(A2D(nn_hls-1)) ) ) *   &
+            &             phml(A2D(nn_hls-1)) / ( svstr(A2D(nn_hls-1))**3 + 0.5_wp * swstrc(A2D(nn_hls-1))**3 + epsln )
+         zsc_ws_1(:,:)  = swbav(A2D(nn_hls-1)) * sws0(A2D(nn_hls-1))  * ( 1.0_wp + EXP( 0.2_wp * shol(A2D(nn_hls-1)) ) ) *   &
+            &             phml(A2D(nn_hls-1)) / ( svstr(A2D(nn_hls-1))**3 + 0.5_wp * swstrc(A2D(nn_hls-1))**3 + epsln )
       ELSEWHERE
          zsc_wth_1(:,:) = 0.0_wp
@@ -2379,9 +2379,9 @@
       !
       zsc_vw_1(:,:) = 0.0_wp
-      WHERE ( l_conv(A2D((nn_hls-1))) )
-         zsc_uw_1(:,:) = -1.0_wp * swb0(A2D((nn_hls-1))) * sustar(A2D((nn_hls-1)))**2 * phml(A2D((nn_hls-1))) /   &
-            &            ( svstr(A2D((nn_hls-1)))**3 + 0.5_wp * swstrc(A2D((nn_hls-1)))**3 + epsln )
-         zsc_uw_2(:,:) =           swb0(A2D((nn_hls-1))) * sustke(A2D((nn_hls-1)))    * phml(A2D((nn_hls-1))) /   &
-            &            ( svstr(A2D((nn_hls-1)))**3 + 0.5_wp * swstrc(A2D((nn_hls-1)))**3 + epsln )**( 2.0_wp / 3.0_wp )
+      WHERE ( l_conv(A2D(nn_hls-1)) )
+         zsc_uw_1(:,:) = -1.0_wp * swb0(A2D(nn_hls-1)) * sustar(A2D(nn_hls-1))**2 * phml(A2D(nn_hls-1)) /   &
+            &            ( svstr(A2D(nn_hls-1))**3 + 0.5_wp * swstrc(A2D(nn_hls-1))**3 + epsln )
+         zsc_uw_2(:,:) =           swb0(A2D(nn_hls-1)) * sustke(A2D(nn_hls-1))    * phml(A2D(nn_hls-1)) /   &
+            &            ( svstr(A2D(nn_hls-1))**3 + 0.5_wp * swstrc(A2D(nn_hls-1))**3 + epsln )**( 2.0_wp / 3.0_wp )
       ELSEWHERE
          zsc_uw_1(:,:) = 0.0_wp
@@ -2443,16 +2443,16 @@
       ! (X0.3) ]
       ! -----------------------------------------------------------------------
-      WHERE ( l_conv(A2D((nn_hls-1))) )
-         zsc_wth_1(:,:) = swth0(A2D((nn_hls-1))) / ( 1.0_wp - 0.56_wp * EXP( shol(A2D((nn_hls-1))) ) )
-         zsc_ws_1(:,:)  = sws0(A2D((nn_hls-1)))  / ( 1.0_wp - 0.56_wp * EXP( shol(A2D((nn_hls-1))) ) )
-         WHERE ( l_pyc(A2D((nn_hls-1))) )   ! Pycnocline scales
-            zsc_wth_pyc(:,:) = -0.003_wp * swstrc(A2D((nn_hls-1))) * ( 1.0_wp - pdh(A2D((nn_hls-1))) / phbl(A2D((nn_hls-1))) ) *   &
-               &               av_dt_ml(A2D((nn_hls-1)))
-            zsc_ws_pyc(:,:)  = -0.003_wp * swstrc(A2D((nn_hls-1))) * ( 1.0_wp - pdh(A2D((nn_hls-1))) / phbl(A2D((nn_hls-1))) ) *   &
-               &               av_ds_ml(A2D((nn_hls-1)))
+      WHERE ( l_conv(A2D(nn_hls-1)) )
+         zsc_wth_1(:,:) = swth0(A2D(nn_hls-1)) / ( 1.0_wp - 0.56_wp * EXP( shol(A2D(nn_hls-1)) ) )
+         zsc_ws_1(:,:)  = sws0(A2D(nn_hls-1))  / ( 1.0_wp - 0.56_wp * EXP( shol(A2D(nn_hls-1)) ) )
+         WHERE ( l_pyc(A2D(nn_hls-1)) )   ! Pycnocline scales
+            zsc_wth_pyc(:,:) = -0.003_wp * swstrc(A2D(nn_hls-1)) * ( 1.0_wp - pdh(A2D(nn_hls-1)) / phbl(A2D(nn_hls-1)) ) *   &
+               &               av_dt_ml(A2D(nn_hls-1))
+            zsc_ws_pyc(:,:)  = -0.003_wp * swstrc(A2D(nn_hls-1)) * ( 1.0_wp - pdh(A2D(nn_hls-1)) / phbl(A2D(nn_hls-1)) ) *   &
+               &               av_ds_ml(A2D(nn_hls-1))
          END WHERE
       ELSEWHERE
-         zsc_wth_1(:,:) = 2.0_wp * swthav(A2D((nn_hls-1)))
-         zsc_ws_1(:,:)  =          sws0(A2D((nn_hls-1)))
+         zsc_wth_1(:,:) = 2.0_wp * swthav(A2D(nn_hls-1))
+         zsc_ws_1(:,:)  =          sws0(A2D(nn_hls-1))
       END WHERE
       DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, MAX( jkm_mld, jkm_bld ) )
@@ -2491,12 +2491,12 @@
       END_3D
       !
-      WHERE ( l_conv(A2D((nn_hls-1))) )
-         zsc_uw_1(:,:) = sustar(A2D((nn_hls-1)))**2
-         zsc_vw_1(:,:) = ff_t(A2D((nn_hls-1))) * sustke(A2D((nn_hls-1))) * phml(A2D((nn_hls-1)))
+      WHERE ( l_conv(A2D(nn_hls-1)) )
+         zsc_uw_1(:,:) = sustar(A2D(nn_hls-1))**2
+         zsc_vw_1(:,:) = ff_t(A2D(nn_hls-1)) * sustke(A2D(nn_hls-1)) * phml(A2D(nn_hls-1))
       ELSEWHERE
-         zsc_uw_1(:,:) = sustar(A2D((nn_hls-1)))**2
+         zsc_uw_1(:,:) = sustar(A2D(nn_hls-1))**2
          zsc_uw_2(:,:) = ( 2.25_wp - 3.0_wp * ( 1.0_wp - EXP( -1.25_wp * 2.0_wp ) ) ) * ( 1.0_wp - EXP( -4.0_wp * 2.0_wp ) ) *   &
             &            zsc_uw_1(:,:)
-         zsc_vw_1(:,:) = ff_t(A2D((nn_hls-1))) * sustke(A2D((nn_hls-1))) * phbl(A2D((nn_hls-1)))
+         zsc_vw_1(:,:) = ff_t(A2D(nn_hls-1)) * sustke(A2D(nn_hls-1)) * phbl(A2D(nn_hls-1))
          zsc_vw_2(:,:) = -0.11_wp * SIN( 3.14159_wp * ( 2.0_wp + 0.4_wp ) ) * EXP( -1.0_wp * ( 1.5_wp + 2.0_wp )**2 ) *   &
             &            zsc_vw_1(:,:)
@@ -2702,11 +2702,11 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER,                              INTENT(in   ) ::   Kmm          ! Time-level index
-      REAL(wp), DIMENSION(A2D(nn_hls)),     INTENT(  out) ::   pmld         ! == Estimated FK BLD used for MLE horizontal gradients == !
-      REAL(wp), DIMENSION(A2D(nn_hls)),     INTENT(inout) ::   pdtdx        ! Horizontal gradient for Fox-Kemper parametrization
-      REAL(wp), DIMENSION(A2D(nn_hls)),     INTENT(inout) ::   pdtdy        ! Horizontal gradient for Fox-Kemper parametrization
-      REAL(wp), DIMENSION(A2D(nn_hls)),     INTENT(inout) ::   pdsdx        ! Horizontal gradient for Fox-Kemper parametrization
-      REAL(wp), DIMENSION(A2D(nn_hls)),     INTENT(inout) ::   pdsdy        ! Horizontal gradient for Fox-Kemper parametrization
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) ::   pdbds_mle    ! Magnitude of horizontal buoyancy gradient
+      INTEGER,                            INTENT(in   ) ::   Kmm          ! Time-level index
+      REAL(wp), DIMENSION(A2D(nn_hls)),   INTENT(  out) ::   pmld         ! == Estimated FK BLD used for MLE horizontal gradients == !
+      REAL(wp), DIMENSION(A2D(nn_hls)),   INTENT(inout) ::   pdtdx        ! Horizontal gradient for Fox-Kemper parametrization
+      REAL(wp), DIMENSION(A2D(nn_hls)),   INTENT(inout) ::   pdtdy        ! Horizontal gradient for Fox-Kemper parametrization
+      REAL(wp), DIMENSION(A2D(nn_hls)),   INTENT(inout) ::   pdsdx        ! Horizontal gradient for Fox-Kemper parametrization
+      REAL(wp), DIMENSION(A2D(nn_hls)),   INTENT(inout) ::   pdsdy        ! Horizontal gradient for Fox-Kemper parametrization
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(inout) ::   pdbds_mle    ! Magnitude of horizontal buoyancy gradient
       !!
       INTEGER                               ::   ji, jj, jk   ! Dummy loop indices
@@ -2809,19 +2809,19 @@
       !!
       !!----------------------------------------------------------------------      
-      INTEGER,                              INTENT(in   ) ::   Kmm         ! Time-level index
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) ::   pwb_fk
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   phbl        ! BL depth
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   phmle       ! MLE depth
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pwb_ent     ! Buoyancy entrainment flux
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pdbds_mle   ! Magnitude of horizontal buoyancy gradient
-      !!
-      INTEGER                              ::   ji, jj, jk        ! Dummy loop indices
-      REAL(wp), DIMENSION(A2D((nn_hls-1))) ::   znd_param
-      REAL(wp)                             ::   zthermal, zbeta
-      REAL(wp)                             ::   zbuoy
-      REAL(wp)                             ::   ztmp
-      REAL(wp)                             ::   zpe_mle_layer
-      REAL(wp)                             ::   zpe_mle_ref
-      REAL(wp)                             ::   zdbdz_mle_int
+      INTEGER,                            INTENT(in   ) ::   Kmm         ! Time-level index
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(inout) ::   pwb_fk
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   phbl        ! BL depth
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   phmle       ! MLE depth
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   pwb_ent     ! Buoyancy entrainment flux
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   pdbds_mle   ! Magnitude of horizontal buoyancy gradient
+      !!
+      INTEGER                            ::   ji, jj, jk        ! Dummy loop indices
+      REAL(wp), DIMENSION(A2D(nn_hls-1)) ::   znd_param
+      REAL(wp)                           ::   zthermal, zbeta
+      REAL(wp)                           ::   zbuoy
+      REAL(wp)                           ::   ztmp
+      REAL(wp)                           ::   zpe_mle_layer
+      REAL(wp)                           ::   zpe_mle_ref
+      REAL(wp)                           ::   zdbdz_mle_int
       !!----------------------------------------------------------------------      
       !
@@ -2921,12 +2921,12 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER,                              INTENT(in   ) ::   Kmm         ! Time-level index
-      REAL(wp), DIMENSION(A2D((nn_hls))),   INTENT(in   ) ::   pmld        ! == Estimated FK BLD used for MLE horiz gradients == !
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) ::   phmle       ! MLE depth
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) ::   pvel_mle    ! Velocity scale for dhdt with stable ML and FK
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(inout) ::   pdiff_mle   ! Extra MLE vertical diff
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pdbds_mle   ! Magnitude of horizontal buoyancy gradient
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   phbl        ! BL depth
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(in   ) ::   pwb0tot     ! Total surface buoyancy flux including insolation
+      INTEGER,                            INTENT(in   ) ::   Kmm         ! Time-level index
+      REAL(wp), DIMENSION(A2D(nn_hls)),   INTENT(in   ) ::   pmld        ! == Estimated FK BLD used for MLE horiz gradients == !
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(inout) ::   phmle       ! MLE depth
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(inout) ::   pvel_mle    ! Velocity scale for dhdt with stable ML and FK
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(inout) ::   pdiff_mle   ! Extra MLE vertical diff
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   pdbds_mle   ! Magnitude of horizontal buoyancy gradient
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   phbl        ! BL depth
+      REAL(wp), DIMENSION(A2D(nn_hls-1)), INTENT(in   ) ::   pwb0tot     ! Total surface buoyancy flux including insolation
       !!
       INTEGER  ::   ji, jj, jk   ! Dummy loop indices
Index: /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/do_loop_substitute.h90
===================================================================
--- /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/do_loop_substitute.h90	(revision 14912)
+++ /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/do_loop_substitute.h90	(revision 14913)
@@ -61,6 +61,6 @@
 #define DO_2D(L, R, B, T) DO jj = ntsj-(B), ntej+(T) ; DO ji = ntsi-(L), ntei+(R)
 #define DO_2D_OVR(L, R, B, T) DO_2D(L-(L+R)*nthl, R-(R+L)*nthr, B-(B+T)*nthb, T-(T+B)*ntht)
-#define A1Di(H) ntsi-H:ntei+H
-#define A1Dj(H) ntsj-H:ntej+H
+#define A1Di(H) ntsi-(H):ntei+(H)
+#define A1Dj(H) ntsj-(H):ntej+(H)
 #define A2D(H) A1Di(H),A1Dj(H)
 #define A1Di_T(T) (ntsi-nn_hls-1)*T+1:
Index: /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/TOP/C14/trcnam_c14.F90
===================================================================
--- /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/TOP/C14/trcnam_c14.F90	(revision 14912)
+++ /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/TOP/C14/trcnam_c14.F90	(revision 14913)
@@ -60,4 +60,5 @@
       ln_trc_cbc(jp_c14) = .false.
       ln_trc_obc(jp_c14) = .false.
+      ln_trc_ais(jp_c14) = .false.
       !
       READ  ( numtrc_ref, namc14_typ, IOSTAT = ios, ERR = 901)
Index: /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/TOP/CFC/trcnam_cfc.F90
===================================================================
--- /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/TOP/CFC/trcnam_cfc.F90	(revision 14912)
+++ /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/TOP/CFC/trcnam_cfc.F90	(revision 14913)
@@ -77,4 +77,5 @@
          ln_trc_cbc(jn) = .false.
          ln_trc_obc(jn) = .false.
+         ln_trc_ais(jn) = .false.
       ENDIF
       !
