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 14899)
+++ NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90	(revision 14900)
@@ -61,4 +61,7 @@
    !!   trc_osm        : compute and add to the passive tracer trend the non-local flux (TBD)
    !!   dyn_osm        : compute and add to u & v trensd the non-local flux
+   !!   zdf_osm_iomput : iom_put wrapper that accepts arrays without halo
+   !!      zdf_osm_iomput_2d                    : iom_put wrapper for 2D fields
+   !!      zdf_osm_iomput_3d                    : iom_put wrapper for 3D fields
    !!----------------------------------------------------------------------
    USE oce                                        ! Ocean dynamics and active tracers
@@ -114,4 +117,12 @@
       MODULE PROCEDURE zdf_osm_velocity_rotation_2d
       MODULE PROCEDURE zdf_osm_velocity_rotation_3d
+   END INTERFACE
+   !
+   INTERFACE zdf_osm_iomput
+      !!---------------------------------------------------------------------
+      !!                 ***  INTERFACE zdf_osm_iomput  ***
+      !!---------------------------------------------------------------------
+      MODULE PROCEDURE zdf_osm_iomput_2d
+      MODULE PROCEDURE zdf_osm_iomput_3d
    END INTERFACE
    !
@@ -370,10 +381,10 @@
       !
       REAL(wp), DIMENSION(A2D((nn_hls-1)))  ::   zhmle   ! MLE depth - grid
-      REAL(wp), DIMENSION(jpi,jpj) ::   zmld   ! ML depth on 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, zdbdz_bl_ext   ! External temperature/salinity and buoyancy gradients
-      REAL(wp), DIMENSION(jpi,jpj) ::   zdtdx, zdtdy, zdsdx, zdsdy   ! Horizontal gradients for Fox-Kemper parametrization
+      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
@@ -384,7 +395,7 @@
       !
       ! For calculating Ri#-dependent mixing
-      REAL(wp), DIMENSION(jpi,jpj) ::   z2du     ! u-shear^2
-      REAL(wp), DIMENSION(jpi,jpj) ::   z2dv     ! v-shear^2
-      REAL(wp)                     ::   zrimix   ! Spatial form of ri#-induced diffusion
+      REAL(wp), DIMENSION(A2D(nn_hls)) ::   z2du     ! u-shear^2
+      REAL(wp), DIMENSION(A2D(nn_hls)) ::   z2dv     ! v-shear^2
+      REAL(wp)                         ::   zrimix   ! Spatial form of ri#-induced diffusion
       !
       ! Temporary variables
@@ -401,10 +412,11 @@
       REAL(wp), PARAMETER ::   pp_large = -1e10_wp
       !
-      nbld(:,:)   = 0
-      nmld(:,:)   = 0
-      sustke(:,:) = pp_large
-      l_pyc(:,:)  = .FALSE.
-      l_flux(:,:) = .FALSE.
-      l_mle(:,:)  = .FALSE.
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         nmld(ji,jj)   = 0
+         sustke(ji,jj) = pp_large
+         l_pyc(ji,jj)  = .FALSE.
+         l_flux(ji,jj) = .FALSE.
+         l_mle(ji,jj)  = .FALSE.
+      END_2D
       ! Mixed layer
       ! No initialization of zhbl or zhml (or zdh?)
@@ -413,7 +425,11 @@
       IF ( ln_osm_mle ) THEN   ! Only initialise arrays if needed
          zdtdx(:,:)  = pp_large ; zdtdy(:,:)    = pp_large ; zdsdx(:,:)     = pp_large
-         zdsdy(:,:)  = pp_large ; dbdx_mle(:,:) = pp_large ; dbdy_mle(:,:)  = pp_large
+         zdsdy(:,:)  = pp_large
          zwb_fk(:,:) = pp_large ; zvel_mle(:,:) = pp_large
          zhmle(:,:)  = pp_large ; zmld(:,:)     = pp_large
+         DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
+            dbdx_mle(ji,jj) = pp_large
+            dbdy_mle(ji,jj) = pp_large
+         END_2D
       ENDIF
       zhbl_t(:,:)   = pp_large
@@ -422,8 +438,16 @@
       zviscos(:,:,:) = 0.0_wp
       !
-      ghamt(:,:,:)      = pp_large ; ghams(:,:,:)      = pp_large
-      ghamt(A2D((nn_hls-1)),:)   = 0.0_wp   ; ghams(A2D((nn_hls-1)),:)   = 0.0_wp
-      ghamu(:,:,:)      = pp_large ; ghamv(:,:,:)      = pp_large
-      ghamu(A2D((nn_hls-1)),:)   = 0.0_wp   ; ghamv(A2D((nn_hls-1)),:)   = 0.0_wp
+      DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 1, jpk )
+         ghamt(ji,jj,jk) = pp_large
+         ghams(ji,jj,jk) = pp_large
+         ghamu(ji,jj,jk) = pp_large
+         ghamv(ji,jj,jk) = pp_large
+      END_3D
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, jpk )
+         ghamt(ji,jj,jk) = 0.0_wp
+         ghams(ji,jj,jk) = 0.0_wp
+         ghamu(ji,jj,jk) = 0.0_wp
+         ghamv(ji,jj,jk) = 0.0_wp
+      END_3D
       !
       zdiff_mle(:,:) = 0.0_wp
@@ -431,5 +455,5 @@
       ! Ensure only positive hbl values are accessed when using extended halo
       ! (nn_hls==2)
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          hbl(ji,jj) = MAX( hbl(ji,jj), epsln )
       END_2D
@@ -442,5 +466,5 @@
       zz0 =           rn_abs   ! Assume two-band radiation model for depth of OSBL - surface equi-partition in 2-bands
       zz1 =  1.0_wp - rn_abs
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          zrad0(ji,jj)  = qsr(ji,jj) * r1_rho0_rcp   ! Surface downward irradiance (so always +ve)
          zradh(ji,jj)  = zrad0(ji,jj) *                                &   ! Downwards irradiance at base of boundary layer
@@ -453,5 +477,5 @@
             &                                                 zradav )                              !    over depth of OSBL
       END_2D
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          sws0(ji,jj)    = -1.0_wp * ( ( emp(ji,jj) - rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm) +   &   ! Upwards surface salinity flux
             &                         sfx(ji,jj) ) * r1_rho0 * tmask(ji,jj,1)                      !    for non-local term
@@ -464,5 +488,5 @@
             &             grav  * zbeta * swsav(ji,jj)                      ! OBSBL
       END_2D
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          suw0(ji,jj)    = -0.5_wp * (utau(ji-1,jj) + utau(ji,jj)) * r1_rho0 * tmask(ji,jj,1)   ! Surface upward velocity fluxes
          zvw0           = -0.5_wp * (vtau(ji,jj-1) + vtau(ji,jj)) * r1_rho0 * tmask(ji,jj,1)
@@ -476,5 +500,5 @@
          ! Assume constant La#=0.3
       CASE(0)
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             zus_x = scos_wind(ji,jj) * sustar(ji,jj) / 0.3_wp**2
             zus_y = ssin_wind(ji,jj) * sustar(ji,jj) / 0.3_wp**2
@@ -485,5 +509,5 @@
          ! Assume Pierson-Moskovitz wind-wave spectrum
       CASE(1)
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             ! Use wind speed wndm included in sbc_oce module
             sustke(ji,jj)  = MAX ( 0.016_wp * wndm(ji,jj), 1e-8_wp )
@@ -494,5 +518,5 @@
          zfac =  2.0_wp * rpi / 16.0_wp
          !
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             IF ( hsw(ji,jj) > 1e-4_wp ) THEN
                ! Use  wave fields
@@ -511,5 +535,5 @@
       IF (ln_zdfosm_ice_shelter) THEN
          ! Reduce both Stokes drift and its depth scale by ocean fraction to represent sheltering by ice
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             sustke(ji,jj)  = sustke(ji,jj)  * ( 1.0_wp - fr_i(ji,jj) )
             dstokes(ji,jj) = dstokes(ji,jj) * ( 1.0_wp - fr_i(ji,jj) )
@@ -524,5 +548,7 @@
          ! It could represent the effects of the spread of wave directions around the mean wind. The effect of this adjustment needs to be tested.
          IF(nn_osm_wave > 0) THEN
-            sustke(A2D((nn_hls-1))) = rn_zdfosm_adjust_sd * sustke(A2D((nn_hls-1)))
+            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+               sustke(ji,jj) = rn_zdfosm_adjust_sd * sustke(ji,jj)
+            END_2D
          END IF
       CASE(1)
@@ -531,5 +557,5 @@
          zsqrtpi = SQRT(rpi)
          z_two_thirds = 2.0_wp / 3.0_wp
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             zthickness = rn_osm_hblfrac*hbl(ji,jj)
             z2k_times_thickness =  zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 1e-7_wp )
@@ -545,5 +571,5 @@
          ! Assumes approximate depth profile of SD from Breivik (2016)
          zsqrtpi = SQRT(rpi)
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             zthickness = rn_osm_hblfrac*hbl(ji,jj)
             z2k_times_thickness =  zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 1e-7_wp )
@@ -567,5 +593,5 @@
       ! Langmuir velocity scale (swstrl), La # (sla)
       ! Mixed scale (svstr), convective velocity scale (swstrc)
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          ! Langmuir velocity scale (swstrl), at T-point
          swstrl(ji,jj) = ( sustar(ji,jj) * sustar(ji,jj) * sustke(ji,jj) )**pthird
@@ -594,13 +620,16 @@
       ! BL must be always 4 levels deep.
       ! For calculation of lateral buoyancy gradients for FK in
-      ! zdf_osm_zmld_horizontal_gradients need halo values for nbld, so must
-      ! previously exist for hbl also.
+      ! zdf_osm_zmld_horizontal_gradients need halo values for nbld
       !
       ! agn 23/6/20: not clear all this is needed, as hbl checked after it is re-calculated anyway
       ! ##########################################################################
