Index: /NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/ZDF/zdfosm.F90
===================================================================
--- /NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/ZDF/zdfosm.F90	(revision 12217)
+++ /NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/ZDF/zdfosm.F90	(revision 12218)
@@ -76,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
@@ -84,9 +86,19 @@
    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
@@ -102,4 +114,20 @@
    LOGICAL  ::   ln_convmix  = .true.   ! Convective instability mixing
    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
+   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 /rau0 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_mle_tau             ! Adjustment timescale for MLE.
+
 
    !                                    !!! ** General constants  **
@@ -122,7 +150,20 @@
       !!                 ***  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 )
+
+!     ALLOCATE( ghamu(jpi,jpj,jpk), ghamv(jpi,jpj,jpk), ghamt(jpi,jpj,jpk),ghams(jpi,jpj,jpk), &    ! would ths be better ?
+!          &       hbl(jpi,jpj), dh(jpi,jpj), hml(jpi,jpj), dstokes(jpi, jpj), &
+!          &       etmean(jpi,jpj,jpk), STAT= zdf_osm_alloc )
+!     IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays')
+!     IF ( ln_osm_mle ) THEN
+!        Allocate( hmle(jpi,jpj), r1_ft(jpi,jpj), STAT= zdf_osm_alloc )
+!        IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm mle arrays')
+!     ENDIF
+
      IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays')
      CALL mpp_sum ( 'zdfosm', zdf_osm_alloc )
@@ -169,4 +210,7 @@
       !!
       INTEGER ::   ji, jj, jk                   ! dummy loop indices
+
+      INTEGER ::   jl                   ! dummy loop indices
+
       INTEGER ::   ikbot, jkmax, jkm1, jkp2     !
 
@@ -200,4 +244,11 @@
       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_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
@@ -217,6 +268,15 @@
       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) :: zdhdt_2                                    ! correction to dhdt due to internal structure.
+      REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_ext,zdsdz_ext,zdbdz_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,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
@@ -289,4 +349,14 @@
       zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp
       !
+      zdtdz_ext(:,:) = 0._wp ; zdsdz_ext(:,:) = 0._wp ; zdbdz_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
@@ -432,12 +502,30 @@
          END DO
       END DO
-
+      ! Averages over well-mixed and boundary layer
+      ! Alan: do we need zb_nl?, zb_ml?
+      CALL zdf_osm_vertical_average(ibld, zt_bl, zs_bl, zu_bl, zv_bl, zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl)
+      CALL zdf_osm_vertical_average(imld, zt_ml, zs_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml)
+! External gradient
+      CALL zdf_osm_external_gradients( zdtdz_ext, zdsdz_ext, zdbdz_ext )
+
+
+      IF ( ln_osm_mle ) THEN
+         CALL zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle )
+         CALL zdf_osm_mle_parameters( hmle, zwb_fk, zvel_mle, zdiff_mle )
+       ENDIF
+
+! Rate of change of hbl
+      CALL zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 )
       ! Calculate averages over depth of boundary layer
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               zhbl_t(ji,jj) = hbl(ji,jj) + (zdhdt(ji,jj) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wn here, so subtract it
+               ! adjustment to represent limiting by ocean bottom
+               zhbl_t(ji,jj) = MIN(zhbl_t(ji,jj), gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)! ht_n(:,:))
+            END DO
+         END DO
+
       imld = ibld           ! use imld to hold previous blayer index
-      ibld(:,:) = 3
-
-      zhbl_t(:,:) = hbl(:,:) + (zdhdt(:,:) - wn(ji,jj,ibld(ji,jj)))* rn_rdt ! certainly need wb here, so subtract it
-      zhbl_t(:,:) = MIN(zhbl_t(:,:), ht_n(:,:))
-      zdhdt(:,:) = MIN(zdhdt(:,:), (zhbl_t(:,:) - hbl(:,:))/rn_rdt + wn(ji,jj,ibld(ji,jj))) ! adjustment to represent limiting by ocean bottom
+      ibld(:,:) = 4
 
       DO jk = 4, jpkm1
@@ -550,6 +638,7 @@
                    ENDIF
                 ENDIF
