Index: NEMO/branches/2020/dev_r13787_OSMOSIS_IMMERSE/src/OCE/DIA/diawri.F90
===================================================================
--- NEMO/branches/2020/dev_r13787_OSMOSIS_IMMERSE/src/OCE/DIA/diawri.F90	(revision 13860)
+++ NEMO/branches/2020/dev_r13787_OSMOSIS_IMMERSE/src/OCE/DIA/diawri.F90	(revision 13867)
@@ -46,4 +46,5 @@
    USE zdfdrg         ! ocean vertical physics: top/bottom friction
    USE zdfmxl         ! mixed layer
+   USE zdfosm         ! mixed layer
    !
    USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
@@ -1076,4 +1077,17 @@
          CALL iom_rstput ( 0, 0, inum, "qz1_abl",  tq_abl(:,:,2,nt_a,2) )   ! now first level humidity
       ENDIF
+      IF( ln_zdfosm ) THEN
+         CALL iom_rstput( 0, 0, inum, 'hbl', hbl*tmask(:,:,1)  )      ! now boundary-layer depth
+         CALL iom_rstput( 0, 0, inum, 'hml', hml*tmask(:,:,1)  )      ! now mixed-layer depth
+         CALL iom_rstput( 0, 0, inum, 'avt_k', avt_k*wmask     )      ! w-level diffusion
+         CALL iom_rstput( 0, 0, inum, 'avm_k', avm_k*wmask     )      ! now w-level viscosity
+         CALL iom_rstput( 0, 0, inum, 'ghamt', ghamt*wmask     )      ! non-local t forcing
+         CALL iom_rstput( 0, 0, inum, 'ghams', ghams*wmask     )      ! non-local s forcing
+         CALL iom_rstput( 0, 0, inum, 'ghamu', ghamu*umask     )      ! non-local u forcing
+         CALL iom_rstput( 0, 0, inum, 'ghamv', ghamv*vmask     )      ! non-local v forcing
+         IF( ln_osm_mle ) THEN
+            CALL iom_rstput( 0, 0, inum, 'hmle', hmle*tmask(:,:,1)  ) ! now transition-layer depth
+         END IF
+      ENDIF
       !
       CALL iom_close( inum )
Index: NEMO/branches/2020/dev_r13787_OSMOSIS_IMMERSE/src/OCE/TRA/traatf.F90
===================================================================
--- NEMO/branches/2020/dev_r13787_OSMOSIS_IMMERSE/src/OCE/TRA/traatf.F90	(revision 13860)
+++ NEMO/branches/2020/dev_r13787_OSMOSIS_IMMERSE/src/OCE/TRA/traatf.F90	(revision 13867)
@@ -117,6 +117,6 @@
       IF( l_trdtra )   THEN                    
          ALLOCATE( ztrdt(jpi,jpj,jpk) , ztrds(jpi,jpj,jpk) )
-         ztrdt(:,:,jpk) = 0._wp
-         ztrds(:,:,jpk) = 0._wp
+         ztrdt(:,:,:) = 0._wp
+         ztrds(:,:,:) = 0._wp
          IF( ln_traldf_iso ) THEN              ! diagnose the "pure" Kz diffusive trend 
             CALL trd_tra( kt, Kmm, Kaa, 'TRA', jp_tem, jptra_zdfp, ztrdt )
Index: NEMO/branches/2020/dev_r13787_OSMOSIS_IMMERSE/src/OCE/TRA/tramle.F90
===================================================================
--- NEMO/branches/2020/dev_r13787_OSMOSIS_IMMERSE/src/OCE/TRA/tramle.F90	(revision 13860)
+++ NEMO/branches/2020/dev_r13787_OSMOSIS_IMMERSE/src/OCE/TRA/tramle.F90	(revision 13867)
@@ -20,4 +20,8 @@
    USE lib_mpp        ! MPP library
    USE lbclnk         ! lateral boundary condition / mpp link
+
+   ! where OSMOSIS_OBL is used with integrated FK
+   USE zdf_oce, ONLY : ln_zdfosm
+   USE zdfosm, ONLY  : ln_osm_mle, hmle, dbdx_mle, dbdy_mle, mld_prof
 
    IMPLICIT NONE
@@ -96,90 +100,135 @@
       !!----------------------------------------------------------------------
       !
-      !                                      !==  MLD used for MLE  ==!
-      !                                                ! compute from the 10m density to deal with the diurnal cycle
-      inml_mle(:,:) = mbkt(:,:) + 1                    ! init. to number of ocean w-level (T-level + 1)
-      IF ( nla10 > 0 ) THEN                            ! avoid case where first level is thicker than 10m
-         DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 )        ! from the bottom to nlb10 (10m)
-            IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle )   inml_mle(ji,jj) = jk      ! Mixed layer
+      !
+      IF(ln_osm_mle.and.ln_zdfosm) THEN
+         ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 )                  ! max level of the computation
+         !
+         !
+         SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts
+         CASE ( 0 )                                               != min of the 2 neighbour MLDs
+            DO_2D( 1, 0, 1, 0 )
+               zhu(ji,jj) = MIN( hmle(ji+1,jj), hmle(ji,jj) )
+               zhv(ji,jj) = MIN( hmle(ji,jj+1), hmle(ji,jj) )
+            END_2D
+         CASE ( 1 )                                               != average of the 2 neighbour MLDs
+            DO_2D( 1, 0, 1, 0 )
+               zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) )
+               zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) )
+            END_2D
+         CASE ( 2 )                                               != max of the 2 neighbour MLDs
+            DO_2D( 1, 0, 1, 0 )
+               zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) )
+               zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) )
+            END_2D
+         END SELECT
+         IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation
+            DO_2D( 1, 0, 1, 0 )
+               zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2u(ji,jj)                                            &
+                    &           * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) )   &
+                    &           / (  MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) )   )
+               !
+               zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj)  * e1v(ji,jj)                                            &
+                    &           * dbdy_mle(ji,jj)  * MIN( 111.e3_wp , e2v(ji,jj) )   &
+                    &           / (  MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) )   )
+            END_2D
+            !
+         ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat)
+            DO_2D( 1, 0, 1, 0 )
+               zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2u(ji,jj)               &
+                    &                  * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) )
+               !
+               zpsim_v(ji,jj) = rc_f *   zhv(ji,jj)   * zhv(ji,jj)   * e1v(ji,jj)               &
+                    &                  * dbdy_mle(ji,jj) * MIN( 111.e3_wp , e2v(ji,jj) )
+            END_2D
+         ENDIF
+
+      ELSE !do not use osn_mle
+         !                                      !==  MLD used for MLE  ==!
+         !                                                ! compute from the 10m density to deal with the diurnal cycle
+         inml_mle(:,:) = mbkt(:,:) + 1                    ! init. to number of ocean w-level (T-level + 1)
+         IF ( nla10 > 0 ) THEN                            ! avoid case where first level is thicker than 10m
+           DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 )        ! from the bottom to nlb10 (10m)
+              IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle )   inml_mle(ji,jj) = jk      ! Mixed layer
+           END_3D
+         ENDIF
+         ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 )                  ! max level of the computation
+         !
+         !
+         zmld(:,:) = 0._wp                      !==   Horizontal shape of the MLE  ==!
+         zbm (:,:) = 0._wp
+         zn2 (:,:) = 0._wp
+         DO_3D( 1, 1, 1, 1, 1, ikmax )                    ! MLD and mean buoyancy and N2 over the mixed layer
+            zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points
+            zmld(ji,jj) = zmld(ji,jj) + zc
+            zbm (ji,jj) = zbm (ji,jj) + zc * (rho0 - rhop(ji,jj,jk) ) * r1_rho0
+            zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp
          END_3D
-      ENDIF
-      ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 )                  ! max level of the computation
-      !
-      !
-      zmld(:,:) = 0._wp                      !==   Horizontal shape of the MLE  ==!
-      zbm (:,:) = 0._wp
-      zn2 (:,:) = 0._wp
-      DO_3D( 1, 1, 1, 1, 1, ikmax )                    ! MLD and mean buoyancy and N2 over the mixed layer
-         zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points
-         zmld(ji,jj) = zmld(ji,jj) + zc
-         zbm (ji,jj) = zbm (ji,jj) + zc * (rho0 - rhop(ji,jj,jk) ) * r1_rho0
-         zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp
-      END_3D
-
-      SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts
-      CASE ( 0 )                                               != min of the 2 neighbour MLDs
-         DO_2D( 1, 0, 1, 0 )
-            zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) )
-            zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) )
-         END_2D
-      CASE ( 1 )                                               != average of the 2 neighbour MLDs
-         DO_2D( 1, 0, 1, 0 )
-            zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp
-            zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp
-         END_2D
-      CASE ( 2 )                                               != max of the 2 neighbour MLDs
-         DO_2D( 1, 0, 1, 0 )
-            zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) )
-            zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) )
-         END_2D
-      END SELECT
-      !                                                ! convert density into buoyancy
-      zbm(:,:) = + grav * zbm(:,:) / MAX( e3t(:,:,1,Kmm), zmld(:,:) )
-      !
-      !
-      !                                      !==  Magnitude of the MLE stream function  ==!
-      !
-      !                 di[bm]  Ds
-      ! Psi = Ce  H^2 ---------------- e2u  mu(z)   where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 )
-      !                  e1u   Lf fu                      and the e2u for the "transport"
-      !                                                      (not *e3u as divided by e3u at the end)
-      !
-      IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation
-         DO_2D( 1, 0, 1, 0 )
-            zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2_e1u(ji,jj)                                            &
-               &           * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )   &
-               &           / (  MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) )   )
+   
+         SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts
+         CASE ( 0 )                                               != min of the 2 neighbour MLDs
+            DO_2D( 1, 0, 1, 0 )
+               zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) )
+               zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) )
+            END_2D
+         CASE ( 1 )                                               != average of the 2 neighbour MLDs
+            DO_2D( 1, 0, 1, 0 )
+               zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp
+               zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp
+            END_2D
+         CASE ( 2 )                                               != max of the 2 neighbour MLDs
+            DO_2D( 1, 0, 1, 0 )
+               zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) )
+               zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) )
+            END_2D
+         END SELECT
+         !                                                ! convert density into buoyancy
+         zbm(:,:) = + grav * zbm(:,:) / MAX( e3t(:,:,1,Kmm), zmld(:,:) )
+         !
+         !
+         !                                      !==  Magnitude of the MLE stream function  ==!
+         !
+         !                 di[bm]  Ds
+         ! Psi = Ce  H^2 ---------------- e2u  mu(z)   where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 )
+         !                  e1u   Lf fu                      and the e2u for the "transport"
+         !                                                      (not *e3u as divided by e3u at the end)
+         !
+         IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation
+            DO_2D( 1, 0, 1, 0 )
+               zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2_e1u(ji,jj)                                            &
+                    &           * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )   &
+                    &           / (  MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) )   )
                !
-            zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj)  * e1_e2v(ji,jj)                                            &
-               &           * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) )   &
-               &           / (  MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) )   )
-         END_2D
-         !
-      ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat)
-         DO_2D( 1, 0, 1, 0 )
-            zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2_e1u(ji,jj)               &
-               &                  * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )
+               zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj)  * e1_e2v(ji,jj)                                            &
+                    &           * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) )   &
+                    &           / (  MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) )   )
+            END_2D
+            !
+         ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat)
+            DO_2D( 1, 0, 1, 0 )
+               zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2_e1u(ji,jj)               &
+                    &                  * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )
                !
-            zpsim_v(ji,jj) = rc_f *   zhv(ji,jj)   * zhv(ji,jj)   * e1_e2v(ji,jj)               &
-               &                  * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) )
-         END_2D
-      ENDIF
-      !
-      IF( nn_conv == 1 ) THEN              ! No MLE in case of convection
-         DO_2D( 1, 0, 1, 0 )
-            IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp )   zpsim_u(ji,jj) = 0._wp
-            IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp )   zpsim_v(ji,jj) = 0._wp
-         END_2D
-      ENDIF
-      !
-      !                                      !==  structure function value at uw- and vw-points  ==!
-      DO_2D( 1, 0, 1, 0 )
-         zhu(ji,jj) = 1._wp / zhu(ji,jj)                   ! hu --> 1/hu
-         zhv(ji,jj) = 1._wp / zhv(ji,jj)
-      END_2D
-      !
-      zpsi_uw(:,:,:) = 0._wp
-      zpsi_vw(:,:,:) = 0._wp
-      !
+               zpsim_v(ji,jj) = rc_f *   zhv(ji,jj)   * zhv(ji,jj)   * e1_e2v(ji,jj)               &
+                    &                  * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) )
+            END_2D
+         ENDIF
+         !
+         IF( nn_conv == 1 ) THEN              ! No MLE in case of convection
+            DO_2D( 1, 0, 1, 0 )
+               IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp )   zpsim_u(ji,jj) = 0._wp
+               IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp )   zpsim_v(ji,jj) = 0._wp
+            END_2D
+         ENDIF
+         !
+      ENDIF  ! end of ln_osm_mle conditional
+    !                                      !==  structure function value at uw- and vw-points  ==!
+    DO_2D( 1, 0, 1, 0 )
+       zhu(ji,jj) = 1._wp / MAX(zhu(ji,jj), rsmall)                   ! hu --> 1/hu
+       zhv(ji,jj) = 1._wp / MAX(zhv(ji,jj), rsmall) 
+    END_2D
+    !
+    zpsi_uw(:,:,:) = 0._wp
+    zpsi_vw(:,:,:) = 0._wp
+    !
       DO_3D( 1, 0, 1, 0, 2, ikmax )                ! start from 2 : surface value = 0
          zcuw = 1._wp - ( gdepw(ji+1,jj,jk,Kmm) + gdepw(ji,jj,jk,Kmm) ) * zhu(ji,jj)
@@ -206,19 +255,23 @@
       END DO
 
-      IF( cdtype == 'TRA') THEN              !==  outputs  ==!
-         !
-         zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:)      ! Lf = N H / f
-         CALL iom_put( "Lf_NHpf" , zLf_NH  )    ! Lf = N H / f
-         !
-         ! divide by cross distance to give streamfunction with dimensions m^2/s
-         DO jk = 1, ikmax+1
-            zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk) * r1_e2u(:,:)
-            zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk) * r1_e1v(:,:)
-         END DO
-         CALL iom_put( "psiu_mle", zpsi_uw )    ! i-mle streamfunction
-         CALL iom_put( "psiv_mle", zpsi_vw )    ! j-mle streamfunction
-      ENDIF
-      !
-   END SUBROUTINE tra_mle_trp
+    IF( cdtype == 'TRA') THEN              !==  outputs  ==!
+       !
+       IF (ln_osm_mle.and.ln_zdfosm) THEN
+          zLf_NH(:,:) = SQRT( rb_c * hmle(:,:) ) * r1_ft(:,:)      ! Lf = N H / f
+       ELSE
+          zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:)      ! Lf = N H / f
+       END IF
+       CALL iom_put( "Lf_NHpf" , zLf_NH  )    ! Lf = N H / f
+       !
+       ! divide by cross distance to give streamfunction with dimensions m^2/s
+       DO jk = 1, ikmax+1
+          zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk) * r1_e2u(:,:)
+          zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk) * r1_e1v(:,:)
+       END DO
+       CALL iom_put( "psiu_mle", zpsi_uw )    ! i-mle streamfunction
+       CALL iom_put( "psiv_mle", zpsi_vw )    ! j-mle streamfunction
+    ENDIF
+    !
+  END SUBROUTINE tra_mle_trp
 
 
Index: NEMO/branches/2020/dev_r13787_OSMOSIS_IMMERSE/src/OCE/TRD/trd_oce.F90
===================================================================
--- NEMO/branches/2020/dev_r13787_OSMOSIS_IMMERSE/src/OCE/TRD/trd_oce.F90	(revision 13860)
+++ NEMO/branches/2020/dev_r13787_OSMOSIS_IMMERSE/src/OCE/TRD/trd_oce.F90	(revision 13867)
@@ -33,5 +33,5 @@
 # endif
    !                                                  !!!* Active tracers trends indexes
-   INTEGER, PUBLIC, PARAMETER ::   jptot_tra  = 20     !: Total trend nb: change it when adding/removing one indice below
+   INTEGER, PUBLIC, PARAMETER ::   jptot_tra  = 21     !: Total trend nb: change it when adding/removing one indice below
    !                               ===============     !  
    INTEGER, PUBLIC, PARAMETER ::   jptra_xad  =  1     !: x- horizontal advection
@@ -46,4 +46,5 @@
    INTEGER, PUBLIC, PARAMETER ::   jptra_bbc  = 10     !: Bottom Boundary Condition (geoth. heating) 
    INTEGER, PUBLIC, PARAMETER ::   jptra_bbl  = 11     !: Bottom Boundary Layer (diffusive and/or advective)
+   INTEGER, PUBLIC, PARAMETER ::   jptra_osm  = 21     !: Non-local terms from OSMOSIS OBL model
    INTEGER, PUBLIC, PARAMETER ::   jptra_npc  = 12     !: non-penetrative convection treatment
    INTEGER, PUBLIC, PARAMETER ::   jptra_dmp  = 13     !: internal restoring (damping)
Index: NEMO/branches/2020/dev_r13787_OSMOSIS_IMMERSE/src/OCE/ZDF/zdfosm.F90
===================================================================
--- NEMO/branches/2020/dev_r13787_OSMOSIS_IMMERSE/src/OCE/ZDF/zdfosm.F90	(revision 13860)
+++ NEMO/branches/2020/dev_r13787_OSMOSIS_IMMERSE/src/OCE/ZDF/zdfosm.F90	(revision 13867)
@@ -25,8 +25,13 @@
    !!            (12) Replace zwstrl with zvstr in calculation of eddy viscosity.
    !! 27/09/2017 (13) Calculate Stokes drift and Stokes penetration depth from wave information
-   !!            (14) Bouyancy flux due to entrainment changed to include contribution from shear turbulence (for testing commented out).
+   !!            (14) Buoyancy flux due to entrainment changed to include contribution from shear turbulence.
    !! 28/09/2017 (15) Calculation of Stokes drift moved into separate do-loops to allow for different options for the determining the Stokes drift to be added.
    !!            (16) Calculation of Stokes drift from windspeed for PM spectrum (for testing, commented out)
    !!            (17) Modification to Langmuir velocity scale to include effects due to the Stokes penetration depth (for testing, commented out)
+   !! ??/??/2018 (18) Revision to code structure, selected using key_osmldpth1. Inline code moved into subroutines. Changes to physics made,
+   !!                  (a) Pycnocline temperature and salinity profies changed for unstable layers
+   !!                  (b) The stable OSBL depth parametrization changed.
+   !! 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
    !!----------------------------------------------------------------------
 
@@ -40,4 +45,6 @@
    !!   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
+   !!
+   !! Subroutines in revised code.
    !!----------------------------------------------------------------------
    USE oce            ! ocean dynamics and active tracers
@@ -69,5 +76,7 @@
    PUBLIC   tra_osm       ! routine called by step.F90
    PUBLIC   trc_osm       ! routine called by trcstp.F90
-   PUBLIC   dyn_osm       ! routine called by 'step.F90'
+   PUBLIC   dyn_osm       ! routine called by step.F90
+
+   PUBLIC   ln_osm_mle    ! logical needed by tra_mle_init in tramle.F90
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamu    !: non-local u-momentum flux
@@ -77,14 +86,28 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   etmean   !: averaging operator for avt
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hbl      !: boundary layer depth
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hbli     !: intial boundary layer depth for stable blayer
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dh       ! depth of pycnocline
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hml      ! ML depth
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dstokes  !: penetration depth of the Stokes drift.
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)           ::   r1_ft    ! inverse of the modified Coriolis parameter at t-pts
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hmle     ! Depth of layer affexted by mixed layer eddies in Fox-Kemper parametrization
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dbdx_mle ! zonal buoyancy gradient in ML
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dbdy_mle ! meridional buoyancy gradient in ML
+   INTEGER,  PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   mld_prof ! level of base of MLE layer.
 
    !                      !!** Namelist  namzdf_osm  **
    LOGICAL  ::   ln_use_osm_la      ! Use namelist  rn_osm_la
+
+   LOGICAL  ::   ln_osm_mle           !: flag to activate the Mixed Layer Eddy (MLE) parameterisation
+
    REAL(wp) ::   rn_osm_la          ! Turbulent Langmuir number
    REAL(wp) ::   rn_osm_dstokes     ! Depth scale of Stokes drift
+   REAL(wp) ::   rn_zdfosm_adjust_sd = 1.0 ! factor to reduce Stokes drift by
+   REAL(wp) ::   rn_osm_hblfrac = 0.1! for nn_osm_wave = 3/4 specify fraction in top of hbl
+   LOGICAL  ::   ln_zdfosm_ice_shelter      ! flag to activate ice sheltering
    REAL(wp) ::   rn_osm_hbl0 = 10._wp       ! Initial value of hbl for 1D runs
    INTEGER  ::   nn_ave             ! = 0/1 flag for horizontal average on avt
    INTEGER  ::   nn_osm_wave = 0    ! = 0/1/2 flag for getting stokes drift from La# / PM wind-waves/Inputs into sbcwave
+   INTEGER  ::   nn_osm_SD_reduce   ! = 0/1/2 flag for getting effective stokes drift from surface value
    LOGICAL  ::   ln_dia_osm         ! Use namelist  rn_osm_la
 
@@ -96,6 +119,26 @@
    REAL(wp) ::   rn_difconv = 1._wp     ! diffusivity when unstable below BL  (m2/s)
 
+! OSMOSIS mixed layer eddy parametrization constants
+   INTEGER  ::   nn_osm_mle             ! = 0/1 flag for horizontal average on avt
+   REAL(wp) ::   rn_osm_mle_ce           ! MLE coefficient
+   !                                        ! parameters used in nn_osm_mle = 0 case
+   REAL(wp) ::   rn_osm_mle_lf               ! typical scale of mixed layer front
+   REAL(wp) ::   rn_osm_mle_time             ! time scale for mixing momentum across the mixed layer
+   !                                        ! parameters used in nn_osm_mle = 1 case
+   REAL(wp) ::   rn_osm_mle_lat              ! reference latitude for a 5 km scale of ML front
+   LOGICAL  ::   ln_osm_hmle_limit           ! If true arbitrarily restrict hmle to rn_osm_hmle_limit*zmld
+   REAL(wp) ::   rn_osm_hmle_limit           ! If ln_osm_hmle_limit true arbitrarily restrict hmle to rn_osm_hmle_limit*zmld
+   REAL(wp) ::   rn_osm_mle_rho_c        ! Density criterion for definition of MLD used by FK
+   REAL(wp) ::   r5_21 = 5.e0 / 21.e0   ! factor used in mle streamfunction computation
+   REAL(wp) ::   rb_c                   ! ML buoyancy criteria = g rho_c /rho0 where rho_c is defined in zdfmld
+   REAL(wp) ::   rc_f                   ! MLE coefficient (= rn_ce / (5 km * fo) ) in nn_osm_mle=1 case
+   REAL(wp) ::   rn_osm_mle_thresh          ! Threshold buoyancy for deepening of MLE layer below OSBL base.
+   REAL(wp) ::   rn_osm_bl_thresh          ! Threshold buoyancy for deepening of OSBL base.
+   REAL(wp) ::   rn_osm_mle_tau             ! Adjustment timescale for MLE.
+
+
    !                                    !!! ** General constants  **