-      hbl(:,:) = MAX(hbl(:,:), gdepw(:,:,4,Kmm) )
-      nbld(:,:) = 4
-      DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 5, jpkm1 )
-         IF ( hbl(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         hbl(ji,jj) = MAX(hbl(ji,jj), gdepw(ji,jj,4,Kmm) )
+      END_2D
+      DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
+         nbld(ji,jj) = 4
+      END_2D
+      DO_3D_OVR( nn_hls, nn_hls, nn_hls, nn_hls, 5, jpkm1 )
+         IF ( MAX( hbl(ji,jj), gdepw(ji,jj,4,Kmm) ) >= gdepw(ji,jj,jk,Kmm) ) THEN
             nbld(ji,jj) = MIN(mbkt(ji,jj)-2, jk)
          ENDIF
@@ -608,5 +637,5 @@
       ! ##########################################################################
       !
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          zhbl(ji,jj) = gdepw(ji,jj,nbld(ji,jj),Kmm)
          nmld(ji,jj) = MAX( 3, nbld(ji,jj) - MAX( INT( dh(ji,jj) / e3t(ji,jj,nbld(ji,jj)-1,Kmm) ), 1 ) )
@@ -617,16 +646,24 @@
       ! Averages over well-mixed and boundary layer, note BL averages use jk_ext=2 everywhere
       jk_ext(:,:) = 1   ! ag 19/03
-      CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D((nn_hls-1))), 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_ext(:,:) = nbld(A2D((nn_hls-1))) - nmld(:,:) + jk_ext(:,:) + 1   ! ag 19/03
-      CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, 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 )
+      CALL zdf_osm_vertical_average( Kbb,                       Kmm,                         &
+         &                           nbld(A2D((nn_hls-1))),     av_t_bl(A2D((nn_hls-1))),    &
+         &                           av_s_bl(A2D((nn_hls-1))),  av_b_bl(A2D((nn_hls-1))),    &
+         &                           av_u_bl(A2D((nn_hls-1))),  av_v_bl(A2D((nn_hls-1))),    &
+         &                           jk_ext,                    av_dt_bl(A2D((nn_hls-1))),   &
+         &                           av_ds_bl(A2D((nn_hls-1))), av_db_bl(A2D((nn_hls-1))),   &
+         &                           av_du_bl(A2D((nn_hls-1))), av_dv_bl(A2D((nn_hls-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,                         &
+         &                           nmld - 1,                  av_t_ml(A2D((nn_hls-1))),    &
+         &                           av_s_ml(A2D((nn_hls-1))),  av_b_ml(A2D((nn_hls-1))),    &
+         &                           av_u_ml(A2D((nn_hls-1))),  av_v_ml(A2D((nn_hls-1))),    &
+         &                           jk_ext,                    av_dt_ml(A2D((nn_hls-1))),   &
+         &                           av_ds_ml(A2D((nn_hls-1))), av_db_ml(A2D((nn_hls-1))),   &
+         &                           av_du_ml(A2D((nn_hls-1))), av_dv_ml(A2D((nn_hls-1))) )
       ! Velocity components in frame aligned with surface stress
-      CALL zdf_osm_velocity_rotation( av_u_ml, av_v_ml  )
-      CALL zdf_osm_velocity_rotation( av_du_ml, av_dv_ml )
-      CALL zdf_osm_velocity_rotation( av_u_bl, av_v_bl  )
-      CALL zdf_osm_velocity_rotation( av_du_bl, av_dv_bl )
+      CALL zdf_osm_velocity_rotation( av_u_ml(A2D((nn_hls-1))),  av_v_ml(A2D((nn_hls-1)))  )
+      CALL zdf_osm_velocity_rotation( av_du_ml(A2D((nn_hls-1))), av_dv_ml(A2D((nn_hls-1))) )
+      CALL zdf_osm_velocity_rotation( av_u_bl(A2D((nn_hls-1))),  av_v_bl(A2D((nn_hls-1)))  )
+      CALL zdf_osm_velocity_rotation( av_du_bl(A2D((nn_hls-1))), av_dv_bl(A2D((nn_hls-1))) )
       !
       ! Determine the state of the OSBL, stable/unstable, shear/no shear
@@ -636,12 +673,16 @@
       IF ( ln_osm_mle ) THEN
          ! Fox-Kemper Scheme
-         mld_prof = 4
-         DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 5, jpkm1 )
+         DO_2D_OVR( nn_hls, nn_hls, nn_hls, nn_hls )
+            mld_prof(ji,jj) = 4
+         END_2D
+         DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 5, jpkm1 )
             IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk)
          END_3D
-         CALL zdf_osm_vertical_average( Kbb, Kmm, mld_prof(A2D((nn_hls-1))), av_t_mle, av_s_mle,   &
-            &                           av_b_mle, av_u_mle, av_v_mle )
+         CALL zdf_osm_vertical_average( Kbb,                       Kmm,                         &
+            &                           mld_prof(A2D((nn_hls-1))), av_t_mle(A2D((nn_hls-1))),   &
+            &                           av_s_mle(A2D((nn_hls-1))), av_b_mle(A2D((nn_hls-1))),   &
+            &                           av_u_mle(A2D((nn_hls-1))), av_v_mle(A2D((nn_hls-1))) )
          !
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             zhmle(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm)
          END_2D
@@ -649,17 +690,17 @@
          ! Calculate fairly-well-mixed depth zmld & its index mld_prof + lateral zmld-averaged gradients
          CALL zdf_osm_zmld_horizontal_gradients( Kmm, zmld, zdtdx, zdtdy, zdsdx,   &
-            &                                    zdsdy, dbdx_mle, dbdy_mle, zdbds_mle )
+            &                                    zdsdy, zdbds_mle )
          ! Calculate max vertical FK flux zwb_fk & set logical descriptors
          CALL zdf_osm_osbl_state_fk( Kmm, zwb_fk, zhbl, zhmle, zwb_ent,   &
             &                        zdbds_mle )
          ! Recalculate hmle, zmle, zvel_mle, zdiff_mle & redefine mld_proc to be index for new hmle
-         CALL zdf_osm_mle_parameters( Kmm, zmld, zhmle, zvel_mle, zdiff_mle,   &
+         CALL zdf_osm_mle_parameters( Kmm, zmld(A2D((nn_hls-1))), zhmle, zvel_mle, zdiff_mle,   &
             &                         zdbds_mle, zhbl, zwb0tot )
       ELSE    ! ln_osm_mle
          ! FK not selected, Boundary Layer only.
-         l_pyc(:,:)  = .TRUE.
-         l_flux(:,:) = .FALSE.
-         l_mle(:,:)  = .FALSE.
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+            l_pyc(ji,jj)  = .TRUE.
+            l_flux(ji,jj) = .FALSE.
+            l_mle(ji,jj)  = .FALSE.
             IF ( l_conv(ji,jj) .AND. av_db_bl(ji,jj) < rn_osm_bl_thresh ) l_pyc(ji,jj) = .FALSE.
          END_2D
@@ -688,11 +729,19 @@
       ! Recalculate bl averages using jk_ext & ml averages .... note no rotation of u & v here..
       jk_ext(:,:) = 1   ! ag 19/03
-      CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D((nn_hls-1))), 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_ext(:,:) = nbld(A2D((nn_hls-1))) - nmld(:,:) + jk_ext(:,:) + 1   ! ag 19/03
-      CALL zdf_osm_vertical_average( Kbb, Kmm, nmld - 1, 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 )   ! ag 19/03
+      CALL zdf_osm_vertical_average( Kbb,                       Kmm,                         &
+         &                           nbld(A2D((nn_hls-1))),     av_t_bl(A2D((nn_hls-1))),    &
+         &                           av_s_bl(A2D((nn_hls-1))),  av_b_bl(A2D((nn_hls-1))),    &
+         &                           av_u_bl(A2D((nn_hls-1))),  av_v_bl(A2D((nn_hls-1))),    &
+         &                           jk_ext,                    av_dt_bl(A2D((nn_hls-1))),   &
+         &                           av_ds_bl(A2D((nn_hls-1))), av_db_bl(A2D((nn_hls-1))),   &
+         &                           av_du_bl(A2D((nn_hls-1))), av_dv_bl(A2D((nn_hls-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,                         &
+         &                           nmld(A2D((nn_hls-1))) - 1, av_t_ml(A2D((nn_hls-1))),    &
+         &                           av_s_ml(A2D((nn_hls-1))),  av_b_ml(A2D((nn_hls-1))),    &
+         &                           av_u_ml(A2D((nn_hls-1))),  av_v_ml(A2D((nn_hls-1))),    &
+         &                           jk_ext,                    av_dt_ml(A2D((nn_hls-1))),   &
+         &                           av_ds_ml(A2D((nn_hls-1))), av_db_ml(A2D((nn_hls-1))),   &
+         &                           av_du_ml(A2D((nn_hls-1))), av_dv_ml(A2D((nn_hls-1))) )   ! ag 19/03
       !
       ! Rate of change of hbl
@@ -700,6 +749,6 @@
          &                         zdbdz_bl_ext, zwb_fk_b, zwb_fk, zvel_mle )
       ! Test if surface boundary layer coupled to bottom
-      l_coup(:,:) = .FALSE.   ! ag 19/03
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         l_coup(ji,jj) = .FALSE.   ! ag 19/03
          zhbl_t(ji,jj) = hbl(ji,jj) + ( zdhdt(ji,jj) - ww(ji,jj,nbld(ji,jj)) ) * rn_Dt   ! Certainly need ww here, so subtract it
          ! Adjustment to represent limiting by ocean bottom
@@ -713,8 +762,10 @@
       END_2D
       !
-      nmld(:,:) = nbld(A2D((nn_hls-1)))           ! use nmld to hold previous blayer index
-      nbld(:,:) = 4
-      !
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 4, jpkm1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         nmld(ji,jj) = nbld(ji,jj)           ! use nmld to hold previous blayer index
+         nbld(ji,jj) = 4
+      END_2D
+      !
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 4, jpkm1 )
          IF ( zhbl_t(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN
             nbld(ji,jj) = jk
@@ -731,7 +782,11 @@
       ! Recalculate BL averages and differences using new BL depth
       jk_ext(:,:) = 1   ! ag 19/03
-      CALL zdf_osm_vertical_average( Kbb, Kmm, nbld(A2D((nn_hls-1))), 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 )
+      CALL zdf_osm_vertical_average( Kbb,                       Kmm,                         &
+         &                           nbld(A2D((nn_hls-1))),     av_t_bl(A2D((nn_hls-1))),    &
+         &                           av_s_bl(A2D((nn_hls-1))),  av_b_bl(A2D((nn_hls-1))),    &
+         &                           av_u_bl(A2D((nn_hls-1))),  av_v_bl(A2D((nn_hls-1))),    &
+         &                           jk_ext,                    av_dt_bl(A2D((nn_hls-1))),   &
+         &                           av_ds_bl(A2D((nn_hls-1))), av_db_bl(A2D((nn_hls-1))),   &
+         &                           av_du_bl(A2D((nn_hls-1))), av_dv_bl(A2D((nn_hls-1))) )
       !
       CALL zdf_osm_pycnocline_thickness( Kmm, zdh, zhml, zdhdt, zhbl,   &
@@ -739,5 +794,5 @@
       !
       ! Reset l_pyc before calculating terms in the flux-gradient relationship
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          IF ( av_db_bl(ji,jj) < rn_osm_bl_thresh .OR. nbld(ji,jj) >= mbkt(ji,jj) - 2 .OR.   &
             & nbld(ji,jj) - nmld(ji,jj) == 1   .OR. zdhdt(ji,jj) < 0.0_wp ) THEN   ! ag 19/03
@@ -753,6 +808,8 @@
       END_2D
       !
-      dstokes(:,:) = MIN ( dstokes(:,:), hbl(A2D((nn_hls-1)))/ 3.0_wp )   ! Limit delta for shallow boundary layers for calculating
-      !                                                       !    flux-gradient terms
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )               ! Limit delta for shallow boundary layers for calculating
+         dstokes(ji,jj) = MIN ( dstokes(ji,jj), hbl(ji,jj) / 3.0_wp )   !    flux-gradient terms
+      END_2D
+      !                                                       
       !
       ! Average over the depth of the mixed layer in the convective boundary layer
@@ -760,14 +817,18 @@
       ! Recalculate ML averages and differences using new ML depth
       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, nmld - 1, 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 )
+      CALL zdf_osm_vertical_average( Kbb,                        Kmm,                         &
+         &                           nmld(A2D((nn_hls-1))) - 1,  av_t_ml(A2D((nn_hls-1))),    &
+         &                           av_s_ml(A2D((nn_hls-1))),   av_b_ml(A2D((nn_hls-1))),    &
+         &                           av_u_ml(A2D((nn_hls-1))),   av_v_ml(A2D((nn_hls-1))),    &
+         &                           jk_ext,                     av_dt_ml(A2D((nn_hls-1))),   &
+         &                           av_ds_ml(A2D((nn_hls-1))),  av_db_ml(A2D((nn_hls-1))),   &
+         &                           av_du_ml(A2D((nn_hls-1))),  av_dv_ml(A2D((nn_hls-1))) )
       !
       CALL zdf_osm_external_gradients( Kmm, nbld(A2D((nn_hls-1))) + 1, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext )
       ! Rotate mean currents and changes onto wind aligned co-ordinates
-      CALL zdf_osm_velocity_rotation( av_u_ml, av_v_ml )
-      CALL zdf_osm_velocity_rotation( av_du_ml, av_dv_ml )
-      CALL zdf_osm_velocity_rotation( av_u_bl, av_v_bl )
-      CALL zdf_osm_velocity_rotation( av_du_bl, av_dv_bl )
+      CALL zdf_osm_velocity_rotation( av_u_ml(A2D((nn_hls-1))),  av_v_ml(A2D((nn_hls-1)))  )
+      CALL zdf_osm_velocity_rotation( av_du_ml(A2D((nn_hls-1))), av_dv_ml(A2D((nn_hls-1))) )
+      CALL zdf_osm_velocity_rotation( av_u_bl(A2D((nn_hls-1))),  av_v_bl(A2D((nn_hls-1)))  )
+      CALL zdf_osm_velocity_rotation( av_du_bl(A2D((nn_hls-1))), av_dv_bl(A2D((nn_hls-1))) )
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship
@@ -791,10 +852,12 @@
       !
       ! Rotate non-gradient velocity terms back to model reference frame
-      CALL zdf_osm_velocity_rotation( ghamu(A2D((nn_hls-1)),:), ghamv(A2D((nn_hls-1)),:), .FALSE., 2, nbld(A2D((nn_hls-1))) )
+      CALL zdf_osm_velocity_rotation( ghamu(A2D((nn_hls-1)),:), ghamv(A2D((nn_hls-1)),:),   &
+         &                            .FALSE.,                  2,                          &
+         &                            nbld(A2D((nn_hls-1))) )
       !
       ! KPP-style Ri# mixing
       IF ( ln_kpprimix ) THEN
          jkflt = jpk
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             IF ( nbld(ji,jj) < jkflt ) jkflt = nbld(ji,jj)
          END_2D
@@ -807,5 +870,5 @@
                   &          wvmask(ji,jj,jk) / ( e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) )
             END_2D
-            DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+            DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
                IF ( jk > nbld(ji,jj) ) THEN
                   ! Shear prod. at w-point weightened by mask
@@ -826,5 +889,5 @@
       ! KPP-style set diffusivity large if unstable below BL
       IF ( ln_convmix) THEN
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             DO jk = nbld(ji,jj) + 1, jpkm1
                IF ( MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1e-12_wp ) zdiffut(ji,jj,jk) = MAX( rn_difconv, zdiffut(ji,jj,jk) )
@@ -834,5 +897,5 @@
       !
       IF ( ln_osm_mle ) THEN   ! Set up diffusivity and non-gradient mixing
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             IF ( l_flux(ji,jj) ) THEN   ! MLE mixing extends below boundary layer
                ! Calculate MLE flux contribution from surface fluxes
@@ -867,23 +930,8 @@
       ! CALL lbc_lnk( 'zdfosm', zviscos(:,:,:), 'W', 1.0_wp )
       ! GN 25/8: need to change tmask --> wmask
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jpkm1 )
          p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk)
          p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk)
       END_3D