-                ! Temporay fix to ensure zdiffut is +ve; won't be necessary with wn taken out
-                zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj)* e3t_n(ji,jj,ibld(ji,jj))
+                ! Temporary fix to ensure zdiffut is +ve; won't be necessary with wn taken out
+                zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)), 1.e-6)
+                zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj) * e3t_n(ji,jj,ibld(ji,jj)), 1.e-5)
                 ! could be taken out, take account of entrainment represents as a diffusivity
                 ! should remove w from here, represents entrainment
@@ -919,8 +1008,4 @@
        END DO
 
-       IF(ln_dia_osm) THEN
-          IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc )
-       END IF
-
 ! KPP-style Ri# mixing
        IF( ln_kpprimix) THEN
@@ -974,4 +1059,53 @@
           END DO
        END IF ! ln_convmix = .true.
+
+
+ 
+       IF ( ln_osm_mle ) THEN  ! set up diffusivity and non-gradient mixing
+          DO jj = 2 , jpjm1
+             DO ji = 2, jpim1
+                IF ( lconv(ji,jj) .AND. mld_prof(ji,jj) - ibld(ji,jj) > 1 ) THEN ! MLE mmixing extends below the OSBL.
+            ! Calculate MLE flux profiles
+                   ! DO jk = 1, mld_prof(ji,jj)
+                   !    znd = - gdepw_n(ji,jj,jk) / MAX(zhmle(ji,jj),epsln)
+                   !    ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + &
+                   !         & zwt_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 )
+                   !    ghams(ji,jj,jk) = ghams(ji,jj,jk) + &
+                   !         & zws_fk(ji,jj) * ( 1.0 - ( 2.0 * znd + 1.0 )**2 ) * ( 1.0 + r5_21 * ( 2.0 * znd + 1.0 )**2 )
+                   ! END DO
+            ! Calculate MLE flux contribution from surface fluxes
+                   DO jk = 1, ibld(ji,jj)
+                     znd = gdepw_n(ji,jj,jk) / 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_n(ji,jj,jk) / 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 = ibld(ji,jj), mld_prof(ji,jj)
+                      zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) )
+                    END DO
+            ! Iterate to find approx vertical index for depth 1.1*zhmle(ji,jj)
+                    jl = MIN(mld_prof(ji,jj) + 2, mbkt(ji,jj))
+                    jl = MIN( MAX(INT( 0.1 * zhmle(ji,jj) / e3t_n(ji,jj,jl)), 2 ) + mld_prof(ji,jj), mbkt(ji,jj))
+                    DO jk = mld_prof(ji,jj),  jl
+                       zdiffut(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), zdiff_mle(ji,jj) * &
+                            & ( gdepw_n(ji,jj,jk) - gdepw_n(ji,jj,jl) ) / &
+                            & ( gdepw_n(ji,jj,mld_prof(ji,jj)) - gdepw_n(ji,jj,jl) - epsln))
+                    END DO
+                ENDIF
+            END DO
+          END DO
+       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
@@ -1058,12 +1192,771 @@
          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. )