-   REAL(wp) ::   epsln   = 1.0e-20_wp   ! a small positive number
+   REAL(wp) ::   epsln   = 1.0e-20_wp   ! a small positive number to ensure no div by zero
+   REAL(wp) ::   depth_tol = 1.0e-6_wp  ! a small-ish positive number to give a hbl slightly shallower than gdepw
    REAL(wp) ::   pthird  = 1._wp/3._wp  ! 1/3
    REAL(wp) ::   p2third = 2._wp/3._wp  ! 2/3
@@ -118,9 +161,14 @@
       !!                 ***  FUNCTION zdf_osm_alloc  ***
       !!----------------------------------------------------------------------
-     ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk), ghams(jpi,jpj,jpk), &
-          &      hbl(jpi,jpj)    ,  hbli(jpi,jpj)    , dstokes(jpi, jpj) ,                     &
-          &   etmean(jpi,jpj,jpk),  STAT= zdf_osm_alloc )
+     ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), &
+          &       hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), &
+          &       etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc )
+
+     ALLOCATE(  hmle(jpi,jpj), r1_ft(jpi,jpj), dbdx_mle(jpi,jpj), dbdy_mle(jpi,jpj), &
+          &       mld_prof(jpi,jpj), STAT= zdf_osm_alloc )
+
+     CALL mpp_sum ( 'zdfosm', zdf_osm_alloc )
      IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays')
-     CALL mpp_sum ( 'zdfosm', zdf_osm_alloc )
+
    END FUNCTION zdf_osm_alloc
 
@@ -166,4 +214,7 @@
       !!
       INTEGER ::   ji, jj, jk                   ! dummy loop indices
+
+      INTEGER ::   jl                   ! dummy loop indices
+
       INTEGER ::   ikbot, jkmax, jkm1, jkp2     !
 
@@ -196,4 +247,12 @@
       REAL(wp), DIMENSION(jpi,jpj) :: zwbav     ! Buoyancy flux - bl average
       REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent   ! Buoyancy entrainment flux
+      REAL(wp), DIMENSION(jpi,jpj) :: zwb_min
+
+
+      REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk_b  ! MLE buoyancy flux averaged over OSBL
+      REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk    ! max MLE buoyancy flux
+      REAL(wp), DIMENSION(jpi,jpj) :: zdiff_mle ! extra MLE vertical diff
+      REAL(wp), DIMENSION(jpi,jpj) :: zvel_mle  ! velocity scale for dhdt with stable ML and FK
+
       REAL(wp), DIMENSION(jpi,jpj) :: zustke    ! Surface Stokes drift
       REAL(wp), DIMENSION(jpi,jpj) :: zla       ! Trubulent Langmuir number
@@ -201,5 +260,9 @@
       REAL(wp), DIMENSION(jpi,jpj) :: zsin_wind ! Sin angle of surface stress
       REAL(wp), DIMENSION(jpi,jpj) :: zhol      ! Stability parameter for boundary layer
-      LOGICAL, DIMENSION(:,:), ALLOCATABLE :: lconv ! unstable/stable bl
+      LOGICAL, DIMENSION(jpi,jpj)  :: lconv     ! unstable/stable bl
+      LOGICAL, DIMENSION(jpi,jpj)  :: lshear    ! Shear layers
+      LOGICAL, DIMENSION(jpi,jpj)  :: lpyc      ! OSBL pycnocline present
+      LOGICAL, DIMENSION(jpi,jpj)  :: lflux     ! surface flux extends below OSBL into MLE layer.
+      LOGICAL, DIMENSION(jpi,jpj)  :: lmle      ! MLE layer increases in hickness.
 
       ! mixed-layer variables
@@ -207,4 +270,6 @@
       INTEGER, DIMENSION(jpi,jpj) :: ibld ! level of boundary layer base
       INTEGER, DIMENSION(jpi,jpj) :: imld ! level of mixed-layer depth (pycnocline top)
+      INTEGER, DIMENSION(jpi,jpj) :: jp_ext, jp_ext_mle ! offset for external level
+      INTEGER, DIMENSION(jpi, jpj) :: j_ddh ! Type of shear layer
 
       REAL(wp) :: ztgrad,zsgrad,zbgrad ! Temporary variables used to calculate pycnocline gradients
@@ -213,12 +278,24 @@
       REAL(wp), DIMENSION(jpi,jpj) :: zhbl  ! bl depth - grid
       REAL(wp), DIMENSION(jpi,jpj) :: zhml  ! ml depth - grid
+
+      REAL(wp), DIMENSION(jpi,jpj) :: zhmle ! MLE depth - grid
+      REAL(wp), DIMENSION(jpi,jpj) :: zmld  ! ML depth on grid
+
       REAL(wp), DIMENSION(jpi,jpj) :: zdh   ! pycnocline depth - grid
       REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency
-      REAL(wp), DIMENSION(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zrh_bl  ! averages over the depth of the blayer
-      REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zrh_ml  ! averages over the depth of the mixed layer
-      REAL(wp), DIMENSION(jpi,jpj) :: zdt_bl,zds_bl,zdu_bl,zdv_bl,zdrh_bl,zdb_bl ! difference between blayer average and parameter at base of blayer
-      REAL(wp), DIMENSION(jpi,jpj) :: zdt_ml,zds_ml,zdu_ml,zdv_ml,zdrh_ml,zdb_ml ! difference between mixed layer average and parameter at base of blayer
-      REAL(wp), DIMENSION(jpi,jpj) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline
-      REAL(wp), DIMENSION(jpi,jpj) :: zuw_bse,zvw_bse  ! momentum fluxes at the top of the pycnocline
+      REAL(wp), DIMENSION(jpi,jpj) :: zddhdt                                    ! correction to dhdt due to internal structure.
+      REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_bl_ext,zdsdz_bl_ext,zdbdz_bl_ext              ! external temperature/salinity and buoyancy gradients
+      REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_mle_ext,zdsdz_mle_ext,zdbdz_mle_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(jpi,jpj) :: zt_bl,zs_bl,zu_bl,zv_bl,zb_bl  ! averages over the depth of the blayer
+      REAL(wp), DIMENSION(jpi,jpj) :: zt_ml,zs_ml,zu_ml,zv_ml,zb_ml  ! averages over the depth of the mixed layer
+      REAL(wp), DIMENSION(jpi,jpj) :: zt_mle,zs_mle,zu_mle,zv_mle,zb_mle  ! averages over the depth of the MLE layer
+      REAL(wp), DIMENSION(jpi,jpj) :: zdt_bl,zds_bl,zdu_bl,zdv_bl,zdb_bl ! difference between blayer average and parameter at base of blayer
+      REAL(wp), DIMENSION(jpi,jpj) :: zdt_ml,zds_ml,zdu_ml,zdv_ml,zdb_ml ! difference between mixed layer average and parameter at base of blayer
+      REAL(wp), DIMENSION(jpi,jpj) :: zdt_mle,zds_mle,zdu_mle,zdv_mle,zdb_mle ! difference between MLE layer average and parameter at base of blayer
+!      REAL(wp), DIMENSION(jpi,jpj) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline
+      REAL(wp) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline
+      REAL(wp) :: zuw_bse,zvw_bse  ! momentum fluxes at the top of the pycnocline
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz_pyc    ! parametrized gradient of temperature in pycnocline
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdsdz_pyc    ! parametrised gradient of salinity in pycnocline
@@ -226,11 +303,13 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz_pyc    ! u-shear across the pycnocline
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdvdz_pyc    ! v-shear across the pycnocline
-
+      REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle    ! Magnitude of horizontal buoyancy gradient.
       ! Flux-gradient relationship variables
+      REAL(wp), DIMENSION(jpi, jpj) :: zshear, zri_i ! Shear production and interfacial richardon number.
 
       REAL(wp) :: zl_c,zl_l,zl_eps  ! Used to calculate turbulence length scale.
 
-      REAL(wp), DIMENSION(jpi,jpj) :: zdifml_sc,zvisml_sc,zdifpyc_sc,zvispyc_sc,zbeta_d_sc,zbeta_v_sc ! Scales for eddy diffusivity/viscosity
+      REAL(wp) :: za_cubic, zb_cubic, zc_cubic, zd_cubic ! coefficients in cubic polynomial specifying diffusivity in pycnocline.  
       REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_1,zsc_ws_1 ! Temporary scales used to calculate scalar non-gradient terms.
+      REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_pyc, zsc_ws_pyc ! Scales for pycnocline transport term/
       REAL(wp), DIMENSION(jpi,jpj) :: zsc_uw_1,zsc_uw_2,zsc_vw_1,zsc_vw_2 ! Temporary scales for non-gradient momentum flux terms.
       REAL(wp), DIMENSION(jpi,jpj) :: zhbl_t ! holds boundary layer depth updated by full timestep
@@ -243,13 +322,23 @@
       ! Temporary variables
       INTEGER :: inhml
-      INTEGER :: i_lconv_alloc
       REAL(wp) :: znd,znd_d,zznd_ml,zznd_pyc,zznd_d ! temporary non-dimensional depths used in various routines
       REAL(wp) :: ztemp, zari, zpert, zzdhdt, zdb   ! temporary variables
       REAL(wp) :: zthick, zz0, zz1 ! temporary variables
       REAL(wp) :: zvel_max, zhbl_s ! temporary variables
-      REAL(wp) :: zfac             ! temporary variable
+      REAL(wp) :: zfac, ztmp       ! temporary variable
       REAL(wp) :: zus_x, zus_y     ! temporary Stokes drift
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zviscos ! viscosity
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity
+      REAL(wp), DIMENSION(jpi,jpj) :: zalpha_pyc
+      REAL(wp), DIMENSION(jpi,jpj) :: ztau_sc_u ! dissipation timescale at baes of WML.
+      REAL(wp) :: zdelta_pyc, zwt_pyc_sc_1, zws_pyc_sc_1, zzeta_pyc
+      REAL(wp) :: zbuoy_pyc_sc, zomega, zvw_max
+      INTEGER :: ibld_ext=0                          ! does not have to be zero for modified scheme
+      REAL(wp) :: zgamma_b_nd, zgamma_b, zdhoh, ztau
+      REAL(wp) :: zzeta_s = 0._wp
+      REAL(wp) :: zzeta_v = 0.46
+      REAL(wp) :: zabsstke
+      REAL(wp) :: zsqrtpi, z_two_thirds, zproportion, ztransp, zthickness
+      REAL(wp) :: z2k_times_thickness, zsqrt_depth, zexp_depth, zdstokes0, zf, zexperfc
 
       ! For debugging
@@ -257,7 +346,4 @@
       !!--------------------------------------------------------------------
       !
-      ALLOCATE( lconv(jpi,jpj),  STAT= i_lconv_alloc )
-      IF( i_lconv_alloc /= 0 )   CALL ctl_warn('zdf_osm: failed to allocate lconv')
-
       ibld(:,:)   = 0     ; imld(:,:)  = 0
       zrad0(:,:)  = 0._wp ; zradh(:,:) = 0._wp ; zradav(:,:)    = 0._wp ; zustar(:,:)    = 0._wp
@@ -267,22 +353,35 @@
       zustke(:,:) = 0._wp ; zla(:,:)   = 0._wp ; zcos_wind(:,:) = 0._wp ; zsin_wind(:,:) = 0._wp
       zhol(:,:)   = 0._wp
-      lconv(:,:)  = .FALSE.
+      lconv(:,:)  = .FALSE.; lpyc(:,:) = .FALSE. ; lflux(:,:) = .FALSE. ;  lmle(:,:) = .FALSE.
       ! mixed layer
       ! no initialization of zhbl or zhml (or zdh?)
       zhbl(:,:)    = 1._wp ; zhml(:,:)    = 1._wp ; zdh(:,:)      = 1._wp ; zdhdt(:,:)   = 0._wp
-      zt_bl(:,:)   = 0._wp ; zs_bl(:,:)   = 0._wp ; zu_bl(:,:)    = 0._wp ; zv_bl(:,:)   = 0._wp
-      zrh_bl(:,:)  = 0._wp ; zt_ml(:,:)   = 0._wp ; zs_ml(:,:)    = 0._wp ; zu_ml(:,:)   = 0._wp
-      zv_ml(:,:)   = 0._wp ; zrh_ml(:,:)  = 0._wp ; zdt_bl(:,:)   = 0._wp ; zds_bl(:,:)  = 0._wp
-      zdu_bl(:,:)  = 0._wp ; zdv_bl(:,:)  = 0._wp ; zdrh_bl(:,:)  = 0._wp ; zdb_bl(:,:)  = 0._wp
+      zt_bl(:,:)   = 0._wp ; zs_bl(:,:)   = 0._wp ; zu_bl(:,:)    = 0._wp
+      zv_bl(:,:)   = 0._wp ; zb_bl(:,:)  = 0._wp
+      zt_ml(:,:)   = 0._wp ; zs_ml(:,:)    = 0._wp ; zu_ml(:,:)   = 0._wp
+      zt_mle(:,:)   = 0._wp ; zs_mle(:,:)    = 0._wp ; zu_mle(:,:)   = 0._wp
+      zb_mle(:,:) = 0._wp
+      zv_ml(:,:)   = 0._wp ; zdt_bl(:,:)   = 0._wp ; zds_bl(:,:)  = 0._wp
+      zdu_bl(:,:)  = 0._wp ; zdv_bl(:,:)  = 0._wp ; zdb_bl(:,:)  = 0._wp
       zdt_ml(:,:)  = 0._wp ; zds_ml(:,:)  = 0._wp ; zdu_ml(:,:)   = 0._wp ; zdv_ml(:,:)  = 0._wp
-      zdrh_ml(:,:) = 0._wp ; zdb_ml(:,:)  = 0._wp ; zwth_ent(:,:) = 0._wp ; zws_ent(:,:) = 0._wp
-      zuw_bse(:,:) = 0._wp ; zvw_bse(:,:) = 0._wp
+      zdb_ml(:,:)  = 0._wp
+      zdt_mle(:,:)  = 0._wp ; zds_mle(:,:)  = 0._wp ; zdu_mle(:,:)   = 0._wp
+      zdv_mle(:,:)  = 0._wp ; zdb_mle(:,:)  = 0._wp
+      zwth_ent = 0._wp ; zws_ent = 0._wp
       !
       zdtdz_pyc(:,:,:) = 0._wp ; zdsdz_pyc(:,:,:) = 0._wp ; zdbdz_pyc(:,:,:) = 0._wp
       zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp
       !
+      zdtdz_bl_ext(:,:) = 0._wp ; zdsdz_bl_ext(:,:) = 0._wp ; zdbdz_bl_ext(:,:) = 0._wp
+
+      IF ( ln_osm_mle ) THEN  ! only initialise arrays if needed
+         zdtdx(:,:) = 0._wp ; zdtdy(:,:) = 0._wp ; zdsdx(:,:) = 0._wp
+         zdsdy(:,:) = 0._wp ; dbdx_mle(:,:) = 0._wp ; dbdy_mle(:,:) = 0._wp
+         zwb_fk(:,:) = 0._wp ; zvel_mle(:,:) = 0._wp; zdiff_mle(:,:) = 0._wp
+         zhmle(:,:) = 0._wp  ; zmld(:,:) = 0._wp
+      ENDIF
+      zwb_fk_b(:,:) = 0._wp   ! must be initialised even with ln_osm_mle=F as used in zdf_osm_calculate_dhdt
+
       ! Flux-Gradient arrays.
-      zdifml_sc(:,:)  = 0._wp ; zvisml_sc(:,:)  = 0._wp ; zdifpyc_sc(:,:) = 0._wp
-      zvispyc_sc(:,:) = 0._wp ; zbeta_d_sc(:,:) = 0._wp ; zbeta_v_sc(:,:) = 0._wp
       zsc_wth_1(:,:)  = 0._wp ; zsc_ws_1(:,:)   = 0._wp ; zsc_uw_1(:,:)   = 0._wp
       zsc_uw_2(:,:)   = 0._wp ; zsc_vw_1(:,:)   = 0._wp ; zsc_vw_2(:,:)   = 0._wp
@@ -292,4 +391,5 @@
       ghams(:,:,:)   = 0._wp ; ghamu(:,:,:)   = 0._wp ; ghamv(:,:,:) = 0._wp
 
+      zddhdt(:,:) = 0._wp
       ! hbl = MAX(hbl,epsln)
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -326,6 +426,6 @@
         zwbav(ji,jj) = grav  * zthermal * zwthav(ji,jj) - grav  * zbeta * zwsav(ji,jj)
         ! Surface upward velocity fluxes
-        zuw0(ji,jj) = -utau(ji,jj) * r1_rho0 * tmask(ji,jj,1)
-        zvw0(ji,jj) = -vtau(ji,jj) * r1_rho0 * tmask(ji,jj,1)
+        zuw0(ji,jj) = - 0.5 * (utau(ji-1,jj) + utau(ji,jj)) * r1_rho0 * tmask(ji,jj,1)
+        zvw0(ji,jj) = - 0.5 * (vtau(ji,jj-1) + vtau(ji,jj)) * r1_rho0 * tmask(ji,jj,1)
         ! Friction velocity (zustar), at T-point : LMD94 eq. 2
         zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * zuw0(ji,jj) + zvw0(ji,jj) * zvw0(ji,jj) ) ), 1.0e-8 )
@@ -340,6 +440,7 @@
            zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2
            zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2
+           ! Linearly
            zustke(ji,jj) = MAX ( SQRT( zus_x*zus_x + zus_y*zus_y), 1.0e-8 )
-           ! dstokes(ji,jj) set to constant value rn_osm_dstokes from namelist in zdf_osm_init
+           dstokes(ji,jj) = rn_osm_dstokes
         END_2D
      ! Assume Pierson-Moskovitz wind-wave spectrum
@@ -347,18 +448,81 @@
         DO_2D( 0, 0, 0, 0 )
            ! Use wind speed wndm included in sbc_oce module
-           zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 )
-           dstokes(ji,jj) = 0.12 * wndm(ji,jj)**2 / grav
+           zustke(ji,jj) =  MAX ( 0.016 * wndm(ji,jj), 1.0e-8 )
+           dstokes(ji,jj) = MAX ( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1)
         END_2D
      ! Use ECMWF wave fields as output from SBCWAVE
      CASE(2)
         zfac =  2.0_wp * rpi / 16.0_wp
+
         DO_2D( 0, 0, 0, 0 )
-           ! The Langmur number from the ECMWF model appears to give La<0.3 for wind-driven seas.
-           !    The coefficient 0.8 gives La=0.3  in this situation.
-           ! It could represent the effects of the spread of wave directions
-           ! around the mean wind. The effect of this adjustment needs to be tested.
-           zustke(ji,jj) = MAX ( 1.0 * ( zcos_wind(ji,jj) * ut0sd(ji,jj ) + zsin_wind(ji,jj)  * vt0sd(ji,jj) ), &
-                &                zustar(ji,jj) / ( 0.45 * 0.45 )                                                  )
-           dstokes(ji,jj) = MAX(zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zustke(ji,jj)*wmp(ji,jj), 1.0e-7 ) ), 5.0e-1) !rn_osm_dstokes !
+           IF (hsw(ji,jj) > 1.e-4) THEN
+              ! Use  wave fields
+              zabsstke = SQRT(ut0sd(ji,jj)**2 + vt0sd(ji,jj)**2)
+              zustke(ji,jj) = MAX ( ( zcos_wind(ji,jj) * ut0sd(ji,jj) + zsin_wind(ji,jj)  * vt0sd(ji,jj) ), 1.0e-8)
+              dstokes(ji,jj) = MAX (zfac * hsw(ji,jj)*hsw(ji,jj) / ( MAX(zabsstke * wmp(ji,jj), 1.0e-7 ) ), 5.0e-1)
+           ELSE
+              ! Assume masking issue (e.g. ice in ECMWF reanalysis but not in model run)
+              ! .. so default to Pierson-Moskowitz
+              zustke(ji,jj) = MAX ( 0.016 * wndm(ji,jj), 1.0e-8 )
+              dstokes(ji,jj) = MAX ( 0.12 * wndm(ji,jj)**2 / grav, 5.e-1)
+           END IF
+        END_2D
+     END SELECT
+
+     IF (ln_zdfosm_ice_shelter) THEN
+        ! Reduce both Stokes drift and its depth scale by ocean fraction to represent sheltering by ice
+        DO_2D( 0, 0, 0, 0 )
+           zustke(ji,jj) = zustke(ji,jj) * (1.0_wp - fr_i(ji,jj))
+           dstokes(ji,jj) = dstokes(ji,jj) * (1.0_wp - fr_i(ji,jj))
+        END_2D
+     END IF
+
+     SELECT CASE (nn_osm_SD_reduce)
+     ! Reduce surface Stokes drift by a constant factor or following Breivik (2016) + van  Roekel (2012) or Grant (2020).
+     CASE(0)
+        ! The Langmur number from the ECMWF model (or from PM)  appears to give La<0.3 for wind-driven seas.
+        !    The coefficient rn_zdfosm_adjust_sd = 0.8 gives La=0.3  in this situation.
+        ! 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
+           zustke(2:jpim1,2:jpjm1) = rn_zdfosm_adjust_sd * zustke(2:jpim1,2:jpjm1)
+        END IF
+     CASE(1)
+        ! van  Roekel (2012): consider average SD over top 10% of boundary layer
+        ! assumes approximate depth profile of SD from Breivik (2016)
+        zsqrtpi = SQRT(rpi)
+        z_two_thirds = 2.0_wp / 3.0_wp
+
+        DO_2D( 0, 0, 0, 0 )
+           zthickness = rn_osm_hblfrac*hbl(ji,jj)
+           z2k_times_thickness =  zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp )
+           zsqrt_depth = SQRT(z2k_times_thickness)
+           zexp_depth  = EXP(-z2k_times_thickness)
+           zustke(ji,jj) = zustke(ji,jj) * (1.0_wp - zexp_depth  &
+                &              - z_two_thirds * ( zsqrtpi*zsqrt_depth*z2k_times_thickness * ERFC(zsqrt_depth) &
+                &              + 1.0_wp - (1.0_wp + z2k_times_thickness)*zexp_depth ) ) / z2k_times_thickness
+
+        END_2D
+     CASE(2)
+        ! Grant (2020): Match to exponential with same SD and d/dz(Sd) at depth 10% of boundary layer
+        ! assumes approximate depth profile of SD from Breivik (2016)
+        zsqrtpi = SQRT(rpi)
+
+        DO_2D( 0, 0, 0, 0 )
+           zthickness = rn_osm_hblfrac*hbl(ji,jj)
+           z2k_times_thickness =  zthickness * 2.0_wp / MAX( ABS( 5.97_wp * dstokes(ji,jj) ), 0.0000001_wp )
+
+           IF(z2k_times_thickness < 50._wp) THEN
+              zsqrt_depth = SQRT(z2k_times_thickness)
+              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
+              ! 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)
+              zexperfc = ((- 1.875_wp/z2k_times_thickness + 0.75_wp)/z2k_times_thickness - 0.5_wp)/z2k_times_thickness + 1.0_wp
+           END IF
+           zf = z2k_times_thickness*(1.0_wp/zexperfc - 1.0_wp)
+           dstokes(ji,jj) = 5.97 * zf * dstokes(ji,jj)
+           zustke(ji,jj) = zustke(ji,jj) * EXP(z2k_times_thickness * ( 1.0_wp / (2. * zf) - 1.0_wp )) * ( 1.0_wp - zexperfc)
         END_2D
      END SELECT
