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 14901)
+++ NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90	(revision 14909)
@@ -34,4 +34,5 @@
    !! 16/05/2019 (19) Fox-Kemper parametrization of restratification through mixed layer eddies added to revised code.
    !! 23/05/19   (20) Old code where key_osmldpth1` is *not* set removed, together with the key key_osmldpth1
+   !!             4.2  !  2021-05  (S. Mueller)  Efficiency and source-code clarity enhancements, adaptation for tiling
    !!----------------------------------------------------------------------
 
@@ -282,25 +283,33 @@
       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
       !
-      ALLOCATE( av_t_bl(A2D((nn_hls-1))), av_s_bl(A2D((nn_hls-1))), av_u_bl(A2D((nn_hls-1))), av_v_bl(A2D((nn_hls-1))), av_b_bl(A2D((nn_hls-1))), STAT=ierr)
+      ALLOCATE( av_t_bl(A2D((nn_hls-1))), av_s_bl(A2D((nn_hls-1))), av_u_bl(A2D((nn_hls-1))), av_v_bl(A2D((nn_hls-1))),   &
+         &      av_b_bl(A2D((nn_hls-1))), STAT=ierr)
       zdf_osm_alloc = zdf_osm_alloc + ierr
       !
-      ALLOCATE( av_dt_bl(A2D((nn_hls-1))), av_ds_bl(A2D((nn_hls-1))), av_du_bl(A2D((nn_hls-1))), av_dv_bl(A2D((nn_hls-1))), av_db_bl(A2D((nn_hls-1))), STAT=ierr)
+      ALLOCATE( av_dt_bl(A2D((nn_hls-1))), av_ds_bl(A2D((nn_hls-1))), av_du_bl(A2D((nn_hls-1))), av_dv_bl(A2D((nn_hls-1))),   &
+         &      av_db_bl(A2D((nn_hls-1))), STAT=ierr)
       zdf_osm_alloc = zdf_osm_alloc + ierr
       !
-      ALLOCATE( av_t_ml(A2D((nn_hls-1))), av_s_ml(A2D((nn_hls-1))), av_u_ml(A2D((nn_hls-1))), av_v_ml(A2D((nn_hls-1))), av_b_ml(A2D((nn_hls-1))), STAT=ierr)
+      ALLOCATE( av_t_ml(A2D((nn_hls-1))), av_s_ml(A2D((nn_hls-1))), av_u_ml(A2D((nn_hls-1))), av_v_ml(A2D((nn_hls-1))),   &
+         &      av_b_ml(A2D((nn_hls-1))), STAT=ierr)
       zdf_osm_alloc = zdf_osm_alloc + ierr
       !
-      ALLOCATE( av_dt_ml(A2D((nn_hls-1))), av_ds_ml(A2D((nn_hls-1))), av_du_ml(A2D((nn_hls-1))), av_dv_ml(A2D((nn_hls-1))), av_db_ml(A2D((nn_hls-1))), STAT=ierr)
+      ALLOCATE( av_dt_ml(A2D((nn_hls-1))), av_ds_ml(A2D((nn_hls-1))), av_du_ml(A2D((nn_hls-1))), av_dv_ml(A2D((nn_hls-1))),   &
+         &      av_db_ml(A2D((nn_hls-1))), STAT=ierr)
       zdf_osm_alloc = zdf_osm_alloc + ierr
       !
-      ALLOCATE( av_t_mle(A2D((nn_hls-1))), av_s_mle(A2D((nn_hls-1))), av_u_mle(A2D((nn_hls-1))), av_v_mle(A2D((nn_hls-1))), av_b_mle(A2D((nn_hls-1))), STAT=ierr)
+      ALLOCATE( av_t_mle(A2D((nn_hls-1))), av_s_mle(A2D((nn_hls-1))), av_u_mle(A2D((nn_hls-1))), av_v_mle(A2D((nn_hls-1))),   &
+         &      av_b_mle(A2D((nn_hls-1))), STAT=ierr)
       zdf_osm_alloc = zdf_osm_alloc + ierr
       !
@@ -364,6 +373,6 @@
       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)                              ::   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
@@ -380,11 +389,12 @@
       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, zdbdz_bl_ext   ! External temperature/salinity and 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))) ::   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
@@ -400,12 +410,12 @@
       !
       ! Temporary variables
-      REAL(wp)                        ::   znd   ! Temporary non-dimensional depth
-      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)                        ::   zabsstke
-      REAL(wp)                        ::   zsqrtpi, z_two_thirds, zthickness
-      REAL(wp)                        ::   z2k_times_thickness, zsqrt_depth, zexp_depth, zf, zexperfc
+      REAL(wp)                                 ::   znd              ! Temporary non-dimensional depth
+      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)                                 ::   zabsstke
+      REAL(wp)                                 ::   zsqrtpi, z_two_thirds, zthickness
+      REAL(wp)                                 ::   z2k_times_thickness, zsqrt_depth, zexp_depth, zf, zexperfc
       !
       ! For debugging
@@ -421,5 +431,7 @@
       ! Mixed layer
       ! No initialization of zhbl or zhml (or zdh?)
-      zhbl(:,:)     = pp_large ; zhml(:,:)     = pp_large ; zdh(:,:)      = pp_large
+      zhbl(:,:) = pp_large
+      zhml(:,:) = pp_large
+      zdh(:,:)  = pp_large
       !
       IF ( ln_osm_mle ) THEN   ! Only initialise arrays if needed
@@ -521,5 +533,5 @@
             IF ( hsw(ji,jj) > 1e-4_wp ) THEN
                ! Use  wave fields
-               zabsstke = SQRT( ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2 )
+               zabsstke       = SQRT( ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2 )
                sustke(ji,jj)  = MAX( ( scos_wind(ji,jj) * ut0sd(ji,jj) + ssin_wind(ji,jj)  * vt0sd(ji,jj) ), 1e-8_wp )
                dstokes(ji,jj) = MAX( zfac * hsw(ji,jj) * hsw(ji,jj) / ( MAX( zabsstke * wmp(ji,jj), 1e-7 ) ), 5e-1_wp )
@@ -578,5 +590,6 @@
                zexperfc    = zsqrtpi * zsqrt_depth * ERFC(zsqrt_depth) * EXP( z2k_times_thickness )
             ELSE
-               ! Asymptotic expansion of sqrt(pi)*zsqrt_depth*EXP(z2k_times_thickness)*ERFC(zsqrt_depth) for large z2k_times_thickness
+               ! Asymptotic expansion of sqrt(pi)*zsqrt_depth*EXP(z2k_times_thickness)*ERFC(zsqrt_depth) for large
+               !    z2k_times_thickness
                ! See Abramowitz and Stegun, Eq. 7.1.23
                ! zexperfc = 1._wp - (1/2)/(z2k_times_thickness) + (3/4)/(z2k_times_thickness**2) - (15/8)/(z2k_times_thickness**3)
@@ -599,6 +612,6 @@
          IF ( sla(ji,jj) > 0.45_wp ) dstokes(ji,jj) = MIN( dstokes(ji,jj), 0.5_wp * hbl(ji,jj) )
          ! Velocity scale that tends to sustar for large Langmuir numbers
-         svstr(ji,jj) = ( swstrl(ji,jj)**3 + ( 1.0_wp - EXP( -0.5_wp * sla(ji,jj)**2 ) ) * sustar(ji,jj) * sustar(ji,jj) *   &
-            &                                sustar(ji,jj) )**pthird
+         svstr(ji,jj)  = ( swstrl(ji,jj)**3 + ( 1.0_wp - EXP( -0.5_wp * sla(ji,jj)**2 ) ) * sustar(ji,jj) * sustar(ji,jj) *   &
+            &                                 sustar(ji,jj) )**pthird
          !
          ! Limit maximum value of Langmuir number as approximate treatment for shear turbulence
@@ -677,5 +690,5 @@
          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)
+            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,                         &
@@ -984,9 +997,14 @@
          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( "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
@@ -1008,5 +1026,6 @@
       ! 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 ( .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 )
@@ -1014,8 +1033,8 @@
             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)
+                  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)
@@ -1027,10 +1046,10 @@
             &                    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
+         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
       !
@@ -1051,5 +1070,5 @@
       !!              knlev+kp_ext
       !!----------------------------------------------------------------------
-      INTEGER,                     INTENT(in   )           ::   Kbb, Kmm   ! Ocean time-level indices
+      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(A2D((nn_hls-1))), INTENT(  out)           ::   pt, ps     ! Average temperature and salinity
@@ -1057,12 +1076,14 @@
       REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(  out)           ::   pu, pv     ! Average current components
       INTEGER,  DIMENSION(A2D((nn_hls-1))), INTENT(in   ), OPTIONAL ::   kp_ext     ! External-level offsets
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(  out), OPTIONAL ::   pdt, pds   ! Difference between average temperature, salinity,
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(  out), OPTIONAL ::   pdb        ! buoyancy,
-      REAL(wp), DIMENSION(A2D((nn_hls-1))), 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))), INTENT(  out), OPTIONAL ::   pdt        ! Difference between average temperature,
+      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(  out), OPTIONAL ::   pds        !    salinity,
+      REAL(wp), DIMENSION(A2D((nn_hls-1))), INTENT(  out), OPTIONAL ::   pdb        !    buoyancy, and
+      REAL(wp), DIMENSION(A2D((nn_hls-1))), 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)                    ::   zthermal, zbeta            ! Thermal/haline expansion/contraction coefficients
+      REAL(wp)                             ::   zthermal                   ! Thermal expansion coefficient
+      REAL(wp)                             ::   zbeta                      ! Haline contraction coefficient
       !!----------------------------------------------------------------------
       !
@@ -1154,6 +1175,6 @@
       !!
       !!----------------------------------------------------------------------      
-      REAL(wp),           INTENT(inout), DIMENSION(A2D((nn_hls-1))) ::   pu, pv           ! Components of current
-      LOGICAL,  OPTIONAL, INTENT(in   )                    ::   fwd              ! Forward (default) or reverse rotation
+      REAL(wp),           INTENT(inout), DIMENSION(A2D((nn_hls-1))) ::   pu, pv   ! Components of current
+      LOGICAL,  OPTIONAL, INTENT(in   )                             ::   fwd      ! Forward (default) or reverse rotation
       !
       INTEGER  ::   ji, jj       ! Loop indices
@@ -1184,8 +1205,8 @@
       !!
       !!----------------------------------------------------------------------      
-      REAL(wp),           INTENT(inout), DIMENSION(A2D((nn_hls-1)),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(A2D((nn_hls-1)),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
@@ -1230,5 +1251,5 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER,                     INTENT(in   ) ::   Kmm       ! Ocean time-level index
+      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
@@ -1243,6 +1264,6 @@
       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)                             ::   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
@@ -1278,10 +1299,11 @@
          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 )
          IF ( l_conv(ji,jj) ) THEN
             IF ( av_db_bl(ji,jj) > 0.0_wp ) THEN
-               zri_p(ji,jj) = MAX (  SQRT( av_db_bl(ji,jj) * pdh(ji,jj) / MAX( av_du_bl(ji,jj)**2 + av_dv_bl(ji,jj)**2,       &
+               zri_p(ji,jj) = MAX (  SQRT( av_db_bl(ji,jj) * pdh(ji,jj) / MAX( av_du_bl(ji,jj)**2 + av_dv_bl(ji,jj)**2,     &
                   &                                                          1e-8_wp ) ) * ( phbl(ji,jj) / pdh(ji,jj) ) *   &
                   &                  ( svstr(ji,jj) / MAX( sustar(ji,jj), 1e-6_wp ) )**2 /                                  &
@@ -1305,5 +1327,6 @@
                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                IF ( pshear(ji,jj) > 1e-10 ) THEN
-                  IF ( zri_p(ji,jj) < pp_ri_p_thresh .AND. MIN( hu(ji,jj,Kmm), hu(ji-1,jj,Kmm), hv(ji,jj,Kmm), hv(ji,jj-1,Kmm) ) > 100.0_wp ) THEN
+                  IF ( zri_p(ji,jj) < pp_ri_p_thresh .AND.   &
+                     & MIN( hu(ji,jj,Kmm), hu(ji-1,jj,Kmm), hv(ji,jj,Kmm), hv(ji,jj-1,Kmm) ) > 100.0_wp ) THEN
                      ! Growing shear layer
                      n_ddh(ji,jj) = 0
@@ -1361,10 +1384,10 @@
                   ! Developing shear layer, additional shear production possible.
 
-                  !              pshear_u = MAX( zustar(ji,jj)**2 * MAX( av_du_ml(ji,jj), 0._wp ) /  phbl(ji,jj), 0._wp )
-                  !              pshear(ji,jj) = pshear(ji,jj) + pshear_u * ( 1.0 - MIN( zri_p(ji,jj) / pp_ri_p_thresh, 1.d0 )**2 )
-                  !              pshear(ji,jj) = MIN( pshear(ji,jj), pshear_u )
+                  !    pshear_u = MAX( zustar(ji,jj)**2 * MAX( av_du_ml(ji,jj), 0._wp ) /  phbl(ji,jj), 0._wp )
+                  !    pshear(ji,jj) = pshear(ji,jj) + pshear_u * ( 1.0 - MIN( zri_p(ji,jj) / pp_ri_p_thresh, 1.d0 )**2 )
+                  !    pshear(ji,jj) = MIN( pshear(ji,jj), pshear_u )
 
-                  !              zwb_shr = zwb_shr - 0.25 * MAX ( pshear_u, 0._wp) * ( 1.0 - MIN( zri_p(ji,jj) / pp_ri_p_thresh, 1._wp )**2 )
-                  !              zwb_shr = MAX( zwb_shr, -0.25 * pshear_u )
+                  !    zwb_shr = zwb_shr - 0.25 * MAX ( pshear_u, 0._wp) * ( 1.0 - MIN( zri_p(ji,jj) / pp_ri_p_thresh, 1._wp )**2 )
+                  !    zwb_shr = MAX( zwb_shr, -0.25 * pshear_u )
                ENDIF
                pwb_ent(ji,jj) = pwb_ent(ji,jj) + zwb_shr
@@ -1391,7 +1414,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, pdbdz   ! External gradients of temperature, salinity 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
       !
       ! Local variables
@@ -1476,11 +1500,14 @@
                      !                                                                 !    entrainment > restratification
                      IF ( av_db_bl(ji,jj) > 1e-15_wp ) THEN
-                        zgamma_b_nd = MAX( pdbdz_bl_ext(ji,jj), 0.0_wp ) * pdh(ji,jj) / ( zvel_max + MAX( av_db_bl(ji,jj), 1e-15_wp ) )
-                        zpsi = ( 1.0_wp - 0.5_wp * pdh(ji,jj) / phbl(ji,jj) ) *   &
-                           &   ( swb0(ji,jj) - MIN( ( pwb_min(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ), 0.0_wp ) ) * pdh(ji,jj) / phbl(ji,jj)
+                        zgamma_b_nd = MAX( pdbdz_bl_ext(ji,jj), 0.0_wp ) * pdh(ji,jj) /   &
+                           &          ( zvel_max + MAX( av_db_bl(ji,jj), 1e-15_wp ) )
+                        zpsi = ( 1.0_wp - 0.5_wp * pdh(ji,jj) / phbl(ji,jj) ) *                                                &
+                           &   ( swb0(ji,jj) - MIN( ( pwb_min(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ), 0.0_wp ) ) * pdh(ji,jj) /   &
+                           &   phbl(ji,jj)
                         zpsi = zpsi + 1.75_wp * ( 1.0_wp - 0.5_wp * pdh(ji,jj) / phbl(ji,jj) ) *   &
-                           &   ( pdh(ji,jj) / phbl(ji,jj) + zgamma_b_nd ) * MIN( ( pwb_min(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ), 0.0_wp )
+                           &          ( pdh(ji,jj) / phbl(ji,jj) + zgamma_b_nd ) *   &
+                           &          MIN( ( pwb_min(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ), 0.0_wp )
                         zpsi = pp_alpha_b * MAX( zpsi, 0.0_wp )
-                        pdhdt(ji,jj) = -1.0_wp * ( pwb_ent(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ) /   &
+                        pdhdt(ji,jj) = -1.0_wp * ( pwb_ent(ji,jj) + 2.0_wp * pwb_fk_b(ji,jj) ) /      &
                            &                      ( zvel_max + MAX( av_db_bl(ji,jj), 1e-15_wp ) ) +   &
                            &            zpsi / ( zvel_max + MAX( av_db_bl(ji,jj), 1e-15_wp ) )
@@ -1488,10 +1515,10 @@
                            IF ( ( swstrc(ji,jj) / svstr(ji,jj) )**3 <= 0.5_wp ) THEN
                               zari = MIN( 1.5_wp * av_db_bl(ji,jj) /                                                   &
-                                 &        ( phbl(ji,jj) * ( MAX( pdbdz_bl_ext(ji,jj), 0.0_wp ) +                     &
+                                 &        ( phbl(ji,jj) * ( MAX( pdbdz_bl_ext(ji,jj), 0.0_wp ) +                       &
                                  &                               av_db_bl(ji,jj)**2 / MAX( 4.5_wp * svstr(ji,jj)**2,   &
                                  &                                                       1e-12_wp ) ) ), 0.2_wp )
                            ELSE
                               zari = MIN( 1.5_wp * av_db_bl(ji,jj) /                                                    &
-                                 &        ( phbl(ji,jj) * ( MAX( pdbdz_bl_ext(ji,jj), 0.0_wp ) +                      &
+                                 &        ( phbl(ji,jj) * ( MAX( pdbdz_bl_ext(ji,jj), 0.0_wp ) +                        &
                                  &                               av_db_bl(ji,jj)**2 / MAX( 4.5_wp * swstrc(ji,jj)**2,   &
                                  &                                                       1e-12_wp ) ) ), 0.2_wp )
@@ -1601,5 +1628,5 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER,                     INTENT(in   ) ::   Kmm        ! Ocean time-level index
+      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
@@ -1701,5 +1728,5 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER,                     INTENT(in   ) ::   Kmm            ! Ocean time-level index
+      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
@@ -1789,7 +1816,7 @@
                            ztmp = swstrc(ji,jj)
                         END IF
-                        zari = MIN( 1.5_wp * av_db_bl(ji,jj) / ( phbl(ji,jj) *                           &
-                           &                                   ( MAX( pdbdz_bl_ext(ji,jj), 0.0_wp) +   &
-                           &                                     av_db_bl(ji,jj)**2 / MAX(4.5_wp * ztmp**2 , 1e-12_wp ) ) ), 0.2_wp )
+                        zari = MIN( 1.5_wp * av_db_bl(ji,jj) /                               &
+                           &        ( phbl(ji,jj) * ( MAX( pdbdz_bl_ext(ji,jj), 0.0_wp ) +   &
+                           &                          av_db_bl(ji,jj)**2 / MAX( 4.5_wp * ztmp**2 , 1e-12_wp ) ) ), 0.2_wp )
                      ELSE
                         zari = 0.2_wp
@@ -1807,7 +1834,7 @@
                         ztmp = swstrc(ji,jj)
                      END IF
-                     zari    = MIN( 1.5_wp * av_db_bl(ji,jj) / ( phbl(ji,jj) *                            &
-                        &                                      ( MAX( pdbdz_bl_ext(ji,jj), 0.0_wp ) +   &
-                        &                                        av_db_bl(ji,jj)**2 / MAX( 4.5_wp * ztmp**2 , 1e-12_wp ) ) ), 0.2_wp )
+                     zari    = MIN( 1.5_wp * av_db_bl(ji,jj) /                               &
+                        &           ( phbl(ji,jj) * ( MAX( pdbdz_bl_ext(ji,jj), 0.0_wp ) +   &
+                        &                             av_db_bl(ji,jj)**2 / MAX( 4.5_wp * ztmp**2 , 1e-12_wp ) ) ), 0.2_wp )
                   ELSE
                      zari    = 0.2_wp
@@ -1862,5 +1889,5 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER,                          INTENT(in   ) ::   Kmm            ! Ocean time-level index
+      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
@@ -1946,5 +1973,6 @@
                      END IF   ! IF (shol >=0.5)
                   END IF      ! IF (av_db_bl> 0.)
-               END IF         ! IF (pdhdt >= 0) pdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero
+               END IF         ! IF (pdhdt >= 0) pdhdt < 0 not considered since pycnocline profile is zero and profile arrays are
+               !              !    intialized to zero
                !
             END IF            ! IF (l_conv)
@@ -1972,5 +2000,5 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER,                          INTENT(in   ) ::   Kbb, Kmm       ! Ocean time-level indices
+      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
@@ -2168,5 +2196,5 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER,                          INTENT(in   ) ::   Kmm            ! Time-level index
+      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
@@ -2185,29 +2213,31 @@
       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
+      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)                             ::   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))) ::   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)                             ::   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))) ::   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
-      !!----------------------------------------------------------------------
+      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
       !
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -2231,6 +2261,8 @@
       ! ------------------------------------------------------
       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 )
+         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)))
@@ -2266,14 +2298,17 @@
       ! ---------------------------------------------------------------------
       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 ) /   &
+         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_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 ) )
@@ -2300,8 +2335,8 @@
       ! ----------------------------------------------------------------------
       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 )
+         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
@@ -2443,6 +2478,8 @@
          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)))
+            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
@@ -2454,8 +2491,11 @@
             IF ( ( jk > 1 ) .AND. ( jk <= nmld(ji,jj) ) ) THEN
                zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj)
-               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3_wp * zsc_wth_1(ji,jj) * ( -2.0_wp + 2.75_wp * ( ( 1.0_wp + 0.6_wp * zznd_ml**4 ) - EXP( -6.0_wp * zznd_ml ) ) ) *   &
-                  &      ( 1.0_wp - EXP( -15.0_wp * ( 1.0_wp - zznd_ml ) ) )
-               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3_wp * zsc_ws_1(ji,jj)  * ( -2.0_wp + 2.75_wp * ( ( 1.0_wp + 0.6_wp * zznd_ml**4 ) - EXP( -6.0_wp * zznd_ml ) ) ) *   &
-                  &      ( 1.0_wp - EXP( -15.0_wp * ( 1.0_wp - zznd_ml ) ) )
+               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3_wp * zsc_wth_1(ji,jj) *                                  &
+                  &                                ( -2.0_wp + 2.75_wp * ( ( 1.0_wp + 0.6_wp * zznd_ml**4 ) -   &
+                  &                                                        EXP( -6.0_wp * zznd_ml ) ) ) *       &
+                  &                                ( 1.0_wp - EXP( -15.0_wp * ( 1.0_wp - zznd_ml ) ) )
+               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3_wp * zsc_ws_1(ji,jj) *                                   &
+                  &                                ( -2.0_wp + 2.75_wp * ( ( 1.0_wp + 0.6_wp * zznd_ml**4 ) -   &
+                  &                                EXP( -6.0_wp * zznd_ml ) ) ) * ( 1.0_wp - EXP( -15.0_wp * ( 1.0_wp - zznd_ml ) ) )
             END IF
             !
@@ -2463,6 +2503,8 @@
             IF ( l_pyc(ji,jj) .AND. ( jk >= nmld(ji,jj) ) .AND. ( jk <= nbld(ji,jj) ) ) THEN   ! Pycnocline
                zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj)
-               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 4.0_wp * zsc_wth_pyc(ji,jj) * ( 0.48_wp - EXP( -1.5_wp * ( zznd_pyc - 0.3_wp )**2 ) )
-               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 4.0_wp * zsc_ws_pyc(ji,jj)  * ( 0.48_wp - EXP( -1.5_wp * ( zznd_pyc - 0.3_wp )**2 ) )
+               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 4.0_wp * zsc_wth_pyc(ji,jj) *   &
+                  &                                ( 0.48_wp - EXP( -1.5_wp * ( zznd_pyc - 0.3_wp )**2 ) )
+               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 4.0_wp * zsc_ws_pyc(ji,jj)  *   &
+                  &                                ( 0.48_wp - EXP( -1.5_wp * ( zznd_pyc - 0.3_wp )**2 ) )
             END IF
          ELSE
@@ -2517,5 +2559,6 @@
                ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0.3_wp * 0.15_wp * SIN( 3.14159_wp * ( 0.65_wp * zznd_d ) ) *   &
                   &                                EXP( -0.25_wp * zznd_d**2 ) * zsc_vw_1(ji,jj)
-               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0.3_wp * 0.15_wp * EXP( -5.0 * ( 1.0 - znd ) ) * ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj)
+               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0.3_wp * 0.15_wp * EXP( -5.0 * ( 1.0 - znd ) ) *   &
+                  &                                ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj)
             END IF
          END IF
@@ -2617,5 +2660,6 @@
                      ENDIF   ! IF (shol >=0.5)
                   ENDIF      ! IF (av_db_bl> 0.)
-               ENDIF         ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero
+               ENDIF         ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are
+               !             !    intialized to zero
             END IF
          END IF
@@ -2689,5 +2733,5 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER,                      INTENT(in   ) ::   Kmm          ! Time-level index
+      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
@@ -2698,9 +2742,9 @@
       !
       ! 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
-      REAL(wp)                         ::   zN2_c        ! Local buoyancy difference from 10m value
+      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
+      REAL(wp)                              ::   zN2_c        ! Local buoyancy difference from 10m value
       REAL(wp), DIMENSION(A2D(nn_hls))      ::   ztm
       REAL(wp), DIMENSION(A2D(nn_hls))      ::   zsm
@@ -2735,5 +2779,6 @@
       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, jk_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)
@@ -2796,5 +2841,5 @@
       !!----------------------------------------------------------------------      
       ! Outputs
-      INTEGER,                     INTENT(in   ) ::   Kmm         ! Time-level index
+      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
@@ -2804,12 +2849,12 @@
       !
       ! Local variables
-      INTEGER                     ::   ji, jj, jk        ! Dummy loop indices
+      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
+      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
       !
       znd_param(:,:) = 0.0_wp
@@ -2908,5 +2953,5 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER,                     INTENT(in   ) ::   Kmm         ! Time-level index
+      INTEGER,                              INTENT(in   ) ::   Kmm         ! Time-level index
       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