+
+CONTAINS
+
+
+   ! Alan: do we need zb?
+   SUBROUTINE zdf_osm_vertical_average( jnlev_av, zt, zs, 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.
+
+        ! Alan: do we need zb?
+        REAL(wp), DIMENSION(jpi,jpj) :: zt, zs        ! 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
+        REAL(wp) :: zthick, zthermal, zbeta
+
+
+        zt   = 0._wp
+        zs   = 0._wp
+        zu   = 0._wp
+        zv   = 0._wp
+        DO jj = 2, jpjm1                                 !  Vertical slab
+         DO ji = 2, jpim1
+            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_n(ji,jj,jk)
+               zt(ji,jj)   = zt(ji,jj)  + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem)
+               zs(ji,jj)   = zs(ji,jj)  + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
+               zu(ji,jj)   = zu(ji,jj)  + e3t_n(ji,jj,jk) &
+                     &            * ( ub(ji,jj,jk) + ub(ji - 1,jj,jk) ) &
+                     &            / MAX( 1. , umask(ji,jj,jk) + umask(ji - 1,jj,jk) )
+               zv(ji,jj)   = zv(ji,jj)  + e3t_n(ji,jj,jk) &
+                     &            * ( vb(ji,jj,jk) + vb(ji,jj - 1,jk) ) &
+                     &            / 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
+            ! Alan: do we need zb?
+            zdt(ji,jj) = zt(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_tem)
+            zds(ji,jj) = zs(ji,jj) - tsn(ji,jj,ibld(ji,jj)+ibld_ext,jp_sal)
+            zdu(ji,jj) = zu(ji,jj) - ( ub(ji,jj,ibld(ji,jj)+ibld_ext) + ub(ji-1,jj,ibld(ji,jj)+ibld_ext ) ) &
+                     &    / MAX(1. , umask(ji,jj,ibld(ji,jj)+ibld_ext ) + umask(ji-1,jj,ibld(ji,jj)+ibld_ext ) )
+            zdv(ji,jj) = zv(ji,jj) - ( vb(ji,jj,ibld(ji,jj)+ibld_ext) + vb(ji,jj-1,ibld(ji,jj)+ibld_ext ) ) &
+                     &   / MAX(1. , vmask(ji,jj,ibld(ji,jj)+ibld_ext ) + vmask(ji,jj-1,ibld(ji,jj)+ibld_ext ) )
+            zdb(ji,jj) = grav * zthermal * zdt(ji,jj) - grav * zbeta * zds(ji,jj)
+         END DO
+        END DO
+   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 jj = 2, jpjm1
+           DO ji = 2, jpim1
+              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 DO
+        END DO
+    END SUBROUTINE zdf_osm_velocity_rotation
+
+    SUBROUTINE zdf_osm_external_gradients( 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.
+     !!
+     !!----------------------------------------------------------------------
+
+     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 jj = 2, jpjm1
+        DO ji = 2, jpim1
+           IF ( ibld(ji,jj) + ibld_ext < 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 = ibld(ji,jj) + ibld_ext
+              jkb1 = MIN(jkb + 1, mbkt(ji,jj))
+              zdtdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_tem) - tsn(ji,jj,jkb,jp_tem ) ) &
+                   &   / e3t_n(ji,jj,ibld(ji,jj))
+              zdsdz(ji,jj) = - ( tsn(ji,jj,jkb1,jp_sal) - tsn(ji,jj,jkb,jp_sal ) ) &
+                   &   / e3t_n(ji,jj,ibld(ji,jj))
+              zdbdz(ji,jj) = grav * zthermal * zdtdz(ji,jj) - grav * zbeta * zdsdz(ji,jj)
+           END IF
+        END DO
+     END DO
+    END SUBROUTINE zdf_osm_external_gradients
+
+    SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz )
+
+     REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz      ! gradients in the pycnocline
+
+     INTEGER :: jk, jj, ji
+     REAL(wp) :: ztgrad, zsgrad, zbgrad
+     REAL(wp) :: zgamma_b_nd, zgamma_c, znd
+     REAL(wp) :: zzeta_s=0.3
+
+     DO jj = 2, jpjm1
+        DO ji = 2, jpim1
+           IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN
+              IF ( lconv(ji,jj) ) THEN  ! convective conditions
+                    IF ( zdbdz_ext(ji,jj) > 0._wp .AND. &
+                         & (zdhdt(ji,jj) > 0._wp .AND. ln_osm_mle .AND. zdb_bl(ji,jj) > rn_osm_mle_thresh &
+                         &  .OR.  zdb_bl(ji,jj) > 0._wp)) THEN  ! zdhdt could be <0 due to FK, hence check zdhdt>0
+                       ztgrad = 0.5 * zdt_ml(ji,jj) / zdh(ji,jj) + zdtdz_ext(ji,jj)
+                       zsgrad = 0.5 * zds_ml(ji,jj) / zdh(ji,jj) + zdsdz_ext(ji,jj)
+                       zbgrad = 0.5 * zdb_ml(ji,jj) / zdh(ji,jj) + zdbdz_ext(ji,jj)
+                       zgamma_b_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj)
+                       zgamma_c = ( 3.14159 / 4.0 ) * ( 0.5 + zgamma_b_nd ) /&
+                            & ( 1.0 - 0.25 * SQRT( 3.14159 / 6.0 ) - 2.0 * zgamma_b_nd * zzeta_s )**2 ! check
+                       DO jk = 2, ibld(ji,jj)+ibld_ext
+                          znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj)
+                          IF ( znd <= zzeta_s ) THEN
+                             zdtdz(ji,jj,jk) = zdtdz_ext(ji,jj) + 0.5 * zdt_ml(ji,jj) / zdh(ji,jj) * &
+                                  & EXP( -6.0 * ( znd -zzeta_s )**2 )
+                             zdsdz(ji,jj,jk) = zdsdz_ext(ji,jj) + 0.5 * zds_ml(ji,jj) / zdh(ji,jj) * &
+                                  & EXP( -6.0 * ( znd -zzeta_s )**2 )
+                             zdbdz(ji,jj,jk) = zdbdz_ext(ji,jj) + 0.5 * zdb_ml(ji,jj) / zdh(ji,jj) * &
+                                  & EXP( -6.0 * ( znd -zzeta_s )**2 )
+                          ELSE
+                             zdtdz(ji,jj,jk) =  ztgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 )
+                             zdsdz(ji,jj,jk) =  zsgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 )
+                             zdbdz(ji,jj,jk) =  zbgrad * EXP( -zgamma_c * ( znd - zzeta_s )**2 )
+                          ENDIF
+                       END DO
+                    ENDIF ! If condition not satisfied, no pycnocline present. Gradients have been initialised to zero, so do nothing
+              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_n(ji,jj,jk) / zhbl(ji,jj)
+                             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.
+                          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_n(ji,jj,jk) - zhml(ji,jj) ) / zdh(ji,jj)
+                             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)
+           END IF      ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) )
+        END DO
+     END DO
+
+    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
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1
+            !
+            IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) ) THEN
+               IF ( lconv (ji,jj) ) THEN
+                  ! Unstable conditions
+                  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_n(ji,jj,jk) - 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_n(ji,jj,jk) / 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 DO
+      END DO
+    END SUBROUTINE zdf_osm_pycnocline_shear_profiles
+
+    SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zdhdt_2 )
+     !!---------------------------------------------------------------------
+     !!                   ***  ROUTINE zdf_osm_calculate_dhdt  ***
+     !!
+     !! ** Purpose : Calculates the rate at which hbl changes.
+     !!
+     !! ** Method  :
+     !!
+     !!----------------------------------------------------------------------
+
+    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zdhdt_2        ! Rate of change of hbl
+
+    INTEGER :: jj, ji
+    REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert
+    REAL(wp) :: zvel_max, zwb_min
+    REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes
+    REAL(wp) :: zzeta_m = 0.3
+    REAL(wp) :: zgamma_c = 2.0
+    REAL(wp) :: zdhoh = 0.1
+    REAL(wp) :: alpha_bc = 0.5
+
+    DO jj = 2, jpjm1
+       DO ji = 2, jpim1
+          IF ( lconv(ji,jj) ) THEN    ! Convective
+             ! Alan is this right?  Yes, it's a bit different from the previous relationship
+             ! zwb_ent(ji,jj) = - 2.0 * 0.2 * zwbav(ji,jj) &
+             !   &            - ( 0.15 * ( 1.0 - EXP( -1.5 * zla(ji,jj) ) ) * zustar(ji,jj)**3 + 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj)
+             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 )
+             zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) &
+                  &                  * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) )
+
+             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)
+             !
+             zwb_min = dh(ji,jj) * zwb0(ji,jj) / zhml(ji,jj) + zwb_ent(ji,jj)
+
+             IF ( ln_osm_mle ) THEN
+                  !  zwb_fk_b(ji,jj) = zwb_fk(ji,jj) *  hmle(ji,jj) / ( 6.0 * hbl(ji,jj) ) * ( 6.0 * hbl(ji,jj) / hmle(ji,jj) - 1.0 + &
+                ! &            ( 1.0 - 2.0 * hbl(ji,jj) / hmle(ji,jj))**3 )          ! Fox-Kemper buoyancy flux average over OSBL
+                IF ( hmle(ji,jj) > hbl(ji,jj) ) THEN
+                   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_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 = - ( 1.0 + 1.0 * ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * rn_rdt / hbl(ji,jj) ) * zwb_ent(ji,jj) / &
+                  &        ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
+                  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
+
+             zdhdt_2(ji,jj) = 0._wp
+
+                ! commented out ajgn 23 July as temporay fix
+!                 IF ( zdb_ml(ji,jj) > 0.0 .and. zdbdz_ext(ji,jj) > 0.0 ) THEN
+! !additional term to account for thickness of pycnocline on dhdt. Small correction, so could get rid of this if necessary.
+!                     zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
+!                     zgamma_b_nd = zdbdz_ext(ji,jj) * zhml(ji,jj) / zdb_ml(ji,jj)
+!                     zgamma_dh_nd = zdbdz_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj)
+!                     zdhdt_2(ji,jj) = ( 1.0 - SQRT( 3.1415 / ( 4.0 * zgamma_c) ) * zdhoh ) * zdh(ji,jj) / zhml(ji,jj)
+!                     zdhdt_2(ji,jj) = zdhdt_2(ji,jj) * ( zwb0(ji,jj) - (1.0 + zgamma_b_nd / alpha_bc ) * zwb_min )
+!                     ! Alan no idea what this should be?
+!                     zdhdt_2(ji,jj) = alpha_bc / ( 4.0 * zgamma_c ) * zdhdt_2(ji,jj) &
+!                        &        + (alpha_bc + zgamma_dh_nd ) * ( 1.0 + SQRT( 3.1414 / ( 4.0 * zgamma_c ) ) * zdh(ji,jj) / zhbl(ji,jj) ) &
+!                        &        * (1.0 / ( 4.0 * zgamma_c * alpha_bc ) ) * zwb_min * zdh(ji,jj) / zhbl(ji,jj)
+!                     zdhdt_2(ji,jj) = zdhdt_2(ji,jj) / ( zvel_max + MAX( zdb_bl(ji,jj), 1.0e-15 ) )
+!                     IF ( zdhdt_2(ji,jj) <= 0.2 * zdhdt(ji,jj) ) THEN
+!                         zdhdt(ji,jj) = zdhdt(ji,jj) + zdhdt_2(ji,jj)
+!                 ENDIF
+            ELSE                        ! Stable
+                zdhdt(ji,jj) = ( 0.06 + 0.52 * zhol(ji,jj) / 2.0 ) * zvstr(ji,jj)**3 / hbl(ji,jj) + zwbav(ji,jj)
+                zdhdt_2(ji,jj) = 0._wp
+                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_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj)
+                ELSE
+                    zpert = MAX( 2.0 * ( 1.0 + 0.0 * 2.0 * zvstr(ji,jj) * rn_rdt / hbl(ji,jj) ) * zvstr(ji,jj)**2 / hbl(ji,jj), zdb_bl(ji,jj) )
+                ENDIF
+                zdhdt(ji,jj) = 2.0 * zdhdt(ji,jj) / zpert
+            ENDIF
+         END DO
+      END DO
+    END SUBROUTINE zdf_osm_calculate_dhdt
+
+    SUBROUTINE zdf_osm_timestep_hbl( zdhdt, zdhdt_2 )
+     !!---------------------------------------------------------------------
+     !!                   ***  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, zdhdt_2   ! rates of change of hbl.
+
+    INTEGER :: jk, jj, ji, jm
+    REAL(wp) :: zhbl_s, zvel_max, zdb
+    REAL(wp) :: zthermal, zbeta
+
+     DO jj = 2, jpjm1
+         DO ji = 2, jpim1
+           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_rdt / 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) - tsn(ji,jj,jm,jp_tem) ) &
+                         & - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), &
+                         &  0.0 ) + zvel_max
+
+
+                    IF ( ln_osm_mle ) THEN
+                       zhbl_s = zhbl_s + MIN( &
+                        & rn_rdt * ( ( -zwb_ent(ji,jj) - 2.0 * zwb_fk_b(ji,jj) )/ zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), &
+                        & e3w_n(ji,jj,jm) )
+                    ELSE
+                      zhbl_s = zhbl_s + MIN( &
+                        & rn_rdt * (  -zwb_ent(ji,jj) / zdb + zdhdt_2(ji,jj) ) / FLOAT(ibld(ji,jj) - imld(ji,jj) ), &
+                        & e3w_n(ji,jj,jm) )
+                    ENDIF
+
+                    zhbl_s = MIN(zhbl_s,  gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)
+
+                    IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) 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) - tsn(ji,jj,jm,jp_tem) )&
+                         &           - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ),&
+                         & 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_rdt / FLOAT( ibld(ji,jj) - imld(ji,jj) ), e3w_n(ji,jj,jm) )
+
+                    zhbl_s = MIN(zhbl_s, gdepw_n(ji,jj, mbkt(ji,jj) + 1) - depth_tol)
+                    IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1
+                 END DO
+             ENDIF   ! IF ( lconv )
+             hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,4) )
+             ibld(ji,jj) = MAX(jm, 4 )
+           ELSE
+! change zero or one model level.
+             hbl(ji,jj) = MAX(zhbl_t(ji,jj), gdepw_n(ji,jj,4) )
+           ENDIF
+           zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj))
+         END DO
+      END DO
+
+    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, zddhdt
+
+
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1
+            IF ( lconv(ji,jj) ) THEN
+
+               IF( ln_osm_mle ) THEN
+                  IF ( ( zwb_ent(ji,jj) + 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_ext(ji,jj) > 0._wp)THEN
+                        IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN  ! near neutral stability
+                           zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / &
+                                (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 )
+                        ELSE                                                     ! unstable
+                           zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / &
+                                (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 )
+                        ENDIF
+                        ztau = 0.2 * hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird
+                        zddhdt = zari * hbl(ji,jj)
+                     ELSE
+                        ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
+                        zddhdt = 0.2 * hbl(ji,jj)
+                     ENDIF
+                  ELSE
+                     ztau = 0.2 * hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
+                     zddhdt = 0.2 * hbl(ji,jj)
+                  ENDIF
+               ELSE ! ln_osm_mle
+                  IF ( zdb_bl(ji,jj) > 0._wp .and. zdbdz_ext(ji,jj) > 0._wp)THEN
+                     IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN  ! near neutral stability
+                        zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / &
+                             (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zvstr(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 )
+                     ELSE                                                     ! unstable
+                        zari = MIN( 1.5 * ( zdb_bl(ji,jj) / ( zdbdz_ext(ji,jj) * zhbl(ji,jj) ) ) / &
+                             (1.0 + zdb_bl(ji,jj)**2 / ( 4.5 * zwstrc(ji,jj)**2 * zdbdz_ext(ji,jj) ) ), 0.2 )
+                     ENDIF
+                     ztau = hbl(ji,jj) / (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird
+                     zddhdt = zari * hbl(ji,jj)
+                  ELSE
+                     ztau = hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird
+                     zddhdt = 0.2 * hbl(ji,jj)
+                  ENDIF
+
+               END IF  ! ln_osm_mle
+
+               dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) )
+               ! Alan: this hml is never defined or used
+            ELSE   ! IF (lconv)
+               ztau = hbl(ji,jj) / zvstr(ji,jj)
+               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 )
+                     zddhdt = MIN( zari, 0.2 ) * hbl(ji,jj)
+                  ELSE
+                     zddhdt = 0.2 * hbl(ji,jj)
+                  ENDIF
+               ELSE     ! IF(dhdt < 0)
+                  zddhdt = 0.2 * hbl(ji,jj)
+               ENDIF    ! IF (dhdt >= 0)
+               dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau )+ zddhdt * ( 1.0 - EXP( -rn_rdt / ztau ) )
+               IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zddhdt       ! can be a problem with dh>hbl for rapid collapse
+               ! Alan: this hml is never defined or used -- do we need it?
+            ENDIF       ! IF (lconv)
+
+            hml(ji,jj) = hbl(ji,jj) - dh(ji,jj)
+            inhml = MAX( INT( dh(ji,jj) / e3t_n(ji,jj,ibld(ji,jj)) ) , 1 )
+            imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 3)
+            zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))
+            zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
+         END DO
+      END DO
+
+    END SUBROUTINE zdf_osm_pycnocline_thickness
+
+
+   SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_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)     :: 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_rau0   ! convert density criteria into N^2 criteria
+      DO jk = nlb10, jpkm1
+         DO jj = 1, jpj                ! Mixed layer level: w-level 
+            DO ji = 1, jpi
+               ikt = mbkt(ji,jj)
+               zmld(ji,jj) = zmld(ji,jj) + MAX( rn2b(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk)
+               IF( zmld(ji,jj) < zN2_c )   mld_prof(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level
+            END DO
+         END DO
+      END DO
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            mld_prof(ji,jj) = MAX(mld_prof(ji,jj),ibld(ji,jj))
+            zmld(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj))
+         END DO
+      END DO
+      ! ensure mld_prof .ge. ibld
+      !
+      ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 )                  ! max level of the computation
+      !
+      ztm(:,:) = 0._wp
+      zsm(:,:) = 0._wp
+      DO jk = 1, ikmax                                 ! MLD and mean buoyancy and N2 over the mixed layer
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zc = e3t_n(ji,jj,jk) * 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 * tsn(ji,jj,jk,jp_tem)
+               zsm(ji,jj) = zsm(ji,jj) + zc * tsn(ji,jj,jk,jp_sal)
+            END DO
+         END DO
+      END DO
+      ! average temperature and salinity.
+      ztm(:,:) = ztm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) )
+      zsm(:,:) = zsm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) )
+      ! calculate horizontal gradients at u & v points
+
+      DO jj = 2, jpjm1
+         DO ji = 1, jpim1
+            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 DO
+      END DO
+
+      DO jj = 1, jpjm1
+         DO ji = 2, jpim1
+            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 DO
+      END DO
+
+      CALL eos_rab(ztsm_midu, zmld_midu, zabu)
+      CALL eos_rab(ztsm_midv, zmld_midv, zabv)
+
+      DO jj = 2, jpjm1
+         DO ji = 1, jpim1
+            dbdx_mle(ji,jj) = grav*(zdtdx(ji,jj)*zabu(ji,jj,jp_tem) - zdsdx(ji,jj)*zabu(ji,jj,jp_sal))
+         END DO
+      END DO
+      DO jj = 1, jpjm1
+         DO ji = 2, jpim1
+            dbdy_mle(ji,jj) = grav*(zdtdy(ji,jj)*zabv(ji,jj,jp_tem) - zdsdy(ji,jj)*zabv(ji,jj,jp_sal))
+         END DO
+      END DO
+
+ END SUBROUTINE zdf_osm_zmld_horizontal_gradients
+  SUBROUTINE zdf_osm_mle_parameters( hmle, zwb_fk, 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
+
+      REAL(wp), DIMENSION(jpi,jpj)     :: hmle, zwb_fk, zvel_mle, zdiff_mle
+      INTEGER  ::   ji, jj, jk          ! dummy loop indices
+      INTEGER  ::   ii, ij, ik  ! local integers
+      INTEGER , DIMENSION(jpi,jpj)     :: inml_mle
+      REAL(wp) ::   zdb_mle, ztmp, zdbds_mle
+
+      mld_prof(:,:) = 4
+      DO jk = 5, jpkm1
+        DO jj = 2, jpjm1
+          DO ji = 2, jpim1
+            IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk)
+          END DO
+        END DO
+      END DO
+      ! DO jj = 2, jpjm1
+      !    DO ji = 1, jpim1
+      !       zhmle(ji,jj) = gdepw_n(ji,jj,mld_prof(ji,jj))
+      !    END DO
+      !  END DO
+   ! Timestep mixed layer eddy depth.
+      DO jj = 2, jpjm1
+        DO ji = 2, jpim1
+          zdb_mle = grav * (rhop(ji,jj,mld_prof(ji,jj)) - rhop(ji,jj,ibld(ji,jj) )) * r1_rau0 ! check ALMG
+          IF ( lconv(ji,jj) .and. ( zdb_bl(ji,jj) < rn_osm_mle_thresh .and. mld_prof(ji,jj) > ibld(ji,jj) .and. zdb_mle > 0.0 ) ) THEN
+             hmle(ji,jj) = hmle(ji,jj) + zwb0(ji,jj) * rn_rdt / MAX( zdb_mle, rn_osm_mle_thresh )  ! MLE layer deepening through encroachment. Don't have a good maximum value for deepening, so use threshold buoyancy.
+          ELSE
+            ! MLE layer relaxes towards mixed layer depth on timescale tau_mle, or tau_mle/10
+             IF ( hmle(ji,jj) > zmld(ji,jj) ) THEN
+               hmle(ji,jj) = hmle(ji,jj) - ( hmle(ji,jj) - zmld(ji,jj) ) * rn_rdt / rn_osm_mle_tau
+             ELSE
+               hmle(ji,jj) = hmle(ji,jj) - 10.0 * ( hmle(ji,jj) - zmld(ji,jj) ) * rn_rdt / rn_osm_mle_tau ! fast relaxation if MLE layer shallower than MLD
+             ENDIF
+          ENDIF
+          hmle(ji,jj) = MIN(hmle(ji,jj), ht_n(ji,jj))
+        END DO
+      END DO
+
+      mld_prof = 4
+      DO jk = 5, jpkm1
+        DO jj = 2, jpjm1
+          DO ji = 2, jpim1
+            IF ( hmle(ji,jj) >= gdepw_n(ji,jj,jk) ) mld_prof(ji,jj) = MIN(mbkt(ji,jj), jk)
+          END DO
+        END DO
+      END DO
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1
+            zhmle(ji,jj) = gdepw_n(ji,jj, mld_prof(ji,jj))
+         END DO
+       END DO
+   ! Calculate vertical buoyancy, heat and salinity fluxes due to MLE.
+
+      DO jj = 2, jpjm1
+        DO ji = 2, jpim1
+          IF ( lconv(ji,jj) ) THEN
+             ztmp =  r1_ft(ji,jj) *  MIN( 111.e3_wp , e1u(ji,jj) ) / rn_osm_mle_lf
+             ! zwt_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdtdx(ji,jj) + dbdy_mle(ji,jj) * zdtdy(ji,jj) &
+             !      & + dbdx_mle(ji-1,jj) * zdtdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdtdy(ji,jj-1) )
+             ! zws_fk(ji,jj) = 0.5_wp * ztmp * ( dbdx_mle(ji,jj) * zdsdx(ji,jj) + dbdy_mle(ji,jj) * zdsdy(ji,jj) &
+             !      & + dbdx_mle(ji-1,jj) * zdsdx(ji-1,jj) + dbdy_mle(ji,jj-1) * zdsdy(ji,jj-1) )
+             zdbds_mle = 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) ) )
+             zwb_fk(ji,jj) = rn_osm_mle_ce * hmle(ji,jj) * hmle(ji,jj) *ztmp * zdbds_mle * zdbds_mle
+      ! This vbelocity scale, defined in Fox-Kemper et al (2008), is needed for calculating dhdt.
+             zvel_mle(ji,jj) = zdbds_mle * ztmp * hmle(ji,jj) * tmask(ji,jj,1) 
+             zdiff_mle(ji,jj) = 1.e-4_wp * ztmp * zdbds_mle * zhmle(ji,jj)**3  / rn_osm_mle_lf
+          ENDIF
+        END DO
+      END DO
+END SUBROUTINE zdf_osm_mle_parameters
+
+END SUBROUTINE zdf_osm
 
 
@@ -1086,5 +1979,10 @@
      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