@@ -369,12 +533,6 @@
         ! Langmuir velocity scale (zwstrl), at T-point
         zwstrl(ji,jj) = ( zustar(ji,jj) * zustar(ji,jj) * zustke(ji,jj) )**pthird
-        ! Modify zwstrl to allow for small and large values of dstokes/hbl.
-        ! Intended as a possible test. Doesn't affect LES results for entrainment,
-        !  but hasn't been shown to be correct as dstokes/h becomes large or small.
-        zwstrl(ji,jj) = zwstrl(ji,jj) *  &
-             & (1.12 * ( 1.0 - ( 1.0 - EXP( -hbl(ji,jj) / dstokes(ji,jj) ) ) * dstokes(ji,jj) / hbl(ji,jj) ))**pthird * &
-             & ( 1.0 - EXP( -15.0 * dstokes(ji,jj) / hbl(ji,jj) ))
-        ! define La this way so effects of Stokes penetration depth on velocity scale are included
-        zla(ji,jj) = SQRT ( zustar(ji,jj) / zwstrl(ji,jj) )**3
+        zla(ji,jj) = MAX(MIN(SQRT ( zustar(ji,jj) / ( zwstrl(ji,jj) + epsln ) )**3, 4.0), 0.2)
+        IF(zla(ji,jj) > 0.45) dstokes(ji,jj) = MIN(dstokes(ji,jj), 0.5_wp*hbl(ji,jj))
         ! Velocity scale that tends to zustar for large Langmuir numbers
         zvstr(ji,jj) = ( zwstrl(ji,jj)**3  + &
@@ -383,5 +541,4 @@
         ! limit maximum value of Langmuir number as approximate treatment for shear turbulence.
         ! Note zustke and zwstrl are not amended.
-        IF ( zla(ji,jj) >= 0.45 ) zla(ji,jj) = 0.45
         !
         ! get convective velocity (zwstrc), stabilty scale (zhol) and logical conection flag lconv
@@ -389,8 +546,6 @@
            zwstrc(ji,jj) = ( 2.0 * zwbav(ji,jj) * 0.9 * hbl(ji,jj) )**pthird
            zhol(ji,jj) = -0.9 * hbl(ji,jj) * 2.0 * zwbav(ji,jj) / (zvstr(ji,jj)**3 + epsln )
-           lconv(ji,jj) = .TRUE.
-        ELSE
+         ELSE
            zhol(ji,jj) = -hbl(ji,jj) *  2.0 * zwbav(ji,jj)/ (zvstr(ji,jj)**3  + epsln )
-           lconv(ji,jj) = .FALSE.
         ENDIF
      END_2D
@@ -399,86 +554,103 @@
      ! Mixed-layer model - calculate averages over the boundary layer, and the change in the boundary layer depth
      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-     ! BL must be always 2 levels deep.
-      hbl(:,:) = MAX(hbl(:,:), gdepw(:,:,3,Kmm) )
-      ibld(:,:) = 3
-      DO_3D( 0, 0, 0, 0, 4, jpkm1 )
+     ! 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 ibld, so must
+     ! previously exist for hbl also.
+
+     ! 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) )
+      ibld(:,:) = 4
+      DO_3D( 1, 1, 1, 1, 5, jpkm1 )
          IF ( hbl(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN
             ibld(ji,jj) = MIN(mbkt(ji,jj), jk)
          ENDIF
       END_3D
+     ! ##########################################################################
 
       DO_2D( 0, 0, 0, 0 )
-            zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
-            zbeta    = rab_n(ji,jj,1,jp_sal)
-            zt   = 0._wp
-            zs   = 0._wp
-            zu   = 0._wp
-            zv   = 0._wp
-            ! average over depth of boundary layer
-            zthick=0._wp
-            DO jm = 2, ibld(ji,jj)
-               zthick=zthick+e3t(ji,jj,jm,Kmm)
-               zt   = zt  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm)
-               zs   = zs  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm)
-               zu   = zu  + e3t(ji,jj,jm,Kmm) &
-                  &            * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) &
-                  &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) )
-               zv   = zv  + e3t(ji,jj,jm,Kmm) &
-                  &            * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) &
-                  &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) )
-            END DO
-            zt_bl(ji,jj) = zt / zthick
-            zs_bl(ji,jj) = zs / zthick
-            zu_bl(ji,jj) = zu / zthick
-            zv_bl(ji,jj) = zv / zthick
-            zdt_bl(ji,jj) = zt_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm)
-            zds_bl(ji,jj) = zs_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm)
-            zdu_bl(ji,jj) = zu_bl(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) &
-                  &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) )
-            zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) &
-                  &   / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) )
-            zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj)
-            IF ( lconv(ji,jj) ) THEN    ! Convective
-                   zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) &
-                        &            + 0.135 * zla(ji,jj) * zwstrl(ji,jj)**3/hbl(ji,jj) )
-
-                   zvel_max =  - ( 1.0 + 1.0 * ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) &
-                        &   * zwb_ent(ji,jj) / ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
-! Entrainment including component due to shear turbulence. Modified Langmuir component, but gives same result for La=0.3 For testing uncomment.
-!                      zwb_ent(ji,jj) = -( 2.0 * 0.2 * zwbav(ji,jj) &
-!                           &            + ( 0.15 * ( 1.0 - EXP( -0.5 * zla(ji,jj) ) ) + 0.03 / zla(ji,jj)**2 ) * zustar(ji,jj)**3/hbl(ji,jj) )
-
-!                      zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / zhbl(ji,jj) ) * zwb_ent(ji,jj) / &
-!                           &       ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
-                   zzdhdt = - zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj),0.0) )
-            ELSE                        ! Stable
-                   zzdhdt = 0.32 * ( hbli(ji,jj) / hbl(ji,jj) -1.0 ) * zwstrl(ji,jj)**3 / hbli(ji,jj) &
-                        &   + ( ( 0.32 / 3.0 ) * exp ( -2.5 * ( hbli(ji,jj) / hbl(ji,jj) - 1.0 ) ) &
-                        & - ( 0.32 / 3.0 - 0.135 * zla(ji,jj) ) * exp ( -12.5 * ( hbli(ji,jj) / hbl(ji,jj) ) ) ) &
-                        &  * zwstrl(ji,jj)**3 / hbli(ji,jj)
-                   zzdhdt = zzdhdt + zwbav(ji,jj)
-                   IF ( zzdhdt < 0._wp ) THEN
-                   ! For long timsteps factor in brackets slows the rapid collapse of the OSBL
-                      zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_Dt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj)
-                   ELSE
-                      zpert   = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_Dt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj) &
-                           &  + MAX( zdb_bl(ji,jj), 0.0 )
-                   ENDIF
-                   zzdhdt = 2.0 * zzdhdt / zpert
-            ENDIF
-            zdhdt(ji,jj) = zzdhdt
+         zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm)
+         imld(ji,jj) = MAX(3,ibld(ji,jj) - MAX( INT( dh(ji,jj) / e3t(ji, jj, ibld(ji,jj), Kmm )) , 1 ))
+         zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm)
+         zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
       END_2D
-
-      ! Calculate averages over depth of boundary layer
-      imld = ibld           ! use imld to hold previous blayer index
-      ibld(:,:) = 3
-
-      zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - ww(ji,jj,ibld(ji,jj)))* rn_Dt ! certainly need wb here, so subtract it
-      zhbl_t(:,:) = MIN(zhbl_t(:,:), ht(:,:))
-      zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_Dt + ww(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom
+      ! Averages over well-mixed and boundary layer
+      jp_ext(:,:) = 2
+      CALL zdf_osm_vertical_average(ibld, jp_ext, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl)
+!      jp_ext(:,:) = ibld(:,:) - imld(:,:) + 1
+      CALL zdf_osm_vertical_average(ibld, ibld-imld+1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml)
+! Velocity components in frame aligned with surface stress.
+      CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml )
+      CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl )
+! Determine the state of the OSBL, stable/unstable, shear/no shear
+      CALL zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear, zri_i )
+
+      IF ( ln_osm_mle ) THEN
+! Fox-Kemper Scheme
+         mld_prof = 4
+         DO_3D( 0, 0, 0, 0, 5, jpkm1 )
+         IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk)
+         END_3D
+         jp_ext_mle(:,:) = 2
+        CALL zdf_osm_vertical_average(mld_prof, jp_ext_mle, zt_mle, zs_mle, zb_mle, zu_mle, zv_mle, zdt_mle, zds_mle, zdb_mle, zdu_mle, zdv_mle)
+
+         DO_2D( 0, 0, 0, 0 )
+           zhmle(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm)
+         END_2D
+
+!! External gradient
+         CALL zdf_osm_external_gradients( ibld+2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext )
+         CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle )
+         CALL zdf_osm_external_gradients( mld_prof, zdtdz_mle_ext, zdsdz_mle_ext, zdbdz_mle_ext )
+         CALL zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk )
+         CALL zdf_osm_mle_parameters( mld_prof, hmle, zhmle, zvel_mle, zdiff_mle )
+      ELSE    ! ln_osm_mle
+! FK not selected, Boundary Layer only.
+         lpyc(:,:) = .TRUE.
+         lflux(:,:) = .FALSE.
+         lmle(:,:) = .FALSE.
+         DO_2D( 0, 0, 0, 0 )
+          IF ( lconv(ji,jj) .AND. zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE.
+         END_2D
+      ENDIF   ! ln_osm_mle
+
+! Test if pycnocline well resolved
+      DO_2D( 0, 0, 0, 0 )
+       IF (lconv(ji,jj) ) THEN
+          ztmp = 0.2 * zhbl(ji,jj) / e3w(ji,jj,ibld(ji,jj),Kmm)
+          IF ( ztmp > 6 ) THEN
+   ! pycnocline well resolved
+            jp_ext(ji,jj) = 1
+          ELSE
+   ! pycnocline poorly resolved
+            jp_ext(ji,jj) = 0
+          ENDIF
+       ELSE
+   ! Stable conditions
+         jp_ext(ji,jj) = 0
+       ENDIF
+      END_2D
+
+      CALL zdf_osm_vertical_average(ibld, jp_ext, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl )
+!      jp_ext = ibld-imld+1
+      CALL zdf_osm_vertical_average(imld-1, ibld-imld+1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml)
+! Rate of change of hbl
+      CALL zdf_osm_calculate_dhdt( zdhdt, zddhdt )
+      DO_2D( 0, 0, 0, 0 )
+       zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - ww(ji,jj,ibld(ji,jj)))* rn_Dt ! certainly need wn here, so subtract it
+            ! adjustment to represent limiting by ocean bottom
+       IF ( zhbl_t(ji,jj) >= gdepw(ji, jj, mbkt(ji,jj) + 1, Kmm ) ) THEN
+          zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw(ji,jj, mbkt(ji,jj) + 1, Kmm) - depth_tol)! ht(:,:))
+          lpyc(ji,jj) = .FALSE.
+       ENDIF
+      END_2D
+
+      imld(:,:) = ibld(:,:)           ! use imld to hold previous blayer index
+      ibld(:,:) = 4
 
       DO_3D( 0, 0, 0, 0, 4, jpkm1 )
          IF ( zhbl_t(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) THEN
-            ibld(ji,jj) =  MIN(mbkt(ji,jj), jk)
+            ibld(ji,jj) = jk
          ENDIF
       END_3D
@@ -487,397 +659,38 @@
 ! Step through model levels taking account of buoyancy change to determine the effect on dhdt
 !
+      CALL zdf_osm_timestep_hbl( zdhdt )
+! is external level in bounds?
+
+      CALL zdf_osm_vertical_average( ibld, jp_ext, zt_bl, zs_bl, zb_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl )
+!
+!
+! Check to see if lpyc needs to be changed 
+
+      CALL zdf_osm_pycnocline_thickness( dh, zdh )
+
       DO_2D( 0, 0, 0, 0 )
-         IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN
+       IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh .or. ibld(ji,jj) + jp_ext(ji,jj) >= mbkt(ji,jj) .or. ibld(ji,jj)-imld(jj,ji) == 1 ) lpyc(ji,jj) = .FALSE. 
+      END_2D
+
+      dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. )  !  Limit delta for shallow boundary layers for calculating flux-gradient terms.
 !
-! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths.
-!
-            zhbl_s = hbl(ji,jj)
-            jm = imld(ji,jj)
-            zthermal = rab_n(ji,jj,1,jp_tem)
-            zbeta = rab_n(ji,jj,1,jp_sal)
-            IF ( lconv(ji,jj) ) THEN
-!unstable
-               zvel_max =  - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) &
-                    &   * zwb_ent(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
-
-               DO jk = imld(ji,jj), ibld(ji,jj)
-                  zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) &
-                       & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0 ) + zvel_max
-
-                  zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_Dt / FLOAT(ibld(ji,jj)-imld(ji,jj) ),   &
-                     &                     e3w(ji,jj,jk,Kmm) )
-                  zhbl_s = MIN(zhbl_s, ht(ji,jj))
-
-                  IF ( zhbl_s >= gdepw(ji,jj,jm+1,Kmm) ) jm = jm + 1
-               END DO
-               hbl(ji,jj) = zhbl_s
-               ibld(ji,jj) = jm
-               hbli(ji,jj) = hbl(ji,jj)
-            ELSE
-! stable
-               DO jk = imld(ji,jj), ibld(ji,jj)
-                  zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) )          &
-                       &               - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), 0.0 ) &
-                       & + 2.0 * zwstrl(ji,jj)**2 / zhbl_s
-
-                  zhbl_s = zhbl_s +  (                                                                                &
-                       &                     0.32         *                         ( hbli(ji,jj) / zhbl_s -1.0 )     &
-                       &               * zwstrl(ji,jj)**3 / hbli(ji,jj)                                               &
-                       &               + ( ( 0.32 / 3.0 )           * EXP( -  2.5 * ( hbli(ji,jj) / zhbl_s -1.0 ) )   &
-                       &               -   ( 0.32 / 3.0  - 0.0485 ) * EXP( - 12.5 * ( hbli(ji,jj) / zhbl_s      ) ) ) &
-                       &          * zwstrl(ji,jj)**3 / hbli(ji,jj) ) / zdb * e3w(ji,jj,jk,Kmm) / zdhdt(ji,jj)  ! ALMG to investigate whether need to include ww here
-
-                  zhbl_s = MIN(zhbl_s, ht(ji,jj))
-                  IF ( zhbl_s >= gdepw(ji,jj,jm,Kmm) ) jm = jm + 1
-               END DO
-               hbl(ji,jj) = MAX(zhbl_s, gdepw(ji,jj,3,Kmm) )
-               ibld(ji,jj) = MAX(jm, 3 )
-               IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj)
-            ENDIF   ! IF ( lconv )
-         ELSE
-! change zero or one model level.
-            hbl(ji,jj) = zhbl_t(ji,jj)
-            IF ( lconv(ji,jj) ) THEN
-               hbli(ji,jj) = hbl(ji,jj)
-            ELSE
-               hbl(ji,jj) = MAX(hbl(ji,jj), gdepw(ji,jj,3,Kmm) )
-               IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj)
-            ENDIF
-         ENDIF
-         zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm)
-      END_2D
-      dstokes(:,:) = MIN ( dstokes(:,:), hbl(:,:)/3. )  !  Limit delta for shallow boundary layers for calculating flux-gradient terms.
-
-! Recalculate averages over boundary layer after depth updated
-     ! Consider later  combining this into the loop above and looking for columns
-     ! where the index for base of the boundary layer have changed
-      DO_2D( 0, 0, 0, 0 )
-            zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
-            zbeta    = rab_n(ji,jj,1,jp_sal)
-            zt   = 0._wp
-            zs   = 0._wp
-            zu   = 0._wp
-            zv   = 0._wp
-            ! average over depth of boundary layer
-            zthick=0._wp
-            DO jm = 2, ibld(ji,jj)
-               zthick=zthick+e3t(ji,jj,jm,Kmm)
-               zt   = zt  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm)
-               zs   = zs  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm)
-               zu   = zu  + e3t(ji,jj,jm,Kmm) &
-                  &            * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) &
-                  &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) )
-               zv   = zv  + e3t(ji,jj,jm,Kmm) &
-                  &            * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) &
-                  &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) )
-            END DO
-            zt_bl(ji,jj) = zt / zthick
-            zs_bl(ji,jj) = zs / zthick
-            zu_bl(ji,jj) = zu / zthick
-            zv_bl(ji,jj) = zv / zthick
-            zdt_bl(ji,jj) = zt_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm)
-            zds_bl(ji,jj) = zs_bl(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm)
-            zdu_bl(ji,jj) = zu_bl(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) &
-                   &   / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) )
-            zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) &
-                   &  / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) )
-            zdb_bl(ji,jj) = grav * zthermal * zdt_bl(ji,jj) - grav * zbeta * zds_bl(ji,jj)
-            zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm)
-            IF ( lconv(ji,jj) ) THEN
-               IF ( zdb_bl(ji,jj) > 0._wp )THEN
-                  IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN  ! near neutral stability
-                        zari = 4.5 * ( zvstr(ji,jj)**2 ) &
-                          & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01
-                  ELSE                                                     ! unstable
-                        zari = 4.5 * ( zwstrc(ji,jj)**2 ) &
-                          & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01
-                  ENDIF
-                  IF ( zari > 0.2 ) THEN                                                ! This test checks for weakly stratified pycnocline
-                     zari = 0.2
-                     zwb_ent(ji,jj) = 0._wp
-                  ENDIF
-                  inhml = MAX( INT( zari * zhbl(ji,jj)   &
-                     &              / e3t(ji,jj,ibld(ji,jj),Kmm) ), 1 )
-                  imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1)
-                  zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm)
-                  zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
-               ELSE  ! IF (zdb_bl)
-                  imld(ji,jj) = ibld(ji,jj) - 1
-                  zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm)
-                  zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
-               ENDIF
-            ELSE   ! IF (lconv)
-               IF ( zdhdt(ji,jj) >= 0.0 ) THEN    ! probably shouldn't include wm here
-               ! boundary layer deepening
-                  IF ( zdb_bl(ji,jj) > 0._wp ) THEN
-               ! pycnocline thickness set by stratification - use same relationship as for neutral conditions.
-                     zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) &
-                       & / ( zdb_bl(ji,jj) * zhbl(ji,jj) ) + 0.01  , 0.2 )
-                     inhml = MAX( INT( zari * zhbl(ji,jj)   &
-                        &             / e3t(ji,jj,ibld(ji,jj),Kmm) ), 1 )
-                     imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1)
-                     zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm)
-                     zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
-                  ELSE
-                     imld(ji,jj) = ibld(ji,jj) - 1
-                     zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm)
-                     zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
-                  ENDIF ! IF (zdb_bl > 0.0)
-               ELSE     ! IF(dhdt >= 0)
-               ! boundary layer collapsing.
-                  imld(ji,jj) = ibld(ji,jj)
-                  zhml(ji,jj) = zhbl(ji,jj)
-                  zdh(ji,jj) = 0._wp
-               ENDIF    ! IF (dhdt >= 0)
-            ENDIF       ! IF (lconv)
-      END_2D
-
-      ! Average over the depth of the mixed layer in the convective boundary layer
-      ! Also calculate entrainment fluxes for temperature and salinity
-      DO_2D( 0, 0, 0, 0 )
-         zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
-         zbeta    = rab_n(ji,jj,1,jp_sal)
-         IF ( lconv(ji,jj) ) THEN
-            zt   = 0._wp
-            zs   = 0._wp
-            zu   = 0._wp
-            zv   = 0._wp
-            ! average over depth of boundary layer
-            zthick=0._wp
-            DO jm = 2, imld(ji,jj)
-               zthick=zthick+e3t(ji,jj,jm,Kmm)
-               zt   = zt  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm)
-               zs   = zs  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm)
-               zu   = zu  + e3t(ji,jj,jm,Kmm) &
-                  &            * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) &
-                  &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) )
-               zv   = zv  + e3t(ji,jj,jm,Kmm) &
-                  &            * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) &
-                  &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) )
-            END DO
-            zt_ml(ji,jj) = zt / zthick
-            zs_ml(ji,jj) = zs / zthick
-            zu_ml(ji,jj) = zu / zthick
-            zv_ml(ji,jj) = zv / zthick
-            zdt_ml(ji,jj) = zt_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm)
-            zds_ml(ji,jj) = zs_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm)
-            zdu_ml(ji,jj) = zu_ml(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) &
-                  &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) )
-            zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) &
-                  &    / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) )
-            zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj)
-         ELSE
-         ! stable, if entraining calulate average below interface layer.
-            IF ( zdhdt(ji,jj) >= 0._wp ) THEN
-               zt   = 0._wp
-               zs   = 0._wp
-               zu   = 0._wp
-               zv   = 0._wp
-               ! average over depth of boundary layer
-               zthick=0._wp
-               DO jm = 2, imld(ji,jj)
-                  zthick=zthick+e3t(ji,jj,jm,Kmm)
-                  zt   = zt  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_tem,Kmm)
-                  zs   = zs  + e3t(ji,jj,jm,Kmm) * ts(ji,jj,jm,jp_sal,Kmm)
-                  zu   = zu  + e3t(ji,jj,jm,Kmm) &
-                     &            * ( uu(ji,jj,jm,Kbb) + uu(ji - 1,jj,jm,Kbb) ) &
-                     &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) )
-                  zv   = zv  + e3t(ji,jj,jm,Kmm) &
-                     &            * ( vv(ji,jj,jm,Kbb) + vv(ji,jj - 1,jm,Kbb) ) &
-                     &            / MAX( 1. , vmask(ji,jj,jm) + vmask(ji,jj - 1,jm) )
-               END DO
-               zt_ml(ji,jj) = zt / zthick
-               zs_ml(ji,jj) = zs / zthick
-               zu_ml(ji,jj) = zu / zthick
-               zv_ml(ji,jj) = zv / zthick
-               zdt_ml(ji,jj) = zt_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_tem,Kmm)
-               zds_ml(ji,jj) = zs_ml(ji,jj) - ts(ji,jj,ibld(ji,jj),jp_sal,Kmm)
-               zdu_ml(ji,jj) = zu_ml(ji,jj) - ( uu(ji,jj,ibld(ji,jj),Kbb) + uu(ji-1,jj,ibld(ji,jj) ,Kbb) ) &
-                     &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) )
-               zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vv(ji,jj,ibld(ji,jj),Kbb) + vv(ji,jj-1,ibld(ji,jj) ,Kbb) ) &
-                     &    / MAX(1. , vmask(ji,jj,ibld(ji,jj) ) + vmask(ji,jj-1,ibld(ji,jj) ) )
-               zdb_ml(ji,jj) = grav * zthermal * zdt_ml(ji,jj) - grav * zbeta * zds_ml(ji,jj)
-            ENDIF
-         ENDIF
-      END_2D
-    !
+    ! Average over the depth of the mixed layer in the convective boundary layer
+!      jp_ext = ibld - imld +1
+      CALL zdf_osm_vertical_average( imld-1, ibld-imld+1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml )
     ! rotate mean currents and changes onto wind align co-ordinates
     !
