Index: /NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/GYRE_PISCES/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/GYRE_PISCES/EXPREF/namelist_cfg	(revision 12321)
+++ /NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/GYRE_PISCES/EXPREF/namelist_cfg	(revision 12322)
@@ -142,4 +142,8 @@
       rn_Ld        = 100.e+3        !  lateral diffusive length   [m]
 /
+!-----------------------------------------------------------------------
+&namtra_mle    !   mixed layer eddy parametrisation (Fox-Kemper)       (default: OFF)
+!-----------------------------------------------------------------------
+/
 !!======================================================================
 !!                      ***  Dynamics namelists  ***                  !!
@@ -212,4 +216,12 @@
 !-----------------------------------------------------------------------
    nn_etau     =   0       !  penetration of tke below the mixed layer (ML) due to internal & intertial waves
+/
+!-----------------------------------------------------------------------
+&namzdf_osm    !   OSM vertical diffusion                               (ln_zdfosm =T)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namosm_mle    !   mixed layer eddy parametrisation (Fox-Kemper)       (default: OFF)
+!-----------------------------------------------------------------------
 /
 !!======================================================================
Index: /NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg	(revision 12321)
+++ /NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg	(revision 12322)
@@ -135,4 +135,8 @@
 /
 !-----------------------------------------------------------------------
+&namtra_mle    !   mixed layer eddy parametrisation (Fox-Kemper)       (default: OFF)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
 &namsbc_ssr    !   surface boundary condition : sea surface restoring   (ln_ssr =T)
 !-----------------------------------------------------------------------
@@ -377,4 +381,12 @@
 /
 !-----------------------------------------------------------------------
+&namzdf_osm    !   OSM vertical diffusion                               (ln_zdfosm =T)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namosm_mle    !   mixed layer eddy parametrisation (Fox-Kemper)       (default: OFF)
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
 &namzdf_iwm    !    internal wave-driven mixing parameterization        (ln_zdfiwm =T)
 !-----------------------------------------------------------------------
Index: /NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/SHARED/field_def_nemo-oce.xml
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/SHARED/field_def_nemo-oce.xml	(revision 12321)
+++ /NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/SHARED/field_def_nemo-oce.xml	(revision 12322)
@@ -300,6 +300,6 @@
           <field id="emp_oce"      long_name="Evap minus Precip over ocean"         standard_name="evap_minus_precip_over_sea_water"                                     unit="kg/m2/s"   />
           <field id="emp_ice"      long_name="Evap minus Precip over ice"           standard_name="evap_minus_precip_over_sea_ice"                                       unit="kg/m2/s"   />
-          <field id="saltflx"      long_name="Downward salt flux"                                                                                                        unit="1e-3/m2/s" />
-          <field id="fmmflx"       long_name="Water flux due to freezing/melting"                                                                                        unit="kg/m2/s"   />
+          <field id="saltflx"      long_name="Total downward salinity flux"             standard_name="total_downward_salinity_flux"                                   unit="1e-3/m2/s" />
+          <field id="fmmflx"       long_name="Water flux due to freezing/melting"     standard_name="water_flux_due_to_freezing/melting"                  unit="kg/m2/s"   />
           <field id="snowpre"      long_name="Snow precipitation"                   standard_name="snowfall_flux"                                                        unit="kg/m2/s"   />
           <field id="runoffs"      long_name="River Runoffs"                        standard_name="water_flux_into_sea_water_from_rivers"                                unit="kg/m2/s"   />
@@ -776,6 +776,6 @@
 
           <!-- ln_zdfosm=T only (OSMOSIS-OBL) -->
-     <field id="ttrd_osm"      long_name="temperature-trend: OSM-OSBL non-local forcing"                             unit="degC/s * m" >  ttrd_osm * e3t </field>
-     <field id="strd_osm"      long_name="salinity   -trend: OSM-OSBL non-local forcing"                             unit="1e-3/s * m" >  strd_osm * e3t </field>
+     <field id="ttrd_osm_e3t"      long_name="temperature-trend: OSM-OSBL non-local forcing"                             unit="degC/s * m" >  ttrd_osm * e3t </field>
+     <field id="strd_osm_e3t"      long_name="salinity   -trend: OSM-OSBL non-local forcing"                             unit="1e-3/s * m" >  strd_osm * e3t </field>
      
      <!-- -->
@@ -787,13 +787,20 @@
      <field id="strd_npc_e3t"      unit="1e-3/s * m"  >  strd_npc * e3t </field>
      <field id="ttrd_qns_e3t"      unit="degC/s * m"  >  ttrd_qns * e3t_surf </field>