+          & ,ln_kpprimix, rn_riinfty, rn_difri, ln_convmix, rn_difconv, nn_osm_wave &
+          & ,ln_osm_mle
+! 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
+
      !!----------------------------------------------------------------------
      !
@@ -1102,6 +2000,7 @@
         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,*) '     Initial hbl for 1D runs                       rn_osm_hbl0   = ', rn_osm_hbl0
@@ -1128,5 +2027,70 @@
      IF( zdf_osm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' )
 
-     call osm_rst( nit000, 'READ' ) !* read or initialize hbl
+
+     IF( ln_osm_mle ) THEN
+! Initialise Fox-Kemper parametrization
+         REWIND( numnam_ref )              ! Namelist namosm_mle in reference namelist : Tracer advection scheme
+         READ  ( numnam_ref, namosm_mle, IOSTAT = ios, ERR = 903)
+903      IF( ios /= 0 )   CALL ctl_nam ( ios , 'namosm_mle in reference namelist')
+
+         REWIND( numnam_cfg )              ! Namelist namosm_mle in configuration namelist : Tracer advection scheme
+         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 ML 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'
+         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 /rau0        ! 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 jj=1,jpj
+            do ji = 1,jpi
+               r1_ft(ji,jj) = MIN(1./( ABS(ff_t(ji,jj)) + epsln ), ABS(ff_t(ji,jj))/z1_t2**2)
+            end do
+         end do
+         ! 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, 'READ' ) !* read or initialize hbl, dh, hmle
+
 
      IF( ln_zdfddm) THEN
@@ -1309,7 +2273,7 @@
         END DO
      END DO
-     hbl = MAX(hbl,epsln)
-     hbli(:,:) = hbl(:,:)
-     DEALLOCATE( imld_rst )
+
+     IF( ln_osm_mle ) hmle(:,:) = hbl(:,:)            ! Initialise MLE depth.
+
      WRITE(numout,*) ' ===>>>> : hbl computed from stratification'
      wn(:,:,:) = 0._wp