-
-      DO_2D( 0, 0, 0, 0 )
-         ztemp = zu_ml(ji,jj)
-         zu_ml(ji,jj) = zu_ml(ji,jj) * zcos_wind(ji,jj) + zv_ml(ji,jj) * zsin_wind(ji,jj)
-         zv_ml(ji,jj) = zv_ml(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj)
-         ztemp = zdu_ml(ji,jj)
-         zdu_ml(ji,jj) = zdu_ml(ji,jj) * zcos_wind(ji,jj) + zdv_ml(ji,jj) * zsin_wind(ji,jj)
-         zdv_ml(ji,jj) = zdv_ml(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj)
- !
-         ztemp = zu_bl(ji,jj)
-         zu_bl = zu_bl(ji,jj) * zcos_wind(ji,jj) + zv_bl(ji,jj) * zsin_wind(ji,jj)
-         zv_bl(ji,jj) = zv_bl(ji,jj) * zcos_wind(ji,jj) - ztemp * zsin_wind(ji,jj)
-         ztemp = zdu_bl(ji,jj)
-         zdu_bl(ji,jj) = zdu_bl(ji,jj) * zcos_wind(ji,jj) + zdv_bl(ji,jj) * zsin_wind(ji,jj)
-         zdv_bl(ji,jj) = zdv_bl(ji,jj) * zsin_wind(ji,jj) - ztemp * zsin_wind(ji,jj)
-      END_2D
-
-     zuw_bse = 0._wp
-     zvw_bse = 0._wp
-     DO_2D( 0, 0, 0, 0 )
-
-        IF ( lconv(ji,jj) ) THEN
-           IF ( zdb_bl(ji,jj) > 0._wp ) THEN
-              zwth_ent(ji,jj) = zwb_ent(ji,jj) * zdt_ml(ji,jj) / (zdb_ml(ji,jj) + epsln)
-              zws_ent(ji,jj) = zwb_ent(ji,jj) * zds_ml(ji,jj) / (zdb_ml(ji,jj) + epsln)
-           ENDIF
-        ELSE
-           zwth_ent(ji,jj) = -2.0 * zwthav(ji,jj) * ( (1.0 - 0.8) - ( 1.0 - 0.8)**(3.0/2.0) )
-           zws_ent(ji,jj) = -2.0 * zwsav(ji,jj) * ( (1.0 - 0.8 ) - ( 1.0 - 0.8 )**(3.0/2.0) )
-        ENDIF
-     END_2D
-
+     CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_ml, zv_ml, zdu_ml, zdv_ml )
+     CALL zdf_osm_velocity_rotation( zcos_wind, zsin_wind, zu_bl, zv_bl, zdu_bl, zdv_bl )
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       !  Pycnocline gradients for scalars and velocity
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 
-       DO_2D( 0, 0, 0, 0 )
-       !
-          IF ( lconv (ji,jj) ) THEN
-          ! Unstable conditions
-             IF( zdb_bl(ji,jj) > 0._wp ) THEN
-             ! calculate pycnocline profiles, no need if zdb_bl <= 0. since profile is zero and arrays have been initialized to zero
-                ztgrad = ( zdt_ml(ji,jj) / zdh(ji,jj) )
-                zsgrad = ( zds_ml(ji,jj) / zdh(ji,jj) )
-                zbgrad = ( zdb_ml(ji,jj) / zdh(ji,jj) )
-                DO jk = 2 , ibld(ji,jj)
-                   znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj)
-                   zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
-                   zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
-                   zdsdz_pyc(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
-                END DO
-             ENDIF
-          ELSE
-          ! stable conditions
-          ! if pycnocline profile only defined when depth steady of increasing.
-             IF ( zdhdt(ji,jj) >= 0.0 ) THEN        ! Depth increasing, or steady.
-                IF ( zdb_bl(ji,jj) > 0._wp ) THEN
-                  IF ( zhol(ji,jj) >= 0.5 ) THEN      ! Very stable - 'thick' pycnocline
-                      ztgrad = zdt_bl(ji,jj) / zhbl(ji,jj)
-                      zsgrad = zds_bl(ji,jj) / zhbl(ji,jj)
-                      zbgrad = zdb_bl(ji,jj) / zhbl(ji,jj)
-                      DO jk = 2, ibld(ji,jj)
-                         znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
-                         zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
-                         zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
-                         zdsdz_pyc(ji,jj,jk) =  zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
-                      END DO
-                  ELSE                                   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form.
-                      ztgrad = zdt_bl(ji,jj) / zdh(ji,jj)
-                      zsgrad = zds_bl(ji,jj) / zdh(ji,jj)
-                      zbgrad = zdb_bl(ji,jj) / zdh(ji,jj)
-                      DO jk = 2, ibld(ji,jj)
-                         znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj)
-                         zdtdz_pyc(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
-                         zdbdz_pyc(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
-                         zdsdz_pyc(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
-                      END DO
-                   ENDIF ! IF (zhol >=0.5)
-                ENDIF    ! IF (zdb_bl> 0.)
-             ENDIF       ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero, profile arrays are intialized to zero
-          ENDIF          ! IF (lconv)
-         !
-       END_2D
-!
-       DO_2D( 0, 0, 0, 0 )
-       !
-          IF ( lconv (ji,jj) ) THEN
-          ! Unstable conditions
-              zugrad = ( zdu_ml(ji,jj) / zdh(ji,jj) ) + 0.275 * zustar(ji,jj)*zustar(ji,jj) / &
-            & (( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) / zla(ji,jj)**(8.0/3.0)
-             zvgrad = ( zdv_ml(ji,jj) / zdh(ji,jj) ) + 3.5 * ff_t(ji,jj) * zustke(ji,jj) / &
-           & ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
-             DO jk = 2 , ibld(ji,jj)-1
-                znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj)
-                zdudz_pyc(ji,jj,jk) =  zugrad * EXP( -1.75 * ( znd + 0.75 )**2 )
-                zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
-             END DO
-          ELSE
-          ! stable conditions
-             zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj)
-             zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj)
-             DO jk = 2, ibld(ji,jj)
-                znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
-                IF ( znd < 1.0 ) THEN
-                   zdudz_pyc(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 )
-                ELSE
-                   zdudz_pyc(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 )
-                ENDIF
-                zdvdz_pyc(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 )
-             END DO
-          ENDIF
-         !
-       END_2D
+      CALL zdf_osm_external_gradients( ibld+2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext )
+      CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc, zalpha_pyc )
+      CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc )
        !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
        ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship
        !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-
-      ! WHERE ( lconv )
-      !     zdifml_sc = zhml * ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird
-      !     zvisml_sc = zdifml_sc
-      !     zdifpyc_sc = 0.165 * ( zwstrl**3 + zwstrc**3 )**pthird * ( zhbl - zhml )
-      !     zvispyc_sc = 0.142 * ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * ( zhbl - zhml )
-      !     zbeta_d_sc = 1.0 - (0.165 / 0.8 * ( zhbl - zhml ) / zhbl )**p2third
-      !     zbeta_v_sc = 1.0 -  2.0 * (0.142 /0.375) * (zhbl - zhml ) / zhml
-      !  ELSEWHERE
-      !     zdifml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 )
-      !     zvisml_sc = zwstrl * zhbl * EXP ( -( zhol / 0.183_wp )**2 )
-      !  ENDWHERE
-       DO_2D( 0, 0, 0, 0 )
-          IF ( lconv(ji,jj) ) THEN
-            zdifml_sc(ji,jj) = zhml(ji,jj) * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
-            zvisml_sc(ji,jj) = zdifml_sc(ji,jj)
-            zdifpyc_sc(ji,jj) = 0.165 * ( zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj)
-            zvispyc_sc(ji,jj) = 0.142 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zdh(ji,jj)
-            zbeta_d_sc(ji,jj) = 1.0 - (0.165 / 0.8 * zdh(ji,jj) / zhbl(ji,jj) )**p2third
-            zbeta_v_sc(ji,jj) = 1.0 -  2.0 * (0.142 /0.375) * zdh(ji,jj) / zhml(ji,jj)
-          ELSE
-            zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 )
-            zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * EXP ( -( zhol(ji,jj) / 0.6_wp )**2 )
-         END IF
-       END_2D
-!
-       DO_2D( 0, 0, 0, 0 )
-          IF ( lconv(ji,jj) ) THEN
-             DO jk = 2, imld(ji,jj)   ! mixed layer diffusivity
-                 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj)
-                 !
-                 zdiffut(ji,jj,jk) = 0.8   * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml    )**1.5
-                 !
-                 zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml    ) &
-                      &            *                                      ( 1.0 -               0.5 * zznd_ml**2 )
-             END DO
-             ! pycnocline - if present linear profile
-             IF ( zdh(ji,jj) > 0._wp ) THEN
-                DO jk = imld(ji,jj)+1 , ibld(ji,jj)
-                    zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj)
-                    !
-                    zdiffut(ji,jj,jk) = zdifpyc_sc(ji,jj) * ( 1.0 + zznd_pyc )
-                    !
-                    zviscos(ji,jj,jk) = zvispyc_sc(ji,jj) * ( 1.0 + zznd_pyc )
-                END DO
-             ENDIF
-             ! Temporay fix to ensure zdiffut is +ve; won't be necessary with ww taken out
-             zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj)* e3t(ji,jj,ibld(ji,jj),Kmm)
-             ! could be taken out, take account of entrainment represents as a diffusivity
-             ! should remove w from here, represents entrainment
-          ELSE
-          ! stable conditions
-             DO jk = 2, ibld(ji,jj)
-                zznd_ml = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
-                zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5
-                zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 )
-             END DO
-          ENDIF   ! end if ( lconv )
-!
-       END_2D
+       CALL zdf_osm_diffusivity_viscosity( zdiffut, zviscos )
 
        !
@@ -918,15 +731,17 @@
        END_2D
 
-
 ! Stokes term in flux-gradient relationship (note in zsc_uw_n don't use zvstr since term needs to go to zero as zwstrl goes to zero)
        WHERE ( lconv )
-          zsc_uw_1 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke /( 1.0 - 1.0 * 6.5 * zla**(8.0/3.0) )
-          zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / ( zla**(8.0/3.0) + epsln )
-          zsc_vw_1 = ff_t * zhml * zustke**3 * zla**(8.0/3.0) / ( ( zvstr**3 + 0.5 * zwstrc**3 )**(2.0/3.0) + epsln )
+          zsc_uw_1 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MAX( ( 1.0 - 1.0 * 6.5 * zla**(8.0/3.0) ), 0.2 )
+          zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MIN( zla**(8.0/3.0) + epsln, 0.12 )
+          zsc_vw_1 = ff_t * zhml * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / ( ( zvstr**3 + 0.5 * zwstrc**3 )**(2.0/3.0) + epsln )
        ELSEWHERE
           zsc_uw_1 = zustar**2
-          zsc_vw_1 = ff_t * zhbl * zustke**3 * zla**(8.0/3.0) / (zvstr**2 + epsln)
+          zsc_vw_1 = ff_t * zhbl * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / (zvstr**2 + epsln)
        ENDWHERE
-
+       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 )
+       END IF
        DO_2D( 0, 0, 0, 0 )
           IF ( lconv(ji,jj) ) THEN
@@ -970,10 +785,41 @@
                 zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           &
                      &     * ( 1.0 - EXP ( - 5.0 * (     1.0 - zznd_ml          ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) )
-                zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( 3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0/2.0)
+                zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0 / 2.0)
                 ! non-gradient buoyancy terms
                 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * 0.5 * zsc_wth_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml )
                 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * 0.5 *  zsc_ws_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml )
              END DO
-          ELSE
+             
+             IF ( lpyc(ji,jj) ) THEN
+               ztau_sc_u(ji,jj) = zhml(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird
+               ztau_sc_u(ji,jj) = ztau_sc_u(ji,jj) * ( 1.4 -0.4 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) )**1.5 )
+               zwth_ent =  -0.003 * ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird * ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zdt_ml(ji,jj)                  
+               zws_ent =  -0.003 * ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird * ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zds_ml(ji,jj)
+! Cubic profile used for buoyancy term
+               za_cubic = 0.755 * ztau_sc_u(ji,jj)
+               zb_cubic = 0.25 * ztau_sc_u(ji,jj)
+               DO jk = 2, ibld(ji,jj)
+                 zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj)
+                 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - 0.045 * ( ( zwth_ent * zdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) * MAX( ( 1.75 * zznd_pyc -0.15 * zznd_pyc**2 - 0.2 * zznd_pyc**3 ), 0.0 )
+
+                 ghams(ji,jj,jk) = ghams(ji,jj,jk) - 0.045 * ( ( zws_ent * zdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) * MAX( ( 1.75 * zznd_pyc -0.15 * zznd_pyc**2 - 0.2 * zznd_pyc**3 ), 0.0 )
+               END DO
+!
+               zbuoy_pyc_sc = zalpha_pyc(ji,jj) * zdb_ml(ji,jj) / zdh(ji,jj) + zdbdz_bl_ext(ji,jj)
+               zdelta_pyc = ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird / SQRT( MAX( zbuoy_pyc_sc, ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**p2third / zdh(ji,jj)**2 ) )
+!
+               zwt_pyc_sc_1 = 0.325 * ( zalpha_pyc(ji,jj) * zdt_ml(ji,jj) / zdh(ji,jj) + zdtdz_bl_ext(ji,jj) ) * zdelta_pyc**2 / zdh(ji,jj)
+!
+               zws_pyc_sc_1 = 0.325 * ( zalpha_pyc(ji,jj) * zds_ml(ji,jj) / zdh(ji,jj) + zdsdz_bl_ext(ji,jj) ) * zdelta_pyc**2 / zdh(ji,jj)
+!
+               zzeta_pyc = 0.15 - 0.175 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) ) 
+               DO jk = 2, ibld(ji,jj)
+                 zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj)
+                 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.05 * zwt_pyc_sc_1 * EXP( -0.25 * ( zznd_pyc / zzeta_pyc )**2 ) * zdh(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird
+!
+                 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.05 * zws_pyc_sc_1 * EXP( -0.25 * ( zznd_pyc / zzeta_pyc )**2 ) * zdh(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird
+               END DO
+            ENDIF ! End of pycnocline                  
+          ELSE ! lconv test - stable conditions
              DO jk = 2, ibld(ji,jj)
                 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj)
@@ -982,5 +828,4 @@
           ENDIF
        END_2D
-
 
        WHERE ( lconv )
@@ -1011,13 +856,54 @@
        END_2D
 
+       DO_2D( 0, 0, 0, 0 )
+        IF ( lpyc(ji,jj) ) THEN
+          IF ( j_ddh(ji,jj) == 0 ) THEN
+! Place holding code. Parametrization needs checking for these conditions.
+            zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) )**pthird )**pthird
+            zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj)
+            zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj)
+          ELSE
+            zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) )**pthird )**pthird
+            zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj)
+            zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj)
+          ENDIF
+          zd_cubic = zdh(ji,jj) / zhbl(ji,jj) * zuw0(ji,jj) - ( 2.0 + zdh(ji,jj) /zhml(ji,jj) ) * zuw_bse
+          zc_cubic = zuw_bse - zd_cubic
+! need ztau_sc_u to be available. Change to array. 
+          DO jk = imld(ji,jj), ibld(ji,jj)
+             zznd_pyc = - ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj)
+             ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045 * ztau_sc_u(ji,jj)**2 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk)
+          END DO
+          zvw_max = 0.7 * ff_t(ji,jj) * ( zustke(ji,jj) * dstokes(ji,jj) + 0.75 * zustar(ji,jj) * zhml(ji,jj) )
+          zd_cubic = zvw_max * zdh(ji,jj) / zhml(ji,jj) - ( 2.0 + zdh(ji,jj) /zhml(ji,jj) ) * zvw_bse
+          zc_cubic = zvw_bse - zd_cubic
+          DO jk = imld(ji,jj), ibld(ji,jj)
+            zznd_pyc = -( gdepw(ji,jj,jk,Kmm) -zhbl(ji,jj) ) / zdh(ji,jj)
+            ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.045 * ztau_sc_u(ji,jj)**2 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk)
+          END DO
+        ENDIF  ! lpyc
+       END_2D
+
+       IF(ln_dia_osm) THEN
+          IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu )
+          IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 )
+       END IF
 ! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ]
 
-       WHERE ( lconv )
-          zsc_wth_1 = zwth0
-          zsc_ws_1 = zws0
-       ELSEWHERE
-          zsc_wth_1 = 2.0 * zwthav
-          zsc_ws_1 = zws0
-       ENDWHERE
+       DO_2D( 1, 0, 1, 0 )
+       
+         IF ( lconv(ji,jj) ) THEN
+           zsc_wth_1(ji,jj) = zwth0(ji,jj) / ( 1.0 - 0.56 * EXP( zhol(ji,jj) ) )
+           zsc_ws_1(ji,jj) = zws0(ji,jj) / (1.0 - 0.56 *EXP( zhol(ji,jj) ) )
+           IF ( lpyc(ji,jj) ) THEN
+! Pycnocline scales
+              zsc_wth_pyc(ji,jj) = -0.2 * zwb0(ji,jj) * zdt_bl(ji,jj) / zdb_bl(ji,jj)
+              zsc_ws_pyc(ji,jj) = -0.2 * zwb0(ji,jj) * zds_bl(ji,jj) / zdb_bl(ji,jj)
+            ENDIF
+         ELSE
+           zsc_wth_1(ji,jj) = 2.0 * zwthav(ji,jj)
+           zsc_ws_1(ji,jj) = zws0(ji,jj)
+         ENDIF
+       END_2D
 
        DO_2D( 0, 0, 0, 0 )
@@ -1035,16 +921,26 @@
                     &          * ( 1.0 - EXP ( -15.0 * (         1.0 - zznd_ml    ) ) )
             END DO
+!
+            IF ( lpyc(ji,jj) ) THEN
+! pycnocline
+              DO jk = imld(ji,jj), ibld(ji,jj)
+                zznd_pyc = - ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj)
+                ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 4.0 * zsc_wth_pyc(ji,jj) * ( 0.48 - EXP( -1.5 * ( zznd_pyc -0.3)**2 ) ) 
+                ghams(ji,jj,jk) = ghams(ji,jj,jk) + 4.0 * zsc_ws_pyc(ji,jj) * ( 0.48 - EXP( -1.5 * ( zznd_pyc -0.3)**2 ) ) 
+              END DO
+           ENDIF
          ELSE
-            DO jk = 2, ibld(ji,jj)
-               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
-               znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
-               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + &
-            &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj)
-               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + &
-            &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_ws_1(ji,jj)
-            END DO
+            IF( zdhdt(ji,jj) > 0. ) THEN
+              DO jk = 2, ibld(ji,jj)
+                 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
+                 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
+                 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + &
+              &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj)
+                 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + &
+               &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_ws_1(ji,jj)
+              END DO
+            ENDIF
          ENDIF
        END_2D
-
 
        WHERE ( lconv )
@@ -1090,22 +986,23 @@
           ENDIF
        END_2D
+
+       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") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 )
+          IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 )
+          IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 )
+          IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 )
+       END IF
 !
 ! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer.
 
+
+ ! Make surface forced velocity non-gradient terms go to zero at the base of the boundary layer.
+
       DO_2D( 0, 0, 0, 0 )
-         IF ( lconv(ji,jj) ) THEN
+         IF ( .not. lconv(ji,jj) ) THEN
             DO jk = 2, ibld(ji,jj)
-               znd = ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about
-               IF ( znd >= 0.0 ) THEN
-                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) )
-                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -30.0 * znd**2 ) )
-               ELSE
-                  ghamu(ji,jj,jk) = 0._wp
-                  ghamv(ji,jj,jk) = 0._wp
-               ENDIF
-            END DO
-         ELSE
-            DO jk = 2, ibld(ji,jj)
-               znd = ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zhml(ji,jj) !ALMG to think about
+               znd = ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zhbl(ji,jj) !ALMG to think about
                IF ( znd >= 0.0 ) THEN
                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) )
@@ -1120,43 +1017,36 @@
 
       ! pynocline contributions
-       ! Temporary fix to avoid instabilities when zdb_bl becomes very very small
-       zsc_uw_1 = 0._wp ! 50.0 * zla**(8.0/3.0) * zustar**2 * zhbl / ( zdb_bl + epsln )
        DO_2D( 0, 0, 0, 0 )
-          DO jk= 2, ibld(ji,jj)
-             znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
-             ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk)
-             ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk)
-             ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk)
-             ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj) * ( 1.0 - znd )**(7.0/4.0) * zdbdz_pyc(ji,jj,jk)
-             ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk)
-          END DO
+         IF ( .not. lconv(ji,jj) ) THEN
+          IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
+             DO jk= 2, ibld(ji,jj)
+                znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
+                ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk)
+                ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk)
+                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk)
+                ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk)
+             END DO
+          END IF
+         END IF
        END_2D