-      ! Lateral boundary conditions on ghamu and ghamv, currently on W-grid (sign unchanged), needed to caclulate gham[uv] on u and
-      !    v grids
-      IF ( nn_hls == 1 ) CALL lbc_lnk( 'zdfosm', ghamu, 'W', 1.0_wp,   &
-         &                                       ghamv, 'W', 1.0_wp )
-      DO_3D( 0, 0, 0, 0, 2, jpkm1 )
-         ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) / MAX( 1.0_wp, tmask(ji,jj,jk) + tmask (ji+1,jj,jk) ) *   &
-            &              umask(ji,jj,jk)
-         ghamv(ji,jj,jk) = ( ghamv(ji,jj,jk) + ghamv(ji,jj+1,jk) ) / MAX( 1.0_wp, tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) *   &
-            &              vmask(ji,jj,jk)
-         ghamt(ji,jj,jk) =  ghamt(ji,jj,jk) * tmask(ji,jj,jk)
-         ghams(ji,jj,jk) =  ghams(ji,jj,jk) * tmask(ji,jj,jk)
-      END_3D
-      ! Lateral boundary conditions on final outputs for hbl, on T-grid (sign unchanged)
-      CALL lbc_lnk( 'zdfosm', hbl,  'T', 1.0_wp,   &
-         &                    hmle, 'T', 1.0_wp )
       !
       IF ( ln_dia_osm ) THEN
@@ -891,157 +939,98 @@
             ! Stokes drift set by assumimg onstant La#=0.3 (=0) or Pierson-Moskovitz spectrum (=1)
          CASE(0:1)
-            IF ( iom_use("us_x") ) THEN                                                           ! x surface Stokes drift
-               osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * sustke * scos_wind
-               CALL iom_put( "us_x", osmdia2d )
-            END IF
-            IF ( iom_use("us_y") ) THEN                                                           ! y surface Stokes drift
-               osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * sustke * ssin_wind
-               CALL iom_put( "us_y", osmdia2d )
-            END IF
-            IF ( iom_use("wind_wave_abs_power") ) THEN
-               osmdia2d(A2D((nn_hls-1))) = 1000.0_wp * rho0 * tmask(A2D((nn_hls-1)),1) * sustar**2 * sustke
-               CALL iom_put( "wind_wave_abs_power", osmdia2d )
-            END IF
+            CALL zdf_osm_iomput( "us_x", tmask(A2D(0),1) * sustke(A2D(0)) * scos_wind(A2D(0)) )   ! x surface Stokes drift
+            CALL zdf_osm_iomput( "us_y", tmask(A2D(0),1) * sustke(A2D(0)) * scos_wind(A2D(0)) )   ! y surface Stokes drift
+            CALL zdf_osm_iomput( "wind_wave_abs_power", 1000.0_wp * rho0 * tmask(A2D(0),1) * sustar(A2D(0))**2 * sustke(A2D(0)) )
             ! Stokes drift read in from sbcwave  (=2).
          CASE(2:3)
-            IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd*umask(:,:,1) )                     ! x surface Stokes drift
-            IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd*vmask(:,:,1) )                     ! y surface Stokes drift
-            IF ( iom_use("wmp") ) CALL iom_put( "wmp", wmp*tmask(:,:,1) )                         ! Wave mean period
-            IF ( iom_use("hsw") ) CALL iom_put( "hsw", hsw*tmask(:,:,1) )                         ! Significant wave height
-            IF ( iom_use("wmp_NP") ) CALL iom_put( "wmp_NP", ( 2.0_wp * rpi * 1.026_wp /      &   ! Wave mean period from NP
-               &                                               ( 0.877_wp * grav ) ) *        &   !    spectrum
-               &                                               wndm * tmask(:,:,1) )
-            IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", ( 0.22_wp / grav ) * wndm**2 *   &   ! Significant wave
-               &                                             tmask(:,:,1) )                       !    height from NP spectrum
-            IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) )                      ! U_10
-            IF ( iom_use("wind_wave_abs_power") ) THEN
-               osmdia2d(A2D((nn_hls-1))) = 1000.0_wp * rho0 * tmask(A2D((nn_hls-1)),1) * sustar**2 * SQRT( ut0sd(A2D((nn_hls-1)))**2 + vt0sd(A2D((nn_hls-1)))**2 )
-               CALL iom_put( "wind_wave_abs_power", osmdia2d )
-            END IF
+            CALL zdf_osm_iomput( "us_x",   ut0sd(A2D(0)) * umask(A2D(0),1) )                         ! x surface Stokes drift
+            CALL zdf_osm_iomput( "us_y",   vt0sd(A2D(0)) * vmask(A2D(0),1) )                         ! y surface Stokes drift
+            CALL zdf_osm_iomput( "wmp",    wmp(A2D(0)) * tmask(A2D(0),1) )                           ! Wave mean period
+            CALL zdf_osm_iomput( "hsw",    hsw(A2D(0)) * tmask(A2D(0),1) )                           ! Significant wave height
+            CALL zdf_osm_iomput( "wmp_NP", ( 2.0_wp * rpi * 1.026_wp / ( 0.877_wp * grav ) ) *   &   ! Wave mean period from NP
+               &                           wndm(A2D(0)) * tmask(A2D(0),1) )                          !    spectrum
+            CALL zdf_osm_iomput( "hsw_NP", ( 0.22_wp / grav ) * wndm**2 * tmask(A2D(0),1) )          ! Significant wave height from
+            !                                                                                        !    NP spectrum
+            CALL zdf_osm_iomput( "wndm",   wndm(A2D(0)) * tmask(A2D(0),1) )                          ! U_10
+            CALL zdf_osm_iomput( "wind_wave_abs_power", 1000.0_wp * rho0 * tmask(A2D(0),1) * sustar(A2D(0))**2 *   &
+               &                                        SQRT( ut0sd(A2D(0))**2 + vt0sd(A2D(0))**2 ) )
          END SELECT