-     <field id="strd_cdt_e3t"      unit="degC/s * m"  >  strd_cdt * e3t_surf </field>
+     <!-- <field id="strd_cdt_e3t"      unit="degC/s * m"  >  strd_cdt * e3t_surf </field> -->
      <field id="ttrd_qsr_e3t"      unit="degC/s * m"  >  ttrd_qsr * e3t </field>
      <field id="ttrd_bbc_e3t"      unit="degC/s * m"  >  ttrd_bbc * e3t </field>
 
      <!-- OMIP  layer-integrated trends -->
+     <field id="sfd"      long_name="Total downward salt flux"                                                     unit="kg/(m^2 s)" > saltflux * 1026.0 * 0.001 </field>
+     <field id="wfo"      long_name="Total downward FW mass flux"                                                   unit="kg/(m^2 s)" > -empmr </field>
+
      <field id="ttrd_totad_li"    long_name="layer integrated heat-trend: total advection"         unit="W/m^2"     > ttrd_totad_e3t * 1026.0 * 3991.86795711963 </field>
      <field id="strd_totad_li"    long_name="layer integrated salt-trend: total advection"         unit="kg/(m^2 s)"    > strd_totad_e3t * 1026.0 * 0.001  </field>
+     <field id="ttrd_osm_li"    long_name="layer integrated heat-trend: non-local OSM"         unit="W/m^2"     > ttrd_osm_e3t * 1026.0 * 3991.86795711963 </field>
+     <field id="strd_osm_li"    long_name="layer integrated salt-trend: non-local OSM"         unit="kg/(m^2 s)"    > strd_osm_e3t * 1026.0 * 0.001  </field>
      <field id="ttrd_evd_li"      long_name="layer integrated heat-trend: EVD convection"          unit="W/m^2"    > ttrd_evd_e3t * 1026.0 * 3991.86795711963 </field>
      <field id="strd_evd_li"      long_name="layer integrated salt-trend: EVD convection"          unit="kg/(m^2 s)"  > strd_evd_e3t * 1026.0 * 0.001  </field>
+     <field id="ttrd_npc_li"      long_name="layer integrated heat-trend: NPC convection"          unit="W/m^2"    > ttrd_npc_e3t * 1026.0 * 3991.86795711963 </field>
+     <field id="strd_npc_li"      long_name="layer integrated salt-trend: NPC convection"          unit="kg/(m^2 s)"  > strd_npc_e3t * 1026.0 * 0.001  </field>
      <field id="ttrd_iso_li"      long_name="layer integrated heat-trend: isopycnal diffusion"     unit="W/m^2" > ttrd_iso_e3t * 1026.0 * 3991.86795711963 </field>
      <field id="strd_iso_li"      long_name="layer integrated salt-trend: isopycnal diffusion"     unit="kg/(m^2 s)" > strd_iso_e3t * 1026.0 * 0.001  </field>
@@ -802,4 +809,5 @@
      <field id="ttrd_qns_li"      long_name="layer integrated heat-trend: non-solar flux + runoff" unit="W/m^2" grid_ref="grid_T_2D"> ttrd_qns_e3t * 1026.0 * 3991.86795711963 </field>
      <field id="ttrd_qsr_li"      long_name="layer integrated heat-trend: solar flux"              unit="W/m^2"  grid_ref="grid_T_3D"> ttrd_qsr_e3t * 1026.0 * 3991.86795711963 </field>
+     <field id="ttrd_bbc_li"      long_name="layer integrated heat-trend: geothermal heating "  unit="W/m^2" > ttrd_bbc_e3t * 1026.0 * 3991.86795711963 </field>
      <field id="ttrd_bbl_li"      long_name="layer integrated heat-trend: bottom boundary layer "  unit="W/m^2" > ttrd_bbl_e3t * 1026.0 * 3991.86795711963 </field>
      <field id="strd_bbl_li"      long_name="layer integrated salt-trend: bottom boundary layer "  unit="kg/(m^2 s)" > strd_bbl_e3t * 1026.0 * 0.001  </field>
Index: /NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/SHARED/grid_def_nemo.xml
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/SHARED/grid_def_nemo.xml	(revision 12321)
+++ /NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/SHARED/grid_def_nemo.xml	(revision 12322)
@@ -128,4 +128,11 @@
          <axis axis_ref="depthw" />
        </grid>