-
-! Entrainment contribution.
+      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 )
+       END IF
 
        DO_2D( 0, 0, 0, 0 )
-          IF ( lconv(ji,jj) ) THEN
-            DO jk = 1, imld(ji,jj) - 1
-               znd=gdepw(ji,jj,jk,Kmm) / zhml(ji,jj)
-               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * znd
-               ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * znd
-               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * znd
-               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * znd
-            END DO
-            DO jk = imld(ji,jj), ibld(ji,jj)
-               znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) / zdh(ji,jj)
-               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth_ent(ji,jj) * ( 1.0 + znd )
-               ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws_ent(ji,jj) * ( 1.0 + znd )
-               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zuw_bse(ji,jj) * ( 1.0 + znd )
-               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zvw_bse(ji,jj) * ( 1.0 + znd )
-             END DO
-          ENDIF
-          ghamt(ji,jj,ibld(ji,jj)) = 0._wp
-          ghams(ji,jj,ibld(ji,jj)) = 0._wp
-          ghamu(ji,jj,ibld(ji,jj)) = 0._wp
-          ghamv(ji,jj,ibld(ji,jj)) = 0._wp
+          ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
+          ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
+          ghamu(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
+          ghamv(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
        END_2D
 
-
+       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("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc )
+          IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc )
+          IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos )
+       END IF
        !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
        ! Need to put in code for contributions that are applied explicitly to
@@ -1180,4 +1070,6 @@
        IF(ln_dia_osm) THEN
           IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc )
+          IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc )
+          IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc )
        END IF
 
@@ -1222,6 +1114,45 @@
        END IF ! ln_convmix = .true.
 
+
+
+       IF ( ln_osm_mle ) THEN  ! set up diffusivity and non-gradient mixing
+          DO_2D( 0, 0, 0, 0 )
+              IF ( lflux(ji,jj) ) THEN ! MLE mixing extends below boundary layer
+             ! Calculate MLE flux contribution from surface fluxes
+                DO jk = 1, ibld(ji,jj)
+                  znd = gdepw(ji,jj,jk,Kmm) / MAX(zhbl(ji,jj),epsln)
+                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - zwth0(ji,jj) * ( 1.0 - znd )
+                  ghams(ji,jj,jk) = ghams(ji,jj,jk) - zws0(ji,jj) * ( 1.0 - znd )
+                 END DO
+                 DO jk = 1, mld_prof(ji,jj)
+                   znd = gdepw(ji,jj,jk,Kmm) / MAX(zhmle(ji,jj),epsln)
+                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zwth0(ji,jj) * ( 1.0 - znd )
+                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + zws0(ji,jj) * ( 1.0 -znd )
+                 END DO
+         ! Viscosity for MLEs
+                 DO jk = 1, mld_prof(ji,jj)
+                   znd = -gdepw(ji,jj,jk,Kmm) / MAX(zhmle(ji,jj),epsln)
+                   zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdiff_mle(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + 5.0 / 21.0 * ( 2.0 * znd + 1.0 )** 2 )
+                 END DO
+              ELSE
+! Surface transports limited to OSBL.                 
+         ! Viscosity for MLEs
+                 DO jk = 1, mld_prof(ji,jj)
+                   znd = -gdepw(ji,jj,jk,Kmm) / MAX(zhmle(ji,jj),epsln)
+                   zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdiff_mle(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + 5.0 / 21.0 * ( 2.0 * znd + 1.0 )** 2 )
+                 END DO
+              ENDIF
+          END_2D
+       ENDIF
+
+       IF(ln_dia_osm) THEN
+          IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc )
+          IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc )
+          IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc )
+       END IF
+
+
        ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids
-       CALL lbc_lnk( 'zdfosm', zviscos(:,:,:), 'W', 1.0_wp )
+       !CALL lbc_lnk( zviscos(:,:,:), 'W', 1. )
 
        ! GN 25/8: need to change tmask --> wmask
@@ -1244,4 +1175,6 @@
             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_multi( 'zdfosm', hbl, 'T', 1., dh, 'T', 1., hmle, 'T', 1. )
         ! Lateral boundary conditions on final outputs for gham[ts],  on W-grid  (sign unchanged)
         ! Lateral boundary conditions on final outputs for gham[uv],  on [UV]-grid  (sign unchanged)
@@ -1249,5 +1182,5 @@
          &                  ghamu, 'U', 1.0_wp , ghamv, 'V', 1.0_wp )
 
-       IF(ln_dia_osm) THEN
+      IF(ln_dia_osm) THEN
          SELECT CASE (nn_osm_wave)
          ! Stokes drift set by assumimg onstant La#=0.3(=0)  or Pierson-Moskovitz spectrum (=1).
@@ -1257,7 +1190,12 @@
             IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*zustar**2*zustke )
          ! Stokes drift read in from sbcwave  (=2).
-         CASE(2)
-            IF ( iom_use("us_x") ) CALL iom_put( "us_x", ut0sd )               ! x surface Stokes drift
-            IF ( iom_use("us_y") ) CALL iom_put( "us_y", vt0sd )               ! y surface Stokes drift
+         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.*rpi*1.026/(0.877*grav) )*wndm*tmask(:,:,1) )                  ! wave mean period from NP spectrum
+            IF ( iom_use("hsw_NP") ) CALL iom_put( "hsw_NP", (0.22/grav)*wndm**2*tmask(:,:,1) )                   ! significant wave height from NP spectrum
+            IF ( iom_use("wndm") ) CALL iom_put( "wndm", wndm*tmask(:,:,1) )                   ! U_10
             IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rho0*tmask(:,:,1)*zustar**2* &
                  & SQRT(ut0sd**2 + vt0sd**2 ) )
@@ -1270,5 +1208,12 @@
          IF ( iom_use("zws0") ) CALL iom_put( "zws0", tmask(:,:,1)*zws0 )                ! <Sw_0>
          IF ( iom_use("hbl") ) CALL iom_put( "hbl", tmask(:,:,1)*hbl )                  ! boundary-layer depth
-         IF ( iom_use("hbli") ) CALL iom_put( "hbli", tmask(:,:,1)*hbli )               ! Initial boundary-layer depth
+         IF ( iom_use("ibld") ) CALL iom_put( "ibld", tmask(:,:,1)*ibld )               ! boundary-layer max k
+         IF ( iom_use("zdt_bl") ) CALL iom_put( "zdt_bl", tmask(:,:,1)*zdt_bl )           ! dt at ml base
+         IF ( iom_use("zds_bl") ) CALL iom_put( "zds_bl", tmask(:,:,1)*zds_bl )           ! ds at ml base
+         IF ( iom_use("zdb_bl") ) CALL iom_put( "zdb_bl", tmask(:,:,1)*zdb_bl )           ! db at ml base
+         IF ( iom_use("zdu_bl") ) CALL iom_put( "zdu_bl", tmask(:,:,1)*zdu_bl )           ! du at ml base
+         IF ( iom_use("zdv_bl") ) CALL iom_put( "zdv_bl", tmask(:,:,1)*zdv_bl )           ! dv at ml base
+         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("dstokes") ) CALL iom_put( "dstokes", tmask(:,:,1)*dstokes )      ! Stokes drift penetration depth
          IF ( iom_use("zustke") ) CALL iom_put( "zustke", tmask(:,:,1)*zustke )            ! Stokes drift magnitude at T-points
@@ -1276,21 +1221,1228 @@
          IF ( iom_use("zwstrl") ) CALL iom_put( "zwstrl", tmask(:,:,1)*zwstrl )         ! Langmuir velocity scale
          IF ( iom_use("zustar") ) CALL iom_put( "zustar", tmask(:,:,1)*zustar )         ! friction velocity scale
+         IF ( iom_use("zvstr") ) CALL iom_put( "zvstr", tmask(:,:,1)*zvstr )         ! mixed velocity scale
+         IF ( iom_use("zla") ) CALL iom_put( "zla", tmask(:,:,1)*zla )         ! langmuir #
          IF ( iom_use("wind_power") ) CALL iom_put( "wind_power", 1000.*rho0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine
          IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rho0*tmask(:,:,1)*zustar**2*zustke )
          IF ( iom_use("zhbl") ) CALL iom_put( "zhbl", tmask(:,:,1)*zhbl )               ! BL depth internal to zdf_osm routine
          IF ( iom_use("zhml") ) CALL iom_put( "zhml", tmask(:,:,1)*zhml )               ! ML depth internal to zdf_osm routine
-         IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh )               ! ML depth internal to zdf_osm routine
+         IF ( iom_use("imld") ) CALL iom_put( "imld", tmask(:,:,1)*imld )               ! index for ML depth internal to zdf_osm routine
+         IF ( iom_use("zdh") ) CALL iom_put( "zdh", tmask(:,:,1)*zdh )                  ! pyc thicknessh internal to zdf_osm routine
          IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol )               ! ML depth internal to zdf_osm routine
-         IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav )               ! ML depth internal to zdf_osm routine
-         IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent )               ! ML depth internal to zdf_osm routine
-         IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml )               ! average T in ML
+         IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav )         ! upward BL-avged turb temp flux
+         IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent )   ! upward turb temp entrainment flux
+         IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent )      ! upward turb buoyancy entrainment flux
+         IF ( iom_use("zws_ent") ) CALL iom_put( "zws_ent", tmask(:,:,1)*zws_ent )      ! upward turb salinity entrainment flux
+         IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml )            ! average T in ML
+
+         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") ) CALL iom_put( "zwb_fk", tmask(:,:,1)*zwb_fk )         ! FK b flux
+         IF ( iom_use("zwb_fk_b") ) CALL iom_put( "zwb_fk_b", tmask(:,:,1)*zwb_fk_b )   ! FK b flux averaged over ML
+         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") ) CALL iom_put( "zdiff_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt
+         IF ( iom_use("zvel_mle") ) CALL iom_put( "zvel_mle", tmask(:,:,1)*zdiff_mle )! FK diff in MLE at t-pt
+
       END IF