-         IF ( iom_use("ghamt") ) CALL iom_put( "ghamt", tmask*ghamt )                      ! <Tw_NL>
-         IF ( iom_use("ghams") ) CALL iom_put( "ghams", tmask*ghams )                      ! <Sw_NL>
-         IF ( iom_use("ghamu") ) CALL iom_put( "ghamu", umask*ghamu )                      ! <uw_NL>
-         IF ( iom_use("ghamv") ) CALL iom_put( "ghamv", vmask*ghamv )                      ! <vw_NL>
-         IF ( iom_use("zwth0") ) THEN                                                      ! <Tw_0>
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * swth0;     CALL iom_put( "zwth0",     osmdia2d )
-         END IF
-         IF ( iom_use("zws0") ) THEN                                                       ! <Sw_0>
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * sws0;      CALL iom_put( "zws0",      osmdia2d )
-         END IF
-         IF ( iom_use("zwb0") ) THEN                                                       ! <Sw_0>
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * swb0;      CALL iom_put( "zwb0",      osmdia2d )
-         END IF
-         IF ( iom_use("zwbav") ) THEN                                                      ! Upward BL-avged turb buoyancy flux
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * swth0;     CALL iom_put( "zwbav",     osmdia2d )
-         END IF
-         IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl )                     ! Boundary-layer depth
-         IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*nbld )                  ! Boundary-layer max k
-         IF ( iom_use("zdt_bl") ) THEN                                                     ! dt at ml base
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_dt_bl;  CALL iom_put( "zdt_bl", osmdia2d )
-         END IF
-         IF ( iom_use("zds_bl") ) THEN                                                     ! ds at ml base
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_ds_bl;  CALL iom_put( "zds_bl", osmdia2d )
-         END IF
-         IF ( iom_use("zdb_bl") ) THEN                                                     ! db at ml base
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_db_bl;  CALL iom_put( "zdb_bl", osmdia2d )
-         END IF
-         IF ( iom_use("zdu_bl") ) THEN                                                     ! du at ml base
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_du_bl;  CALL iom_put( "zdu_bl", osmdia2d )
-         END IF
-         IF ( iom_use("zdv_bl") ) THEN                                                     ! dv at ml base
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_dv_bl;  CALL iom_put( "zdv_bl", osmdia2d )
-         END IF
-         IF ( iom_use("dh") ) CALL iom_put( "dh", tmask(:,:,1)*dh )                        ! Initial boundary-layer depth
-         IF ( iom_use("hml") ) CALL iom_put( "hml", tmask(:,:,1)*hml )                     ! Initial boundary-layer depth
-         IF ( iom_use("zdt_ml") ) THEN                                                     ! dt at ml base
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_dt_ml;  CALL iom_put( "zdt_ml", osmdia2d )
-         END IF
-         IF ( iom_use("zds_ml") ) THEN                                                     ! ds at ml base
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_ds_ml;  CALL iom_put( "zds_ml", osmdia2d )
-         END IF
-         IF ( iom_use("zdb_ml") ) THEN                                                     ! db at ml base
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_db_ml;  CALL iom_put( "zdb_ml", osmdia2d )
-         END IF
-         IF ( iom_use("dstokes") ) THEN                                                    ! Stokes drift penetration depth
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * dstokes;   CALL iom_put( "dstokes",   osmdia2d )
-         END IF
-         IF ( iom_use("zustke") ) THEN                                                     ! Stokes drift magnitude at T-points
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * sustke;    CALL iom_put( "zustke",    osmdia2d )
-         END IF
-         IF ( iom_use("zwstrc") ) THEN                                                     ! Convective velocity scale
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * swstrc;    CALL iom_put( "zwstrc",    osmdia2d )
-         END IF
-         IF ( iom_use("zwstrl") ) THEN                                                     ! Langmuir velocity scale
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * swstrl;    CALL iom_put( "zwstrl",    osmdia2d )
-         END IF
-         IF ( iom_use("zustar") ) THEN                                                     ! Friction velocity scale
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * sustar;    CALL iom_put( "zustar",    osmdia2d )
-         END IF
-         IF ( iom_use("zvstr") ) THEN                                                      ! Mixed velocity scale
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * svstr;     CALL iom_put( "zvstr",     osmdia2d )
-         END IF
-         IF ( iom_use("zla") ) THEN                                                        ! Langmuir #
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * sla;       CALL iom_put( "zla",       osmdia2d )
-         END IF
-         IF ( iom_use("wind_power") ) THEN                                                 ! BL depth internal to zdf_osm routine
-            osmdia2d(A2D((nn_hls-1))) = 1000.0_wp * rho0 * tmask(A2D((nn_hls-1)),1) * sustar**3
-            CALL iom_put( "wind_power", osmdia2d )
-         END IF
-         IF ( iom_use("wind_wave_power") ) THEN
-            osmdia2d(A2D((nn_hls-1))) = 1000.0_wp * rho0 * tmask(A2D((nn_hls-1)),1) * sustar**2 * sustke
-            CALL iom_put( "wind_wave_power", osmdia2d )
-         END IF
-         IF ( iom_use("zhbl") ) THEN                                                       ! BL depth internal to zdf_osm routine
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zhbl;      CALL iom_put( "zhbl",      osmdia2d )
-         END IF
-         IF ( iom_use("zhml") ) THEN                                                       ! ML depth internal to zdf_osm routine
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zhml;      CALL iom_put( "zhml",      osmdia2d )
-         END IF
-         IF ( iom_use("imld") ) THEN                                                       ! Index for ML depth internal to zdf_osm routine
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * nmld;      CALL iom_put( "imld",      osmdia2d )
-         END IF
-         IF ( iom_use("jp_ext") ) THEN                                                     ! =1 if pycnocline resolved internal to zdf_osm routine
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * jk_ext;    CALL iom_put( "jp_ext",    osmdia2d )
-         END IF
-         IF ( iom_use("j_ddh") ) THEN                                                      ! Index forpyc thicknessh internal to zdf_osm routine
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * n_ddh;     CALL iom_put( "j_ddh",     osmdia2d )
-         END IF
-         IF ( iom_use("zshear") ) THEN                                                     ! Shear production of TKE internal to zdf_osm routine
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zshear;    CALL iom_put( "zshear",    osmdia2d )
-         END IF
-         IF ( iom_use("zdh") ) THEN                                                        ! Pyc thicknessh internal to zdf_osm routine
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zdh;       CALL iom_put( "zdh",       osmdia2d )
-         END IF
-         IF ( iom_use("zhol") ) THEN                                                       ! ML depth internal to zdf_osm routine
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * shol;      CALL iom_put( "zhol",      osmdia2d )
-         END IF
-         IF ( iom_use("zwb_ent") ) THEN                                                    ! Upward turb buoyancy entrainment flux
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zwb_ent;   CALL iom_put( "zwb_ent",   osmdia2d )
-         END IF
-         IF ( iom_use("zt_ml") ) THEN                                                      ! Average T in ML
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * av_t_ml;   CALL iom_put( "zt_ml",     osmdia2d )
-         END IF
-         IF ( iom_use("hmle") ) CALL iom_put( "hmle", tmask(:,:,1)*hmle )                  ! FK layer depth
-         IF ( iom_use("zmld") ) CALL iom_put( "zmld", tmask(:,:,1)*zmld )                  ! FK target layer depth
-         IF ( iom_use("zwb_fk") ) THEN                                                     ! FK b flux
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zwb_fk;    CALL iom_put( "zwb_fk",    osmdia2d )
-         END IF
-         IF ( iom_use("zwb_fk_b") ) THEN                                                   ! FK b flux averaged over ML
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zwb_fk_b;  CALL iom_put( "zwb_fk_b",  osmdia2d )
-         END IF
-         IF ( iom_use("mld_prof") ) CALL iom_put( "mld_prof", tmask(:,:,1)*mld_prof )      ! FK layer max k
-         IF ( iom_use("zdtdx") ) CALL iom_put( "zdtdx", umask(:,:,1)*zdtdx )               ! FK dtdx at u-pt
-         IF ( iom_use("zdtdy") ) CALL iom_put( "zdtdy", vmask(:,:,1)*zdtdy )               ! FK dtdy at v-pt
-         IF ( iom_use("zdsdx") ) CALL iom_put( "zdsdx", umask(:,:,1)*zdsdx )               ! FK dtdx at u-pt
-         IF ( iom_use("zdsdy") ) CALL iom_put( "zdsdy", vmask(:,:,1)*zdsdy )               ! FK dsdy at v-pt
-         IF ( iom_use("dbdx_mle") ) CALL iom_put( "dbdx_mle", umask(:,:,1)*dbdx_mle )      ! FK dbdx at u-pt
-         IF ( iom_use("dbdy_mle") ) CALL iom_put( "dbdy_mle", vmask(:,:,1)*dbdy_mle )      ! FK dbdy at v-pt
-         IF ( iom_use("zdiff_mle") ) THEN                                                  ! FK diff in MLE at t-pt
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zdiff_mle; CALL iom_put( "zdiff_mle", osmdia2d )
-         END IF
-         IF ( iom_use("zvel_mle") ) THEN                                                   ! FK diff in MLE at t-pt
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zdiff_mle; CALL iom_put( "zvel_mle",  osmdia2d )
-         END IF
+         CALL zdf_osm_iomput( "zwth0",           tmask(A2D(0),1) * swth0(A2D(0))     )      ! <Tw_0>
+         CALL zdf_osm_iomput( "zws0",            tmask(A2D(0),1) * sws0(A2D(0))      )      ! <Sw_0>
+         CALL zdf_osm_iomput( "zwb0",            tmask(A2D(0),1) * swb0(A2D(0))      )      ! <Sw_0>
+         CALL zdf_osm_iomput( "zwbav",           tmask(A2D(0),1) * swth0(A2D(0))     )      ! Upward BL-avged turb buoyancy flux
+         CALL zdf_osm_iomput( "ibld",            tmask(A2D(0),1) * nbld(A2D(0))      )      ! Boundary-layer max k
+         CALL zdf_osm_iomput( "zdt_bl",          tmask(A2D(0),1) * av_dt_bl(A2D(0))  )      ! dt at ml base
+         CALL zdf_osm_iomput( "zds_bl",          tmask(A2D(0),1) * av_ds_bl(A2D(0))  )      ! ds at ml base
+         CALL zdf_osm_iomput( "zdb_bl",          tmask(A2D(0),1) * av_db_bl(A2D(0))  )      ! db at ml base
+         CALL zdf_osm_iomput( "zdu_bl",          tmask(A2D(0),1) * av_du_bl(A2D(0))  )      ! du at ml base
+         CALL zdf_osm_iomput( "zdv_bl",          tmask(A2D(0),1) * av_dv_bl(A2D(0))  )      ! dv at ml base
+         CALL zdf_osm_iomput( "dh",              tmask(A2D(0),1) * dh(A2D(0))        )      ! Initial boundary-layer depth
+         CALL zdf_osm_iomput( "hml",             tmask(A2D(0),1) * hml(A2D(0))       )      ! Initial boundary-layer depth
+         CALL zdf_osm_iomput( "zdt_ml",          tmask(A2D(0),1) * av_dt_ml(A2D(0))  )      ! dt at ml base
+         CALL zdf_osm_iomput( "zds_ml",          tmask(A2D(0),1) * av_ds_ml(A2D(0))  )      ! ds at ml base
+         CALL zdf_osm_iomput( "zdb_ml",          tmask(A2D(0),1) * av_db_ml(A2D(0))  )      ! db at ml base
+         CALL zdf_osm_iomput( "dstokes",         tmask(A2D(0),1) * dstokes(A2D(0))   )      ! Stokes drift penetration depth
+         CALL zdf_osm_iomput( "zustke",          tmask(A2D(0),1) * sustke(A2D(0))    )      ! Stokes drift magnitude at T-points
+         CALL zdf_osm_iomput( "zwstrc",          tmask(A2D(0),1) * swstrc(A2D(0))    )      ! Convective velocity scale
+         CALL zdf_osm_iomput( "zwstrl",          tmask(A2D(0),1) * swstrl(A2D(0))    )      ! Langmuir velocity scale
+         CALL zdf_osm_iomput( "zustar",          tmask(A2D(0),1) * sustar(A2D(0))    )      ! Friction velocity scale
+         CALL zdf_osm_iomput( "zvstr",           tmask(A2D(0),1) * svstr(A2D(0))     )      ! Mixed velocity scale
+         CALL zdf_osm_iomput( "zla",             tmask(A2D(0),1) * sla(A2D(0))       )      ! Langmuir #
+         CALL zdf_osm_iomput( "wind_power",      1000.0_wp * rho0 * tmask(A2D(0),1) *   &   ! BL depth internal to zdf_osm routine
+            &                                    sustar(A2D(0))**3 )
+         CALL zdf_osm_iomput( "wind_wave_power", 1000.0_wp * rho0 * tmask(A2D(0),1) *   &
+            &                                    sustar(A2D(0))**2 * sustke(A2D(0))  )
+         CALL zdf_osm_iomput( "zhbl",            tmask(A2D(0),1) * zhbl(A2D(0))      )      ! BL depth internal to zdf_osm routine
+         CALL zdf_osm_iomput( "zhml",            tmask(A2D(0),1) * zhml(A2D(0))      )      ! ML depth internal to zdf_osm routine
+         CALL zdf_osm_iomput( "imld",            tmask(A2D(0),1) * nmld(A2D(0))      )      ! Index for ML depth internal to zdf_osm routine
+         CALL zdf_osm_iomput( "jp_ext",          tmask(A2D(0),1) * jk_ext(A2D(0))    )      ! =1 if pycnocline resolved internal to zdf_osm routine
+         CALL zdf_osm_iomput( "j_ddh",           tmask(A2D(0),1) * n_ddh(A2D(0))     )      ! Index forpyc thicknessh internal to zdf_osm routine
+         CALL zdf_osm_iomput( "zshear",          tmask(A2D(0),1) * zshear(A2D(0))    )      ! Shear production of TKE internal to zdf_osm routine
+         CALL zdf_osm_iomput( "zdh",             tmask(A2D(0),1) * zdh(A2D(0))       )      ! Pyc thicknessh internal to zdf_osm routine
+         CALL zdf_osm_iomput( "zhol",            tmask(A2D(0),1) * shol(A2D(0))      )      ! ML depth internal to zdf_osm routine
+         CALL zdf_osm_iomput( "zwb_ent",         tmask(A2D(0),1) * zwb_ent(A2D(0))   )      ! Upward turb buoyancy entrainment flux
+         CALL zdf_osm_iomput( "zt_ml",           tmask(A2D(0),1) * av_t_ml(A2D(0))   )      ! Average T in ML
+         CALL zdf_osm_iomput( "zmld",            tmask(A2D(0),1) * zmld(A2D(0))      )      ! FK target layer depth
+         CALL zdf_osm_iomput( "zwb_fk",          tmask(A2D(0),1) * zwb_fk(A2D(0))    )      ! FK b flux
+         CALL zdf_osm_iomput( "zwb_fk_b",        tmask(A2D(0),1) * zwb_fk_b(A2D(0))  )      ! FK b flux averaged over ML
+         CALL zdf_osm_iomput( "mld_prof",        tmask(A2D(0),1) * mld_prof(A2D(0))  )      ! FK layer max k
+         CALL zdf_osm_iomput( "zdtdx",           umask(A2D(0),1) * zdtdx(A2D(0))     )      ! FK dtdx at u-pt
+         CALL zdf_osm_iomput( "zdtdy",           vmask(A2D(0),1) * zdtdy(A2D(0))     )      ! FK dtdy at v-pt
+         CALL zdf_osm_iomput( "zdsdx",           umask(A2D(0),1) * zdsdx(A2D(0))     )      ! FK dtdx at u-pt
+         CALL zdf_osm_iomput( "zdsdy",           vmask(A2D(0),1) * zdsdy(A2D(0))     )      ! FK dsdy at v-pt
+         CALL zdf_osm_iomput( "dbdx_mle",        umask(A2D(0),1) * dbdx_mle(A2D(0))  )      ! FK dbdx at u-pt
+         CALL zdf_osm_iomput( "dbdy_mle",        vmask(A2D(0),1) * dbdy_mle(A2D(0))  )      ! FK dbdy at v-pt
+         CALL zdf_osm_iomput( "zdiff_mle",       tmask(A2D(0),1) * zdiff_mle(A2D(0)) )      ! FK diff in MLE at t-pt
+         CALL zdf_osm_iomput( "zvel_mle",        tmask(A2D(0),1) * zdiff_mle(A2D(0)) )      ! FK diff in MLE at t-pt
+      END IF
+      !
+      ! Lateral boundary conditions on ghamu and ghamv, currently on W-grid (sign unchanged), needed to caclulate gham[uv] on u and
+      !    v grids
+      IF ( .NOT. l_istiled .OR. ntile == nijtile ) THEN   ! Finalise ghamu, ghamv, hbl, and hmle only after full domain has been processed
+         IF ( nn_hls == 1 ) CALL lbc_lnk( 'zdfosm', ghamu, 'W', 1.0_wp,   &
+            &                                       ghamv, 'W', 1.0_wp )
+         DO jk = 2, jpkm1
+            DO jj = Njs0, Nje0
+               DO ji = Nis0, Nie0
+                  ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) / MAX( 1.0_wp, tmask(ji,jj,jk) + tmask (ji+1,jj,jk) ) *   &
+                     &              umask(ji,jj,jk)
+                  ghamv(ji,jj,jk) = ( ghamv(ji,jj,jk) + ghamv(ji,jj+1,jk) ) / MAX( 1.0_wp, tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) *   &
+                     &              vmask(ji,jj,jk)
+                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) * tmask(ji,jj,jk)
+                  ghams(ji,jj,jk) = ghams(ji,jj,jk) * tmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         ! Lateral boundary conditions on final outputs for hbl, on T-grid (sign unchanged)
+         CALL lbc_lnk( 'zdfosm', hbl,  'T', 1.0_wp,   &
+            &                    hmle, 'T', 1.0_wp )
+         !
+         CALL zdf_osm_iomput( "ghamt", tmask*ghamt )          ! <Tw_NL>
+         CALL zdf_osm_iomput( "ghams", tmask*ghams )          ! <Sw_NL>
+         CALL zdf_osm_iomput( "ghamu", umask*ghamu )          ! <uw_NL>
+         CALL zdf_osm_iomput( "ghamv", vmask*ghamv )          ! <vw_NL>
+         CALL zdf_osm_iomput( "hbl",   tmask(:,:,1) * hbl )   ! Boundary-layer depth
+         CALL zdf_osm_iomput( "hmle",  tmask(:,:,1)*hmle )    ! FK layer depth
       END IF
       !