-    </grid_definition>
+        <!--  -->
+       <grid id="grid_T_SFC">
+         <domain domain_ref="grid_T" />
+         <scalar>
+           <extract_axis position="0" />
+         </scalar>
+       </grid>
+     </grid_definition>
 
Index: /NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/SHARED/namelist_ref
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/SHARED/namelist_ref	(revision 12321)
+++ /NEMO/branches/UKMO/NEMO_4.0_FKOSM/cfgs/SHARED/namelist_ref	(revision 12322)
@@ -253,4 +253,7 @@
    ln_COARE_3p5 = .false.   ! "COARE 3.5" algorithm   (Edson et al. 2013)
    ln_ECMWF    = .false.   ! "ECMWF"     algorithm   (IFS cycle 31)
+
+   ln_humi_dpt    = .false.   !  Supply dewpoint tempearture instead of specific humidity (is true for ERA5)
+
       !
       rn_zqt      = 10.       !  Air temperature & humidity reference height (m)
Index: /NEMO/branches/UKMO/NEMO_4.0_FKOSM/src/OCE/SBC/sbcblk.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_FKOSM/src/OCE/SBC/sbcblk.F90	(revision 12321)
+++ /NEMO/branches/UKMO/NEMO_4.0_FKOSM/src/OCE/SBC/sbcblk.F90	(revision 12322)
@@ -124,4 +124,7 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cdn_oce, chn_oce, cen_oce ! needed by Lupkes 2015 bulk scheme
 
+   LOGICAL  ::   ln_humi_dpt = .FALSE.                                        ! calculate specific hunidity from dewpoint
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qair                      ! specific humidity of air at input height
+
    INTEGER  ::   nblk           ! choice of the bulk algorithm
    !                            ! associated indices:
@@ -145,5 +148,5 @@
       !!-------------------------------------------------------------------
       ALLOCATE( Cd_atm (jpi,jpj), Ch_atm (jpi,jpj), Ce_atm (jpi,jpj), t_zu(jpi,jpj), q_zu(jpi,jpj), &
-         &      cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), STAT=sbc_blk_alloc )
+         &      cdn_oce(jpi,jpj), chn_oce(jpi,jpj), cen_oce(jpi,jpj), qair(jpi,jpj), STAT=sbc_blk_alloc )
       !
       CALL mpp_sum ( 'sbcblk', sbc_blk_alloc )
@@ -171,5 +174,5 @@
       NAMELIST/namsbc_blk/ sn_wndi, sn_wndj, sn_humi, sn_qsr, sn_qlw ,                &   ! input fields
          &                 sn_tair, sn_prec, sn_snow, sn_slp, sn_tdif,                &
-         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF,             &   ! bulk algorithm
+         &                 ln_NCAR, ln_COARE_3p0, ln_COARE_3p5, ln_ECMWF, ln_humi_dpt,&   ! bulk algorithm
          &                 cn_dir , ln_taudif, rn_zqt, rn_zu,                         & 
          &                 rn_pfac, rn_efac, rn_vfac, ln_Cd_L12, ln_Cd_L15
@@ -323,4 +326,11 @@
       !
       !                                            ! compute the surface ocean fluxes using bulk formulea
+      ! ..... if dewpoint supplied instead of specific humidaity, calculate specific humidity
+      IF(ln_humi_dpt) THEN
+         qair(:,:) = q_sat( sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
+      ELSE
+         qair(:,:) = sf(jp_humi)%fnow(:,:,1)
+      END IF
+      
       IF( MOD( kt - 1, nn_fsbc ) == 0 )   CALL blk_oce( kt, sf, sst_m, ssu_m, ssv_m )
 
@@ -332,5 +342,5 @@
          ENDIF 
          tatm_ice(:,:)    = sf(jp_tair)%fnow(:,:,1)
-         qatm_ice(:,:)    = sf(jp_humi)%fnow(:,:,1)
+         qatm_ice(:,:)    = qair(:,:)
          tprecip(:,:)     = sf(jp_prec)%fnow(:,:,1) * rn_pfac
          sprecip(:,:)     = sf(jp_snow)%fnow(:,:,1) * rn_pfac
@@ -434,15 +444,15 @@
       !!    (see Josey, Gulev & Yu, 2013) / doi=10.1016/B978-0-12-391851-2.00005-2
       !!    (since reanalysis products provide T at z, not theta !)
-      ztpot = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1) ) * rn_zqt
+      ztpot(:,:) = sf(jp_tair)%fnow(:,:,1) + gamma_moist( sf(jp_tair)%fnow(:,:,1), qair(:,:) ) * rn_zqt
 
       SELECT CASE( nblk )        !==  transfer coefficients  ==!   Cd, Ch, Ce at T-point
       !