-      ! Lateral boundary conditions on p_avt  (sign unchanged)
-      CALL lbc_lnk( 'zdfosm', p_avt(:,:,:), 'W', 1.0_wp )
+
+CONTAINS
+! subroutine code changed, needs syntax checking.
+  SUBROUTINE zdf_osm_diffusivity_viscosity( zdiffut, zviscos )
+
+!!---------------------------------------------------------------------
+     !!                   ***  ROUTINE zdf_osm_diffusivity_viscosity  ***
+     !!
+     !! ** Purpose : Determines the eddy diffusivity and eddy viscosity profiles in the mixed layer and the pycnocline.
+     !!
+     !! ** Method  : 
+     !!
+     !! !!----------------------------------------------------------------------
+     REAL(wp), DIMENSION(:,:,:) :: zdiffut
+     REAL(wp), DIMENSION(:,:,:) :: zviscos
+! local
+
+! Scales used to calculate eddy diffusivity and viscosity profiles
+      REAL(wp), DIMENSION(jpi,jpj) :: zdifml_sc, zvisml_sc
+      REAL(wp), DIMENSION(jpi,jpj) :: zdifpyc_n_sc, zdifpyc_s_sc, zdifpyc_shr
+      REAL(wp), DIMENSION(jpi,jpj) :: zvispyc_n_sc, zvispyc_s_sc,zvispyc_shr
+      REAL(wp), DIMENSION(jpi,jpj) :: zbeta_d_sc, zbeta_v_sc
+!
+      REAL(wp) :: zvel_sc_pyc, zvel_sc_ml, zstab_fac
+      
+      REAL(wp), PARAMETER :: rn_dif_ml = 0.8, rn_vis_ml = 0.375
+      REAL(wp), PARAMETER :: rn_dif_pyc = 0.15, rn_vis_pyc = 0.142
+      REAL(wp), PARAMETER :: rn_vispyc_shr = 0.15
+      
+      DO_2D( 0, 0, 0, 0 )
+          IF ( lconv(ji,jj) ) THEN
+          
+            zvel_sc_pyc = ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.25 * zshear(ji,jj) * zhbl(ji,jj) )**pthird
+            zvel_sc_ml = ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
+            zstab_fac = ( zhml(ji,jj) / zvel_sc_ml * ( 1.4 - 0.4 / ( 1.0 + EXP(-3.5 * LOG10(-zhol(ji,jj) ) ) )**1.25 ) )**2
+
+            zdifml_sc(ji,jj) = rn_dif_ml * zhml(ji,jj) * zvel_sc_ml
+            zvisml_sc(ji,jj) = rn_vis_ml * zdifml_sc(ji,jj)
+
+            IF ( lpyc(ji,jj) ) THEN
+              zdifpyc_n_sc(ji,jj) =  rn_dif_pyc * zvel_sc_ml * zdh(ji,jj)
+
+              IF ( lshear(ji,jj) .and. j_ddh(ji,jj) == 1 ) THEN
+                zdifpyc_n_sc(ji,jj) = zdifpyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj) )**pthird * zhbl(ji,jj)
+              ENDIF
+            
+              zdifpyc_s_sc(ji,jj) = zwb_ent(ji,jj) + 0.0025 * zvel_sc_pyc * ( zhbl(ji,jj) / zdh(ji,jj) - 1.0 ) * ( zb_ml(ji,jj) - zb_bl(ji,jj) )
+              zdifpyc_s_sc(ji,jj) = 0.09 * zdifpyc_s_sc(ji,jj) * zstab_fac
+              zdifpyc_s_sc(ji,jj) = MAX( zdifpyc_s_sc(ji,jj), -0.5 * zdifpyc_n_sc(ji,jj) )
+              
+              zvispyc_n_sc(ji,jj) = 0.09 * zvel_sc_pyc * ( 1.0 - zhbl(ji,jj) / zdh(ji,jj) )**2 * ( 0.005 * ( zu_ml(ji,jj)-zu_bl(ji,jj) )**2 + 0.0075 * ( zv_ml(ji,jj)-zv_bl(ji,jj) )**2 ) / zdh(ji,jj)
+              zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * zdh(ji,jj) + zvispyc_n_sc(ji,jj) * zstab_fac
+              IF ( lshear(ji,jj) .and. j_ddh(ji,jj) == 1 ) THEN
+                zvispyc_n_sc(ji,jj) = zvispyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj ) )**pthird * zhbl(ji,jj)
+              ENDIF
+             
+              zvispyc_s_sc(ji,jj) = 0.09 * ( zwb_min(ji,jj) + 0.0025 * zvel_sc_pyc * ( zhbl(ji,jj) / zdh(ji,jj) - 1.0 ) * ( zb_ml(ji,jj) - zb_bl(ji,jj) ) )
+              zvispyc_s_sc(ji,jj) = zvispyc_s_sc(ji,jj) * zstab_fac
+              zvispyc_s_sc(ji,jj) = MAX( zvispyc_s_sc(ji,jj), -0.5 * zvispyc_s_sc(ji,jj) )
+
+              zbeta_d_sc(ji,jj) = 1.0 - ( ( zdifpyc_n_sc(ji,jj) + 1.4 * zdifpyc_s_sc(ji,jj) ) / ( zdifml_sc(ji,jj) + epsln ) )**p2third
+              zbeta_v_sc(ji,jj) = 1.0 -  2.0 * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) / ( zvisml_sc(ji,jj) + epsln )
+            ELSE
+              zbeta_d_sc = 1.0
+              zbeta_v_sc = 1.0
+            ENDIF
+          ELSE
+            zdifml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp)
+            zvisml_sc(ji,jj) = zvstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( zhol(ji,jj) / 0.6_wp )**2 ), 0.2_wp)
+          END IF
+      END_2D
+!
+       DO_2D( 0, 0, 0, 0 )
+          IF ( lconv(ji,jj) ) THEN
+             DO jk = 2, imld(ji,jj)   ! mixed layer diffusivity
+                 zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj)
+                 !
+                 zdiffut(ji,jj,jk) = zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml )**1.5
+                 !
+                 zviscos(ji,jj,jk) = zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml ) &
+   &            *                                      ( 1.0 - 0.5 * zznd_ml**2 )
+             END DO
+! pycnocline
+             IF ( lpyc(ji,jj) ) THEN
+! Diffusivity profile in the pycnocline given by cubic polynomial.
+                za_cubic = 0.5
+                zb_cubic = -1.75 * zdifpyc_s_sc(ji,jj) / zdifpyc_n_sc(ji,jj)
+                zd_cubic = ( zdh(ji,jj) * zdifml_sc(ji,jj) / zhml(ji,jj) * SQRT( 1.0 - zbeta_d_sc(ji,jj) ) * ( 2.5 * zbeta_d_sc(ji,jj) - 1.0 ) &
+                     & - 0.85 * zdifpyc_s_sc(ji,jj) ) / MAX(zdifpyc_n_sc(ji,jj), 1.e-8)
+                zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic  - zb_cubic )
+                zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic
+                DO jk = imld(ji,jj) , ibld(ji,jj)
+                  zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6)
+                      !
+                  zdiffut(ji,jj,jk) = zdifpyc_n_sc(ji,jj) * ( za_cubic + zb_cubic * zznd_pyc + zc_cubic * zznd_pyc**2 +   zd_cubic * zznd_pyc**3 )
+
+                  zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdifpyc_s_sc(ji,jj) * ( 1.75 * zznd_pyc - 0.15 * zznd_pyc**2 - 0.2 * zznd_pyc**3 )
+                END DO
+ ! viscosity profiles.
+                za_cubic = 0.5
+                zb_cubic = -1.75 * zvispyc_s_sc(ji,jj) / zvispyc_n_sc(ji,jj)
+                zd_cubic = ( 0.5 * zvisml_sc(ji,jj) * zdh(ji,jj) / zhml(ji,jj) - 0.85 * zvispyc_s_sc(ji,jj)  )  / MAX(zvispyc_n_sc(ji,jj), 1.e-8)
+                zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic - zd_cubic )
+                zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic
+                DO jk = imld(ji,jj) , ibld(ji,jj)
+                   zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6)
+                    zviscos(ji,jj,jk) = zvispyc_n_sc(ji,jj) * ( za_cubic + zb_cubic * zznd_pyc + zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 )
+                    zviscos(ji,jj,jk) = zviscos(ji,jj,jk) + zvispyc_s_sc(ji,jj) * ( 1.75 * zznd_pyc - 0.15 * zznd_pyc**2 -0.2 * zznd_pyc**3 )
+                END DO
+                IF ( zdhdt(ji,jj) > 0._wp ) THEN
+                 zdiffut(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w(ji,jj,ibld(ji,jj)+1,Kmm), 1.0e-6 )
+                 zviscos(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w(ji,jj,ibld(ji,jj)+1,Kmm), 1.0e-6 )
+                ELSE
+                  zdiffut(ji,jj,ibld(ji,jj)) = 0._wp
+                  zviscos(ji,jj,ibld(ji,jj)) = 0._wp
+                ENDIF
+             ENDIF
+          ELSE
+          ! stable conditions
+             DO jk = 2, ibld(ji,jj)
+                zznd_ml = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
+                zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5
+                zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 )
+             END DO
+
+             IF ( zdhdt(ji,jj) > 0._wp ) THEN
+                zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w(ji, jj, ibld(ji,jj), Kmm)
+                zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w(ji, jj, ibld(ji,jj), Kmm)
+             ENDIF
+          ENDIF   ! end if ( lconv )
+          !
+       END_2D
+       
+  END SUBROUTINE zdf_osm_diffusivity_viscosity
+  
+  SUBROUTINE zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear, zri_i )
+
+!!---------------------------------------------------------------------
+     !!                   ***  ROUTINE zdf_osm_osbl_state  ***
+     !!
+     !! ** Purpose : Determines the state of the OSBL, stable/unstable, shear/ noshear. Also determines shear production, entrainment buoyancy flux and interfacial Richardson number
+     !!
+     !! ** Method  : 
+     !!
+     !! !!----------------------------------------------------------------------
+
+     INTEGER, DIMENSION(jpi,jpj) :: j_ddh  ! j_ddh = 0, active shear layer; j_ddh=1, shear layer not active; j_ddh=2 shear production low.
+     
+     LOGICAL, DIMENSION(jpi,jpj) :: lconv, lshear
+
+     REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent, zwb_min ! Buoyancy fluxes at base of well-mixed layer.
+     REAL(wp), DIMENSION(jpi,jpj) :: zshear  ! production of TKE due to shear across the pycnocline
+     REAL(wp), DIMENSION(jpi,jpj) :: zri_i  ! Interfacial Richardson Number
+
+! Local Variables
+
+     INTEGER :: jj, ji
+     
+     REAL(wp), DIMENSION(jpi,jpj) :: zekman
+     REAL(wp) :: 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, PARAMETER :: za_shr = 0.4, zb_shr = 6.5, za_wb_s = 0.1
+     REAL, PARAMETER :: rn_ri_thres_a = 0.5, rn_ri_thresh_b = 0.59
+     REAL, PARAMETER :: zalpha_c = 0.2, zalpha_lc = 0.04     
+     REAL, PARAMETER :: zalpha_ls = 0.06, zalpha_s = 0.15
+     REAL, PARAMETER :: rn_ri_p_thresh = 27.0
+     REAL, PARAMETER :: zrot=0._wp  ! dummy rotation rate of surface stress.
+     
+! Determins stability and set flag lconv
+     DO_2D( 0, 0, 0, 0 )
+      IF ( zhol(ji,jj) < 0._wp ) THEN
+         lconv(ji,jj) = .TRUE.
+       ELSE
+          lconv(ji,jj) = .FALSE.
+       ENDIF
+     END_2D
+ 
+     zekman(:,:) = EXP( - 4.0 * ABS( ff_t(:,:) ) * zhbl(:,:) / MAX(zustar(:,:), 1.e-8 ) )
+     
+     WHERE ( lconv )
+       zri_i = zdb_ml * zhml**2 / MAX( ( zvstr**3 + 0.5 * zwstrc**3 )**p2third * zdh, 1.e-12 )
+     END WHERE
+
+     zshear(:,:) = 0._wp
+     j_ddh(:,:) = 1     
+ 
+     DO_2D( 0, 0, 0, 0 )
+      IF ( lconv(ji,jj) ) THEN
+         IF ( zdb_bl(ji,jj) > 0._wp ) THEN
+           zri_p = MAX (  SQRT( zdb_bl(ji,jj) * zdh(ji,jj) / MAX( zdu_bl(ji,jj)**2 + zdv_bl(ji,jj)**2, 1.e-8) )  *  ( zhbl(ji,jj) / zdh(ji,jj) ) * ( zvstr(ji,jj) / MAX( zustar(ji,jj), 1.e-6 ) )**2 &
+                & / MAX( zekman(ji,jj), 1.e-6 )  , 5._wp )
+         
+           zri_b = zdb_ml(ji,jj) * zdh(ji,jj) / MAX( zdu_ml(ji,jj)**2 + zdv_ml(ji,jj)**2, 1.e-8 )
+                     
+           zshear(ji,jj) = za_shr * zekman(ji,jj) * ( MAX( zustar(ji,jj)**2 * zdu_ml(ji,jj) / zhbl(ji,jj), 0._wp ) + zb_shr * MAX( -ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) * zdv_ml(ji,jj) / zhbl(ji,jj), 0._wp ) )
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! Test ensures j_ddh=0 is not selected. Change to zri_p<27 when  !
+! full code available                                          !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+           IF ( zri_p < -rn_ri_p_thresh .and. zshear(ji,jj) > 0._wp ) THEN
+! Growing shear layer
+             j_ddh(ji,jj) = 0
+             lshear(ji,jj) = .TRUE.
+           ELSE
+             j_ddh(ji,jj) = 1
+             IF ( zri_b <= 1.5 .and. zshear(ji,jj) > 0._wp ) THEN
+! shear production large enough to determine layer charcteristics, but can't maintain a shear layer.
+               lshear(ji,jj) = .TRUE.
+             ELSE
+! Shear production may not be zero, but is small and doesn't determine characteristics of pycnocline.
+               zshear(ji,jj) = 0.5 * zshear(ji,jj)
+               lshear(ji,jj) = .FALSE.
+             ENDIF 
+           ENDIF                
+         ELSE                ! zdb_bl test, note zshear set to zero
+           j_ddh(ji,jj) = 2
+           lshear(ji,jj) = .FALSE.
+         ENDIF
+       ENDIF
+     END_2D
+ 
+! Calculate entrainment buoyancy flux due to surface fluxes.
+
+     DO_2D( 0, 0, 0, 0 )
+      IF ( lconv(ji,jj) ) THEN
+        zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln
+        zrf_conv = TANH( ( zwstrc(ji,jj) / zwcor )**0.69 )
+        zrf_shear = TANH( ( zustar(ji,jj) / zwcor )**0.69 )
+        zrf_langmuir = TANH( ( zwstrl(ji,jj) / zwcor )**0.69 )
+        IF (nn_osm_SD_reduce > 0 ) THEN
+        ! Effective Stokes drift already reduced from surface value
+           zr_stokes = 1.0_wp
+        ELSE
+         ! Effective Stokes drift only reduced by factor rn_zdfodm_adjust_sd,
+          ! requires further reduction where BL is deep
+           zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) &
+         &                  * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) )
+        END IF
+        zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) &
+               &                  - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) &
+               &         + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 &
+               &                                         - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj)
+          !
+      ENDIF
+     END_2D
+
+     zwb_min(:,:) = 0._wp
+
+     DO_2D( 0, 0, 0, 0 )
+      IF ( lshear(ji,jj) ) THEN
+        IF ( lconv(ji,jj) ) THEN
+! Unstable OSBL
+           zwb_shr = -za_wb_s * zshear(ji,jj)
+           IF ( j_ddh(ji,jj) == 0 ) THEN
+
+! Developing shear layer, additional shear production possible.
+
+             zshear_u = MAX( zustar(ji,jj)**2 * zdu_ml(ji,jj) /  zhbl(ji,jj), 0._wp )
+             zshear(ji,jj) = zshear(ji,jj) + zshear_u * ( 1.0 - MIN( zri_p / rn_ri_p_thresh, 1.d0 ) )
+             zshear(ji,jj) = MIN( zshear(ji,jj), zshear_u )
+             
+             zwb_shr = -za_wb_s * zshear(ji,jj)
+             
+           ENDIF                
+           zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr
+           zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj)
+        ELSE    ! IF ( lconv ) THEN - ENDIF
+! Stable OSBL  - shear production not coded for first attempt.           
+        ENDIF  ! lconv
+      ELSE  ! lshear
+        IF ( lconv(ji,jj) ) THEN
+! Unstable OSBL
+           zwb_shr = -za_wb_s * zshear(ji,jj)
+           zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr
+           zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj)
+        ENDIF  ! lconv
+      ENDIF    ! lshear
+     END_2D
+   END SUBROUTINE zdf_osm_osbl_state
+     
+     
+   SUBROUTINE zdf_osm_vertical_average( jnlev_av, jp_ext, zt, zs, zb, zu, zv, zdt, zds, zdb, zdu, zdv )
+     !!---------------------------------------------------------------------
+     !!                   ***  ROUTINE zdf_vertical_average  ***
+     !!
+     !! ** Purpose : Determines vertical averages from surface to jnlev.
+     !!
+     !! ** Method  : Averages are calculated from the surface to jnlev.
+     !!              The external level used to calculate differences is ibld+ibld_ext
+     !!
+     !!----------------------------------------------------------------------
+
+        INTEGER, DIMENSION(jpi,jpj) :: jnlev_av  ! Number of levels to average over.
+        INTEGER, DIMENSION(jpi,jpj) :: jp_ext
+
+        ! Alan: do we need zb?
+        REAL(wp), DIMENSION(jpi,jpj) :: zt, zs, zb        ! Average temperature and salinity
+        REAL(wp), DIMENSION(jpi,jpj) :: zu,zv         ! Average current components
+        REAL(wp), DIMENSION(jpi,jpj) :: zdt, zds, zdb ! Difference between average and value at base of OSBL
+        REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv      ! Difference for velocity components.
+
+        INTEGER :: jk, ji, jj, ibld_ext
+        REAL(wp) :: zthick, zthermal, zbeta
+
+
+        zt   = 0._wp
+        zs   = 0._wp
+        zu   = 0._wp
+        zv   = 0._wp
+        DO_2D( 0, 0, 0, 0 )
+         zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
+         zbeta    = rab_n(ji,jj,1,jp_sal)
+            ! average over depth of boundary layer
+         zthick = epsln
+         DO jk = 2, jnlev_av(ji,jj)
+            zthick = zthick + e3t(ji,jj,jk,Kmm)
+            zt(ji,jj)   = zt(ji,jj)  + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm)
+            zs(ji,jj)   = zs(ji,jj)  + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm)
+            zu(ji,jj)   = zu(ji,jj)  + e3t(ji,jj,jk,Kmm) &
+                  &            * ( uu(ji,jj,jk,Kbb) + uu(ji - 1,jj,jk,Kbb) ) &
+                  &            / MAX( 1. , umask(ji,jj,jk) + umask(ji - 1,jj,jk) )
+            zv(ji,jj)   = zv(ji,jj)  + e3t(ji,jj,jk,Kmm) &
+                  &            * ( vv(ji,jj,jk,Kbb) + vv(ji,jj - 1,jk,Kbb) ) &
+                  &            / MAX( 1. , vmask(ji,jj,jk) + vmask(ji,jj - 1,jk) )
+         END DO
+         zt(ji,jj) = zt(ji,jj) / zthick
+         zs(ji,jj) = zs(ji,jj) / zthick
+         zu(ji,jj) = zu(ji,jj) / zthick
+         zv(ji,jj) = zv(ji,jj) / zthick
+         zb(ji,jj) = grav * zthermal * zt(ji,jj) - grav * zbeta * zs(ji,jj)
+         ibld_ext = jnlev_av(ji,jj) + jp_ext(ji,jj)
+         IF ( ibld_ext < mbkt(ji,jj) ) THEN
+           zdt(ji,jj) = zt(ji,jj) - ts(ji,jj,ibld_ext,jp_tem,Kmm)
+           zds(ji,jj) = zs(ji,jj) - ts(ji,jj,ibld_ext,jp_sal,Kmm)
+           zdu(ji,jj) = zu(ji,jj) - ( uu(ji,jj,ibld_ext,Kbb) + uu(ji-1,jj,ibld_ext,Kbb ) ) &
+                  &    / MAX(1. , umask(ji,jj,ibld_ext ) + umask(ji-1,jj,ibld_ext ) )
+           zdv(ji,jj) = zv(ji,jj) - ( vv(ji,jj,ibld_ext,Kbb) + vv(ji,jj-1,ibld_ext,Kbb ) ) &
+                  &   / MAX(1. , vmask(ji,jj,ibld_ext ) + vmask(ji,jj-1,ibld_ext ) )
+           zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj)
+         ELSE
+           zdt(ji,jj) = 0._wp
+           zds(ji,jj) = 0._wp
+           zdu(ji,jj) = 0._wp
+           zdv(ji,jj) = 0._wp
+           zdb(ji,jj) = 0._wp
+         ENDIF
+        END_2D
+   END SUBROUTINE zdf_osm_vertical_average
+
+   SUBROUTINE zdf_osm_velocity_rotation( zcos_w, zsin_w, zu, zv, zdu, zdv )
+     !!---------------------------------------------------------------------
+     !!                   ***  ROUTINE zdf_velocity_rotation  ***
+     !!
+     !! ** Purpose : Rotates frame of reference of averaged velocity components.
+     !!
+     !! ** Method  : The velocity components are rotated into frame specified by zcos_w and zsin_w
+     !!
+     !!----------------------------------------------------------------------
+
+        REAL(wp), DIMENSION(jpi,jpj) :: zcos_w, zsin_w       ! Cos and Sin of rotation angle
+        REAL(wp), DIMENSION(jpi,jpj) :: zu, zv               ! Components of current
+        REAL(wp), DIMENSION(jpi,jpj) :: zdu, zdv             ! Change in velocity components across pycnocline
+
+        INTEGER :: ji, jj
+        REAL(wp) :: ztemp
+
+        DO_2D( 0, 0, 0, 0 )
+           ztemp = zu(ji,jj)
+           zu(ji,jj) = zu(ji,jj) * zcos_w(ji,jj) + zv(ji,jj) * zsin_w(ji,jj)
+           zv(ji,jj) = zv(ji,jj) * zcos_w(ji,jj) - ztemp * zsin_w(ji,jj)
+           ztemp = zdu(ji,jj)
+           zdu(ji,jj) = zdu(ji,jj) * zcos_w(ji,jj) + zdv(ji,jj) * zsin_w(ji,jj)
+           zdv(ji,jj) = zdv(ji,jj) * zsin_w(ji,jj) - ztemp * zsin_w(ji,jj)
+        END_2D
+    END SUBROUTINE zdf_osm_velocity_rotation
+
+    SUBROUTINE zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk )
+     !!---------------------------------------------------------------------
+     !!                   ***  ROUTINE zdf_osm_osbl_state_fk  ***
+     !!
+     !! ** Purpose : Determines the state of the OSBL and MLE layer. Info is returned in the logicals lpyc,lflux and lmle. Used with Fox-Kemper scheme.
+     !!  lpyc :: determines whether pycnocline flux-grad relationship needs to be determined
+     !!  lflux :: determines whether effects of surface flux extend below the base of the OSBL
+     !!  lmle  :: determines whether the layer with MLE is increasing with time or if base is relaxing towards hbl. 
+     !!
+     !! ** Method  : 
+     !!
+     !! 
+     !!----------------------------------------------------------------------
+      
+! Outputs
+      LOGICAL,  DIMENSION(jpi,jpj)  :: lpyc, lflux, lmle
+      REAL(wp), DIMENSION(jpi,jpj)  :: zwb_fk
+!
+      REAL(wp), DIMENSION(jpi,jpj)  :: znd_param
+      REAL(wp)                      :: zbuoy, ztmp, zpe_mle_layer
+      REAL(wp)                      :: zpe_mle_ref, zwb_ent, zdbdz_mle_int
+      
+      znd_param(:,:) = 0._wp
+
+        DO_2D( 0, 0, 0, 0 )
+          ztmp =  r1_ft(ji,jj) *  MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf
+          zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) * ztmp * zdbds_mle(ji,jj) * zdbds_mle(ji,jj)
+        END_2D
+        DO_2D( 0, 0, 0, 0 )
+                 !
+         IF ( lconv(ji,jj) ) THEN
+           IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN
+             zt_mle(ji,jj) = ( zt_mle(ji,jj) * zhmle(ji,jj) - zt_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) )
+             zs_mle(ji,jj) = ( zs_mle(ji,jj) * zhmle(ji,jj) - zs_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) )
+             zb_mle(ji,jj) = ( zb_mle(ji,jj) * zhmle(ji,jj) - zb_bl(ji,jj) * zhbl(ji,jj) ) / ( zhmle(ji,jj) - zhbl(ji,jj) )
+             zdbdz_mle_int = ( zb_bl(ji,jj) - ( 2.0 * zb_mle(ji,jj) -zb_bl(ji,jj) ) ) / ( zhmle(ji,jj) - zhbl(ji,jj) )
+! Calculate potential energies of actual profile and reference profile.
+             zpe_mle_layer = 0._wp
+             zpe_mle_ref = 0._wp
+             DO jk = ibld(ji,jj), mld_prof(ji,jj)
+               zbuoy = grav * ( zthermal * ts(ji,jj,jk,jp_tem,Kmm) - zbeta * ts(ji,jj,jk,jp_sal,Kmm) )
+               zpe_mle_layer = zpe_mle_layer + zbuoy * gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)
+               zpe_mle_ref = zpe_mle_ref + ( zb_bl(ji,jj) - zdbdz_mle_int * ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) ) * gdepw(ji,jj,jk,Kmm) * e3w(ji,jj,jk,Kmm)
+             END DO
+! Non-dimensional parameter to diagnose the presence of thermocline
+                
+             znd_param(ji,jj) = ( zpe_mle_layer - zpe_mle_ref ) * ABS( ff_t(ji,jj) ) / ( MAX( zwb_fk(ji,jj), 1.0e-10 ) * zhmle(ji,jj) )
+           ENDIF
+         ENDIF
+        END_2D
+
+! Diagnosis
+        DO_2D( 0, 0, 0, 0 )
+          IF ( lconv(ji,jj) ) THEN
+            zwb_ent = - 2.0 * 0.2 * zwbav(ji,jj) &
+               &                  - 0.15 * zustar(ji,jj)**3 /zhml(ji,jj) &
+               &         + ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zustar(ji,jj)**3 &
+               &         - 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj)
+            IF ( -2.0 * zwb_fk(ji,jj) / zwb_ent > 0.5 ) THEN
+              IF ( zhmle(ji,jj) > 1.2 * zhbl(ji,jj) ) THEN
+! MLE layer growing
+                IF ( znd_param (ji,jj) > 100. ) THEN
+! Thermocline present
+                  lflux(ji,jj) = .FALSE.
+                  lmle(ji,jj) =.FALSE.
+                ELSE
+! Thermocline not present
+                  lflux(ji,jj) = .TRUE.
+                  lmle(ji,jj) = .TRUE.
+                ENDIF  ! znd_param > 100
+!
+                IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN
+                  lpyc(ji,jj) = .FALSE.
+                ELSE
+                   lpyc = .TRUE.
+                ENDIF
+              ELSE
+! MLE layer restricted to OSBL or just below.
+                IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) THEN
+! Weak stratification MLE layer can grow.
+                  lpyc(ji,jj) = .FALSE.
+                  lflux(ji,jj) = .TRUE.
+                  lmle(ji,jj) = .TRUE.
+                ELSE
+! Strong stratification
+                  lpyc(ji,jj) = .TRUE.
+                  lflux(ji,jj) = .FALSE.
+                  lmle(ji,jj) = .FALSE.
+                ENDIF ! zdb_bl < rn_mle_thresh_bl and 
+              ENDIF  ! zhmle > 1.2 zhbl
+            ELSE
+              lpyc(ji,jj) = .TRUE.
+              lflux(ji,jj) = .FALSE.
+              lmle(ji,jj) = .FALSE.
+              IF ( zdb_bl(ji,jj) < rn_osm_bl_thresh ) lpyc(ji,jj) = .FALSE.
+            ENDIF !  -2.0 * zwb_fk(ji,jj) / zwb_ent > 0.5 
+          ELSE
+! Stable Boundary Layer
+            lpyc(ji,jj) = .FALSE.
+            lflux(ji,jj) = .FALSE.
+            lmle(ji,jj) = .FALSE.
+          ENDIF  ! lconv
+        END_2D
+    END SUBROUTINE zdf_osm_osbl_state_fk
+
+    SUBROUTINE zdf_osm_external_gradients(jbase, zdtdz, zdsdz, zdbdz )
+     !!---------------------------------------------------------------------
+     !!                   ***  ROUTINE zdf_osm_external_gradients  ***
+     !!
+     !! ** Purpose : Calculates the gradients below the OSBL
+     !!
+     !! ** Method  : Uses ibld and ibld_ext to determine levels to calculate the gradient.
+     !!
+     !!----------------------------------------------------------------------
+
+     INTEGER, DIMENSION(jpi,jpj)  :: jbase
+     REAL(wp), DIMENSION(jpi,jpj) :: zdtdz, zdsdz, zdbdz   ! External gradients of temperature, salinity and buoyancy.
+
+     INTEGER :: jj, ji, jkb, jkb1
+     REAL(wp) :: zthermal, zbeta
+
+
+     DO_2D( 0, 0, 0, 0 )
+        IF ( jbase(ji,jj)+1 < mbkt(ji,jj) ) THEN
+           zthermal = rab_n(ji,jj,1,jp_tem) !ideally use ibld not 1??
+           zbeta    = rab_n(ji,jj,1,jp_sal)
+           jkb = jbase(ji,jj)
+           jkb1 = MIN(jkb + 1, mbkt(ji,jj))
+           zdtdz(ji,jj) = - ( ts(ji,jj,jkb1,jp_tem,Kmm) - ts(ji,jj,jkb,jp_tem,Kmm ) ) &
+                &   / e3t(ji,jj,ibld(ji,jj),Kmm)
+           zdsdz(ji,jj) = - ( ts(ji,jj,jkb1,jp_sal,Kmm) - ts(ji,jj,jkb,jp_sal,Kmm ) ) &
+                &   / e3t(ji,jj,ibld(ji,jj),Kmm)
+           zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj)
+        ELSE
+           zdtdz(ji,jj) = 0._wp
+           zdsdz(ji,jj) = 0._wp
+           zdbdz(ji,jj) = 0._wp
+        END IF
+     END_2D
+    END SUBROUTINE zdf_osm_external_gradients
+
+    SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz, zalpha )
+
+     REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz      ! gradients in the pycnocline
+     REAL(wp), DIMENSION(jpi,jpj) :: zalpha
+
+     INTEGER :: jk, jj, ji
+     REAL(wp) :: ztgrad, zsgrad, zbgrad
+     REAL(wp) :: zgamma_b_nd, znd
+     REAL(wp) :: zzeta_m, zzeta_en, zbuoy_pyc_sc
+     REAL(wp), PARAMETER :: zgamma_b = 2.25, zzeta_sh = 0.15
+
+     DO_2D( 0, 0, 0, 0 )
+        IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
+           IF ( lconv(ji,jj) ) THEN  ! convective conditions
+             IF ( lpyc(ji,jj) ) THEN
+                zzeta_m = 0.1 + 0.3 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) )
+                zalpha(ji,jj) = 2.0 * ( 1.0 - ( 0.80 * zzeta_m + 0.5 * SQRT( 3.14159 / zgamma_b ) ) * zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) ) / ( 0.723 + SQRT( 3.14159 / zgamma_b ) )
+                zalpha(ji,jj) = MAX( zalpha(ji,jj), 0._wp )
+
+                ztmp = 1._wp/MAX(zdh(ji,jj), epsln)
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+! Commented lines in this section are not needed in new code, once tested !
+! can be removed                                                          !
+!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+!                   ztgrad = zalpha * zdt_ml(ji,jj) * ztmp + zdtdz_bl_ext(ji,jj)
+!                   zsgrad = zalpha * zds_ml(ji,jj) * ztmp + zdsdz_bl_ext(ji,jj)
+                zbgrad = zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp + zdbdz_bl_ext(ji,jj)
+                zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln)
+                DO jk = 2, ibld(ji,jj)+ibld_ext
+                  znd = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) * ztmp
+                  IF ( znd <= zzeta_m ) THEN
+!                        zdtdz(ji,jj,jk) = zdtdz_bl_ext(ji,jj) + zalpha * zdt_ml(ji,jj) * ztmp * &
+!                &        EXP( -6.0 * ( znd -zzeta_m )**2 )
+!                        zdsdz(ji,jj,jk) = zdsdz_bl_ext(ji,jj) + zalpha * zds_ml(ji,jj) * ztmp * &
+!                                  & EXP( -6.0 * ( znd -zzeta_m )**2 )
+                     zdbdz(ji,jj,jk) = zdbdz_bl_ext(ji,jj) + zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp * &
+                               & EXP( -6.0 * ( znd -zzeta_m )**2 )
+                  ELSE
+!                         zdtdz(ji,jj,jk) =  ztgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 )
+!                         zdsdz(ji,jj,jk) =  zsgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 )
+                      zdbdz(ji,jj,jk) =  zbgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 )
+                  ENDIF
+               END DO
+            ENDIF ! if no pycnocline pycnocline gradients set to zero
+           ELSE
+              ! stable conditions
+              ! if pycnocline profile only defined when depth steady of increasing.
+              IF ( zdhdt(ji,jj) > 0.0 ) THEN        ! Depth increasing, or steady.
+                 IF ( zdb_bl(ji,jj) > 0._wp ) THEN
+                    IF ( zhol(ji,jj) >= 0.5 ) THEN      ! Very stable - 'thick' pycnocline
+                       ztmp = 1._wp/MAX(zhbl(ji,jj), epsln)
+                       ztgrad = zdt_bl(ji,jj) * ztmp
+                       zsgrad = zds_bl(ji,jj) * ztmp
+                       zbgrad = zdb_bl(ji,jj) * ztmp
+                       DO jk = 2, ibld(ji,jj)
+                          znd = gdepw(ji,jj,jk,Kmm) * ztmp
+                          zdtdz(ji,jj,jk) =  ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
+                          zdbdz(ji,jj,jk) =  zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
+                          zdsdz(ji,jj,jk) =  zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
+                       END DO
+                    ELSE                                   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form.
+                       ztmp = 1._wp/MAX(zdh(ji,jj), epsln)
+                       ztgrad = zdt_bl(ji,jj) * ztmp
+                       zsgrad = zds_bl(ji,jj) * ztmp
+                       zbgrad = zdb_bl(ji,jj) * ztmp
+                       DO jk = 2, ibld(ji,jj)
+                          znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) * ztmp
+                          zdtdz(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
+                          zdbdz(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
+                          zdsdz(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
+                       END DO
+                    ENDIF ! IF (zhol >=0.5)
+                 ENDIF    ! IF (zdb_bl> 0.)
+              ENDIF       ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero
+           ENDIF          ! IF (lconv)
+        ENDIF      ! IF ( ibld(ji,jj) < mbkt(ji,jj) )
+     END_2D
+
+    END SUBROUTINE zdf_osm_pycnocline_scalar_profiles
+
+    SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE zdf_osm_pycnocline_shear_profiles  ***
+      !!
+      !! ** Purpose : Calculates velocity shear in the pycnocline
+      !!
+      !! ** Method  :
+      !!
+      !!----------------------------------------------------------------------
+
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz, zdvdz
+
+      INTEGER :: jk, jj, ji
+      REAL(wp) :: zugrad, zvgrad, znd
+      REAL(wp) :: zzeta_v = 0.45
       !
-   END SUBROUTINE zdf_osm
-
-
-   SUBROUTINE zdf_osm_init( Kmm ) 
+      DO_2D( 0, 0, 0, 0 )
+         !
+         IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
+            IF ( lconv (ji,jj) ) THEN
+               ! Unstable conditions. Shouldn;t be needed with no pycnocline code.
+!                  zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / &
+!                       &      ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * &
+!                      &      MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 ))
+               !Alan is this right?
+!                  zvgrad = ( 0.7 * zdv_ml(ji,jj) + &
+!                       &    2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / &
+!                       &          ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird  + epsln ) &
+!                       &      )/ (zdh(ji,jj)  + epsln )
+!                  DO jk = 2, ibld(ji,jj) - 1 + ibld_ext
+!                     znd = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v
+!                     IF ( znd <= 0.0 ) THEN
+!                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd )
+!                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd )
+!                     ELSE
+!                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd )
+!                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd )
+!                     ENDIF
+!                  END DO
+            ELSE
+               ! stable conditions
+               zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj)
+               zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj)
+               DO jk = 2, ibld(ji,jj)
+                  znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
+                  IF ( znd < 1.0 ) THEN
+                     zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 )
+                  ELSE
+                     zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 )
+                  ENDIF
+                  zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 )
+               END DO
+            ENDIF
+            !
+         END IF      ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) )
+      END_2D
+    END SUBROUTINE zdf_osm_pycnocline_shear_profiles
+
+   SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zddhdt )
+     !!---------------------------------------------------------------------
+     !!                   ***  ROUTINE zdf_osm_calculate_dhdt  ***
+     !!
+     !! ** Purpose : Calculates the rate at which hbl changes.
+     !!
+     !! ** Method  :
+     !!
+     !!----------------------------------------------------------------------
+
+    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zddhdt        ! Rate of change of hbl
+
+    INTEGER :: jj, ji
+    REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert, zpsi
+    REAL(wp) :: zvel_max!, zwb_min
+    REAL(wp) :: zzeta_m = 0.3
+    REAL(wp) :: zgamma_c = 2.0
+    REAL(wp) :: zdhoh = 0.1
+    REAL(wp) :: alpha_bc = 0.5
+    REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth
+ 
+  DO_2D( 0, 0, 0, 0 )
+    
+    IF ( lshear(ji,jj) ) THEN
+       IF ( lconv(ji,jj) ) THEN    ! Convective
+
+          IF ( ln_osm_mle ) THEN
+
+             IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN
+       ! Fox-Kemper buoyancy flux average over OSBL
+                zwb_fk_b(ji,jj) = zwb_fk(ji,jj) *  &
+                     (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) )
+             ELSE
+                zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj)
+             ENDIF
+             zvel_max = ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj)
+             IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN
+                ! OSBL is deepening, entrainment > restratification
+                IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_bl_ext(ji,jj) > 0.0 ) THEN
+! *** Used for shear Needs to be changed to work stabily
+!                zgamma_b_nd = zdbdz_bl_ext * dh / zdb_ml
+!                zalpha_b = 6.7 * zgamma_b_nd / ( 1.0 + zgamma_b_nd )
+!                zgamma_b = zgamma_b_nd / ( 0.12 * ( 1.25 + zgamma_b_nd ) )
+!                za_1 = 1.0 / zgamma_b**2 - 0.017
+!                za_2 = 1.0 / zgamma_b**3 - 0.0025
+!                zpsi = zalpha_b * ( 1.0 + zgamma_b_nd ) * ( za_1 - 2.0 * za_2 * dh / hbl )
+                   zpsi = 0._wp
+                   zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) )
+                   zdhdt(ji,jj) = zdhdt(ji,jj)! - zpsi * ( -1.0 / zhml(ji,jj) + 2.4 * zdbdz_bl_ext(ji,jj) / zdb_ml(ji,jj) ) * zwb_min(ji,jj) * zdh(ji,jj) / zdb_bl(ji,jj)
+                   IF ( j_ddh(ji,jj) == 1 ) THEN
+                     IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN
+                        zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )
+                     ELSE
+                        zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )
+                     ENDIF
+! Relaxation to dh_ref = zari * hbl
+                     zddhdt(ji,jj) = -a_ddh_2 * ( 1.0 - zdh(ji,jj) / ( zari * zhbl(ji,jj) ) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj)
+                     
+                   ELSE  ! j_ddh == 0
+! Growing shear layer
+                     zddhdt(ji,jj) = -a_ddh * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj)
+                   ENDIF ! j_ddh
+                     zdhdt(ji,jj) = zdhdt(ji,jj) ! + zpsi * zddhdt(ji,jj)
+                ELSE    ! zdb_bl >0
+                   zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) /  MAX( zvel_max, 1.0e-15)
+                ENDIF
+             ELSE   ! zwb_min + 2*zwb_fk_b < 0
+                ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008)
+                zdhdt(ji,jj) = - zvel_mle(ji,jj)
+
+
+             ENDIF
+
+          ELSE
+             ! Fox-Kemper not used.
+
+               zvel_max = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) * zwb_ent(ji,jj) / &
+               &        MAX((zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird, epsln)
+               zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) )
+             ! added ajgn 23 July as temporay fix
+
+          ENDIF  ! ln_osm_mle
+
+         ELSE    ! lconv - Stable
+             zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj)
+             IF ( zdhdt(ji,jj) < 0._wp ) THEN
+                ! For long timsteps factor in brackets slows the rapid collapse of the OSBL
+                 zpert = 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_Dt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj)
+             ELSE
+                 zpert = MAX( 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) )
+             ENDIF
+             zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln)
+         ENDIF  ! lconv
+    ELSE ! lshear
+      IF ( lconv(ji,jj) ) THEN    ! Convective
+
+          IF ( ln_osm_mle ) THEN
+
+             IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN
+       ! Fox-Kemper buoyancy flux average over OSBL
+                zwb_fk_b(ji,jj) = zwb_fk(ji,jj) *  &
+                     (1.0 + hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * (-1.0 + ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3) )
+             ELSE
+                zwb_fk_b(ji,jj) = 0.5 * zwb_fk(ji,jj) * hmle(ji,jj) / hbl(ji,jj)
+             ENDIF
+             zvel_max = ( zwstrl(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj)
+             IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0.0 ) THEN
+                ! OSBL is deepening, entrainment > restratification
+                IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_bl_ext(ji,jj) > 0.0 ) THEN
+                   zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) )
+                ELSE
+                   zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) /  MAX( zvel_max, 1.0e-15)
+                ENDIF
+             ELSE
+                ! OSBL shoaling due to restratification flux. This is the velocity defined in Fox-Kemper et al (2008)
+                zdhdt(ji,jj) = - zvel_mle(ji,jj)
+
+
+             ENDIF
+
+          ELSE
+             ! Fox-Kemper not used.
+
+               zvel_max = -zwb_ent(ji,jj) / &
+               &        MAX((zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird, epsln)
+               zdhdt(ji,jj) = -zwb_ent(ji,jj) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) )
+             ! added ajgn 23 July as temporay fix
+
+          ENDIF  ! ln_osm_mle
+
+         ELSE                        ! Stable
+             zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj)
+             IF ( zdhdt(ji,jj) < 0._wp ) THEN
+                ! For long timsteps factor in brackets slows the rapid collapse of the OSBL
+                 zpert = 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj)
+             ELSE
+                 zpert = MAX( 2.0 * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) )
+             ENDIF
+             zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / MAX(zpert, epsln)
+         ENDIF  ! lconv
+      ENDIF ! lshear
+  END_2D
+    END SUBROUTINE zdf_osm_calculate_dhdt
+
+    SUBROUTINE zdf_osm_timestep_hbl( zdhdt )
+     !!---------------------------------------------------------------------
+     !!                   ***  ROUTINE zdf_osm_timestep_hbl  ***
+     !!
+     !! ** Purpose : Increments hbl.
+     !!
+     !! ** Method  : If thechange in hbl exceeds one model level the change is
+     !!              is calculated by moving down the grid, changing the buoyancy
+     !!              jump. This is to ensure that the change in hbl does not
+     !!              overshoot a stable layer.
+     !!
+     !!----------------------------------------------------------------------
+
+
+    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt   ! rates of change of hbl.
+
+    INTEGER :: jk, jj, ji, jm
+    REAL(wp) :: zhbl_s, zvel_max, zdb
+    REAL(wp) :: zthermal, zbeta
+
+     DO_2D( 0, 0, 0, 0 )
+        IF ( ibld(ji,jj) - imld(ji,jj) > 1 ) THEN
+!
+! If boundary layer changes by more than one level, need to check for stable layers between initial and final depths.
+!
+           zhbl_s = hbl(ji,jj)
+           jm = imld(ji,jj)
+           zthermal = rab_n(ji,jj,1,jp_tem)
+           zbeta = rab_n(ji,jj,1,jp_sal)
+
+
+           IF ( lconv(ji,jj) ) THEN
+!unstable
+
+              IF( ln_osm_mle ) THEN
+                 zvel_max = ( zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 )**p2third / hbl(ji,jj)
+              ELSE
+
+                 zvel_max = -( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_Dt / hbl(ji,jj) ) * zwb_ent(ji,jj) / &
+                   &      ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
+
+              ENDIF
+
+              DO jk = imld(ji,jj), ibld(ji,jj)
+                 zdb = MAX( grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) ) &
+                      & - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ), &
+                      &  0.0 ) + zvel_max
+
+
+                 IF ( ln_osm_mle ) THEN
+                    zhbl_s = zhbl_s + MIN( &
+                     & rn_Dt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), &
+                     & e3w(ji,jj,jm,Kmm) )
+                 ELSE
+                   zhbl_s = zhbl_s + MIN( &
+                     & rn_Dt * (  -zwb_ent(ji,jj) / zdb ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), &
+                     & e3w(ji,jj,jm,Kmm) )
+                 ENDIF
+
+!                    zhbl_s = MIN(zhbl_s,  gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol)
+                 IF ( zhbl_s >= gdepw(ji,jj,mbkt(ji,jj) + 1,Kmm) ) THEN
+                   zhbl_s = MIN(zhbl_s,  gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol)
+                   lpyc(ji,jj) = .FALSE.
+                 ENDIF
+                 IF ( zhbl_s >= gdepw(ji,jj,jm+1,Kmm) ) jm = jm + 1
+              END DO
+              hbl(ji,jj) = zhbl_s
+              ibld(ji,jj) = jm
+          ELSE
+! stable
+              DO jk = imld(ji,jj), ibld(ji,jj)
+                 zdb = MAX( &
+                      & grav * ( zthermal * ( zt_bl(ji,jj) - ts(ji,jj,jm,jp_tem,Kmm) )&
+                      &           - zbeta * ( zs_bl(ji,jj) - ts(ji,jj,jm,jp_sal,Kmm) ) ),&
+                      & 0.0 ) + &
+          &       2.0 * zvstr(ji,jj)**2 / zhbl_s
+
+                 ! Alan is thuis right? I have simply changed hbli to hbl
+                 zhol(ji,jj) = -zhbl_s / ( ( zvstr(ji,jj)**3 + epsln )/ zwbav(ji,jj) )
+                 zdhdt(ji,jj) = -( zwbav(ji,jj) - 0.04 / 2.0 * zwstrl(ji,jj)**3 / zhbl_s - 0.15 / 2.0 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * &
+            &                  zustar(ji,jj)**3 / zhbl_s ) * ( 0.725 + 0.225 * EXP( -7.5 * zhol(ji,jj) ) )
+                 zdhdt(ji,jj) = zdhdt(ji,jj) + zwbav(ji,jj)
+                 zhbl_s = zhbl_s + MIN( zdhdt(ji,jj) / zdb * rn_Dt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w(ji,jj,jm,Kmm) )
+
+!                    zhbl_s = MIN(zhbl_s, gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol)
+                 IF ( zhbl_s >= mbkt(ji,jj) + 1 ) THEN
+                   zhbl_s = MIN(zhbl_s,  gdepw(ji,jj, mbkt(ji,jj) + 1,Kmm) - depth_tol)
+                   lpyc(ji,jj) = .FALSE.
+                 ENDIF
+                 IF ( zhbl_s >= gdepw(ji,jj,jm,Kmm) ) jm = jm + 1
+              END DO
+          ENDIF   ! IF ( lconv )
+          hbl(ji,jj) = MAX(zhbl_s, gdepw(ji,jj,4,Kmm) )
+          ibld(ji,jj) = MAX(jm, 4 )
+        ELSE
+! change zero or one model level.
+          hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw(ji,jj,4,Kmm) )
+        ENDIF
+        zhbl(ji,jj) = gdepw(ji,jj,ibld(ji,jj),Kmm)
+     END_2D
+
+    END SUBROUTINE zdf_osm_timestep_hbl
+
+    SUBROUTINE zdf_osm_pycnocline_thickness( dh, zdh )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE zdf_osm_pycnocline_thickness  ***
+      !!
+      !! ** Purpose : Calculates thickness of the pycnocline
+      !!
+      !! ** Method  : The thickness is calculated from a prognostic equation
+      !!              that relaxes the pycnocine thickness to a diagnostic
+      !!              value. The time change is calculated assuming the
+      !!              thickness relaxes exponentially. This is done to deal
+      !!              with large timesteps.
+      !!
+      !!----------------------------------------------------------------------
+
+      REAL(wp), DIMENSION(jpi,jpj) :: dh, zdh     ! pycnocline thickness.
+       !
+      INTEGER :: jj, ji
+      INTEGER :: inhml
+      REAL(wp) :: zari, ztau, zdh_ref
+      REAL, PARAMETER :: a_ddh_2 = 3.5 ! also in pycnocline_depth
+
+    DO_2D( 0, 0, 0, 0 )
+
+      IF ( lshear(ji,jj) ) THEN
+         IF ( lconv(ji,jj) ) THEN
+           IF ( j_ddh(ji,jj) == 0 ) THEN
+! ddhdt for pycnocline determined in osm_calculate_dhdt
+             dh(ji,jj) = dh(ji,jj) + zddhdt(ji,jj) * rn_Dt
+           ELSE
+! Temporary (probably) Recalculate dh_ref to ensure dh doesn't go negative. Can't do this using zddhdt from calculate_dhdt 
+             IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN
+               zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )
+             ELSE
+               zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )
+             ENDIF
+             ztau = MAX( zdb_bl(ji,jj) * ( zari * hbl(ji,jj) ) / ( a_ddh_2 * MAX(-zwb_ent(ji,jj), 1.e-12) ), 2.0 * rn_Dt )
+             dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau ) + zari * zhbl(ji,jj) * ( 1.0 - EXP( -rn_Dt / ztau ) )
+             IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zari * zhbl(ji,jj)
+           ENDIF
+            
+         ELSE ! lconv
+! Initially shear only for entraining OSBL. Stable code will be needed if extended to stable OSBL 
+
+            ztau = hbl(ji,jj) / MAX(zvstr(ji,jj), epsln)
+            IF ( zdhdt(ji,jj) >= 0.0 ) THEN    ! probably shouldn't include wm here
+               ! boundary layer deepening
+               IF ( zdb_bl(ji,jj) > 0._wp ) THEN
+                  ! pycnocline thickness set by stratification - use same relationship as for neutral conditions.
+                  zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) &
+                       & /  MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01  , 0.2 )
+                  zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj)
+               ELSE
+                  zdh_ref = 0.2 * hbl(ji,jj)
+               ENDIF
+            ELSE     ! IF(dhdt < 0)
+               zdh_ref = 0.2 * hbl(ji,jj)
+            ENDIF    ! IF (dhdt >= 0)
+            dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau )+ zdh_ref * ( 1.0 - EXP( -rn_Dt / ztau ) )
+            IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref       ! can be a problem with dh>hbl for rapid collapse
+            ! Alan: this hml is never defined or used -- do we need it?
+         ENDIF
+          
+      ELSE   ! lshear  
+! for lshear = .FALSE. calculate ddhdt here
+
+          IF ( lconv(ji,jj) ) THEN
+
+            IF( ln_osm_mle ) THEN
+               IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) < 0._wp ) THEN
+                  ! OSBL is deepening. Note wb_fk_b is zero if ln_osm_mle=F
+                  IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN
+                     IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN  ! near neutral stability
+                        zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )
+                     ELSE                                                     ! unstable
+                        zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )
+                     ENDIF
+                     ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird)
+                     zdh_ref = zari * hbl(ji,jj)
+                  ELSE
+                     ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird)
+                     zdh_ref = 0.2 * hbl(ji,jj)
+                  ENDIF
+               ELSE
+                  ztau = 0.2 * hbl(ji,jj) /  MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird)
+                  zdh_ref = 0.2 * hbl(ji,jj)
+               ENDIF
+            ELSE ! ln_osm_mle
+               IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_bl_ext(ji,jj) > 0._wp)THEN
+                  IF ( ( zwstrc(ji,jj) / MAX(zvstr(ji,jj), epsln) )**3 <= 0.5 ) THEN  ! near neutral stability
+                     zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )
+                  ELSE                                                     ! unstable
+                     zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )
+                  ENDIF
+                  ztau = hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird)
+                  zdh_ref = zari * hbl(ji,jj)
+               ELSE
+                  ztau = hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird)
+                  zdh_ref = 0.2 * hbl(ji,jj)
+               ENDIF
+
+            END IF  ! ln_osm_mle
+
+            dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau ) + zdh_ref * ( 1.0 - EXP( -rn_Dt / ztau ) )
+!               IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref
+            IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref
+            ! Alan: this hml is never defined or used
+         ELSE   ! IF (lconv)
+            ztau = hbl(ji,jj) / MAX(zvstr(ji,jj), epsln)
+            IF ( zdhdt(ji,jj) >= 0.0 ) THEN    ! probably shouldn't include wm here
+               ! boundary layer deepening
+               IF ( zdb_bl(ji,jj) > 0._wp ) THEN
+                  ! pycnocline thickness set by stratification - use same relationship as for neutral conditions.
+                  zari = MIN( 4.5 * ( zvstr(ji,jj)**2 ) &
+                       & /  MAX(zdb_bl(ji,jj) * zhbl(ji,jj), epsln ) + 0.01  , 0.2 )
+                  zdh_ref = MIN( zari, 0.2 ) * hbl(ji,jj)
+               ELSE
+                  zdh_ref = 0.2 * hbl(ji,jj)
+               ENDIF
+            ELSE     ! IF(dhdt < 0)
+               zdh_ref = 0.2 * hbl(ji,jj)
+            ENDIF    ! IF (dhdt >= 0)
+            dh(ji,jj) = dh(ji,jj) * EXP( -rn_Dt / ztau )+ zdh_ref * ( 1.0 - EXP( -rn_Dt / ztau ) )
+            IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zdh_ref       ! can be a problem with dh>hbl for rapid collapse
+         ENDIF       ! IF (lconv)
+      ENDIF  ! lshear
+ 
+      hml(ji,jj) = hbl(ji,jj) - dh(ji,jj)
+      inhml = MAX( INT( dh(ji,jj) / MAX(e3t(ji,jj,ibld(ji,jj),Kmm), 1.e-3) ) , 1 )
+      imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3)
+      zhml(ji,jj) = gdepw(ji,jj,imld(ji,jj),Kmm)
+      zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
+    END_2D
+
+    END SUBROUTINE zdf_osm_pycnocline_thickness
+
+
+   SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_osm_horizontal_gradients  ***
+      !!
+      !! ** Purpose :   Calculates horizontal gradients of buoyancy for use with Fox-Kemper parametrization.
+      !!
+      !! ** Method  :
+      !!
+      !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008
+      !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008
+
+
+      REAL(wp), DIMENSION(jpi,jpj)     :: dbdx_mle, dbdy_mle ! MLE horiz gradients at u & v points
+      REAL(wp), DIMENSION(jpi,jpj)     :: zdbds_mle ! Magnitude of horizontal buoyancy gradient.
+      REAL(wp), DIMENSION(jpi,jpj)     :: zmld ! ==  estimated FK BLD used for MLE horiz gradients  == !
+      REAL(wp), DIMENSION(jpi,jpj)     :: zdtdx, zdtdy, zdsdx, zdsdy
+
+      INTEGER  ::   ji, jj, jk          ! dummy loop indices
+      INTEGER  ::   ii, ij, ik, ikmax   ! local integers
+      REAL(wp)                         :: zc
+      REAL(wp)                         :: zN2_c           ! local buoyancy difference from 10m value
+      REAL(wp), DIMENSION(jpi,jpj)     :: ztm, zsm, zLf_NH, zLf_MH
+      REAL(wp), DIMENSION(jpi,jpj,jpts):: ztsm_midu, ztsm_midv, zabu, zabv
+      REAL(wp), DIMENSION(jpi,jpj)     :: zmld_midu, zmld_midv
+!!----------------------------------------------------------------------
+      !
+      !                                      !==  MLD used for MLE  ==!
+
+      mld_prof(:,:)  = nlb10               ! Initialization to the number of w ocean point
+      zmld(:,:)  = 0._wp               ! here hmlp used as a dummy variable, integrating vertically N^2
+      zN2_c = grav * rn_osm_mle_rho_c * r1_rho0   ! convert density criteria into N^2 criteria
+      DO_3D( 1, 1, 1, 1, nlb10, jpkm1 )
+         ikt = mbkt(ji,jj)
+         zmld(ji,jj) = zmld(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm)
+         IF( zmld(ji,jj) < zN2_c )   mld_prof(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level
+      END_3D
+      DO_2D( 1, 1, 1, 1 )
+         mld_prof(ji,jj) = MAX(mld_prof(ji,jj),ibld(ji,jj))
+         zmld(ji,jj) = gdepw(ji,jj,mld_prof(ji,jj),Kmm)
+      END_2D
+      ! ensure mld_prof .ge. ibld
+      !
+      ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 )                  ! max level of the computation
+      !
+      ztm(:,:) = 0._wp
+      zsm(:,:) = 0._wp
+      DO_3D( 1, 1, 1, 1, 1, ikmax )
+         zc = e3t(ji,jj,jk,Kmm) * REAL( MIN( MAX( 0, mld_prof(ji,jj)-jk ) , 1  )  )    ! 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), zmld(:,:) )
+      zsm(:,:) = zsm(:,:) / MAX( e3t(:,:,1,Kmm), zmld(:,:) )
+      ! calculate horizontal gradients at u & v points
+
+      DO_2D( 0, 0, 1, 0 )
+         zdtdx(ji,jj) = ( ztm(ji+1,jj) - ztm( ji,jj) )  * umask(ji,jj,1) / e1u(ji,jj)
+         zdsdx(ji,jj) = ( zsm(ji+1,jj) - zsm( ji,jj) )  * umask(ji,jj,1) / e1u(ji,jj)
+         zmld_midu(ji,jj) = 0.25_wp * (zmld(ji+1,jj) + zmld( ji,jj))
+         ztsm_midu(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji+1,jj) + ztm( ji,jj) )
+         ztsm_midu(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji+1,jj) + zsm( ji,jj) )
+      END_2D
+
+      DO_2D( 1, 0, 0, 0 )
+         zdtdy(ji,jj) = ( ztm(ji,jj+1) - ztm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj)
+         zdsdy(ji,jj) = ( zsm(ji,jj+1) - zsm( ji,jj) ) * vmask(ji,jj,1) / e1v(ji,jj)
+         zmld_midv(ji,jj) = 0.25_wp * (zmld(ji,jj+1) + zmld( ji,jj))
+         ztsm_midv(ji,jj,jp_tem) = 0.5_wp * ( ztm(ji,jj+1) + ztm( ji,jj) )
+         ztsm_midv(ji,jj,jp_sal) = 0.5_wp * ( zsm(ji,jj+1) + zsm( ji,jj) )
+      END_2D
+
+      CALL eos_rab(ztsm_midu, zmld_midu, zabu, Kmm)
+      CALL eos_rab(ztsm_midv, zmld_midv, zabv, Kmm)
+
+      DO_2D( 0, 0, 1, 0 )
+         dbdx_mle(ji,jj) = grav*(zdtdx(ji,jj)*zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal))
+      END_2D
+      DO_2D( 1, 0, 0, 0 )
+         dbdy_mle(ji,jj) = grav*(zdtdy(ji,jj)*zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal))
+      END_2D
+
+      DO_2D( 0, 0, 0, 0 )
+        ztmp =  r1_ft(ji,jj) *  MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf
+        zdbds_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
+      
+ END SUBROUTINE zdf_osm_zmld_horizontal_gradients
+  SUBROUTINE zdf_osm_mle_parameters( mld_prof, hmle, zhmle, zvel_mle, zdiff_mle )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_osm_mle_parameters  ***
+      !!
+      !! ** Purpose :   Timesteps the mixed layer eddy depth, hmle and calculates the mixed layer eddy fluxes for buoyancy, heat and salinity.
+      !!
+      !! ** Method  :
+      !!
+      !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008
+      !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008
+
+      INTEGER, DIMENSION(jpi,jpj)      :: mld_prof
+      REAL(wp), DIMENSION(jpi,jpj)     :: hmle, zhmle, zwb_fk, zvel_mle, zdiff_mle
+      INTEGER  ::   ji, jj, jk          ! dummy loop indices
+      INTEGER  ::   ii, ij, ik, jkb, jkb1  ! local integers
+      INTEGER , DIMENSION(jpi,jpj)     :: inml_mle
+      REAL(wp) ::  ztmp, zdbdz, zdtdz, zdsdz, zthermal,zbeta, zbuoy, zdb_mle
+
+   ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE.
+
+      DO_2D( 0, 0, 0, 0 )
+       IF ( lconv(ji,jj) ) THEN
+          ztmp =  r1_ft(ji,jj) *  MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf
+   ! This velocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt.
+          zvel_mle(ji,jj) = zdbds_mle(ji,jj) * ztmp * hmle(ji,jj) * tmask(ji,jj,1)
+          zdiff_mle(ji,jj) = 5.e-4_wp * rn_osm_mle_ce * ztmp * zdbds_mle(ji,jj) * zhmle(ji,jj)**2
+       ENDIF
+      END_2D
+   ! Timestep mixed layer eddy depth.
+      DO_2D( 0, 0, 0, 0 )
+        IF ( lmle(ji,jj) ) THEN  ! MLE layer growing.
+! Buoyancy gradient at base of MLE layer.
+           zthermal = rab_n(ji,jj,1,jp_tem)
+           zbeta    = rab_n(ji,jj,1,jp_sal)
+           jkb = mld_prof(ji,jj)
+           jkb1 = MIN(jkb + 1, mbkt(ji,jj))
+!              
+           zbuoy = grav * ( zthermal * ts(ji,jj,mld_prof(ji,jj)+2,jp_tem,Kmm) - zbeta * ts(ji,jj,mld_prof(ji,jj)+2,jp_sal,Kmm) )
+           zdb_mle = zb_bl(ji,jj) - zbuoy 
+! Timestep hmle. 
+           hmle(ji,jj) = hmle(ji,jj) + zwb0(ji,jj) * rn_Dt / zdb_mle
+        ELSE
+           IF ( zhmle(ji,jj) > zhbl(ji,jj) ) THEN
+              hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - hbl(ji,jj) ) * rn_Dt / rn_osm_mle_tau
+           ELSE
+              hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - hbl(ji,jj) ) * rn_Dt /rn_osm_mle_tau
+           ENDIF
+        ENDIF
+        hmle(ji,jj) = MIN(hmle(ji,jj), ht(ji,jj))
+       IF(ln_osm_hmle_limit) hmle(ji,jj) = MIN(hmle(ji,jj), MAX(rn_osm_hmle_limit,1.2*hbl(ji,jj)) )
+      END_2D
+
+      mld_prof = 4
+      DO_3D( 0, 0, 0, 0, 5, jpkm1 )
+      IF ( hmle(ji,jj) >= gdepw(ji,jj,jk,Kmm) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk)
+      END_3D
+      DO_2D( 0, 0, 0, 0 )
+         zhmle(ji,jj) = gdepw(ji,jj, mld_prof(ji,jj),Kmm)
+      END_2D
+END SUBROUTINE zdf_osm_mle_parameters
+
+END SUBROUTINE zdf_osm
+
+
+   SUBROUTINE zdf_osm_init( Kmm )
      !!----------------------------------------------------------------------
      !!                  ***  ROUTINE zdf_osm_init  ***