@@ -1079,16 +1068,18 @@
       !
       ! Averages over depth of boundary layer
-      pt(:,:)     = 0.0_wp
-      ps(:,:)     = 0.0_wp
-      pu(:,:)     = 0.0_wp
-      pv(:,:)     = 0.0_wp
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         pt(ji,jj) = 0.0_wp
+         ps(ji,jj) = 0.0_wp
+         pu(ji,jj) = 0.0_wp
+         pv(ji,jj) = 0.0_wp
+      END_2D
       zthick(:,:) = epsln
       jkflt = jpk
       jkmax = 0
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          IF ( knlev(ji,jj) < jkflt ) jkflt = knlev(ji,jj)
          IF ( knlev(ji,jj) > jkmax ) jkmax = knlev(ji,jj)
       END_2D
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkflt )   ! Upper, flat part of layer
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkflt )   ! Upper, flat part of layer
          zthick(ji,jj) = zthick(ji,jj) + e3t(ji,jj,jk,Kmm)
          pt(ji,jj)     = pt(ji,jj)     + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm)
@@ -1101,5 +1092,5 @@
             &                               MAX( 1.0_wp , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) )         
       END_3D
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jkflt+1, jkmax )   ! Lower, non-flat part of layer
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jkflt+1, jkmax )   ! Lower, non-flat part of layer
          IF ( knlev(ji,jj) >= jk ) THEN
             zthick(ji,jj) = zthick(ji,jj) + e3t(ji,jj,jk,Kmm)
@@ -1114,5 +1105,5 @@
          END IF
       END_3D
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          pt(ji,jj) = pt(ji,jj) / zthick(ji,jj)
          ps(ji,jj) = ps(ji,jj) / zthick(ji,jj)
@@ -1126,5 +1117,5 @@
       ! Differences between vertical averages and values at an external layer
       IF ( PRESENT( kp_ext ) ) THEN
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             ibld_ext = knlev(ji,jj) + kp_ext(ji,jj)
             IF ( ibld_ext <= mbkt(ji,jj)-1 ) THEN   ! ag 09/03
@@ -1171,5 +1162,5 @@
       zfwd = 1.0_wp
       IF( PRESENT(fwd) .AND. ( .NOT. fwd ) ) zfwd = -1.0_wp
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          ztmp      = pu(ji,jj)
          pu(ji,jj) = pu(ji,jj) * scos_wind(ji,jj) + zfwd * pv(ji,jj) * ssin_wind(ji,jj)
@@ -1208,5 +1199,5 @@
       IF( PRESENT(knlev) ) THEN
          jkmax = 0
-         DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
             IF ( knlev(ji,jj) > jkmax ) jkmax = knlev(ji,jj)
          END_2D
@@ -1216,5 +1207,5 @@
          llkbot = .TRUE.
       END IF
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jktop, jkmax )
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jktop, jkmax )
          IF ( llkbot .OR. knlev(ji,jj) >= jk ) THEN
             ztmp         = pu(ji,jj,jk)
@@ -1264,13 +1255,17 @@
       !
       ! Initialise arrays
-      l_conv(:,:)  = .FALSE.
-      l_shear(:,:) = .FALSE.
-      n_ddh(:,:)   = 1
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         l_conv(ji,jj)  = .FALSE.
+         l_shear(ji,jj) = .FALSE.
+         n_ddh(ji,jj)   = 1
+      END_2D
       ! Initialise INTENT(  out) arrays
-      pwb_ent(:,:) = pp_large
-      pwb_min(:,:) = pp_large
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         pwb_ent(ji,jj) = pp_large
+         pwb_min(ji,jj) = pp_large
+      END_2D
       !
       ! Determins stability and set flag l_conv
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          IF ( shol(ji,jj) < 0.0_wp ) THEN
             l_conv(ji,jj) = .TRUE.
@@ -1280,8 +1275,10 @@
       END_2D
       !
-      pshear(:,:) = 0.0_wp
-      zekman(:,:) = EXP( -1.0_wp * pp_ek * ABS( ff_t(A2D((nn_hls-1))) ) * phbl(:,:) / MAX( sustar(A2D((nn_hls-1))), 1.e-8 ) )
-      !
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         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 ) )
+      !
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          IF ( l_conv(ji,jj) ) THEN
             IF ( av_db_bl(ji,jj) > 0.0_wp ) THEN
@@ -1335,5 +1332,5 @@
       !
       ! Calculate entrainment buoyancy flux due to surface fluxes.
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          IF ( l_conv(ji,jj) ) THEN
             zwcor        = ABS( ff_t(ji,jj) ) * phbl(ji,jj) + epsln
@@ -1356,5 +1353,5 @@
       END_2D
       !
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          IF ( l_shear(ji,jj) ) THEN
             IF ( l_conv(ji,jj) ) THEN
@@ -1404,9 +1401,11 @@
       REAL(wp), PARAMETER ::   pp_large = -1e10_wp
       !
