Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8920)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8930)
@@ -162,6 +162,6 @@
    REAL(wp)          , PUBLIC ::   rn_amax_n       !: maximum ice concentration Northern hemisphere
    REAL(wp)          , PUBLIC ::   rn_amax_s       !: maximum ice concentration Southern hemisphere
-   CHARACTER(len=32) , PUBLIC ::   cn_icerst_in    !: suffix of ice restart name (input)
-   CHARACTER(len=32) , PUBLIC ::   cn_icerst_out   !: suffix of ice restart name (output)
+   CHARACTER(len=256), PUBLIC ::   cn_icerst_in    !: suffix of ice restart name (input)
+   CHARACTER(len=256), PUBLIC ::   cn_icerst_out   !: suffix of ice restart name (output)
    CHARACTER(len=256), PUBLIC ::   cn_icerst_indir !: ice restart input directory
    CHARACTER(len=256), PUBLIC ::   cn_icerst_outdir!: ice restart output directory
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/C1D/dyncor_c1d.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/C1D/dyncor_c1d.F90	(revision 8920)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/C1D/dyncor_c1d.F90	(revision 8930)
@@ -21,4 +21,7 @@
    USE prtctl         ! Print control
 
+   USE sbcwave        ! Surface Waves (add Stokes-Coriolis force)
+   USE sbc_oce , ONLY : ln_stcor    ! use Stoke-Coriolis force
+   
    IMPLICIT NONE
    PRIVATE
@@ -71,12 +74,24 @@
       ENDIF
       !
-      DO jk = 1, jpkm1
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               ua(ji,jj,jk) = ua(ji,jj,jk) + ff_t(ji,jj) * vn(ji,jj,jk)
-               va(ji,jj,jk) = va(ji,jj,jk) - ff_t(ji,jj) * un(ji,jj,jk)
+      IF( ln_stcor ) THEN
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ua(ji,jj,jk) = ua(ji,jj,jk) + ff_t(ji,jj) * (vn(ji,jj,jk) + vsd(ji,jj,jk))
+                  va(ji,jj,jk) = va(ji,jj,jk) - ff_t(ji,jj) * (un(ji,jj,jk) + usd(ji,jj,jk))
+               END DO
             END DO
          END DO