@@ -1304,12 +2456,17 @@
      !! ** input   :   Namlist namosm
      !!----------------------------------------------------------------------
-     INTEGER, INTENT(in)    :: Kmm ! time level index (middle)
-     !
+     INTEGER, INTENT(in)   ::   Kmm       ! time level
      INTEGER  ::   ios            ! local integer
      INTEGER  ::   ji, jj, jk     ! dummy loop indices
+     REAL z1_t2
      !!
      NAMELIST/namzdf_osm/ ln_use_osm_la, rn_osm_la, rn_osm_dstokes, nn_ave &
-          & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0 &
-          & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv
+          & ,nn_osm_wave, ln_dia_osm, rn_osm_hbl0, rn_zdfosm_adjust_sd &
+          & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave &
+          & ,nn_osm_SD_reduce, ln_osm_mle, rn_osm_hblfrac, rn_osm_bl_thresh, ln_zdfosm_ice_shelter
+! Namelist for Fox-Kemper parametrization.
+      NAMELIST/namosm_mle/ nn_osm_mle, rn_osm_mle_ce, rn_osm_mle_lf, rn_osm_mle_time, rn_osm_mle_lat,&
+           & rn_osm_mle_rho_c, rn_osm_mle_thresh, rn_osm_mle_tau, ln_osm_hmle_limit, rn_osm_hmle_limit
+
      !!----------------------------------------------------------------------
      !