-      pdtdz(:,:) = pp_large
-      pdsdz(:,:) = pp_large
-      pdbdz(:,:) = pp_large
-      !
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         pdtdz(ji,jj) = pp_large
+         pdsdz(ji,jj) = pp_large
+         pdbdz(ji,jj) = pp_large
+      END_2D
+      !
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          IF ( kbase(ji,jj)+1 < mbkt(ji,jj) ) THEN
             zthermal = rab_n(ji,jj,1,jp_tem)   ! Ideally use nbld not 1??
@@ -1455,8 +1454,10 @@
       REAL(wp), PARAMETER ::   pp_large   = -1e10_wp
       !
-      pdhdt(:,:)    = pp_large
-      pwb_fk_b(:,:) = pp_large
-      !
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         pdhdt(ji,jj)    = pp_large
+         pwb_fk_b(ji,jj) = pp_large
+      END_2D
+      !
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          !
          IF ( l_shear(ji,jj) ) THEN
@@ -1612,5 +1613,5 @@
       REAL(wp) ::   zthermal, zbeta
       !
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          IF ( nbld(ji,jj) - nmld(ji,jj) > 1 ) THEN
             !
@@ -1718,5 +1719,5 @@
       REAL, PARAMETER ::   pp_ddh = 2.5_wp, pp_ddh_2 = 3.5_wp   ! Also in pycnocline_depth
       !
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          !
          IF ( l_shear(ji,jj) ) THEN
@@ -1881,8 +1882,10 @@
       REAL(wp), PARAMETER ::   pp_large   = -1e10_wp
       !
-      pdbdz(:,:,:) = pp_large
-      palpha(:,:)  = pp_large
-      !
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )      
+         pdbdz(ji,jj,:) = pp_large
+         palpha(ji,jj)  = pp_large
+      END_2D
+      !
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          !
          IF ( nbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN
@@ -1952,7 +1955,5 @@
       !
       IF ( ln_dia_pyc_scl ) THEN   ! Output of pycnocline gradient profiles
-         IF ( iom_use("zdbdz_pyc") ) THEN
-            osmdia3d(A2D((nn_hls-1)),:) = wmask(A2D((nn_hls-1)),:) * pdbdz(:,:,:); CALL iom_put( "zdbdz_pyc", osmdia3d )
-         END IF
+         CALL zdf_osm_iomput( "zdbdz_pyc", wmask(A2D(0),:) * pdbdz(A2D(0),:) )
       END IF
       !
@@ -2004,5 +2005,5 @@
       zb_coup(:,:) = 0.0_wp
       !
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          IF ( l_conv(ji,jj) ) THEN
             !
@@ -2083,5 +2084,5 @@
       END_2D
       !
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          IF ( l_conv(ji,jj) ) THEN
             DO jk = 2, nmld(ji,jj)   ! Mixed layer diffusivity
@@ -2152,7 +2153,5 @@
          !
       END_2D
-      IF( iom_use("pb_coup") ) THEN
-         osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zb_coup(:,:); CALL iom_put( "pb_coup", osmdia2d )   ! BBL-coupling velocity scale
-      END IF
+      CALL zdf_osm_iomput( "pb_coup", tmask(A2D(0),1) * zb_coup(A2D(0)) )   ! BBL-coupling velocity scale
       !
    END SUBROUTINE zdf_osm_diffusivity_viscosity
@@ -2223,5 +2222,5 @@
       jkf_mld = jpk
       jkm_mld = 0
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          IF ( nbld(ji,jj) > jkm_bld ) jkm_bld = nbld(ji,jj)
          IF ( nmld(ji,jj) < jkf_mld ) jkf_mld = nmld(ji,jj)
@@ -2238,5 +2237,5 @@
          zsc_ws_1(:,:)  = 2.0_wp * swsav(A2D((nn_hls-1)))
       ENDWHERE
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) )
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) )
          IF ( l_conv(ji,jj) ) THEN
             IF ( jk <= nmld(ji,jj) ) THEN
@@ -2259,6 +2258,6 @@
       !
       IF ( ln_dia_osm ) THEN
-         IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask * ghamu )
-         IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask * ghamv )
+         CALL zdf_osm_iomput( "ghamu_00", wmask(A2D(0),:) * ghamu(A2D(0),:) )
+         CALL zdf_osm_iomput( "ghamv_00", wmask(A2D(0),:) * ghamv(A2D(0),:) )
       END IF
       !
@@ -2278,5 +2277,5 @@
             &            ( svstr(A2D((nn_hls-1)))**2 + epsln )
       ENDWHERE
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) )
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) )
          IF ( l_conv(ji,jj) ) THEN
             IF ( jk <= nmld(ji,jj) ) THEN
@@ -2309,5 +2308,5 @@
          zsc_ws_1(:,:)  = 0.0_wp
       ENDWHERE
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) )
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) )
          IF ( l_conv(ji,jj) ) THEN
             IF ( jk <= nmld(ji,jj) ) THEN
@@ -2330,5 +2329,5 @@
          END IF
       END_3D
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          IF ( l_conv(ji,jj) .AND. l_pyc(ji,jj) ) THEN
             ztau_sc_u(ji,jj)    = phml(ji,jj) / ( svstr(ji,jj)**3 + swstrc(ji,jj)**3 )**pthird *                             &
@@ -2350,5 +2349,5 @@
          END IF
       END_2D
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkm_bld )
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkm_bld )
          IF ( l_conv(ji,jj) .AND. l_pyc(ji,jj) .AND. ( jk <= nbld(ji,jj) ) ) THEN
             zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj)
@@ -2371,10 +2370,6 @@
       !
       IF ( ln_dia_osm ) THEN
-         IF ( iom_use("zwth_ent") ) THEN   ! Upward turb. temperature entrainment flux
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zwth_ent(:,:); CALL iom_put( "zwth_ent", osmdia2d )
-         END IF
-         IF ( iom_use("zws_ent")  ) THEN   ! Upward turb. salinity entrainment flux
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zws_ent(:,:);  CALL iom_put( "zws_ent", osmdia2d )
-         END IF
+         CALL iom_put( "zwth_ent", tmask(A2D(0),1) * zwth_ent(A2D(0)) )   ! Upward turb. temperature entrainment flux
+         CALL iom_put( "zws_ent",  tmask(A2D(0),1) * zws_ent(A2D(0))  )   ! Upward turb. salinity entrainment flux
       END IF
       !
@@ -2388,5 +2383,5 @@
          zsc_uw_1(:,:) = 0.0_wp
       ENDWHERE
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) )
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) )
          IF ( l_conv(ji,jj) ) THEN
             IF ( jk <= nmld(ji,jj) ) THEN
@@ -2405,5 +2400,5 @@
       END_3D
       !
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          IF ( l_conv(ji,jj) .AND. l_pyc(ji,jj) ) THEN
             IF ( n_ddh(ji,jj) == 0 ) THEN
@@ -2424,5 +2419,5 @@
          END IF
       END_2D
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jkf_mld, jkm_bld )   ! Need ztau_sc_u to be available. Change to array.
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, jkf_mld, jkm_bld )   ! Need ztau_sc_u to be available. Change to array.
          IF ( l_conv(ji,jj) .AND. l_pyc(ji,jj) .AND. ( jk >= nmld(ji,jj) ) .AND. ( jk <= nbld(ji,jj) ) ) THEN
             zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj)
@@ -2437,8 +2432,6 @@
       !
       IF ( ln_dia_osm ) THEN
-         IF ( iom_use("ghamu_0") )    CALL iom_put( "ghamu_0", wmask * ghamu )
-         IF ( iom_use("zsc_uw_1_0") ) THEN
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zsc_uw_1(:,:); CALL iom_put( "zsc_uw_1_0", osmdia2d )
-         END IF
+         CALL iom_put( "ghamu_0",    wmask(A2D(0),:) * ghamu(A2D(0),:)  )
+         CALL iom_put( "zsc_uw_1_0", tmask(A2D(0),1) * zsc_uw_1(A2D(0)) )
       END IF
       !
@@ -2457,5 +2450,5 @@
          zsc_ws_1(:,:)  =          sws0(A2D((nn_hls-1)))
       END WHERE
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, MAX( jkm_mld, jkm_bld ) )
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 1, MAX( jkm_mld, jkm_bld ) )
          IF ( l_conv(ji,jj) ) THEN
             IF ( ( jk > 1 ) .AND. ( jk <= nmld(ji,jj) ) ) THEN
@@ -2498,5 +2491,5 @@
             &            zsc_vw_1(:,:)
       ENDWHERE
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) )
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, MAX( jkm_mld, jkm_bld ) )
          IF ( l_conv(ji,jj) ) THEN
             IF ( jk <= nmld(ji,jj) ) THEN
@@ -2530,18 +2523,10 @@
       !
       IF ( ln_dia_osm ) THEN
-         IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask * ghamu )
-         IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask * ghamv )
-         IF ( iom_use("zsc_uw_1_f") ) THEN
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zsc_uw_1(:,:); CALL iom_put( "zsc_uw_1_f", osmdia2d )
-         END IF
-         IF ( iom_use("zsc_vw_1_f") ) THEN
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zsc_vw_1(:,:); CALL iom_put( "zsc_vw_1_f", osmdia2d )
-         END IF
-         IF ( iom_use("zsc_uw_2_f") ) THEN
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zsc_uw_2(:,:); CALL iom_put( "zsc_uw_2_f", osmdia2d )
-         END IF
-         IF ( iom_use("zsc_vw_2_f") ) THEN
-            osmdia2d(A2D((nn_hls-1))) = tmask(A2D((nn_hls-1)),1) * zsc_vw_2(:,:); CALL iom_put( "zsc_vw_2_f", osmdia2d )
-         END IF
+         CALL zdf_osm_iomput( "ghamu_f",    wmask(A2D(0),:) * ghamu(A2D(0),:)  )
+         CALL zdf_osm_iomput( "ghamv_f",    wmask(A2D(0),:) * ghamv(A2D(0),:)  )
+         CALL zdf_osm_iomput( "zsc_uw_1_f", tmask(A2D(0),1) * zsc_uw_1(A2D(0)) )
+         CALL zdf_osm_iomput( "zsc_vw_1_f", tmask(A2D(0),1) * zsc_vw_1(A2D(0)) )
+         CALL zdf_osm_iomput( "zsc_uw_2_f", tmask(A2D(0),1) * zsc_uw_2(A2D(0)) )
+         CALL zdf_osm_iomput( "zsc_vw_2_f", tmask(A2D(0),1) * zsc_vw_2(A2D(0)) )
       END IF
       !
@@ -2551,5 +2536,5 @@
       ! Make surface forced velocity non-gradient terms go to zero at the base
       ! of the boundary layer.
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkm_bld )
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkm_bld )
          IF ( ( .NOT. l_conv(ji,jj) ) .AND. ( jk <= nbld(ji,jj) ) ) THEN
             znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / phbl(ji,jj)   ! ALMG to think about
@@ -2567,10 +2552,10 @@
       !
       IF ( ln_dia_pyc_scl .OR. ln_dia_pyc_shr ) THEN   ! Allocate arrays for output of pycnocline gradient/shear profiles