-      END DO   
+      ELSE
+         DO jk = 1, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1   ! vector opt.
+                  ua(ji,jj,jk) = ua(ji,jj,jk) + ff_t(ji,jj) * vn(ji,jj,jk)
+                  va(ji,jj,jk) = va(ji,jj,jk) - ff_t(ji,jj) * un(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+      END IF
+      
       !
       IF(ln_ctl)   CALL prt_ctl( tab3d_1=ua, clinfo1=' cor  - Ua: ', mask1=umask,  &
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/C1D/step_c1d.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/C1D/step_c1d.F90	(revision 8920)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/C1D/step_c1d.F90	(revision 8930)
@@ -78,4 +78,8 @@
                          CALL zdf_phy( kstp )         ! vertical physics update (bfr, avt, avs, avm + MLD)
 
+      IF(.NOT.ln_linssh )   CALL ssh_nxt       ( kstp )  ! after ssh (includes call to div_hor)
+      IF(.NOT.ln_linssh )   CALL dom_vvl_sf_nxt( kstp )  ! after vertical scale factors 
+
+      IF(.NOT.ln_linssh )   CALL wzv           ( kstp )  ! now cross-level velocity 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! diagnostics and outputs             (ua, va, ta, sa used as workspace)
@@ -100,4 +104,6 @@
       IF( ln_traqsr )   CALL tra_qsr( kstp )       ! penetrative solar radiation qsr
       IF( ln_tradmp )   CALL tra_dmp( kstp )       ! internal damping trends- tracers
+      IF(.NOT.ln_linssh)CALL tra_adv( kstp )       ! horizontal & vertical advection
+      IF( lk_zdfosm  )  CALL tra_osm( kstp )       ! OSMOSIS non-local tracer fluxes
                         CALL tra_zdf( kstp )       ! vertical mixing
                         CALL eos( tsn, rhd, rhop, gdept_0(:,:,:) )   ! now potential density for zdfmxl
@@ -115,6 +121,10 @@
       IF( ln_dyndmp )   CALL dyn_dmp    ( kstp )   ! internal damping trends- momentum
                         CALL dyn_cor_c1d( kstp )   ! vorticity term including Coriolis
+      IF( lk_zdfosm  )  CALL dyn_osm    ( kstp )   ! OSMOSIS non-local velocity fluxes
                         CALL dyn_zdf    ( kstp )   ! vertical diffusion
                         CALL dyn_nxt    ( kstp )   ! lateral velocity at next time step
+      IF(.NOT.ln_linssh)CALL ssh_swp    ( kstp )   ! swap of sea surface height
+
+      IF(.NOT.ln_linssh)CALL dom_vvl_sf_swp( kstp )! swap of vertical scale factors
 
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/DIA/diacfl.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/DIA/diacfl.F90	(revision 8920)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/DIA/diacfl.F90	(revision 8930)
@@ -105,7 +105,7 @@
       !                    ! write out to file
       IF( lwp ) THEN
-         WRITE(numcfl,FMT='(2x,i4,5x,a6,5x,f6.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zCu_max, iloc_u(1), iloc_u(2), iloc_u(3)
-         WRITE(numcfl,FMT='(11x,     a6,5x,f6.4,1x,i4,1x,i4,1x,i4)')     'Max Cv', zCv_max, iloc_v(1), iloc_v(2), iloc_v(3)
-         WRITE(numcfl,FMT='(11x,     a6,5x,f6.4,1x,i4,1x,i4,1x,i4)')     'Max Cw', zCw_max, iloc_w(1), iloc_w(2), iloc_w(3)
+         WRITE(numcfl,FMT='(2x,i4,5x,a6,4x,f7.4,1x,i4,1x,i4,1x,i4)') kt, 'Max Cu', zCu_max, iloc_u(1), iloc_u(2), iloc_u(3)
+         WRITE(numcfl,FMT='(11x,     a6,4x,f7.4,1x,i4,1x,i4,1x,i4)')     'Max Cv', zCv_max, iloc_v(1), iloc_v(2), iloc_v(3)
+         WRITE(numcfl,FMT='(11x,     a6,4x,f7.4,1x,i4,1x,i4,1x,i4)')     'Max Cw', zCw_max, iloc_w(1), iloc_w(2), iloc_w(3)
       ENDIF
       !
@@ -119,11 +119,11 @@
          ! to ascii file
          WRITE(numcfl,*) '******************************************'
-         WRITE(numcfl,FMT='(3x,a12,7x,f6.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cu', rCu_max, nCu_loc(1), nCu_loc(2), nCu_loc(3)
+         WRITE(numcfl,FMT='(3x,a12,6x,f7.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cu', rCu_max, nCu_loc(1), nCu_loc(2), nCu_loc(3)
          WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', z2dt/rCu_max
          WRITE(numcfl,*) '******************************************'
-         WRITE(numcfl,FMT='(3x,a12,7x,f6.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cv', rCv_max, nCv_loc(1), nCv_loc(2), nCv_loc(3)
+         WRITE(numcfl,FMT='(3x,a12,6x,f7.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cv', rCv_max, nCv_loc(1), nCv_loc(2), nCv_loc(3)
          WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', z2dt/rCv_max
          WRITE(numcfl,*) '******************************************'
-         WRITE(numcfl,FMT='(3x,a12,7x,f6.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cw', rCw_max, nCw_loc(1), nCw_loc(2), nCw_loc(3)
+         WRITE(numcfl,FMT='(3x,a12,6x,f7.4,1x,i4,1x,i4,1x,i4)') 'Run Max Cw', rCw_max, nCw_loc(1), nCw_loc(2), nCw_loc(3)
          WRITE(numcfl,FMT='(3x,a8,11x,f15.1)') ' => dt/C', z2dt/rCw_max
          CLOSE( numcfl ) 
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90	(revision 8920)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/DIA/diahth.F90	(revision 8930)
@@ -139,5 +139,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            zztmp = bathy(ji,jj)
+            zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 
             hth     (ji,jj) = zztmp
             zabs2   (ji,jj) = zztmp
@@ -150,5 +150,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zztmp = bathy(ji,jj)
+               zztmp = gdepw_n(ji,jj,mbkt(ji,jj)+1) 
                zrho0_3(ji,jj) = zztmp
                zrho0_1(ji,jj) = zztmp
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 8920)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/DOM/domain.F90	(revision 8930)
@@ -37,5 +37,4 @@
    USE domvvl         ! variable volume
    USE c1d            ! 1D configuration
-   USE domc1d         ! 1D configuration: column location
    USE dyncor_c1d     ! 1D configuration: Coriolis term    (cor_c1d routine)
    USE wet_dry        ! wetting and drying
@@ -117,9 +116,4 @@
          WRITE(numout,*)     '              nn_cfg = ', nn_cfg
       ENDIF
-      !
-      !      
-!!gm  This should be removed with the new configuration interface
-      IF( lk_c1d .AND. ln_c1d_locpt )  CALL dom_c1d( rn_lat1d, rn_lon1d )
-!!gm end
       !
       !           !==  Reference coordinate system  ==!
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 8920)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 8930)
@@ -87,7 +87,10 @@
       INTEGER  ::   jj, ji, jk   ! dummy loop argument
       INTEGER  ::   ik           ! local integer 
-      REAL(wp) ::  ztransp, zfac, ztemp, zsp0
-      REAL(wp) ::  zdep_u, zdep_v, zkh_u, zkh_v, zda_u, zda_v
+      REAL(wp) ::  ztransp, zfac, zsp0
+      REAL(wp) ::  zdepth, zsqrt_depth,  zexp_depth, z_two_thirds, zsqrtpi !sqrt of pi
+      REAL(wp) ::  zbot_u, zbot_v, zkb_u, zkb_v, zke3_u, zke3_v, zda_u, zda_v
+      REAL(wp) ::  zstokes_psi_u_bot, zstokes_psi_v_bot
       REAL(wp), DIMENSION(:,:)  , POINTER ::   zk_t, zk_u, zk_v, zu0_sd, zv0_sd   ! 2D workspace
+      REAL(wp), DIMENSION(:,:)  , POINTER ::   zstokes_psi_u_top, zstokes_psi_v_top   ! 2D workspace
       REAL(wp), DIMENSION(:,:,:), POINTER ::   ze3divh                            ! 3D workspace
       !!---------------------------------------------------------------------
@@ -95,6 +98,9 @@
       CALL wrk_alloc( jpi,jpj,jpk,   ze3divh )
       CALL wrk_alloc( jpi,jpj,       zk_t, zk_u, zk_v, zu0_sd, zv0_sd )
-      !
-      !
+      CALL wrk_alloc( jpi,jpj,       zstokes_psi_u_top, zstokes_psi_v_top)
+      !
+      !
+      zsqrtpi = SQRT(rpi)
+      z_two_thirds = 2.0_wp / 3.0_wp
       zfac =  2.0_wp * rpi / 16.0_wp
       DO jj = 1, jpj                ! exp. wave number at t-point    (Eq. (19) in Breivick et al. (2014) )
@@ -106,7 +112,8 @@
                tsd2d(ji,jj) = zsp0
                ! Wavenumber scale
-               zk_t(ji,jj) = ABS( zsp0 ) / MAX( ABS( 5.97_wp*ztransp ) , 0.0000001_wp )
-         END DO
-      END DO      
+               zk_t(ji,jj) =  (1.0_wp-2.0_wp/3.0_wp)*zsp0/MAX(2.0_wp*ztransp,0.0000001_wp)
+         END DO
+      END DO
+      !
       DO jj = 1, jpjm1              ! exp. wave number & Stokes drift velocity at u- & v-points
          DO ji = 1, jpim1
@@ -120,16 +127,42 @@
       !
       !                       !==  horizontal Stokes Drift 3D velocity  ==!
+      DO jj = 1, jpjm1              ! exp. wave number & Stokes drift velocity at u- & v-points
+         DO ji = 1, jpim1
+            zstokes_psi_u_top(ji,jj) = 0._wp
+            zstokes_psi_v_top(ji,jj) = 0._wp
+         END DO
+      END DO
+
       DO jk = 1, jpkm1
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
-               zdep_u = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji+1,jj,jk) )
-               zdep_v = 0.5_wp * ( gdept_n(ji,jj,jk) + gdept_n(ji,jj+1,jk) )
-               !                          
-               zkh_u = zk_u(ji,jj) * zdep_u     ! k * depth
-               zkh_v = zk_v(ji,jj) * zdep_v
-               !                                ! Depth attenuation
-               zda_u = EXP( -2.0_wp*zkh_u ) / ( 1.0_wp + 8.0_wp*zkh_u )
-               zda_v = EXP( -2.0_wp*zkh_v ) / ( 1.0_wp + 8.0_wp*zkh_v )
+               zbot_u = 0.5_wp * ( gdepw_n(ji,jj,jk+1) + gdepw_n(ji+1,jj,jk+1) )
+               zbot_v = 0.5_wp * ( gdepw_n(ji,jj,jk+1) + gdepw_n(ji,jj+1,jk+1) )
+               zkb_u = 2.0_wp * zk_u(ji,jj) * zbot_u     ! 2k * bottom depth
+               zkb_v = 2.0_wp * zk_v(ji,jj) * zbot_v     ! 2k * bottom depth
                !
+               zke3_u = MAX(1.e-8_wp, 2.0_wp * zk_u(ji,jj) * e3u_n(ji,jj,jk))     ! 2k * thickness
+               zke3_v = MAX(1.e-8_wp, 2.0_wp * zk_v(ji,jj) * e3v_n(ji,jj,jk))     ! 2k * thickness
+
+               ! Depth attenuation .... do u component first..
+               zdepth=zkb_u
+               zsqrt_depth = SQRT(zdepth)
+               zexp_depth = EXP(-zdepth)
+               zstokes_psi_u_bot = 1.0_wp - zexp_depth  &
+                    & - z_two_thirds*(zsqrtpi*zsqrt_depth*zdepth*ERFC(zsqrt_depth) &
+                    &                  + 1.0_wp - (1.0_wp + zdepth)*zexp_depth)
+               zda_u = (zstokes_psi_u_bot - zstokes_psi_u_top(ji,jj))/zke3_u
+               zstokes_psi_u_top(ji,jj) = zstokes_psi_u_bot
+
+               !         ... and then v component
+               zdepth=zkb_v
+               zsqrt_depth = SQRT(zdepth)
+               zexp_depth = EXP(-zdepth)
+               zstokes_psi_v_bot = 1.0_wp - zexp_depth  &
+                    & - z_two_thirds*(zsqrtpi*zsqrt_depth*zdepth*ERFC(zsqrt_depth) &
+                    &                  + 1.0_wp - (1.0_wp + zdepth)*zexp_depth)
+               zda_v = (zstokes_psi_v_bot - zstokes_psi_v_top(ji,jj))/zke3_v
+               zstokes_psi_v_top(ji,jj) = zstokes_psi_v_bot
+              !
                usd(ji,jj,jk) = zda_u * zu0_sd(ji,jj) * umask(ji,jj,jk)
                vsd(ji,jj,jk) = zda_v * zv0_sd(ji,jj) * vmask(ji,jj,jk)
@@ -190,4 +223,5 @@
       CALL wrk_dealloc( jpi,jpj,jpk,   ze3divh )
       CALL wrk_dealloc( jpi,jpj,       zk_t, zk_u, zk_v, zu0_sd, zv0_sd )
+      CALL wrk_dealloc( jpi,jpj,       zstokes_psi_u_top, zstokes_psi_v_top)
       !
    END SUBROUTINE sbc_stokes
@@ -347,4 +381,10 @@
          ALLOCATE( div_sd(jpi,jpj) )
          ALLOCATE( tsd2d (jpi,jpj) )
+
+         ut0sd(:,:) = 0._wp
+         vt0sd(:,:) = 0._wp
+         hsw(:,:) = 0._wp
+         wmp(:,:) = 0._wp
+
          usd(:,:,:) = 0._wp
          vsd(:,:,:) = 0._wp
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 8920)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdf_oce.F90	(revision 8930)
@@ -22,5 +22,6 @@
    LOGICAL , PUBLIC ::   ln_zdfric   !: Richardson depend coefficients
    LOGICAL , PUBLIC ::   ln_zdftke   !: Turbulent Kinetic Energy closure
-   LOGICAL , PUBLIC ::   ln_zdfgls   !: Generic Length Sclare closure
+   LOGICAL , PUBLIC ::   ln_zdfgls   !: Generic Length Scale closure
+   LOGICAL , PUBLIC ::   ln_zdfosm   !: OSMOSIS BL closure
    !                             ! convection
    LOGICAL , PUBLIC ::   ln_zdfevd   !: convection: enhanced vertical diffusion flag
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfosm.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfosm.F90	(revision 8930)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfosm.F90	(revision 8930)
@@ -0,0 +1,1720 @@
+MODULE zdfosm
+   !!======================================================================
+   !!                       ***  MODULE  zdfosm  ***
+   !! Ocean physics:  vertical mixing coefficient compute from the OSMOSIS
+   !!                 turbulent closure parameterization
+   !!=====================================================================
+   !!  History : NEMO 4.0  ! A. Grant, G. Nurser
+
+   !! 15/03/2017  Changed calculation of pycnocline thickness in unstable conditions and stable conditions AG
+   !! 15/03/2017  Calculation of pycnocline gradients for stable conditions changed. Pycnocline gradients now depend on stability of the OSBL. A.G
+   !! 06/06/2017  (1) Checks on sign of buoyancy jump in calculation of  OSBL depth.  A.G.
+   !!             (2) Removed variable zbrad0, zbradh and zbradav since they are not used.
+   !!             (3) Approximate treatment for shear turbulence.
+   !!                        Minimum values for zustar and zustke.
+   !!                        Add velocity scale, zvstr, that tends to zustar for large Langmuir numbers.
+   !!                        Limit maximum value for Langmuir number.
+   !!                        Use zvstr in definition of stability parameter zhol.
+   !!             (4) Modified parametrization of entrainment flux, changing original coefficient 0.0485 for Langmuir contribution to 0.135 * zla
+   !!             (5) For stable boundary layer add factor that depends on length of timestep to 'slow' collapse and growth. Make sure buoyancy jump not negative.
+   !!             (6) For unstable conditions when growth is over multiple levels, limit change to maximum of one level per cycle through loop.
+   !!             (7) Change lower limits for loops that calculate OSBL averages from 1 to 2. Large gradients between levels 1 and 2 can cause problems.
+   !!             (8) Change upper limits from ibld-1 to ibld.
+   !!             (9) Calculation of pycnocline thickness in unstable conditions. Check added to ensure that buoyancy jump is positive before calculating Ri.
+   !!            (10) Thickness of interface layer at base of the stable OSBL set by Richardson number. Gives continuity in transition from unstable OSBL.
+   !!            (11) Checks that buoyancy jump is poitive when calculating pycnocline profiles.
+   !!            (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).
+   !! 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)
+   !!
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   'key_zdfosm'                                             OSMOSIS scheme
+   !!----------------------------------------------------------------------
+   !!   zdf_osm       : update momentum and tracer Kz from osm scheme
+   !!   zdf_osm_init  : initialization, namelist read, and parameters control
+   !!   osm_rst       : read (or initialize) and write osmosis restart fields
+   !!   tra_osm       : compute and add to the T & S trend the non-local flux
+   !!   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
+   !!----------------------------------------------------------------------
+   USE oce !, wb => wn  ! ocean dynamics and active tracers
+                      ! uses wn from previous time step (which is now wb) to calculate hbl
+   USE dom_oce        ! ocean space and time domain
+   USE zdf_oce        ! ocean vertical physics
+   USE sbc_oce        ! surface boundary condition: ocean
+   USE sbcwave        ! surface wave parameters
+   USE phycst         ! physical constants
+   USE eosbn2         ! equation of state
+   USE traqsr         ! details of solar radiation absorption
+   USE zdfddm         ! double diffusion mixing (avs array)
+   USE iom            ! I/O library
+   USE lib_mpp        ! MPP library
+   USE trd_oce        ! ocean trends definition
+   USE trdtra         ! tracers trends
+   !
+   USE in_out_manager ! I/O manager
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE prtctl         ! Print control
+   USE timing         ! Timing
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   zdf_osm       ! routine called by step.F90
+   PUBLIC   zdf_osm_init  ! routine called by nemogcm.F90
+   PUBLIC   osm_rst       ! routine called by step.F90
+   PUBLIC   tra_osm       ! routine called by step.F90
+   PUBLIC   trc_osm       ! routine called by trcstp.F90
+   PUBLIC   dyn_osm       ! routine called by 'step.F90'
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamu    !: non-local u-momentum flux
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamv    !: non-local v-momentum flux
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamt    !: non-local temperature flux (gamma/<ws>o)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghams    !: non-local salinity flux (gamma/<ws>o)
+   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(:,:)   ::   dstokes  ! penetration depth of the Stokes drift.
+
+   !                      !!** Namelist  namzdf_osm  **
+   LOGICAL  ::   ln_use_osm_la      ! Use namelist  rn_osm_la
+   REAL(wp) ::   rn_osm_la          ! Turbulent Langmuir number
+   REAL(wp) ::   rn_osm_dstokes     ! Depth scale of Stokes drift
+   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
+   LOGICAL  ::   ln_dia_osm         ! Use namelist  rn_osm_la
+
+
+   LOGICAL  ::   ln_kpprimix  = .true.  ! Shear instability mixing
+   REAL(wp) ::   rn_riinfty   = 0.7     ! local Richardson Number limit for shear instability
+   REAL(wp) ::   rn_difri    =  0.005   ! maximum shear mixing at Rig = 0    (m2/s)
+   LOGICAL  ::   ln_convmix  = .true.   ! Convective instability mixing
+   REAL(wp) ::   rn_difconv = 1._wp     ! diffusivity when unstable below BL  (m2/s)
+
+   !                                    !!! ** General constants  **
+   REAL(wp) ::   epsln   = 1.0e-20_wp   ! a small positive number
+   REAL(wp) ::   pthird  = 1._wp/3._wp  ! 1/3
+   REAL(wp) ::   p2third = 2._wp/3._wp  ! 2/3
+
+   INTEGER :: idebug = 236
+   INTEGER :: jdebug = 228
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   INTEGER FUNCTION zdf_osm_alloc()
+      !!----------------------------------------------------------------------
+      !!                 ***  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 )
+     IF( zdf_osm_alloc /= 0 )   CALL ctl_warn('zdf_osm_alloc: failed to allocate zdf_osm arrays')
+     IF( lk_mpp             )   CALL mpp_sum ( zdf_osm_alloc )
+   END FUNCTION zdf_osm_alloc
+
+   SUBROUTINE zdf_osm( kt, p_avm, p_avt )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE zdf_osm  ***
+      !!
+      !! ** Purpose :   Compute the vertical eddy viscosity and diffusivity
+      !!      coefficients and non local mixing using the OSMOSIS scheme
+      !!
+      !! ** Method :   The boundary layer depth hosm is diagnosed at tracer points
+      !!      from profiles of buoyancy, and shear, and the surface forcing.
+      !!      Above hbl (sigma=-z/hbl <1) the mixing coefficients are computed from
+      !!
+      !!                      Kx =  hosm  Wx(sigma) G(sigma)
+      !!
+      !!             and the non local term ghamt = Cs / Ws(sigma) / hosm
+      !!      Below hosm  the coefficients are the sum of mixing due to internal waves
+      !!      shear instability and double diffusion.
+      !!
+      !!      -1- Compute the now interior vertical mixing coefficients at all depths.
+      !!      -2- Diagnose the boundary layer depth.
+      !!      -3- Compute the now boundary layer vertical mixing coefficients.
+      !!      -4- Compute the now vertical eddy vicosity and diffusivity.
+      !!      -5- Smoothing
+      !!
+      !!        N.B. The computation is done from jk=2 to jpkm1
+      !!             Surface value of avt are set once a time to zero
+      !!             in routine zdf_osm_init.
+      !!
+      !! ** Action  :   update the non-local terms ghamts
+      !!                update avt (before vertical eddy coef.)
+      !!
+      !! References : Large W.G., Mc Williams J.C. and Doney S.C.
+      !!         Reviews of Geophysics, 32, 4, November 1994
+      !!         Comments in the code refer to this paper, particularly
+      !!         the equation number. (LMD94, here after)
+      !!----------------------------------------------------------------------
+
+      INTEGER, INTENT( in  )                    ::   kt            ! ocean time step
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::  p_avm, p_avt   ! momentum and tracer Kz (w-points)
+      !!
+      INTEGER ::   ji, jj, jk                   ! dummy loop indices
+      INTEGER ::   ikbot, jkmax, jkm1, jkp2     !
+
+      REAL(wp) ::   ztx, zty, zflageos, zstabl, zbuofdep,zucube      !
+      REAL(wp) ::   zalbet, zbeta, zthermal, zatt1                   !
+      REAL(wp) ::   zehat, zeta, zhrib, zsig, zscale, zwst, zws, zwm ! Velocity scales
+      REAL(wp) ::   zwsun, zwmun, zcons, zconm, zwcons, zwconm       !
+      REAL(wp) ::   zsr, zbw, ze, zb, zd, zc, zaw, za, zb1, za1, zkw, zk0, zcomp , zrhd,zrhdr,zbvzed   ! In situ density
+      INTEGER  ::   jm                          ! dummy loop indices
+      REAL(wp) ::   zr1, zr2, zr3, zr4, zrhop   ! Compression terms
+      REAL(wp) ::   zflag, zrn2, zdep21, zdep32, zdep43
+      REAL(wp) ::   zesh2, zri, zfri            ! Interior richardson mixing
+      REAL(wp) ::   zdelta, zdelta2, zdzup, zdzdn, zdzh, zvath, zgat1, zdat1, zkm1m, zkm1t
+      REAL(wp) :: zt,zs,zu,zv,zrh               ! variables used in constructing averages
+! Scales
+      REAL(wp), DIMENSION(jpi,jpj) :: zrad0     ! Surface solar temperature flux (deg m/s)
+      REAL(wp), DIMENSION(jpi,jpj) :: zradh     ! Radiative flux at bl base (Buoyancy units)
+      REAL(wp), DIMENSION(jpi,jpj) :: zradav    ! Radiative flux, bl average (Buoyancy Units)
+      REAL(wp), DIMENSION(jpi,jpj) :: zustar    ! friction velocity
+      REAL(wp), DIMENSION(jpi,jpj) :: zwstrl    ! Langmuir velocity scale
+      REAL(wp), DIMENSION(jpi,jpj) :: zvstr     ! Velocity scale that ends to zustar for large Langmuir number.
+      REAL(wp), DIMENSION(jpi,jpj) :: zwstrc    ! Convective velocity scale
+      REAL(wp), DIMENSION(jpi,jpj) :: zuw0      ! Surface u-momentum flux
+      REAL(wp), DIMENSION(jpi,jpj) :: zvw0      ! Surface v-momentum flux
+      REAL(wp), DIMENSION(jpi,jpj) :: zwth0     ! Surface heat flux (Kinematic)
+      REAL(wp), DIMENSION(jpi,jpj) :: zws0      ! Surface freshwater flux
+      REAL(wp), DIMENSION(jpi,jpj) :: zwb0      ! Surface buoyancy flux
+      REAL(wp), DIMENSION(jpi,jpj) :: zwthav    ! Heat flux - bl average
+      REAL(wp), DIMENSION(jpi,jpj) :: zwsav     ! freshwater flux - bl average
+      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) :: zustke    ! Surface Stokes drift
+      REAL(wp), DIMENSION(jpi,jpj) :: zla       ! Trubulent Langmuir number
+      REAL(wp), DIMENSION(jpi,jpj) :: zcos_wind ! Cos angle of surface stress
+      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
+
+      ! mixed-layer variables
+
+      INTEGER, DIMENSION(jpi,jpj) :: ibld ! level of boundary layer base
+      INTEGER, DIMENSION(jpi,jpj) :: imld ! level of mixed-layer depth (pycnocline top)
+
+      REAL(wp) :: ztgrad,zsgrad,zbgrad ! Temporary variables used to calculate pycnocline gradients
+      REAL(wp) :: zugrad,zvgrad        ! temporary variables for calculating pycnocline shear
+
+      REAL(wp), DIMENSION(jpi,jpj) :: zhbl  ! bl depth - grid
+      REAL(wp), DIMENSION(jpi,jpj) :: zhml  ! ml depth - 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,jpk) :: zdtdz_pyc    ! parametrized gradient of temperature in pycnocline
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdsdz_pyc    ! parametrised gradient of salinity in pycnocline
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdbdz_pyc    ! parametrised gradient of buoyancy in the pycnocline
+      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
+
+      ! Flux-gradient relationship variables
+
+      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), DIMENSION(jpi,jpj) :: zsc_wth_1,zsc_ws_1 ! Temporary scales used to calculate scalar non-gradient terms.
+      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
+
+      ! For calculating Ri#-dependent mixing
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3du   ! u-shear^2
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3dv   ! v-shear^2
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrimix ! spatial form of ri#-induced diffusion
+
+      ! 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) :: zus_x, zus_y     ! temporary Stokes drift
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zviscos ! viscosity
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity
+
+      ! For debugging
+      INTEGER :: ikt
+
+      !!--------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('zdf_osm')
+      !
+      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
+      zwstrl(:,:) = 0._wp ; zvstr(:,:) = 0._wp ; zwstrc(:,:)    = 0._wp ; zuw0(:,:)      = 0._wp
+      zvw0(:,:)   = 0._wp ; zwth0(:,:) = 0._wp ; zws0(:,:)      = 0._wp ; zwb0(:,:)      = 0._wp
+      zwthav(:,:) = 0._wp ; zwsav(:,:) = 0._wp ; zwbav(:,:)     = 0._wp ; zwb_ent(:,:)   = 0._wp
+      zustke(:,:) = 0._wp ; zla(:,:)   = 0._wp ; zcos_wind(:,:) = 0._wp ; zsin_wind(:,:) = 0._wp
+      zhol(:,:)   = 0._wp
+      lconv(:,:)  = .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
+      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
+      !
+      zdtdz_pyc(:,:,:) = 0._wp ; zdsdz_pyc(:,:,:) = 0._wp ; zdbdz_pyc(:,:,:) = 0._wp
+      zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp
+      !
+      ! 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
+      zhbl_t(:,:)     = 0._wp ; zdhdt(:,:)      = 0._wp
+
+      zdiffut(:,:,:) = 0._wp ; zviscos(:,:,:) = 0._wp ; ghamt(:,:,:) = 0._wp
+      ghams(:,:,:)   = 0._wp ; ghamu(:,:,:)   = 0._wp ; ghamv(:,:,:) = 0._wp
+
+      ! hbl = MAX(hbl,epsln)
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      ! Calculate boundary layer scales
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+      ! Assume two-band radiation model for depth of OSBL
+     zz0 =       rn_abs       ! surface equi-partition in 2-bands
+     zz1 =  1. - rn_abs
+     DO jj = 2, jpjm1
+        DO ji = 2, jpim1
+           ! Surface downward irradiance (so always +ve)
+           zrad0(ji,jj) = qsr(ji,jj) * r1_rau0_rcp
+           ! Downwards irradiance at base of boundary layer
+           zradh(ji,jj) = zrad0(ji,jj) * ( zz0 * EXP( -hbl(ji,jj)/rn_si0 ) + zz1 * EXP( -hbl(ji,jj)/rn_si1) )
+           ! Downwards irradiance averaged over depth of the OSBL
+           zradav(ji,jj) = zrad0(ji,jj) * ( zz0 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si0 ) )*rn_si0 &
+                 &                         + zz1 * ( 1.0 - EXP( -hbl(ji,jj)/rn_si1 ) )*rn_si1 ) / hbl(ji,jj)
+        END DO
+     END DO
+     ! Turbulent surface fluxes and fluxes averaged over depth of the OSBL
+     DO jj = 2, jpjm1
+        DO ji = 2, jpim1
+           zthermal = rab_n(ji,jj,1,jp_tem)
+           zbeta    = rab_n(ji,jj,1,jp_sal)
+           ! Upwards surface Temperature flux for non-local term
+           zwth0(ji,jj) = - qns(ji,jj) * r1_rau0_rcp * tmask(ji,jj,1)
+           ! Upwards surface salinity flux for non-local term
+           zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * tsn(ji,jj,1,jp_sal)  + sfx(ji,jj) ) * r1_rau0 * tmask(ji,jj,1)
+           ! Non radiative upwards surface buoyancy flux
+           zwb0(ji,jj) = grav * zthermal * zwth0(ji,jj) -  grav * zbeta * zws0(ji,jj)
+           ! turbulent heat flux averaged over depth of OSBL
+           zwthav(ji,jj) = 0.5 * zwth0(ji,jj) - ( 0.5*( zrad0(ji,jj) + zradh(ji,jj) ) - zradav(ji,jj) )
+           ! turbulent salinity flux averaged over depth of the OBSL
+           zwsav(ji,jj) = 0.5 * zws0(ji,jj)
+           ! turbulent buoyancy flux averaged over the depth of the OBSBL
+           zwbav(ji,jj) = grav  * zthermal * zwthav(ji,jj) - grav  * zbeta * zwsav(ji,jj)
+           ! Surface upward velocity fluxes
+           zuw0(ji,jj) = -utau(ji,jj) * r1_rau0 * tmask(ji,jj,1)
+           zvw0(ji,jj) = -vtau(ji,jj) * r1_rau0 * 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 )
+           zcos_wind(ji,jj) = -zuw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) )
+           zsin_wind(ji,jj) = -zvw0(ji,jj) / ( zustar(ji,jj) * zustar(ji,jj) )
+        END DO
+     END DO
+     ! Calculate Stokes drift in direction of wind (zustke) and Stokes penetration depth (dstokes)
+     SELECT CASE (nn_osm_wave)
+     ! Assume constant La#=0.3
+     CASE(0)
+        DO jj = 2, jpjm1
+           DO ji = 2, jpim1
+              zus_x = zcos_wind(ji,jj) * zustar(ji,jj) / 0.3**2
+              zus_y = zsin_wind(ji,jj) * zustar(ji,jj) / 0.3**2
+              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
+           END DO
+        END DO
+     ! Assume Pierson-Moskovitz wind-wave spectrum
+     CASE(1)
+        DO jj = 2, jpjm1
+           DO ji = 2, jpim1
+              ! 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
+           END DO
+        END DO
+     ! Use ECMWF wave fields as output from SBCWAVE
+     CASE(2)
+        zfac =  2.0_wp * rpi / 16.0_wp
+        DO jj = 2, jpjm1
+           DO ji = 2, jpim1
+              ! 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 !
+           END DO
+        END DO
+     END SELECT
+
+     ! Langmuir velocity scale (zwstrl), La # (zla)
+     ! mixed scale (zvstr), convective velocity scale (zwstrc)
+     DO jj = 2, jpjm1
+        DO ji = 2, jpim1
+           ! 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
+           ! Velocity scale that tends to zustar for large Langmuir numbers
+           zvstr(ji,jj) = ( zwstrl(ji,jj)**3  + &
+                & ( 1.0 - EXP( -0.5 * zla(ji,jj)**2 ) ) * zustar(ji,jj) * zustar(ji,jj) * zustar(ji,jj) )**pthird
+
+           ! 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
+           IF ( zwbav(ji,jj) > 0.0) THEN
+              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
+              zhol(ji,jj) = -hbl(ji,jj) *  2.0 * zwbav(ji,jj)/ (zvstr(ji,jj)**3  + epsln )
+              lconv(ji,jj) = .FALSE.
+           ENDIF
+        END DO
+     END DO
+
+     !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+     ! 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_n(:,:,3) )
+      ibld(:,:) = 3
+      DO jk = 4, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               IF ( hbl(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN
+                  ibld(ji,jj) = MIN(mbkt(ji,jj), jk)
+               ENDIF
+            END DO
+         END DO
+      END DO
+
+      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)
+               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_n(ji,jj,jm)
+                  zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem)
+                  zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal)
+                  zu   = zu  + e3t_n(ji,jj,jm) &
+                     &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) &
+                     &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) )
+                  zv   = zv  + e3t_n(ji,jj,jm) &
+                     &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) &
+                     &            / 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) - tsn(ji,jj,ibld(ji,jj),jp_tem)
+               zds_bl(ji,jj) = zs_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal)
+               zdu_bl(ji,jj) = zu_bl(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) &
+                     &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) )
+               zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) &
+                     &   / 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_rdt / 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_rdt / 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_rdt / hbl(ji,jj) ) * zwstrl(ji,jj)**2 / hbl(ji,jj)
+                      ELSE
+                          zpert = 2.0 * ( 1.0 + 2.0 * zwstrl(ji,jj) * rn_rdt / 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
+           END DO
+      END DO
+
+      ! Calculate averages over depth of boundary layer
+      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
+
+      DO jk = 4, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               IF ( zhbl_t(ji,jj) >= gdepw_n(ji,jj,jk) ) THEN
+                  ibld(ji,jj) =  MIN(mbkt(ji,jj), jk)
+               ENDIF
+            END DO
+         END DO
+      END DO
+
+!
+! Step through model levels taking account of buoyancy change to determine the effect on dhdt
+!
+      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
+                 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
+
+                 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
+
+                    zhbl_s = zhbl_s + MIN( - zwb_ent(ji,jj) / zdb * rn_rdt / FLOAT(ibld(ji,jj)-imld(ji,jj) ), e3w_n(ji,jj,jk) )
+                    zhbl_s = MIN(zhbl_s, ht_n(ji,jj))
+
+                    IF ( zhbl_s >= gdepw_n(ji,jj,jm+1) ) 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) - tsn(ji,jj,jm,jp_tem) ) - zbeta * ( zs_bl(ji,jj) - tsn(ji,jj,jm,jp_sal) ) ), 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_n(ji,jj,jk) / zdhdt(ji,jj)  ! ALMG to investigate whether need to include wn here
+
+                    zhbl_s = MIN(zhbl_s, ht_n(ji,jj))
+                    IF ( zhbl_s >= gdepw_n(ji,jj,jm) ) jm = jm + 1
+                 END DO
+                 hbl(ji,jj) = MAX(zhbl_s, gdepw_n(ji,jj,3) )
+                 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_n(ji,jj,3) )
+                 IF ( hbl(ji,jj) > hbli(ji,jj) ) hbli(ji,jj) = hbl(ji,jj)
+             ENDIF
+           ENDIF
+           zhbl(ji,jj) = gdepw_n(ji,jj,ibld(ji,jj))
+         END DO
+      END DO
+      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 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)
+               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_n(ji,jj,jm)
+                  zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem)
+                  zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal)
+                  zu   = zu  + e3t_n(ji,jj,jm) &
+                     &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) &
+                     &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) )
+                  zv   = zv  + e3t_n(ji,jj,jm) &
+                     &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) &
+                     &            / 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) - tsn(ji,jj,ibld(ji,jj),jp_tem)
+               zds_bl(ji,jj) = zs_bl(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal)
+               zdu_bl(ji,jj) = zu_bl(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) &
+                      &   / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) )
+               zdv_bl(ji,jj) = zv_bl(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) &
+                      &  / 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_n(ji,jj,ibld(ji,jj))
+               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_n(ji,jj,ibld(ji,jj)) ) , 1 )
+                     imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1)
+                     zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))
+                     zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
+                  ELSE  ! IF (zdb_bl)
+                     imld(ji,jj) = ibld(ji,jj) - 1
+                     zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))
+                     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_n(ji,jj,ibld(ji,jj)) ) , 1 )
+                        imld(ji,jj) = MAX( ibld(ji,jj) - inhml, 1)
+                        zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))
+                        zdh(ji,jj) = zhbl(ji,jj) - zhml(ji,jj)
+                     ELSE
+                        imld(ji,jj) = ibld(ji,jj) - 1
+                        zhml(ji,jj) = gdepw_n(ji,jj,imld(ji,jj))
+                        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 DO
+      END DO
+
+      ! Average over the depth of the mixed layer in the convective boundary layer
+      ! Also calculate entrainment fluxes for temperature and salinity
+      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)
+            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_n(ji,jj,jm)
+                  zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem)
+                  zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal)
+                  zu   = zu  + e3t_n(ji,jj,jm) &
+                     &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) &
+                     &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) )
+                  zv   = zv  + e3t_n(ji,jj,jm) &
+                     &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) &
+                     &            / 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) - tsn(ji,jj,ibld(ji,jj),jp_tem)
+               zds_ml(ji,jj) = zs_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal)
+               zdu_ml(ji,jj) = zu_ml(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) &
+                     &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) )
+               zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) &
+                     &    / 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_n(ji,jj,jm)
+                     zt   = zt  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_tem)
+                     zs   = zs  + e3t_n(ji,jj,jm) * tsn(ji,jj,jm,jp_sal)
+                     zu   = zu  + e3t_n(ji,jj,jm) &
+                        &            * ( ub(ji,jj,jm) + ub(ji - 1,jj,jm) ) &
+                        &            / MAX( 1. , umask(ji,jj,jm) + umask(ji - 1,jj,jm) )
+                     zv   = zv  + e3t_n(ji,jj,jm) &
+                        &            * ( vb(ji,jj,jm) + vb(ji,jj - 1,jm) ) &
+                        &            / 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) - tsn(ji,jj,ibld(ji,jj),jp_tem)
+                  zds_ml(ji,jj) = zs_ml(ji,jj) - tsn(ji,jj,ibld(ji,jj),jp_sal)
+                  zdu_ml(ji,jj) = zu_ml(ji,jj) - ( ub(ji,jj,ibld(ji,jj)) + ub(ji-1,jj,ibld(ji,jj) ) ) &
+                        &    / MAX(1. , umask(ji,jj,ibld(ji,jj) ) + umask(ji-1,jj,ibld(ji,jj) ) )
+                  zdv_ml(ji,jj) = zv_ml(ji,jj) - ( vb(ji,jj,ibld(ji,jj)) + vb(ji,jj-1,ibld(ji,jj) ) ) &
+                        &    / 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 DO
+      END DO
+    !
+    ! rotate mean currents and changes onto wind align co-ordinates
+    !
+
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1
+            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 DO
+      END DO
+
+     zuw_bse = 0._wp
+     zvw_bse = 0._wp
+     DO jj = 2, jpjm1
+        DO ji = 2, jpim1
+
+           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 DO
+     END DO
+
+      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+      !  Pycnocline gradients for scalars and velocity
+      !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+       DO jj = 2, jpjm1
+          DO ji = 2, jpim1
+          !
+             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_n(ji,jj,jk) - 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_n(ji,jj,jk) / 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_n(ji,jj,jk) - 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 DO
+       END DO
+!
+       DO jj = 2, jpjm1
+          DO ji = 2, jpim1
+          !
+             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_n(ji,jj,jk) - 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_n(ji,jj,jk) / 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 DO
+       END DO
+       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+       ! 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 jj = 2, jpjm1
+          DO ji = 2, jpim1
+             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 DO
+    END DO
+!
+       DO jj = 2, jpjm1
+          DO ji = 2, jpim1
+             IF ( lconv(ji,jj) ) THEN
+                DO jk = 2, imld(ji,jj)   ! mixed layer diffusivity
+                    zznd_ml = gdepw_n(ji,jj,jk) / 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_n(ji,jj,jk) - 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 wn taken out
+                zdiffut(ji,jj,ibld(ji,jj)) = zdhdt(ji,jj)* e3t_n(ji,jj,ibld(ji,jj))
+                ! 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_n(ji,jj,jk) / 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 DO  ! end of ji loop
+       END DO     ! end of jj loop
+
+       !
+       ! calculate non-gradient components of the flux-gradient relationships
+       !
+! Stokes term in scalar flux, flux-gradient relationship
+       WHERE ( lconv )
+          zsc_wth_1 = zwstrl**3 * zwth0 / ( zvstr**3 + 0.5 * zwstrc**3 + epsln)
+          !
+          zsc_ws_1 = zwstrl**3 * zws0 / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
+       ELSEWHERE
+          zsc_wth_1 = 2.0 * zwthav
+          !
+          zsc_ws_1 = 2.0 * zwsav
+       ENDWHERE
+
+
+       DO jj = 2, jpjm1
+          DO ji = 2, jpim1
+            IF ( lconv(ji,jj) ) THEN
+              DO jk = 2, imld(ji,jj)
+                 zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
+                 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) +  1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_wth_1(ji,jj)
+                 !
+                 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_ws_1(ji,jj)
+              END DO ! end jk loop
+            ELSE     ! else for if (lconv)
+ ! Stable conditions
+               DO jk = 2, ibld(ji,jj)
+                  zznd_d=gdepw_n(ji,jj,jk) / dstokes(ji,jj)
+                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_wth_1(ji,jj)
+                  !
+                  ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_ws_1(ji,jj)
+               END DO
+            ENDIF               ! endif for check on lconv
+
+          END DO  ! end of ji loop
+       END DO     ! end of jj loop
+
+
+! 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 )
+       ELSEWHERE
+          zsc_uw_1 = zustar**2
+          zsc_vw_1 = ff_t * zhbl * zustke**3 * zla**(8.0/3.0) / (zvstr**2 + epsln)
+       ENDWHERE
+
+       DO jj = 2, jpjm1
+          DO ji = 2, jpim1
+             IF ( lconv(ji,jj) ) THEN
+                DO jk = 2, imld(ji,jj)
+                   zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
+                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + ( -0.05 * EXP ( -0.4 * zznd_d ) * zsc_uw_1(ji,jj) + 0.00125 * EXP ( -zznd_d ) * zsc_uw_2(ji,jj) ) * ( 1.0 - EXP ( -2.0 * zznd_d ) )
+!
+                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.65 * 0.15 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_vw_1(ji,jj)
+                END DO   ! end jk loop
+             ELSE
+! Stable conditions
+                DO jk = 2, ibld(ji,jj) ! corrected to ibld
+                   zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
+                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75 * 1.3 * EXP ( -0.5 * zznd_d ) * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_uw_1(ji,jj)
+                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0._wp
+                END DO   ! end jk loop
+             ENDIF
+          END DO  ! ji loop
+       END DO     ! jj loo
+
+! Buoyancy term in flux-gradient relationship [note : includes ROI ratio (X0.3) and pressure (X0.5)]
+
+       WHERE ( lconv )
+          zsc_wth_1 = zwbav * zwth0 * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
+          zsc_ws_1 = zwbav * zws0 * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
+       ELSEWHERE
+          zsc_wth_1 = 0._wp
+          zsc_ws_1 = 0._wp
+       ENDWHERE
+
+       DO jj = 2, jpjm1
+          DO ji = 2, jpim1
+             IF (lconv(ji,jj) ) THEN
+                DO jk = 2, imld(ji,jj)
+                   zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj)
+            ! calculate turbulent length scale
+                   zl_c = 0.9 * ( 1.0 - EXP ( -7.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) ) * ( 1.0 - EXP ( -15.0 * ( 1.1 - zznd_ml ) ) )
+                   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)
+            ! 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
+                DO jk = 2, ibld(ji,jj)
+                   ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj)
+                   ghams(ji,jj,jk) = ghams(ji,jj,jk) + zsc_ws_1(ji,jj)
+                END DO
+             ENDIF
+          END DO   ! ji loop
+       END DO      ! jj loop
+
+
+       WHERE ( lconv )
+          zsc_uw_1 = -zwb0 * zustar**2 * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
+          zsc_uw_2 = zwb0 *  zustke * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )**(2.0/3.0)
+          zsc_vw_1 = 0._wp
+       ELSEWHERE
+         zsc_uw_1 = 0._wp
+         zsc_vw_1 = 0._wp
+       ENDWHERE
+
+       DO jj = 2, jpjm1
+          DO ji = 2, jpim1
+             IF ( lconv(ji,jj) ) THEN
+                DO jk = 2 , imld(ji,jj)
+                   zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
+                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * 0.5 * ( zsc_uw_1(ji,jj) + 0.125 * EXP( -0.5 * zznd_d ) * ( 1.0 - EXP ( -0.5 * zznd_d) ) * zsc_uw_2(ji,jj) )
+                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj)
+                END DO  ! jk loop
+             ELSE
+             ! stable conditions
+                DO jk = 2, ibld(ji,jj)
+                   ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj)
+                   ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj)
+                END DO
+             ENDIF
+          END DO        ! ji loop
+       END DO           ! jj loop
+
+! 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 jj = 2, jpjm1
+          DO ji = 2, jpim1
+            IF ( lconv(ji,jj) ) THEN
+               DO jk = 2, imld(ji,jj)
+                  zznd_ml=gdepw_n(ji,jj,jk) / zhml(ji,jj)
+                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * zsc_wth_1(ji,jj) * ( -2.0 + 2.75 * ( ( 1.0 + 0.6 * zznd_ml**4 ) - EXP ( -6.0 *zznd_ml ) ) )  &
+      &           * ( 1.0 -EXP ( -15.0 * ( 1.0 - zznd_ml ) ) )
+                  !
+                  ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * zsc_ws_1(ji,jj) * ( -2.0 + 2.75 * ( ( 1.0 + 0.6 * zznd_ml**4 ) - EXP ( -6.0 *zznd_ml ) ) )  &
+      &           * ( 1.0 -EXP ( -15.0 * ( 1.0 - zznd_ml ) ) )
+               END DO
+            ELSE
+               DO jk = 2, ibld(ji,jj)
+                  zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
+                  znd = gdepw_n(ji,jj,jk) / 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
+          ENDDO    ! ji loop
+       END DO      ! jj loop
+
+
+       WHERE ( lconv )
+          zsc_uw_1 = zustar**2
+          zsc_vw_1 = ff_t * zustke * zhml
+       ELSEWHERE
+          zsc_uw_1 = zustar**2
+          zsc_uw_2 = (2.25 - 3.0 * ( 1.0 - EXP( -1.25 * 2.0 ) ) ) * ( 1.0 - EXP( -4.0 * 2.0 ) ) * zsc_uw_1
+          zsc_vw_1 = ff_t * zustke * zhbl
+          zsc_vw_2 = -0.11 * SIN( 3.14159 * ( 2.0 + 0.4 ) ) * EXP(-( 1.5 + 2.0 )**2 ) * zsc_vw_1
+       ENDWHERE
+
+       DO jj = 2, jpjm1
+          DO ji = 2, jpim1
+             IF ( lconv(ji,jj) ) THEN
+               DO jk = 2, imld(ji,jj)
+                  zznd_ml = gdepw_n(ji,jj,jk) / zhml(ji,jj)
+                  zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
+                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk)&
+                       & + 0.3 * ( -2.0 + 2.5 * ( 1.0 + 0.1 * zznd_ml**4 ) - EXP ( -8.0 * zznd_ml ) ) * zsc_uw_1(ji,jj)
+                  !
+                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
+                       & + 0.3 * 0.1 * ( EXP( -zznd_d ) + EXP( -5.0 * ( 1.0 - zznd_ml ) ) ) * zsc_vw_1(ji,jj)
+               END DO
+             ELSE
+               DO jk = 2, ibld(ji,jj)
+                  znd = gdepw_n(ji,jj,jk) / zhbl(ji,jj)
+                  zznd_d = gdepw_n(ji,jj,jk) / dstokes(ji,jj)
+                  IF ( zznd_d <= 2.0 ) THEN
+                     ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5 * 0.3 &
+                          &*  ( 2.25 - 3.0  * ( 1.0 - EXP( - 1.25 * zznd_d ) ) * ( 1.0 - EXP( -2.0 * zznd_d ) ) ) * zsc_uw_1(ji,jj)
+                     !
+                  ELSE
+                     ghamu(ji,jj,jk) = ghamu(ji,jj,jk)&
+                          & + 0.5 * 0.3 * ( 1.0 - EXP( -5.0 * ( 1.0 - znd ) ) ) * zsc_uw_2(ji,jj)
+                     !
+                  ENDIF
+
+                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
+                       & + 0.3 * 0.15 * SIN( 3.14159 * ( 0.65 * zznd_d ) ) * EXP( -0.25 * zznd_d**2 ) * zsc_vw_1(ji,jj)
+                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
+                       & + 0.3 * 0.15 * EXP( -5.0 * ( 1.0 - znd ) ) * ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj)
+               END DO
+             ENDIF
+          END DO
+       END DO
+!
+! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer.
+
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1
+            IF ( lconv(ji,jj) ) THEN
+               DO jk = 2, ibld(ji,jj)
+                  znd = ( gdepw_n(ji,jj,jk) - 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_n(ji,jj,jk) - 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( -10.0 * znd**2 ) )
+                     ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) )
+                  ELSE
+                     ghamu(ji,jj,jk) = 0._wp
+                     ghamv(ji,jj,jk) = 0._wp
+                  ENDIF
+               END DO
+            ENDIF
+         END DO
+      END DO
+
+      ! 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 jj = 2, jpjm1
+          DO ji = 2, jpim1
+             DO jk= 2, ibld(ji,jj)
+                znd = gdepw_n(ji,jj,jk) / 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
+           END DO
+       END DO
+
+! Entrainment contribution.
+
+       DO jj=2, jpjm1
+          DO ji = 2, jpim1
+             IF ( lconv(ji,jj) ) THEN
+               DO jk = 1, imld(ji,jj) - 1
+                  znd=gdepw_n(ji,jj,jk) / 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_n(ji,jj,jk) - 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
+          END DO       ! ji loop
+       END DO          ! jj loop
+
+
+       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+       ! Need to put in code for contributions that are applied explicitly to
+       ! the prognostic variables
+       !  1. Entrainment flux
+       !
+       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
+
+
+
+       ! rotate non-gradient velocity terms back to model reference frame
+
+       DO jj = 2, jpjm1
+          DO ji = 2, jpim1
+             DO jk = 2, ibld(ji,jj)
+                ztemp = ghamu(ji,jj,jk)
+                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * zcos_wind(ji,jj) - ghamv(ji,jj,jk) * zsin_wind(ji,jj)
+                ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * zcos_wind(ji,jj) + ztemp * zsin_wind(ji,jj)
+             END DO
+          END DO
+       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
+          DO jk = 2, jpkm1           !* Shear production at uw- and vw-points (energy conserving form)
+             DO jj = 1, jpjm1
+                DO ji = 1, jpim1   ! vector opt.
+                   z3du(ji,jj,jk) = 0.5 * (  un(ji,jj,jk-1) -  un(ji  ,jj,jk) )   &
+                        &                 * (  ub(ji,jj,jk-1) -  ub(ji  ,jj,jk) ) * wumask(ji,jj,jk) &
+                        &                 / (  e3uw_n(ji,jj,jk) * e3uw_b(ji,jj,jk) )
+                   z3dv(ji,jj,jk) = 0.5 * (  vn(ji,jj,jk-1) -  vn(ji,jj  ,jk) )   &
+                        &                 * (  vb(ji,jj,jk-1) -  vb(ji,jj  ,jk) ) * wvmask(ji,jj,jk) &
+                        &                 / (  e3vw_n(ji,jj,jk) * e3vw_b(ji,jj,jk) )
+                END DO
+             END DO
+          END DO
+      !
+         DO jk = 2, jpkm1
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1   ! vector opt.
+                  !                                          ! shear prod. at w-point weightened by mask
+                  zesh2  =  ( z3du(ji-1,jj,jk) + z3du(ji,jj,jk) ) / MAX( 1._wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
+                     &    + ( z3dv(ji,jj-1,jk) + z3dv(ji,jj,jk) ) / MAX( 1._wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
+                  !                                          ! local Richardson number
+                  zri   = MAX( rn2b(ji,jj,jk), 0._wp ) / MAX(zesh2, epsln)
+                  zfri =  MIN( zri / rn_riinfty , 1.0_wp )
+                  zfri  = ( 1.0_wp - zfri * zfri )
+                  zrimix(ji,jj,jk)  =  zfri * zfri  * zfri * wmask(ji, jj, jk)
+                END DO
+             END DO
+          END DO
+
+          DO jj = 2, jpjm1
+             DO ji = 2, jpim1
+                DO jk = ibld(ji,jj) + 1, jpkm1
+                   zdiffut(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
+                   zviscos(ji,jj,jk) = zrimix(ji,jj,jk)*rn_difri
+                END DO
+             END DO
+          END DO
+
+       END IF ! ln_kpprimix = .true.
+
+! KPP-style set diffusivity large if unstable below BL
+       IF( ln_convmix) THEN
+          DO jj = 2, jpjm1
+             DO ji = 2, jpim1
+                DO jk = ibld(ji,jj) + 1, jpkm1
+                  IF(  MIN( rn2(ji,jj,jk), rn2b(ji,jj,jk) ) <= -1.e-12 ) zdiffut(ji,jj,jk) = rn_difconv
+                END DO
+             END DO
+          END DO
+       END IF ! ln_convmix = .true.
+
+       ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids
+       CALL lbc_lnk( zviscos(:,:,:), 'W', 1. )
+
+       ! GN 25/8: need to change tmask --> wmask
+
+     DO jk = 2, jpkm1
+         DO jj = 2, jpjm1
+             DO ji = 2, jpim1
+                p_avt(ji,jj,jk) = MAX( zdiffut(ji,jj,jk), avtb(jk) ) * tmask(ji,jj,jk)
+                p_avm(ji,jj,jk) = MAX( zviscos(ji,jj,jk), avmb(jk) ) * tmask(ji,jj,jk)
+             END DO
+         END DO
+     END DO
+
+
+     CALL lbc_lnk( p_avt, 'W', 1. )                      ! Lateral boundary conditions on p_avt  (sign unchanged)
+     CALL lbc_lnk( p_avm, 'W', 1. )                      ! Lateral boundary conditions on p_avm  (sign unchanged)
+
+      ! Lateral boundary conditions on ghamu and ghamv, currently on W-grid  (sign unchanged), needed to caclulate gham[uv] on u and v grids
+      CALL lbc_lnk( ghamu(:,:,:), 'W', 1. )
+      CALL lbc_lnk( ghamv(:,:,:), 'W', 1. )
+       DO jk = 2, jpkm1
+           DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  ghamu(ji,jj,jk) = ( ghamu(ji,jj,jk) + ghamu(ji+1,jj,jk) ) &
+                     &  / MAX( 1., tmask(ji,jj,jk) + tmask (ji + 1,jj,jk) ) * umask(ji,jj,jk)
+
+                  ghamv(ji,jj,jk) = ( ghamv(ji,jj,jk) + ghamv(ji,jj+1,jk) ) &
+                      &  / MAX( 1., tmask(ji,jj,jk) + tmask (ji,jj+1,jk) ) * vmask(ji,jj,jk)
+
+                  ghamt(ji,jj,jk) =  ghamt(ji,jj,jk) * tmask(ji,jj,jk)
+                  ghams(ji,jj,jk) =  ghams(ji,jj,jk) * tmask(ji,jj,jk)
+               END DO
+           END DO
+        END DO
+        ! Lateral boundary conditions on final outputs for gham[ts],  on W-grid  (sign unchanged)
+       CALL lbc_lnk( ghamt(:,:,:), 'W', 1. )
+       CALL lbc_lnk( ghams(:,:,:), 'W', 1. )
+      ! Lateral boundary conditions on final outputs for gham[uv],  on [UV]-grid  (sign unchanged)
+       CALL lbc_lnk( ghamu(:,:,:), 'U', 1. )
+       CALL lbc_lnk( ghamv(:,:,:), 'V', 1. )
+
+       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).
+         CASE(0:1)
+            IF ( iom_use("us_x") ) CALL iom_put( "us_x", tmask(:,:,1)*zustke*zcos_wind )   ! x surface Stokes drift
+            IF ( iom_use("us_y") ) CALL iom_put( "us_y", tmask(:,:,1)*zustke*zsin_wind )  ! y surface Stokes drift
+            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*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
+            IF ( iom_use("wind_wave_abs_power") ) CALL iom_put( "wind_wave_abs_power", 1000.*rau0*tmask(:,:,1)*zustar**2* &
+                 & SQRT(ut0sd**2 + vt0sd**2 ) )
+         END SELECT
+         IF ( iom_use("ghamt") ) CALL iom_put( "ghamt", tmask*ghamt )            ! <Tw_NL>
+         IF ( iom_use("ghams") ) CALL iom_put( "ghams", tmask*ghams )            ! <Sw_NL>
+         IF ( iom_use("ghamu") ) CALL iom_put( "ghamu", umask*ghamu )            ! <uw_NL>
+         IF ( iom_use("ghamv") ) CALL iom_put( "ghamv", vmask*ghamv )            ! <vw_NL>
+         IF ( iom_use("zwth0") ) CALL iom_put( "zwth0", tmask(:,:,1)*zwth0 )            ! <Tw_0>
+         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("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
+         IF ( iom_use("zwstrc") ) CALL iom_put( "zwstrc", tmask(:,:,1)*zwstrc )         ! convective velocity scale
+         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("wind_power") ) CALL iom_put( "wind_power", 1000.*rau0*tmask(:,:,1)*zustar**3 ) ! BL depth internal to zdf_osm routine
+         IF ( iom_use("wind_wave_power") ) CALL iom_put( "wind_wave_power", 1000.*rau0*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("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
+       END IF
+         ! Lateral boundary conditions on p_avt  (sign unchanged)
+         CALL lbc_lnk( p_avt(:,:,:), 'W', 1. )
+
+   END SUBROUTINE zdf_osm
+
+   SUBROUTINE zdf_osm_init
+     !!----------------------------------------------------------------------
+     !!                  ***  ROUTINE zdf_osm_init  ***
+     !!
+     !! ** Purpose :   Initialization of the vertical eddy diffivity and
+     !!      viscosity when using a osm turbulent closure scheme
+     !!
+     !! ** Method  :   Read the namosm namelist and check the parameters
+     !!      called at the first timestep (nit000)
+     !!
+     !! ** input   :   Namlist namosm
+     !!----------------------------------------------------------------------
+     INTEGER  ::   ios            ! local integer
+     INTEGER  ::   ji, jj, jk     ! dummy loop indices
+     !!
+     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
+     !!----------------------------------------------------------------------
+     !
+     IF( nn_timing == 1 )  CALL timing_start('zdf_osm_init')
+     !
+     REWIND( numnam_ref )              ! Namelist namzdf_osm in reference namelist : Osmosis ML model
+     READ  ( numnam_ref, namzdf_osm, IOSTAT = ios, ERR = 901)
+901    IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_osm in reference namelist', lwp )
+
+     REWIND( numnam_cfg )              ! Namelist namzdf_tke in configuration namelist : Turbulent Kinetic Energy
+     READ  ( numnam_cfg, namzdf_osm, IOSTAT = ios, ERR = 902 )
+902    IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_osm in configuration namelist', lwp )
+     IF(lwm) WRITE ( numond, namzdf_osm )
+
+     IF(lwp) THEN                    ! Control print
+        WRITE(numout,*)
+        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,*) '     Turbulent Langmuir number                     rn_osm_la   = ', rn_osm_la
+        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,*) '     horizontal average flag                       nn_ave      = ', nn_ave
+        WRITE(numout,*) '     Stokes drift                                  nn_osm_wave = ', nn_osm_wave
+        SELECT CASE (nn_osm_wave)
+        CASE(0)
+           WRITE(numout,*) '     calculated assuming constant La#=0.3'
+        CASE(1)
+           WRITE(numout,*) '     calculated from Pierson Moskowitz wind-waves'
+        CASE(2)
+           WRITE(numout,*) '     calculated from ECMWF wave fields'
+        END SELECT
+        WRITE(numout,*) '     Output osm diagnostics                       ln_dia_osm  = ',  ln_dia_osm
+        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
+        WRITE(numout,*) '     maximum shear diffusivity at Rig = 0    (m2/s) rn_difri = ', rn_difri
+        WRITE(numout,*) '     Use large mixing below BL when unstable       ln_convmix = ', ln_convmix
+        WRITE(numout,*) '     diffusivity when unstable below BL     (m2/s) rn_difconv = ', rn_difconv
+     ENDIF
+
+     !                              ! allocate zdfosm arrays
+     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_zdfddm) THEN
+        IF(lwp) THEN
+           WRITE(numout,*)
+           WRITE(numout,*) '    Double diffusion mixing on temperature and salinity '
+           WRITE(numout,*) '    CAUTION : done in routine zdfosm, not in routine zdfddm '
+        ENDIF
+     ENDIF
+
+
+     !set constants not in namelist
+     !-----------------------------
+
+     IF(lwp) THEN
+        WRITE(numout,*)
+     ENDIF
+
+     IF (nn_osm_wave == 0) THEN
+        dstokes(:,:) = rn_osm_dstokes
+     END IF
+
+     ! Horizontal average : initialization of weighting arrays
+     ! -------------------
+
+     SELECT CASE ( nn_ave )
+
+     CASE ( 0 )                ! no horizontal average
+        IF(lwp) WRITE(numout,*) '          no horizontal average on avt'
+        IF(lwp) WRITE(numout,*) '          only in very high horizontal resolution !'
+        ! weighting mean arrays etmean
+        !           ( 1  1 )
+        ! avt = 1/4 ( 1  1 )
+        !
+        etmean(:,:,:) = 0.e0
+
+        DO jk = 1, jpkm1
+           DO jj = 2, jpjm1
+              DO ji = 2, jpim1   ! vector opt.
+                 etmean(ji,jj,jk) = tmask(ji,jj,jk)                     &
+                      &  / MAX( 1.,  umask(ji-1,jj  ,jk) + umask(ji,jj,jk)   &
+                      &            + vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)  )
+              END DO
+           END DO
+        END DO
+
+     CASE ( 1 )                ! horizontal average
+        IF(lwp) WRITE(numout,*) '          horizontal average on avt'
+        ! weighting mean arrays etmean
+        !           ( 1/2  1  1/2 )
+        ! avt = 1/8 ( 1    2  1   )
+        !           ( 1/2  1  1/2 )
+        etmean(:,:,:) = 0.e0
+
+        DO jk = 1, jpkm1
+           DO jj = 2, jpjm1
+              DO ji = 2, jpim1   ! vector opt.
+                 etmean(ji,jj,jk) = tmask(ji, jj,jk)                           &
+                      & / MAX( 1., 2.* tmask(ji,jj,jk)                           &
+                      &      +.5 * ( tmask(ji-1,jj+1,jk) + tmask(ji-1,jj-1,jk)   &
+                      &             +tmask(ji+1,jj+1,jk) + tmask(ji+1,jj-1,jk) ) &
+                      &      +1. * ( tmask(ji-1,jj  ,jk) + tmask(ji  ,jj+1,jk)   &
+                      &             +tmask(ji  ,jj-1,jk) + tmask(ji+1,jj  ,jk) ) )
+              END DO
+           END DO
+        END DO
+
+     CASE DEFAULT
+        WRITE(ctmp1,*) '          bad flag value for nn_ave = ', nn_ave
+        CALL ctl_stop( ctmp1 )
+
+     END SELECT
+
+     ! Initialization of vertical eddy coef. to the background value
+     ! -------------------------------------------------------------
+     DO jk = 1, jpk
+        avt (:,:,jk) = avtb(jk) * tmask(:,:,jk)
+     END DO
+
+     ! zero the surface flux for non local term and osm mixed layer depth
+     ! ------------------------------------------------------------------
+     ghamt(:,:,:) = 0.
+     ghams(:,:,:) = 0.
+     ghamu(:,:,:) = 0.
+     ghamv(:,:,:) = 0.
+
+     !
+     IF( nn_timing == 1 )  CALL timing_stop('zdf_osm_init')
+     !
+   END SUBROUTINE zdf_osm_init
+
+   SUBROUTINE osm_rst( kt, cdrw )
+     !!---------------------------------------------------------------------
+     !!                   ***  ROUTINE osm_rst  ***
+     !!
+     !! ** Purpose :   Read or write BL fields in restart file
+     !!
+     !! ** Method  :   use of IOM library. If the restart does not contain
+     !!                required fields, they are recomputed from stratification
+     !!----------------------------------------------------------------------
+
+     INTEGER, INTENT(in) :: kt
+     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
+
+     INTEGER ::   id1, id2   ! iom enquiry index
+     INTEGER  ::   ji, jj, jk     ! dummy loop indices
+     INTEGER  ::   iiki, ikt ! local integer
+     REAL(wp) ::   zhbf           ! tempory scalars
+     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)
+     !!----------------------------------------------------------------------
+     !
+     !!-----------------------------------------------------------------------------
+     ! If READ/WRITE Flag is 'READ', try to get hbl from restart file. If successful then return
+     !!-----------------------------------------------------------------------------
+     IF( TRIM(cdrw) == 'READ'.AND. ln_rstart) THEN
+        id1 = iom_varid( numror, 'wn'   , ldstop = .FALSE. )
+        IF( id1 > 0 ) THEN                       ! 'wn' exists; read
+           CALL iom_get( numror, jpdom_autoglo, 'wn', wn )
+           WRITE(numout,*) ' ===>>>> :  wn read from restart file'
+        ELSE
+           wn(:,:,:) = 0._wp
+           WRITE(numout,*) ' ===>>>> :  wn not in restart file, set to zero initially'
+        END IF
+        id1 = iom_varid( numror, 'hbl'   , ldstop = .FALSE. )
+        id2 = iom_varid( numror, 'hbli'   , ldstop = .FALSE. )
+        IF( id1 > 0 .AND. id2 > 0) THEN                       ! 'hbl' exists; read and return
+           CALL iom_get( numror, jpdom_autoglo, 'hbl', hbl )
+           CALL iom_get( numror, jpdom_autoglo, 'hbli', hbli )
+           WRITE(numout,*) ' ===>>>> :  hbl & hbli read from restart file'
+           RETURN
+        ELSE                      ! 'hbl' & 'hbli' not in restart file, recalculate
+           WRITE(numout,*) ' ===>>>> : previous run without osmosis scheme, hbl computed from stratification'
+        END IF
+     END IF
+
+     !!-----------------------------------------------------------------------------
+     ! 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(lwp) WRITE(numout,*) '---- osm-rst ----'
+         CALL iom_rstput( kt, nitrst, numrow, 'wn'   , wn     )
+         CALL iom_rstput( kt, nitrst, numrow, 'hbl'   , hbl     )
+         CALL iom_rstput( kt, nitrst, numrow, 'hbli'   , hbli     )
+        RETURN
+     END IF
+
+     !!-----------------------------------------------------------------------------
+     ! Getting hbl, no restart file with hbl, so calculate from surface stratification
+     !!-----------------------------------------------------------------------------
+     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( tsn, rab_n )
+     CALL bn2(tsn, rab_n, rn2)
+     imld_rst(:,:)  = nlb10         ! Initialization to the number of w ocean point
+     hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2
+     zN2_c = grav * rho_c * r1_rau0 ! convert density criteria into N^2 criteria
+     !
+     hbl(:,:)  = 0._wp              ! here hbl used as a dummy variable, integrating vertically N^2
+     DO jk = 1, jpkm1
+        DO jj = 1, jpj              ! Mixed layer level: w-level
+           DO ji = 1, jpi
+              ikt = mbkt(ji,jj)
+              hbl(ji,jj) = hbl(ji,jj) + MAX( rn2(ji,jj,jk) , 0._wp ) * e3w_n(ji,jj,jk)
+              IF( hbl(ji,jj) < zN2_c )   imld_rst(ji,jj) = MIN( jk , ikt ) + 1   ! Mixed layer level
+           END DO
+        END DO
+     END DO
+     !
+     DO jj = 1, jpj
+        DO ji = 1, jpi
+           iiki = imld_rst(ji,jj)
+           hbl (ji,jj) = gdepw_n(ji,jj,iiki  ) * ssmask(ji,jj)    ! Turbocline depth
+        END DO
+     END DO
+     hbl = MAX(hbl,epsln)
+     hbli(:,:) = hbl(:,:)
+     DEALLOCATE( imld_rst )
+     WRITE(numout,*) ' ===>>>> : hbl computed from stratification'
+   END SUBROUTINE osm_rst
+
+   SUBROUTINE tra_osm( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE tra_osm  ***
+      !!
+      !! ** Purpose :   compute and add to the tracer trend the non-local tracer flux
+      !!
+      !! ** Method  :   ???
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   ztrdt, ztrds   ! 3D workspace
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) :: kt
+      INTEGER :: ji, jj, jk
+      !
+      IF( nn_timing == 1 )  CALL timing_start('tra_osm')
+      !
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'tra_osm : OSM non-local tracer fluxes'
+         IF(lwp) WRITE(numout,*) '~~~~~~~   '
+      ENDIF
+
+      IF( l_trdtra )   THEN                    !* Save ta and sa trends
+         ALLOCATE( ztrdt(jpi,jpj,jpk) )   ;    ztrdt(:,:,:) = tsa(:,:,:,jp_tem)
+         ALLOCATE( ztrds(jpi,jpj,jpk) )   ;    ztrds(:,:,:) = tsa(:,:,:,jp_sal)
+      ENDIF
+
+      ! add non-local temperature and salinity flux
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               tsa(ji,jj,jk,jp_tem) =  tsa(ji,jj,jk,jp_tem)                      &
+                  &                 - (  ghamt(ji,jj,jk  )  &
+                  &                    - ghamt(ji,jj,jk+1) ) /e3t_n(ji,jj,jk)
+               tsa(ji,jj,jk,jp_sal) =  tsa(ji,jj,jk,jp_sal)                      &
+                  &                 - (  ghams(ji,jj,jk  )  &
+                  &                    - ghams(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+
+
+      ! save the non-local tracer flux trends for diagnostic
+      IF( l_trdtra )   THEN
+         ztrdt(:,:,:) = tsa(:,:,:,jp_tem) - ztrdt(:,:,:)
+         ztrds(:,:,:) = tsa(:,:,:,jp_sal) - ztrds(:,:,:)
+!!bug gm jpttdzdf ==> jpttosm
+         CALL trd_tra( kt, 'TRA', jp_tem, jptra_zdf, ztrdt )
+         CALL trd_tra( kt, 'TRA', jp_sal, jptra_zdf, ztrds )
+         DEALLOCATE( ztrdt )      ;     DEALLOCATE( ztrds )
+      ENDIF
+
+      IF(ln_ctl) THEN
+         CALL prt_ctl( tab3d_1=tsa(:,:,:,jp_tem), clinfo1=' osm  - Ta: ', mask1=tmask,   &
+         &             tab3d_2=tsa(:,:,:,jp_sal), clinfo2=       ' Sa: ', mask2=tmask, clinfo3='tra' )
+      ENDIF
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('tra_osm')
+      !
+   END SUBROUTINE tra_osm
+
+   SUBROUTINE trc_osm( kt )          ! Dummy routine
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE trc_osm  ***
+      !!
+      !! ** Purpose :   compute and add to the passive tracer trend the non-local
+      !!                 passive tracer flux
+      !!
+      !!
+      !! ** Method  :   ???
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) :: kt
+      WRITE(*,*) 'trc_osm: Not written yet', kt
+   END SUBROUTINE trc_osm
+
+   SUBROUTINE dyn_osm( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dyn_osm  ***
+      !!
+      !! ** Purpose :   compute and add to the velocity trend the non-local flux
+      !! copied/modified from tra_osm
+      !!
+      !! ** Method  :   ???
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) :: kt
+      INTEGER :: ji, jj, jk
+      !
+      IF( nn_timing == 1 )  CALL timing_start('dyn_osm')
+      !
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'dyn_osm : OSM non-local velocity'
+         IF(lwp) WRITE(numout,*) '~~~~~~~   '
+      ENDIF
+      !code saving tracer trends removed, replace with trdmxl_oce
+
+      ! add non-local u and v fluxes
+      DO jk = 1, jpkm1
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               ua(ji,jj,jk) =  ua(ji,jj,jk)                      &
+                  &                 - (  ghamu(ji,jj,jk  )  &
+                  &                    - ghamu(ji,jj,jk+1) ) / e3u_n(ji,jj,jk)
+               va(ji,jj,jk) =  va(ji,jj,jk)                      &
+                  &                 - (  ghamv(ji,jj,jk  )  &
+                  &                    - ghamv(ji,jj,jk+1) ) / e3v_n(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+
+      ! code for saving tracer trends removed
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('dyn_osm')
+      !
+   END SUBROUTINE dyn_osm
+
+END MODULE zdfosm
+
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfphy.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfphy.F90	(revision 8920)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfphy.F90	(revision 8930)
@@ -18,4 +18,5 @@
    USE zdftke         ! vertical physics: TKE vertical mixing
    USE zdfgls         ! vertical physics: GLS vertical mixing
+   USE zdfosm         ! vertical physics: OSMOSIS vertical mixing
    USE zdfddm         ! vertical physics: double diffusion mixing      
    USE zdfevd         ! vertical physics: convection via enhanced vertical diffusion  
@@ -49,4 +50,5 @@
    INTEGER, PARAMETER ::   np_TKE = 3   ! Turbulente Kinetic Eenergy closure scheme for Kz
    INTEGER, PARAMETER ::   np_GLS = 4   ! Generic Length Scale closure scheme for Kz
+   INTEGER, PARAMETER ::   np_OSM = 5   ! OSMOSIS-OBL closure scheme for Kz
 
    LOGICAL ::   l_zdfsh2   ! shear production term flag (=F for CST, =T otherwise (i.e. TKE, GLS, RIC))
@@ -72,4 +74,5 @@
       !!
       NAMELIST/namzdf/ ln_zdfcst, ln_zdfric, ln_zdftke, ln_zdfgls,   &     ! type of closure scheme
+         &             ln_zdfosm,                                    &     ! type of closure scheme
          &             ln_zdfevd, nn_evdm, rn_evd ,                  &     ! convection : evd
          &             ln_zdfnpc, nn_npc , nn_npcp,                  &     ! convection : npc
@@ -102,4 +105,5 @@
          WRITE(numout,*) '         Turbulent Kinetic Energy closure (TKE)  ln_zdftke = ', ln_zdftke
          WRITE(numout,*) '         Generic Length Scale closure (GLS)      ln_zdfgls = ', ln_zdfgls
+         WRITE(numout,*) '         OSMOSIS-OBL closure (OSM)               ln_zdfosm = ', ln_zdfosm
          WRITE(numout,*) '      convection: '
          WRITE(numout,*) '         enhanced vertical diffusion             ln_zdfevd = ', ln_zdfevd
@@ -156,5 +160,7 @@
       !
       IF( ln_zdfnpc .AND. ln_zdfevd )   CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfnpc and ln_zdfevd' )
+      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(lwp) THEN
          WRITE(numout,*)
@@ -178,4 +184,5 @@
       IF( ln_zdftke ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_TKE   ;   CALL zdf_tke_init   ;   ENDIF
       IF( ln_zdfgls ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_GLS   ;   CALL zdf_gls_init   ;   ENDIF
+      IF( ln_zdfosm ) THEN   ;   ioptio = ioptio + 1   ;    nzdf_phy = np_OSM   ;   CALL zdf_osm_init   ;   ENDIF
       !
       IF( ioptio /= 1 )    CALL ctl_stop( 'zdf_phy_init: one and only one vertical diffusion option has to be defined ' )
@@ -247,4 +254,5 @@
       CASE( np_TKE )   ;   CALL zdf_tke( kt         , zsh2, avm_k, avt_k )    ! TKE closure scheme for Kz
       CASE( np_GLS )   ;   CALL zdf_gls( kt         , zsh2, avm_k, avt_k )    ! GLS closure scheme for Kz
+      CASE( np_OSM )   ;   CALL zdf_osm( kt               , avm_k, avt_k )    ! OSMOSIS closure scheme for Kz
 !     CASE( np_CST )                                  ! Constant Kz (reset avt, avm to the background value)
 !         ! avt_k and avm_k set one for all at initialisation phase
@@ -290,5 +298,5 @@
       CALL lbc_lnk( avm  , 'W', 1. )                  ! needed to compute avm at u- and v-points
       CALL lbc_lnk( avt  , 'W', 1. )                  !!gm  a priori only avm_k and avm are required
-      CALL lbc_lnk( avs  , 'W', 1. )                  !!gm  for calculation, keeped here for output only
+      CALL lbc_lnk( avs  , 'W', 1. )                  !!gm  for calculation, kept here for output only
       !
       IF( l_zdfdrg ) THEN     ! drag  have been updated (non-linear cases)
@@ -303,4 +311,5 @@
          IF( ln_zdfgls )   CALL gls_rst( kt, 'WRITE' )
          IF( ln_zdfric )   CALL ric_rst( kt, 'WRITE' ) 
+         ! NB. OSMOSIS restart (osm_rst) will be called in step.F90 after wn has been updated
       ENDIF
       !
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 8920)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 8930)
@@ -186,4 +186,5 @@
                          CALL dyn_vor       ( kstp )  ! vorticity term including Coriolis
                          CALL dyn_ldf       ( kstp )  ! lateral mixing
+      IF( ln_zdfosm  )   CALL dyn_osm       ( kstp )  ! OSMOSIS non-local velocity fluxes
                          CALL dyn_hpg       ( kstp )  ! horizontal gradient of Hydrostatic pressure
                          CALL dyn_spg       ( kstp )  ! surface pressure gradient
@@ -243,4 +244,7 @@
 #endif
                          CALL tra_adv       ( kstp )  ! horizontal & vertical advection
+      IF( ln_zdfosm  )   CALL tra_osm       ( kstp )  ! OSMOSIS non-local tracer fluxes
+      IF( lrst_oce .AND. ln_zdfosm ) &
+           &             CALL osm_rst( kstp, 'WRITE' )! write OSMOSIS outputs + wn (so must do here) to restarts
                          CALL tra_ldf       ( kstp )  ! lateral mixing
 
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 8920)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 8930)
@@ -65,4 +65,5 @@
 
    USE zdfphy          ! vertical physics manager      (zdf_phy_init routine)
+   USE zdfosm  , ONLY : osm_rst, dyn_osm, tra_osm      ! OSMOSIS routines used in step.F90
 
    USE step_diu        ! Time stepping for diurnal sst