@@ -1325,9 +2482,11 @@
         WRITE(numout,*) 'zdf_osm_init : OSMOSIS Parameterisation'
         WRITE(numout,*) '~~~~~~~~~~~~'
-        WRITE(numout,*) '   Namelist namzdf_osm : set tke mixing parameters'
-        WRITE(numout,*) '     Use namelist  rn_osm_la                     ln_use_osm_la = ', ln_use_osm_la
+        WRITE(numout,*) '   Namelist namzdf_osm : set osm mixing parameters'
+        WRITE(numout,*) '     Use  rn_osm_la                                ln_use_osm_la = ', ln_use_osm_la
+        WRITE(numout,*) '     Use  MLE in OBL, i.e. Fox-Kemper param        ln_osm_mle = ', ln_osm_mle
         WRITE(numout,*) '     Turbulent Langmuir number                     rn_osm_la   = ', rn_osm_la
+        WRITE(numout,*) '     Stokes drift reduction factor                 rn_zdfosm_adjust_sd   = ', rn_zdfosm_adjust_sd
         WRITE(numout,*) '     Initial hbl for 1D runs                       rn_osm_hbl0   = ', rn_osm_hbl0
-        WRITE(numout,*) '     Depth scale of Stokes drift                rn_osm_dstokes = ', rn_osm_dstokes
+        WRITE(numout,*) '     Depth scale of Stokes drift                   rn_osm_dstokes = ', rn_osm_dstokes
         WRITE(numout,*) '     horizontal average flag                       nn_ave      = ', nn_ave
         WRITE(numout,*) '     Stokes drift                                  nn_osm_wave = ', nn_osm_wave
@@ -1339,6 +2498,19 @@
         CASE(2)
            WRITE(numout,*) '     calculated from ECMWF wave fields'
+         END SELECT
+        WRITE(numout,*) '     Stokes drift reduction                        nn_osm_SD_reduce', nn_osm_SD_reduce
+        WRITE(numout,*) '     fraction of hbl to average SD over/fit'
+        WRITE(numout,*) '     exponential with nn_osm_SD_reduce = 1 or 2    rn_osm_hblfrac =  ', rn_osm_hblfrac
+        SELECT CASE (nn_osm_SD_reduce)
+        CASE(0)
+           WRITE(numout,*) '     No reduction'
+        CASE(1)
+           WRITE(numout,*) '     Average SD over upper rn_osm_hblfrac of BL'
+        CASE(2)
+           WRITE(numout,*) '     Fit exponential to slope rn_osm_hblfrac of BL'
         END SELECT
+        WRITE(numout,*) '     reduce surface SD and depth scale under ice   ln_zdfosm_ice_shelter=', ln_zdfosm_ice_shelter
         WRITE(numout,*) '     Output osm diagnostics                       ln_dia_osm  = ',  ln_dia_osm
+        WRITE(numout,*) '         Threshold used to define BL              rn_osm_bl_thresh  = ', rn_osm_bl_thresh, 'm^2/s'
         WRITE(numout,*) '     Use KPP-style shear instability mixing       ln_kpprimix = ', ln_kpprimix
         WRITE(numout,*) '     local Richardson Number limit for shear instability rn_riinfty = ', rn_riinfty
@@ -1359,5 +2531,68 @@
      IF( zdf_osm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' )
 
-     call osm_rst( nit000, Kmm, 'READ' ) !* read or initialize hbl
+
+     IF( ln_osm_mle ) THEN
+! Initialise Fox-Kemper parametrization
+         READ  ( numnam_ref, namosm_mle, IOSTAT = ios, ERR = 903)
+903      IF( ios /= 0 )   CALL ctl_nam ( ios , 'namosm_mle in reference namelist')
+
+         READ  ( numnam_cfg, namosm_mle, IOSTAT = ios, ERR = 904 )
+904      IF( ios >  0 )   CALL ctl_nam ( ios , 'namosm_mle in configuration namelist')
+         IF(lwm) WRITE ( numond, namosm_mle )
+
+         IF(lwp) THEN                     ! Namelist print
+            WRITE(numout,*)
+            WRITE(numout,*) 'zdf_osm_init : initialise mixed layer eddy (MLE)'
+            WRITE(numout,*) '~~~~~~~~~~~~~'
+            WRITE(numout,*) '   Namelist namosm_mle : '
+            WRITE(numout,*) '         MLE type: =0 standard Fox-Kemper ; =1 new formulation        nn_osm_mle    = ', nn_osm_mle
+            WRITE(numout,*) '         magnitude of the MLE (typical value: 0.06 to 0.08)           rn_osm_mle_ce    = ', rn_osm_mle_ce
+            WRITE(numout,*) '         scale of ML front (ML radius of deformation) (nn_osm_mle=0)      rn_osm_mle_lf     = ', rn_osm_mle_lf, 'm'
+            WRITE(numout,*) '         maximum time scale of MLE                    (nn_osm_mle=0)      rn_osm_mle_time   = ', rn_osm_mle_time, 's'
+            WRITE(numout,*) '         reference latitude (degrees) of MLE coef.    (nn_osm_mle=1)      rn_osm_mle_lat    = ', rn_osm_mle_lat, 'deg'
+            WRITE(numout,*) '         Density difference used to define ML for FK              rn_osm_mle_rho_c  = ', rn_osm_mle_rho_c
+            WRITE(numout,*) '         Threshold used to define MLE for FK                      rn_osm_mle_thresh  = ', rn_osm_mle_thresh, 'm^2/s'
+            WRITE(numout,*) '         Timescale for OSM-FK                                         rn_osm_mle_tau  = ', rn_osm_mle_tau, 's'
+            WRITE(numout,*) '         switch to limit hmle                                      ln_osm_hmle_limit  = ', ln_osm_hmle_limit
+            WRITE(numout,*) '         fraction of zmld to limit hmle to if ln_osm_hmle_limit =.T.  rn_osm_hmle_limit = ', rn_osm_hmle_limit
+         ENDIF         !
+     ENDIF
+      !
+      IF(lwp) THEN
+         WRITE(numout,*)
+         IF( ln_osm_mle ) THEN
+            WRITE(numout,*) '   ==>>>   Mixed Layer Eddy induced transport added to OSMOSIS BL calculation'
+            IF( nn_osm_mle == 0 )   WRITE(numout,*) '              Fox-Kemper et al 2010 formulation'
+            IF( nn_osm_mle == 1 )   WRITE(numout,*) '              New formulation'
+         ELSE
+            WRITE(numout,*) '   ==>>>   Mixed Layer induced transport NOT added to OSMOSIS BL calculation'
+         ENDIF
+      ENDIF
+      !
+      IF( ln_osm_mle ) THEN                ! MLE initialisation
+         !
+         rb_c = grav * rn_osm_mle_rho_c /rho0        ! Mixed Layer buoyancy criteria
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '      ML buoyancy criteria = ', rb_c, ' m/s2 '
+         IF(lwp) WRITE(numout,*) '      associated ML density criteria defined in zdfmxl = ', rn_osm_mle_rho_c, 'kg/m3'
+         !
+         IF( nn_osm_mle == 0 ) THEN           ! MLE array allocation & initialisation            !
+!
+         ELSEIF( nn_osm_mle == 1 ) THEN           ! MLE array allocation & initialisation
+            rc_f = rn_osm_mle_ce/ (  5.e3_wp * 2._wp * omega * SIN( rad * rn_osm_mle_lat )  )
+            !
+         ENDIF
+         !                                ! 1/(f^2+tau^2)^1/2 at t-point (needed in both nn_osm_mle case)
+         z1_t2 = 2.e-5
+         DO_2D( 1, 1, 1, 1 )
+            r1_ft(ji,jj) = MIN(1./( ABS(ff_t(ji,jj)) + epsln ), ABS(ff_t(ji,jj))/z1_t2**2)
+         END_2D
+         ! z1_t2 = 1._wp / ( rn_osm_mle_time * rn_osm_mle_timeji,jj )
+         ! r1_ft(:,:) = 1._wp / SQRT(  ff_t(:,:) * ff_t(:,:) + z1_t2  )
+         !
+      ENDIF
+
+     call osm_rst( nit000, Kmm,  'READ' ) !* read or initialize hbl, dh, hmle
+
 
      IF( ln_zdfddm) THEN
@@ -1440,5 +2675,8 @@
         CALL iom_set_rstw_var_active('wn')
         CALL iom_set_rstw_var_active('hbl')
-        CALL iom_set_rstw_var_active('hbli')
+        CALL iom_set_rstw_var_active('dh')
+        IF( ln_osm_mle ) THEN
+            CALL iom_set_rstw_var_active('hmle')
+        END IF
      ENDIF
    END SUBROUTINE zdf_osm_init
@@ -1459,5 +2697,5 @@
      CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
 
-     INTEGER ::   id1, id2   ! iom enquiry index
+     INTEGER ::   id1, id2, id3   ! iom enquiry index
      INTEGER  ::   ji, jj, jk     ! dummy loop indices
      INTEGER  ::   iiki, ikt ! local integer
@@ -1465,5 +2703,5 @@
      REAL(wp) ::   zN2_c           ! local scalar
      REAL(wp) ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth
-     INTEGER, DIMENSION(:,:), ALLOCATABLE :: imld_rst ! level of mixed-layer depth (pycnocline top)
+     INTEGER, DIMENSION(jpi,jpj) :: imld_rst ! level of mixed-layer depth (pycnocline top)
      !!----------------------------------------------------------------------
      !
@@ -1480,12 +2718,23 @@
            WRITE(numout,*) ' ===>>>> :  ww not in restart file, set to zero initially'
         END IF
+
         id1 = iom_varid( numror, 'hbl'   , ldstop = .FALSE. )
-        id2 = iom_varid( numror, 'hbli'   , ldstop = .FALSE. )
+        id2 = iom_varid( numror, 'dh'   , ldstop = .FALSE. )
         IF( id1 > 0 .AND. id2 > 0) THEN                       ! 'hbl' exists; read and return
            CALL iom_get( numror, jpdom_auto, 'hbl' , hbl , ldxios = lrxios )
-           CALL iom_get( numror, jpdom_auto, 'hbli', hbli, ldxios = lrxios  )
-           WRITE(numout,*) ' ===>>>> :  hbl & hbli read from restart file'
+           CALL iom_get( numror, jpdom_auto, 'dh', dh, ldxios = lrxios  )
+           WRITE(numout,*) ' ===>>>> :  hbl & dh read from restart file'
+           IF( ln_osm_mle ) THEN
+              id3 = iom_varid( numror, 'hmle'   , ldstop = .FALSE. )
+              IF( id3 > 0) THEN
+                 CALL iom_get( numror, jpdom_auto, 'hmle' , hmle , ldxios = lrxios )
+                 WRITE(numout,*) ' ===>>>> :  hmle read from restart file'
+              ELSE
+                 WRITE(numout,*) ' ===>>>> :  hmle not found, set to hbl'
+                 hmle(:,:) = hbl(:,:)            ! Initialise MLE depth.
+              END IF
+           END IF
            RETURN
-        ELSE                      ! 'hbl' & 'hbli' not in restart file, recalculate
+        ELSE                      ! 'hbl' & 'dh' not in restart file, recalculate
            WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification'
         END IF
@@ -1495,9 +2744,12 @@
      ! If READ/WRITE Flag is 'WRITE', write hbl into the restart file, then return
      !!-----------------------------------------------------------------------------
-     IF( TRIM(cdrw) == 'WRITE') THEN     !* Write hbli into the restart file, then return
+     IF( TRIM(cdrw) == 'WRITE') THEN     !* Write hbl into the restart file, then return
         IF(lwp) WRITE(numout,*) '---- osm-rst ----'
-         CALL iom_rstput( kt, nitrst, numrow, 'wn'     , ww  , ldxios = lwxios )
-         CALL iom_rstput( kt, nitrst, numrow, 'hbl'    , hbl , ldxios = lwxios )
-         CALL iom_rstput( kt, nitrst, numrow, 'hbli'   , hbli, ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'ww'     , ww,   ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'hbl'    , hbl,  ldxios = lwxios )
+         CALL iom_rstput( kt, nitrst, numrow, 'dh'     , dh,   ldxios = lwxios )
+         IF( ln_osm_mle ) THEN
+            CALL iom_rstput( kt, nitrst, numrow, 'hmle', hmle, ldxios = lwxios )
+         END IF
         RETURN
      END IF
@@ -1507,5 +2759,4 @@
      !!-----------------------------------------------------------------------------
      IF( lwp ) WRITE(numout,*) ' ===>>>> : calculating hbl computed from stratification'
-     ALLOCATE( imld_rst(jpi,jpj) )
      ! w-level of the mixing and mixed layers
      CALL eos_rab( ts(:,:,:,:,Kmm), rab_n, Kmm )
@@ -1516,5 +2767,5 @@
      !
      hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2
-     DO_3D( 1, 1, 1, 1, 1, jpkm1 )  ! Mixed layer level: w-level
+     DO_3D( 1, 1, 1, 1, 1, jpkm1 )
         ikt = mbkt(ji,jj)
         hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w(ji,jj,jk,Kmm)
@@ -1523,11 +2774,18 @@
      !
      DO_2D( 1, 1, 1, 1 )
-        iiki = imld_rst(ji,jj)
-        hbl (ji,jj) = gdepw(ji,jj,iiki  ,Kmm) * ssmask(ji,jj)    ! Turbocline depth
+        iiki = MAX(4,imld_rst(ji,jj))
+        hbl (ji,jj) = gdepw(ji,jj,iiki,Kmm  )    ! Turbocline depth
+        dh (ji,jj) = e3t(ji,jj,iiki-1,Kmm  )     ! Turbocline depth
      END_2D
-     hbl = MAX(hbl,epsln)
-     hbli(:,:) = hbl(:,:)
-     DEALLOCATE( imld_rst )
+
      WRITE(numout,*) ' ===>>>> : hbl computed from stratification'
+
+     IF( ln_osm_mle ) THEN
+        hmle(:,:) = hbl(:,:)            ! Initialise MLE depth.
+        WRITE(numout,*) ' ===>>>> : hmle set = to hbl'
+     END IF
+
+     ww(:,:,:) = 0._wp
+     WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially'
    END SUBROUTINE osm_rst
 
@@ -1560,5 +2818,4 @@
       ENDIF
 
-      ! add non-local temperature and salinity flux
       DO_3D( 0, 0, 0, 0, 1, jpkm1 )
          pts(ji,jj,jk,jp_tem,Krhs) =  pts(ji,jj,jk,jp_tem,Krhs)                      &
@@ -1570,12 +2827,11 @@
       END_3D
 
-
-      ! save the non-local tracer flux trends for diagnostic
+      ! save the non-local tracer flux trends for diagnostics
       IF( l_trdtra )   THEN
          ztrdt(:,:,:) = pts(:,:,:,jp_tem,Krhs) - ztrdt(:,:,:)
          ztrds(:,:,:) = pts(:,:,:,jp_sal,Krhs) - ztrds(:,:,:)
-!!bug gm jpttdzdf ==> jpttosm
-         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_zdf, ztrdt )
-         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_zdf, ztrds )
+
+         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_tem, jptra_osm, ztrdt )
+         CALL trd_tra( kt, Kmm, Krhs, 'TRA', jp_sal, jptra_osm, ztrds )
          DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds )
       ENDIF
@@ -1643,3 +2899,4 @@
 
    !!======================================================================
+
 END MODULE zdfosm