-         ALLOCATE( z3ddz_pyc_1(jpi,jpj,jpk), z3ddz_pyc_2(jpi,jpj,jpk), STAT=istat )
+         ALLOCATE( z3ddz_pyc_1(A2D(nn_hls),jpk), z3ddz_pyc_2(A2D(nn_hls),jpk), STAT=istat )
          IF ( istat /= 0 ) CALL ctl_stop( 'zdf_osm: failed to allocate temporary arrays' )
          z3ddz_pyc_1(:,:,:) = 0.0_wp
          z3ddz_pyc_2(:,:,:) = 0.0_wp
       END IF
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkm_bld )
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkm_bld )
          IF ( l_conv (ji,jj) ) THEN
             ! Unstable conditions. Shouldn;t be needed with no pycnocline code.
@@ -2637,8 +2622,8 @@
       END_3D
       IF ( ln_dia_pyc_scl ) THEN   ! Output of pycnocline gradient profiles
-         IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask(:,:,:) * z3ddz_pyc_1(:,:,:) )
-         IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask(:,:,:) * z3ddz_pyc_2(:,:,:) )
+         CALL zdf_osm_iomput( "zdtdz_pyc", wmask(A2D(0),:) * z3ddz_pyc_1(A2D(0),:) )
+         CALL zdf_osm_iomput( "zdsdz_pyc", wmask(A2D(0),:) * z3ddz_pyc_2(A2D(0),:) )
       END IF
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkm_bld )
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 2, jkm_bld )
          IF ( .NOT. l_conv (ji,jj) ) THEN
             IF ( nbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN
@@ -2664,10 +2649,10 @@
       END_3D
       IF ( ln_dia_pyc_shr ) THEN   ! Output of pycnocline shear profiles
-         IF ( iom_use("dudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask(:,:,:) * z3ddz_pyc_1(:,:,:) )
-         IF ( iom_use("dvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask(:,:,:) * z3ddz_pyc_2(:,:,:) )
+         CALL zdf_osm_iomput( "zdudz_pyc", wmask(A2D(0),:) * z3ddz_pyc_1(A2D(0),:) )
+         CALL zdf_osm_iomput( "zdvdz_pyc", wmask(A2D(0),:) * z3ddz_pyc_2(A2D(0),:) )
       END IF
       IF ( ln_dia_osm ) THEN
-         IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask * ghamu )
-         IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask * ghamv )
+         CALL zdf_osm_iomput( "ghamu_b", wmask(A2D(0),:) * ghamu(A2D(0),:) )
+         CALL zdf_osm_iomput( "ghamv_b", wmask(A2D(0),:) * ghamv(A2D(0),:) )
       END IF
       IF ( ln_dia_pyc_scl .OR. ln_dia_pyc_shr ) THEN   ! Deallocate arrays used for output of pycnocline gradient/shear profiles
@@ -2675,5 +2660,5 @@
       END IF
       !
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          ghamt(ji,jj,nbld(ji,jj)) = 0.0_wp
          ghams(ji,jj,nbld(ji,jj)) = 0.0_wp
@@ -2683,9 +2668,7 @@
       !
       IF ( ln_dia_osm ) THEN
-         IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask * ghamu )
-         IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask * ghamv )
-         IF ( iom_use("zviscos") ) THEN
-            osmdia3d(A2D((nn_hls-1)),:) = wmask(A2D((nn_hls-1)),:) * pviscos; CALL iom_put( "zviscos", osmdia3d )
-         END IF
+         CALL zdf_osm_iomput( "ghamu_1", wmask(A2D(0),:) * ghamu(A2D(0),:)   )
+         CALL zdf_osm_iomput( "ghamv_1", wmask(A2D(0),:) * ghamv(A2D(0),:)   )
+         CALL zdf_osm_iomput( "zviscos", wmask(A2D(0),:) * pviscos(A2D(0),:) )
       END IF
       !
@@ -2693,5 +2676,5 @@
 
    SUBROUTINE zdf_osm_zmld_horizontal_gradients( Kmm, pmld, pdtdx, pdtdy, pdsdx,   &
-      &                                          pdsdy, pdbdx_mle, pdbdy_mle, pdbds_mle )
+      &                                          pdsdy, pdbds_mle )
       !!----------------------------------------------------------------------
       !!          ***  ROUTINE zdf_osm_zmld_horizontal_gradients  ***
@@ -2707,15 +2690,14 @@
       !!----------------------------------------------------------------------
       INTEGER,                      INTENT(in   ) ::   Kmm          ! Time-level index
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   pmld         ! == Estimated FK BLD used for MLE horizontal gradients == !
-      REAL(wp), DIMENSION(:,:),     INTENT(inout) ::   pdtdx        ! Horizontal gradient for Fox-Kemper parametrization
-      REAL(wp), DIMENSION(:,:),     INTENT(inout) ::   pdtdy        ! Horizontal gradient for Fox-Kemper parametrization
-      REAL(wp), DIMENSION(:,:),     INTENT(inout) ::   pdsdx        ! Horizontal gradient for Fox-Kemper parametrization
-      REAL(wp), DIMENSION(:,:),     INTENT(inout) ::   pdsdy        ! Horizontal gradient for Fox-Kemper parametrization
-      REAL(wp), DIMENSION(:,:),     INTENT(inout) ::   pdbdx_mle    ! MLE horiz gradients at u points
-      REAL(wp), DIMENSION(:,:),     INTENT(inout) ::   pdbdy_mle    ! MLE horiz gradients at v points
-      REAL(wp), DIMENSION(A2D((nn_hls-1))),  INTENT(inout) ::   pdbds_mle    ! Magnitude of horizontal buoyancy gradient
+      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
       !
       ! Local variables
       INTEGER                          ::   ji, jj, jk   ! Dummy loop indices
+      INTEGER,  DIMENSION(A2D(nn_hls)) ::   jk_mld_prof  ! Base level of MLE layer
       INTEGER                          ::   ikt, ikmax   ! Local integers      
       REAL(wp)                         ::   zc
@@ -2731,28 +2713,33 @@
       !
       ! ==  MLD used for MLE  ==!
-      mld_prof(:,:) = nlb10   ! Initialization to the number of w ocean point
-      pmld(:,:)  = 0.0_wp     ! Here hmlp used as a dummy variable, integrating vertically N^2
+      DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
+         jk_mld_prof(ji,jj) = nlb10    ! Initialization to the number of w ocean point
+         pmld(ji,jj)        = 0.0_wp   ! Here hmlp used as a dummy variable, integrating vertically N^2
+      END_2D
       zN2_c = grav * rn_osm_mle_rho_c * r1_rho0   ! Convert density criteria into N^2 criteria
       DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, nlb10, jpkm1 )
          ikt = mbkt(ji,jj)
          pmld(ji,jj) = pmld(ji,jj) + MAX( rn2b(ji,jj,jk), 0.0_wp ) * e3w(ji,jj,jk,Kmm)
-         IF( pmld(ji,jj) < zN2_c ) mld_prof(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level
+         IF( pmld(ji,jj) < zN2_c ) jk_mld_prof(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level
       END_3D
       DO_2D( nn_hls, nn_hls, nn_hls, nn_hls )
-         mld_prof(ji,jj) = MAX( mld_prof(ji,jj), nbld(ji,jj) )   ! Ensure mld_prof .ge. nbld
-         pmld(ji,jj)     = gdepw(ji,jj,mld_prof(ji,jj),Kmm)
-      END_2D
-      !
-      ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 )   ! Max level of the computation
+         jk_mld_prof(ji,jj) = MAX( jk_mld_prof(ji,jj), nbld(ji,jj) )   ! Ensure jk_mld_prof .ge. nbld
+         pmld(ji,jj)     = gdepw(ji,jj,jk_mld_prof(ji,jj),Kmm)
+      END_2D
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+         mld_prof(ji,jj) = jk_mld_prof(ji,jj)
+      END_2D
+      !
+      ikmax = MIN( MAXVAL( jk_mld_prof(A2D(nn_hls)) ), jpkm1 )   ! Max level of the computation
       ztm(:,:) = 0.0_wp
       zsm(:,:) = 0.0_wp
       DO_3D( nn_hls, nn_hls, nn_hls, nn_hls, 1, ikmax )
-         zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, mld_prof(ji,jj) - jk ), 1  ), KIND=wp )   ! zc being 0 outside the ML t-points
+         zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, jk_mld_prof(ji,jj) - jk ), 1  ), KIND=wp )   ! zc being 0 outside the ML t-points
          ztm(ji,jj) = ztm(ji,jj) + zc * ts(ji,jj,jk,jp_tem,Kmm)
          zsm(ji,jj) = zsm(ji,jj) + zc * ts(ji,jj,jk,jp_sal,Kmm)
       END_3D
       ! Average temperature and salinity
-      ztm(:,:) = ztm(:,:) / MAX( e3t(:,:,1,Kmm), pmld(:,:) )
-      zsm(:,:) = zsm(:,:) / MAX( e3t(:,:,1,Kmm), pmld(:,:) )
+      ztm(:,:) = ztm(:,:) / MAX( e3t(A2D(nn_hls),1,Kmm), pmld(A2D(nn_hls)) )
+      zsm(:,:) = zsm(:,:) / MAX( e3t(A2D(nn_hls),1,Kmm), pmld(A2D(nn_hls)) )
       ! Calculate horizontal gradients at u & v points
       zmld_midu(:,:)   =  0.0_wp
@@ -2776,13 +2763,13 @@
       CALL eos_rab( ztsm_midu, zmld_midu, zabu, Kmm )
       CALL eos_rab( ztsm_midv, zmld_midv, zabv, Kmm )
-      DO_2D( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 )
-         pdbdx_mle(ji,jj) = grav * ( pdtdx(ji,jj) * zabu(ji,jj,jp_tem) - pdsdx(ji,jj) * zabu(ji,jj,jp_sal) )
-      END_2D
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 )
-         pdbdy_mle(ji,jj) = grav * ( pdtdy(ji,jj) * zabv(ji,jj,jp_tem) - pdsdy(ji,jj) * zabv(ji,jj,jp_sal) )
+      DO_2D_OVR( nn_hls, nn_hls-1, nn_hls-1, nn_hls-1 )
+         dbdx_mle(ji,jj) = grav * ( pdtdx(ji,jj) * zabu(ji,jj,jp_tem) - pdsdx(ji,jj) * zabu(ji,jj,jp_sal) )
+      END_2D
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls, nn_hls-1 )
+         dbdy_mle(ji,jj) = grav * ( pdtdy(ji,jj) * zabv(ji,jj,jp_tem) - pdsdy(ji,jj) * zabv(ji,jj,jp_sal) )
       END_2D
       DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
-         pdbds_mle(ji,jj) = SQRT( 0.5_wp * ( pdbdx_mle(ji,  jj) * pdbdx_mle(ji,  jj) + pdbdy_mle(ji,jj  ) * pdbdy_mle(ji,jj  ) +   &
-            &                                pdbdx_mle(ji-1,jj) * pdbdx_mle(ji-1,jj) + pdbdy_mle(ji,jj-1) * pdbdy_mle(ji,jj-1) ) )
+         pdbds_mle(ji,jj) = SQRT( 0.5_wp * ( dbdx_mle(ji,  jj) * dbdx_mle(ji,  jj) + dbdy_mle(ji,jj  ) * dbdy_mle(ji,jj  ) +   &
+            &                                dbdx_mle(ji-1,jj) * dbdx_mle(ji-1,jj) + dbdy_mle(ji,jj-1) * dbdy_mle(ji,jj-1) ) )
       END_2D
       !