-      CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! NCAR-COREv2
+      CASE( np_NCAR      )   ;   CALL turb_ncar    ( rn_zqt, rn_zu, zst, ztpot, zsq, qair, wndm,   &  ! NCAR-COREv2
          &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
-      CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.0
+      CASE( np_COARE_3p0 )   ;   CALL turb_coare   ( rn_zqt, rn_zu, zst, ztpot, zsq, qair, wndm,   &  ! COARE v3.0
          &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
-      CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! COARE v3.5
+      CASE( np_COARE_3p5 )   ;   CALL turb_coare3p5( rn_zqt, rn_zu, zst, ztpot, zsq, qair, wndm,   &  ! COARE v3.5
          &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
-      CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, sf(jp_humi)%fnow, wndm,   &  ! ECMWF
+      CASE( np_ECMWF     )   ;   CALL turb_ecmwf   ( rn_zqt, rn_zu, zst, ztpot, zsq, qair, wndm,   &  ! ECMWF
          &                                           Cd_atm, Ch_atm, Ce_atm, t_zu, q_zu, zU_zu, cdn_oce, chn_oce, cen_oce )
       CASE DEFAULT
@@ -454,5 +464,5 @@
          zrhoa(:,:) = rho_air( t_zu(:,:)              , q_zu(:,:)              , sf(jp_slp)%fnow(:,:,1) )
       ELSE                                      ! At zt:
-         zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
+         zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), qair(:,:), sf(jp_slp)%fnow(:,:,1) )
       END IF
 
@@ -495,11 +505,11 @@
       IF( ABS( rn_zu - rn_zqt) < 0.01_wp ) THEN
          !! q_air and t_air are given at 10m (wind reference height)
-         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - sf(jp_humi)%fnow(:,:,1)) ) ! Evaporation, using bulk wind speed
-         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:)             )   ! Sensible Heat, using bulk wind speed
+         zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - qair(:,:)) ) ! Evaporation, using bulk wind speed
+         zqsb (:,:) = cp_air(qair(:,:))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - ztpot(:,:)             )   ! Sensible Heat, using bulk wind speed
       ELSE
          !! q_air and t_air are not given at 10m (wind reference height)
          ! Values of temp. and hum. adjusted to height of wind during bulk algorithm iteration must be used!!!
          zevap(:,:) = rn_efac*MAX( 0._wp,             zqla(:,:)*Ce_atm(:,:)*(zsq(:,:) - q_zu(:,:) ) ) ! Evaporation, using bulk wind speed
-         zqsb (:,:) = cp_air(sf(jp_humi)%fnow(:,:,1))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) )   ! Sensible Heat, using bulk wind speed
+         zqsb (:,:) = cp_air(qair(:,:))*zqla(:,:)*Ch_atm(:,:)*(zst(:,:) - t_zu(:,:) )   ! Sensible Heat, using bulk wind speed
       ENDIF
 
@@ -742,5 +752,5 @@
       ! local scalars ( place there for vector optimisation purposes)
       ! Computing density of air! Way denser that 1.2 over sea-ice !!!
-      zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1))
+      zrhoa (:,:) =  rho_air(sf(jp_tair)%fnow(:,:,1), qair(:,:), sf(jp_slp)%fnow(:,:,1))
 
       !!gm brutal....
@@ -806,5 +816,5 @@
       zcoef_dqla = -Ls * 11637800. * (-5897.8)
       !
-      zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), sf(jp_humi)%fnow(:,:,1), sf(jp_slp)%fnow(:,:,1) )
+      zrhoa(:,:) = rho_air( sf(jp_tair)%fnow(:,:,1), qair(:,:), sf(jp_slp)%fnow(:,:,1) )
       !
       zztmp = 1. / ( 1. - albo )
@@ -837,5 +847,5 @@
                ! Latent Heat
                qla_ice(ji,jj,jl) = rn_efac * MAX( 0.e0, zrhoa(ji,jj) * Ls  * Ch_atm(ji,jj) * wndm_ice(ji,jj) *  &
-                  &                ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - sf(jp_humi)%fnow(ji,jj,1) ) )
+                  &                ( 11637800. * EXP( -5897.8 * z1_st(ji,jj,jl) ) / zrhoa(ji,jj) - qair(ji,jj) ) )
                ! Latent heat sensitivity for ice (Dqla/Dt)
                IF( qla_ice(ji,jj,jl) > 0._wp ) THEN
Index: /NEMO/branches/UKMO/NEMO_4.0_FKOSM/src/OCE/ZDF/zdfosm.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_FKOSM/src/OCE/ZDF/zdfosm.F90	(revision 12321)
+++ /NEMO/branches/UKMO/NEMO_4.0_FKOSM/src/OCE/ZDF/zdfosm.F90	(revision 12322)
@@ -141,5 +141,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$
+   !! $Id: zdfosm.F90 12317 2020-01-14 12:40:47Z agn $
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
@@ -483,9 +483,12 @@
      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
      ! 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.
       hbl(:,:) = MAX(hbl(:,:), gdepw_n(:,:,4) )
       ibld(:,:) = 4
       DO jk = 5, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
                IF ( hbl(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN
                   ibld(ji,jj) = MIN(mbkt(ji,jj), jk)
@@ -526,5 +529,5 @@
          END DO
 
-      imld = ibld           ! use imld to hold previous blayer index
+      imld(:,:) = ibld(:,:)           ! use imld to hold previous blayer index
       ibld(:,:) = 4
 
@@ -1008,4 +1011,10 @@
        END DO
 
+       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
+
 ! KPP-style Ri# mixing
        IF( ln_kpprimix) THEN
@@ -1139,4 +1148,6 @@
            END DO
         END DO
+        ! 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)
@@ -2186,5 +2197,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
@@ -2204,5 +2218,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
@@ -2232,4 +2246,14 @@
            CALL iom_get( numror, jpdom_autoglo, '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_autoglo, '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' & 'dh' not in restart file, recalculate
@@ -2243,9 +2267,10 @@
      IF( TRIM(cdrw) == 'WRITE') THEN     !* Write hbli into the restart file, then return
         IF(lwp) WRITE(numout,*) '---- osm-rst ----'
-         IF( lwxios ) CALL iom_swap(      cwxios_context          )
-            CALL iom_rstput( kt, nitrst, numrow, 'wn'     , wn  , ldxios = lwxios )
-            CALL iom_rstput( kt, nitrst, numrow, 'hbl'    , hbl , ldxios = lwxios )
-            CALL iom_rstput( kt, nitrst, numrow, 'dh'     , dh, ldxios = lwxios )
-         IF( lwxios ) CALL iom_swap(      cxios_context          )
+         CALL iom_rstput( kt, nitrst, numrow, 'wn'     , wn,   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
@@ -2281,7 +2306,11 @@
      END DO
 
-     IF( ln_osm_mle ) hmle(:,:) = hbl(:,:)            ! Initialise MLE depth.
-
      WRITE(numout,*) ' ===>>>> : hbl computed from stratification'
+
+     IF( ln_osm_mle ) THEN
+        hmle(:,:) = hbl(:,:)            ! Initialise MLE depth.
+        WRITE(numout,*) ' ===>>>> : hmle set = to hbl'
+     END IF
+
      wn(:,:,:) = 0._wp
      WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially'
Index: /NEMO/branches/UKMO/NEMO_4.0_FKOSM/src/OCE/ZDF/zdfphy.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_FKOSM/src/OCE/ZDF/zdfphy.F90	(revision 12321)
+++ /NEMO/branches/UKMO/NEMO_4.0_FKOSM/src/OCE/ZDF/zdfphy.F90	(revision 12322)
@@ -56,5 +56,5 @@
    !!----------------------------------------------------------------------
    !! NEMO/OCE 4.0 , NEMO Consortium (2018)
-   !! $Id$
+   !! $Id: zdfphy.F90 12178 2019-12-11 11:02:38Z agn $
    !! Software governed by the CeCILL license (see ./LICENSE)
    !!----------------------------------------------------------------------
@@ -172,5 +172,5 @@
       IF( ln_zdfosm .AND. ln_zdfevd )   CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfosm and ln_zdfevd' )
       IF( lk_top    .AND. ln_zdfnpc )   CALL ctl_stop( 'zdf_phy_init: npc scheme is not working with key_top' )
-      IF( lk_top    .AND. ln_zdfosm )   CALL ctl_stop( 'zdf_phy_init: osmosis scheme is not working with key_top' )
+      IF( lk_top    .AND. ln_zdfosm )   CALL ctl_warn( 'zdf_phy_init: osmosis gives no non-local fluxes for TOP tracers yet' )
       IF(lwp) THEN
          WRITE(numout,*)