@@ -2826,12 +2813,12 @@
       REAL(wp)                    ::   zdbdz_mle_int
       !
-      znd_param(A2D((nn_hls-1))) = 0.0_wp
-      !
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      znd_param(:,:) = 0.0_wp
+      !
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          ztmp =  r1_ft(ji,jj) *  MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf
          pwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) * ztmp * pdbds_mle(ji,jj) * pdbds_mle(ji,jj)
       END_2D
       !
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          !
          IF ( l_conv(ji,jj) ) THEN
@@ -2861,5 +2848,5 @@
       !
       ! Diagnosis
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          !
          IF ( l_conv(ji,jj) ) THEN
@@ -2922,5 +2909,5 @@
       !!----------------------------------------------------------------------
       INTEGER,                     INTENT(in   ) ::   Kmm         ! Time-level index
-      REAL(wp), DIMENSION(:,:),    INTENT(in   ) ::   pmld        ! == Estimated FK BLD used for MLE horiz gradients == !
+      REAL(wp), DIMENSION(A2D((nn_hls-1))), 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
@@ -2942,5 +2929,5 @@
       !
       ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          IF ( l_conv(ji,jj) ) THEN
             ztmp =  r1_ft(ji,jj) * MIN( 111e3_wp, e1u(ji,jj) ) / rn_osm_mle_lf
@@ -2951,5 +2938,5 @@
       END_2D
       ! Timestep mixed layer eddy depth
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          IF ( l_mle(ji,jj) ) THEN   ! MLE layer growing
             ! Buoyancy gradient at base of MLE layer
@@ -2973,9 +2960,8 @@
       END_2D
       !
-      mld_prof(:,:) = 4
-      DO_3D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 5, jpkm1 )
+      DO_3D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1, 5, jpkm1 )
          IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN( mbkt(ji,jj), jk )
       END_3D
-      DO_2D( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
+      DO_2D_OVR( nn_hls-1, nn_hls-1, nn_hls-1, nn_hls-1 )
          phmle(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm)
       END_2D
@@ -3432,4 +3418,48 @@
    END SUBROUTINE dyn_osm
 
+   SUBROUTINE zdf_osm_iomput_2d( cdname, posmdia2d )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE zdf_osm_iomput_2d  ***
+      !!
+      !! ** Purpose :   Wrapper for subroutine iom_put that accepts 2D arrays
+      !!                with and without halo
+      !!
+      !!----------------------------------------------------------------------
+      CHARACTER(LEN=*),         INTENT(in   ) ::   cdname
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   posmdia2d
+      !
+      IF ( ln_dia_osm .AND. iom_use( cdname ) ) THEN
+         IF ( SIZE( posmdia2d, 1 ) == ntei-ntsi+1 .AND. SIZE( posmdia2d, 2 ) == ntej-ntsj+1 ) THEN   ! Halo absent
+            osmdia2d(A2D(0)) = posmdia2d(:,:)
+            CALL iom_put( cdname, osmdia2d(A2D(nn_hls)) )
+         ELSE   ! Halo present
+            CALL iom_put( cdname, osmdia2d )
+         END IF
+      END IF
+      !
+   END SUBROUTINE zdf_osm_iomput_2d
+
+   SUBROUTINE zdf_osm_iomput_3d( cdname, posmdia3d )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE zdf_osm_iomput_3d  ***
+      !!
+      !! ** Purpose :   Wrapper for subroutine iom_put that accepts 3D arrays
+      !!                with and without halo
+      !!
+      !!----------------------------------------------------------------------
+      CHARACTER(LEN=*),           INTENT(in   ) ::   cdname
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   posmdia3d
+      !
+      IF ( ln_dia_osm .AND. iom_use( cdname ) ) THEN
+         IF ( SIZE( posmdia3d, 1 ) == ntei-ntsi+1 .AND. SIZE( posmdia3d, 2 ) == ntej-ntsj+1 ) THEN   ! Halo absent
+            osmdia3d(A2D(0),:) = posmdia3d(:,:,:)
+            CALL iom_put( cdname, osmdia3d(A2D(nn_hls),:) )
+         ELSE   ! Halo present
+            CALL iom_put( cdname, osmdia3d )
+         END IF
+      END IF
+      !
+   END SUBROUTINE zdf_osm_iomput_3d
+
    !!======================================================================
 
Index: NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfphy.F90
===================================================================
--- NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfphy.F90	(revision 14899)
+++ NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfphy.F90	(revision 14900)
@@ -186,6 +186,4 @@
       IF( lk_top    .AND. ln_zdfnpc )   CALL ctl_stop( 'zdf_phy_init: npc scheme is not working with key_top' )
       IF( lk_top    .AND. ln_zdfosm )   CALL ctl_warn( 'zdf_phy_init: osmosis gives no non-local fluxes for TOP tracers yet' )
-      ! TEMP: [tiling] This change not necessary after finalisation of zdf_osm (not yet tiled)
-      IF( ln_tile   .AND. ln_zdfosm )   CALL ctl_warn( 'zdf_phy_init: osmosis does not yet work with tiling' )
       IF( lk_top    .AND. ln_zdfmfc )   CALL ctl_stop( 'zdf_phy_init: Mass Flux scheme is not working with key_top' )
       IF(lwp) THEN
@@ -256,20 +254,7 @@
       INTEGER ::   ji, jj, jk   ! dummy loop indice
       REAL(wp), DIMENSION(A2D(nn_hls),jpk) ::   zsh2   ! shear production
-      ! TEMP: [tiling] This change not necessary after finalisation of zdf_osm (not yet tiled)
-      LOGICAL :: lskip
       !! ---------------------------------------------------------------------
       !
       IF( ln_timing )   CALL timing_start('zdf_phy')
-
-      ! TEMP: [tiling] These changes not necessary after finalisation of zdf_osm (not yet tiled)
-      lskip = .FALSE.
-
-      IF( ln_tile .AND. nzdf_phy == np_OSM )  THEN
-         IF( ntile == 1 ) THEN
-            CALL dom_tile_stop( ldhold=.TRUE. )
-         ELSE
-            lskip = .TRUE.
-         ENDIF
-      ENDIF
       !
       IF( l_zdfdrg ) THEN     !==  update top/bottom drag  ==!   (non-linear cases)
@@ -301,34 +286,29 @@
       !
       CALL zdf_mxl( kt, Kmm )                        !* mixed layer depth, and level
-
-      ! TEMP: [tiling] These changes not necessary after finalisation of zdf_osm (not yet tiled)
-      IF( .NOT. lskip ) THEN
-         !                       !==  Kz from chosen turbulent closure  ==!   (avm_k, avt_k)
-         !
-         ! NOTE: [tiling] the closure schemes (zdf_tke etc) will update avm_k. With tiling, the calculation of zsh2 on adjacent tiles then uses both updated (next timestep) and non-updated (current timestep) values of avm_k. To preserve results, we save a read-only copy of the "now" avm_k to use in the calculation of zsh2.
-         IF( l_zdfsh2 ) THEN        !* shear production at w-points (energy conserving form)
-            IF( ln_tile ) THEN
-               IF( ntile == 1 ) avm_k_n(:,:,:) = avm_k(:,:,:)     ! Preserve "now" avm_k for calculation of zsh2
-               CALL zdf_sh2( Kbb, Kmm, avm_k_n, &     ! <<== in
-                  &                     zsh2    )     ! ==>> out : shear production
-            ELSE
-               CALL zdf_sh2( Kbb, Kmm, avm_k,   &     ! <<== in
-                  &                     zsh2    )     ! ==>> out : shear production
-            ENDIF
-         ENDIF
-         !
-         SELECT CASE ( nzdf_phy )                  !* Vertical eddy viscosity and diffusivity coefficients at w-points
-         CASE( np_RIC )   ;   CALL zdf_ric( kt,      Kmm, zsh2, avm_k, avt_k )    ! Richardson number dependent Kz
-         CASE( np_TKE )   ;   CALL zdf_tke( kt, Kbb, Kmm, zsh2, avm_k, avt_k )    ! TKE closure scheme for Kz
-         CASE( np_GLS )   ;   CALL zdf_gls( kt, Kbb, Kmm, zsh2, avm_k, avt_k )    ! GLS closure scheme for Kz
-         CASE( np_OSM )   ;   CALL zdf_osm( kt, Kbb, Kmm, Krhs, avm_k, avt_k )    ! OSMOSIS closure scheme for Kz
+      !
+      !                       !==  Kz from chosen turbulent closure  ==!   (avm_k, avt_k)
+      !
+      ! NOTE: [tiling] the closure schemes (zdf_tke etc) will update avm_k. With tiling, the calculation of zsh2 on adjacent tiles then uses both updated (next timestep) and non-updated (current timestep) values of avm_k. To preserve results, we save a read-only copy of the "now" avm_k to use in the calculation of zsh2.
+      IF( l_zdfsh2 ) THEN        !* shear production at w-points (energy conserving form)
+         IF( ln_tile ) THEN
+            IF( ntile == 1 ) avm_k_n(:,:,:) = avm_k(:,:,:)     ! Preserve "now" avm_k for calculation of zsh2
+            CALL zdf_sh2( Kbb, Kmm, avm_k_n, &     ! <<== in
+               &                     zsh2    )     ! ==>> out : shear production
+         ELSE
+            CALL zdf_sh2( Kbb, Kmm, avm_k,   &     ! <<== in
+               &                     zsh2    )     ! ==>> out : shear production
+         ENDIF
+      ENDIF
+      !
+      SELECT CASE ( nzdf_phy )                  !* Vertical eddy viscosity and diffusivity coefficients at w-points
+      CASE( np_RIC )   ;   CALL zdf_ric( kt,      Kmm, zsh2, avm_k, avt_k )    ! Richardson number dependent Kz
+      CASE( np_TKE )   ;   CALL zdf_tke( kt, Kbb, Kmm, zsh2, avm_k, avt_k )    ! TKE closure scheme for Kz
+      CASE( np_GLS )   ;   CALL zdf_gls( kt, Kbb, Kmm, zsh2, avm_k, avt_k )    ! GLS closure scheme for Kz
+      CASE( np_OSM )   ;   CALL zdf_osm( kt, Kbb, Kmm, Krhs, avm_k, avt_k )    ! OSMOSIS closure scheme for Kz
    !     CASE( np_CST )                                  ! Constant Kz (reset avt, avm to the background value)
    !         ! avt_k and avm_k set one for all at initialisation phase
 !!gm         avt(2:jpim1,2:jpjm1,1:jpkm1) = rn_avt0 * wmask(2:jpim1,2:jpjm1,1:jpkm1)
 !!gm         avm(2:jpim1,2:jpjm1,1:jpkm1) = rn_avm0 * wmask(2:jpim1,2:jpjm1,1:jpkm1)
-         END SELECT
-
-         IF( ln_tile .AND. .NOT. l_istiled ) CALL dom_tile_start( ldhold=.TRUE. )
-      ENDIF
+      END SELECT
       !
       !                          !==  ocean Kz  ==!   (avt, avs, avm)
