Index: /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90
===================================================================
--- /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90	(revision 14774)
+++ /NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90	(revision 14775)
@@ -39,13 +39,17 @@
    !!   'ln_zdfosm'                                             OSMOSIS scheme
    !!----------------------------------------------------------------------
-   !!   zdf_osm                  : update momentum and tracer Kz from osm scheme
-   !!   zdf_osm_vertical_average : compute vertical averages over boundary layers
-   !!   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
-   !!
-   !! Subroutines in revised code.
+   !!   zdf_osm        : update momentum and tracer Kz from osm scheme
+   !!      zdf_osm_vertical_average        : compute vertical averages over boundary layers
+   !!      zdf_osm_velocity_rotation       : rotate velocity components
+   !!         zdf_osm_velocity_rotation_2d :    rotation of 2d fields
+   !!         zdf_osm_velocity_rotation_3d :    rotation of 3d fields
+   !!      zdf_osm_osbl_state              : determine the state of the OSBL
+   !!      zdf_osm_diffusivity_viscosity   : compute eddy diffusivity and viscosity profiles
+   !!      zdf_osm_fgr_terms               : compute flux-gradient relationship terms
+   !!   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            ! ocean dynamics and active tracers
@@ -83,23 +87,4 @@
    PUBLIC   ln_osm_mle    ! logical needed by tra_mle_init in tramle.F90
 
-   ! Scales
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: swth0     ! Surface heat flux (Kinematic)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: sws0      ! Surface freshwater flux
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: swb0      ! Surface buoyancy flux
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: suw0      ! Surface u-momentum flux
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: sustar    ! Friction velocity
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: scos_wind ! Cos angle of surface stress
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: ssin_wind ! Sin angle of surface stress
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: swthav    ! Heat flux - bl average
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: swsav     ! Freshwater flux - bl average
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: swbav     ! Buoyancy flux - bl average
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: sustke    ! Surface Stokes drift
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: dstokes   ! Penetration depth of the Stokes drift
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: swstrl    ! Langmuir velocity scale
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: swstrc    ! Convective velocity scale
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: sla       ! Trubulent Langmuir number
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: svstr     ! Velocity scale that tends to sustar for large Langmuir number
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: shol      ! Stability parameter for boundary layer
-
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamu    !: non-local u-momentum flux
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ghamv    !: non-local v-momentum flux
@@ -111,9 +96,37 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hml      ! ML depth
 
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)           ::   r1_ft    ! inverse of the modified Coriolis parameter at t-pts
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hmle     ! Depth of layer affexted by mixed layer eddies in Fox-Kemper parametrization
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dbdx_mle ! zonal buoyancy gradient in ML
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dbdy_mle ! meridional buoyancy gradient in ML
    INTEGER,  PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   mld_prof ! level of base of MLE layer.
+
+   INTERFACE zdf_osm_velocity_rotation
+      !!---------------------------------------------------------------------
+      !!              ***  INTERFACE zdf_velocity_rotation  ***
+      !!---------------------------------------------------------------------
+      MODULE PROCEDURE zdf_osm_velocity_rotation_2d
+      MODULE PROCEDURE zdf_osm_velocity_rotation_3d
+   END INTERFACE
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: r1_ft       ! inverse of the modified Coriolis parameter at t-pts
+
+   ! Scales
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: swth0       ! Surface heat flux (Kinematic)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: sws0        ! Surface freshwater flux
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: swb0        ! Surface buoyancy flux
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: suw0        ! Surface u-momentum flux
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: sustar      ! Friction velocity
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: scos_wind   ! Cos angle of surface stress
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: ssin_wind   ! Sin angle of surface stress
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: swthav      ! Heat flux - bl average
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: swsav       ! Freshwater flux - bl average
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: swbav       ! Buoyancy flux - bl average
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: sustke      ! Surface Stokes drift
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: dstokes     ! Penetration depth of the Stokes drift
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: swstrl      ! Langmuir velocity scale
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: swstrc      ! Convective velocity scale
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: sla         ! Trubulent Langmuir number
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: svstr       ! Velocity scale that tends to sustar for large Langmuir number
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   :: shol        ! Stability parameter for boundary layer
 
    !                      !!** Namelist  namzdf_osm  **
@@ -174,12 +187,4 @@
    INTEGER :: idebug = 236
    INTEGER :: jdebug = 228
-
-   INTERFACE zdf_osm_velocity_rotation
-      !!---------------------------------------------------------------------
-      !!              ***  INTERFACE zdf_velocity_rotation  ***
-      !!---------------------------------------------------------------------
-      MODULE PROCEDURE zdf_osm_velocity_rotation_2d
-      MODULE PROCEDURE zdf_osm_velocity_rotation_3d
-   END INTERFACE
 
    !! * Substitutions
@@ -287,5 +292,4 @@
       REAL(wp), DIMENSION(jpi,jpj) :: zwb_min
 
-
       REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk_b  ! MLE buoyancy flux averaged over OSBL
       REAL(wp), DIMENSION(jpi,jpj) :: zwb_fk    ! max MLE buoyancy flux
@@ -366,8 +370,8 @@
       swstrl(:,:) = zlarge ; svstr(:,:)      = zlarge ; swstrc(:,:)     = zlarge ; suw0(:,:)      = zlarge
       swth0(:,:)  = zlarge ; sws0(:,:)       = zlarge ; swb0(:,:)       = zlarge
-      swthav(:,:) = zlarge ; swsav(:,:)      = zlarge ; swbav(:,:)      = zlarge ; zwb_ent(:,:)   = zlarge
+      swthav(:,:) = zlarge ; swsav(:,:)      = zlarge ; swbav(:,:)      = zlarge
       sustke(:,:) = zlarge ; sla(:,:)        = zlarge ; scos_wind(:,:)  = zlarge ; ssin_wind(:,:) = zlarge
       shol(:,:)   = zlarge ; zalpha_pyc(:,:) = zlarge
-      lconv(:,:)  = .FALSE.; lpyc(:,:) = .FALSE. ; lflux(:,:) = .FALSE. ;  lmle(:,:) = .FALSE.
+      lpyc(:,:) = .FALSE. ; lflux(:,:) = .FALSE. ;  lmle(:,:) = .FALSE.
       ! mixed layer
       ! no initialization of zhbl or zhml (or zdh?)
@@ -705,5 +709,6 @@
 
       ! Determine the state of the OSBL, stable/unstable, shear/no shear
-      CALL zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear )
+      CALL zdf_osm_osbl_state( Kmm, lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear,   &
+         &                     zhbl, zhml, zdh, zdb_bl, zdu_bl, zdv_bl, zdb_ml, zdu_ml, zdv_ml )
 
 #ifdef key_osm_debug
@@ -965,5 +970,7 @@
       ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship
       !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      CALL zdf_osm_diffusivity_viscosity( zdiffut, zviscos )
+      CALL zdf_osm_diffusivity_viscosity( Kbb, Kmm, ibld, imld, lconv, lshear, lpyc, lcoup, j_ddh, zdiffut, zviscos,           &
+         &                                zhbl, zhml, zdh, zdhdt, zshear, zb_bl, zu_bl, zv_bl, zb_ml, zu_ml, zv_ml, zwb_ent,   &
+         &                                zwb_min )
 #ifdef key_osm_debug
       IF(narea==nn_narea_db) THEN
@@ -1212,464 +1219,4 @@
 
    CONTAINS
-      ! subroutine code changed, needs syntax checking.
-      SUBROUTINE zdf_osm_diffusivity_viscosity( zdiffut, zviscos )
-         
-         !!---------------------------------------------------------------------
-         !!                   ***  ROUTINE zdf_osm_diffusivity_viscosity  ***
-         !!
-         !! ** Purpose : Determines the eddy diffusivity and eddy viscosity profiles in the mixed layer and the pycnocline.
-         !!
-         !! ** Method  :
-         !!
-         !! !!----------------------------------------------------------------------
-         REAL(wp), DIMENSION(:,:,:) :: zdiffut
-         REAL(wp), DIMENSION(:,:,:) :: zviscos
-         ! local
-
-         ! Scales used to calculate eddy diffusivity and viscosity profiles
-         REAL(wp), DIMENSION(jpi,jpj) :: zdifml_sc, zvisml_sc
-         REAL(wp), DIMENSION(jpi,jpj) :: zdifpyc_n_sc, zdifpyc_s_sc, zdifpyc_shr
-         REAL(wp), DIMENSION(jpi,jpj) :: zvispyc_n_sc, zvispyc_s_sc,zvispyc_shr
-         REAL(wp), DIMENSION(jpi,jpj) :: zbeta_d_sc, zbeta_v_sc
-         REAL(wp), DIMENSION(jpi,jpj) ::   zb_coup, zc_coup_vis, zc_coup_dif
-         !
-         REAL(wp) :: zvel_sc_pyc, zvel_sc_ml, zstab_fac, zz_b
-         REAL(wp) ::   za_cubic, zb_cubic, zc_cubic, zd_cubic   ! Coefficients in cubic polynomial specifying diffusivity in pycnocline
-         REAL(wp) ::   zznd_ml, zznd_pyc
-         REAL(wp) ::   zmsku, zmskv
-
-         REAL(wp), PARAMETER :: rn_dif_ml = 0.8, rn_vis_ml = 0.375
-         REAL(wp), PARAMETER :: rn_dif_pyc = 0.15, rn_vis_pyc = 0.142
-         REAL(wp), PARAMETER :: rn_vispyc_shr = 0.15
-
-         IF( ln_timing ) CALL timing_start('zdf_osm_dv')
-
-         zb_coup(:,:) = 0.0_wp
-
-         DO_2D( 0, 0, 0, 0 )
-            IF ( lconv(ji,jj) ) THEN
-
-               zvel_sc_pyc = ( 0.15 * svstr(ji,jj)**3 + swstrc(ji,jj)**3 + 4.25 * zshear(ji,jj) * zhbl(ji,jj) )**pthird
-               zvel_sc_ml = ( svstr(ji,jj)**3 + 0.5 * swstrc(ji,jj)**3 )**pthird
-               zstab_fac = ( zhml(ji,jj) / zvel_sc_ml * ( 1.4 - 0.4 / ( 1.0 + EXP(-3.5 * LOG10(-shol(ji,jj) ) ) )**1.25 ) )**2
-
-               zdifml_sc(ji,jj) = rn_dif_ml * zhml(ji,jj) * zvel_sc_ml
-               zvisml_sc(ji,jj) = rn_vis_ml * zdifml_sc(ji,jj)
-#ifdef key_osm_debug
-               IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
-                  WRITE(narea+100,'(2(a,g11.3))')'Start of 1st major loop of osm_diffusivity_viscositys, lconv=T: zdifml_sc=',zdifml_sc(ji,jj),' zvisml_sc=',zvisml_sc(ji,jj)
-                  WRITE(narea+100,'(3(a,g11.3))')'zvel_sc_pyc=',zvel_sc_pyc,' zvel_sc_ml=',zvel_sc_ml,' zstab_fac=',zstab_fac
-                  FLUSH(narea+100)
-               END IF
-#endif
-
-               IF ( lpyc(ji,jj) ) THEN
-                  zdifpyc_n_sc(ji,jj) =  rn_dif_pyc * zvel_sc_ml * zdh(ji,jj)
-                  zvispyc_n_sc(ji,jj) = 0.09 * zvel_sc_pyc * ( 1.0 - zhbl(ji,jj) / zdh(ji,jj) )**2 * ( 0.005 * ( zu_ml(ji,jj)-zu_bl(ji,jj) )**2 + 0.0075 * ( zv_ml(ji,jj)-zv_bl(ji,jj) )**2 ) / zdh(ji,jj)
-                  zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * zdh(ji,jj) + zvispyc_n_sc(ji,jj) * zstab_fac
-#ifdef key_osm_debug
-                  IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
-                     WRITE(narea+100,'(2(a,g11.3))')' lpyc=lconv=T, variables w/o shear contributions: zdifpyc_n_sc',zdifpyc_n_sc(ji,jj) ,' zvispyc_n_sc=',zvispyc_n_sc(ji,jj)
-                     FLUSH(narea+100)
-                  END IF
-#endif
-
-                  IF ( lshear(ji,jj) .AND. j_ddh(ji,jj) /= 2 ) THEN
-                     zdifpyc_n_sc(ji,jj) = zdifpyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj) )**pthird * zhbl(ji,jj)
-                     zvispyc_n_sc(ji,jj) = zvispyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj ) )**pthird * zhbl(ji,jj)
-                  ENDIF
-#ifdef key_osm_debug
-                  IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
-                     WRITE(narea+100,'(2(a,g11.3))')' lpyc=lconv=T, variables w shear contributions: zdifpyc_n_sc',zdifpyc_n_sc(ji,jj) ,' zvispyc_n_sc=',zvispyc_n_sc(ji,jj)
-                     FLUSH(narea+100)
-                  END IF
-#endif
-
-                  zdifpyc_s_sc(ji,jj) = zwb_ent(ji,jj) + 0.0025 * zvel_sc_pyc * ( zhbl(ji,jj) / zdh(ji,jj) - 1.0 ) * ( zb_ml(ji,jj) - zb_bl(ji,jj) )
-                  zvispyc_s_sc(ji,jj) = 0.09 * ( zwb_min(ji,jj) + 0.0025 * zvel_sc_pyc * ( zhbl(ji,jj) / zdh(ji,jj) - 1.0 ) * ( zb_ml(ji,jj) - zb_bl(ji,jj) ) )
-#ifdef key_osm_debug
-                  IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
-                     WRITE(narea+100,'(2(a,g11.3))')' 1st shot at: zdifpyc_s_sc',zdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',zvispyc_s_sc(ji,jj)
-                     FLUSH(narea+100)
-                  END IF
-#endif
-                  zdifpyc_s_sc(ji,jj) = 0.09 * zdifpyc_s_sc(ji,jj) * zstab_fac
-                  zvispyc_s_sc(ji,jj) = zvispyc_s_sc(ji,jj) * zstab_fac
-#ifdef key_osm_debug
-                  IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
-                     WRITE(narea+100,'(2(a,g11.3))')' 2nd shot at: zdifpyc_s_sc',zdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',zvispyc_s_sc(ji,jj)
-                     FLUSH(narea+100)
-                  END IF
-#endif
-
-                  zdifpyc_s_sc(ji,jj) = MAX( zdifpyc_s_sc(ji,jj), -0.5 * zdifpyc_n_sc(ji,jj) )
-                  zvispyc_s_sc(ji,jj) = MAX( zvispyc_s_sc(ji,jj), -0.5_wp * zvispyc_n_sc(ji,jj) )
-
-                  zbeta_d_sc(ji,jj) = 1.0 - ( ( zdifpyc_n_sc(ji,jj) + 1.4 * zdifpyc_s_sc(ji,jj) ) / ( zdifml_sc(ji,jj) + epsln ) )**p2third
-                  zbeta_v_sc(ji,jj) = 1.0 -  2.0 * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) / ( zvisml_sc(ji,jj) + epsln )
-               ELSE
-                  zdifpyc_n_sc(ji,jj) = rn_dif_pyc * zvel_sc_ml * zdh(ji,jj)   ! ag 19/03
-                  zdifpyc_s_sc(ji,jj) = 0.0_wp   ! ag 19/03
-                  zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * zdh(ji,jj)   ! ag 19/03
-                  zvispyc_s_sc(ji,jj) = 0.0_wp   ! ag 19/03
-                  IF(lcoup(ji,jj) ) THEN   ! ag 19/03
-                     ! code from SUBROUTINE tke_tke zdftke.F90; uses bottom drag velocity rCdU_bot(ji,jj) = -Cd|ub|
-                     !     already calculated at T-points in SUBROUTINE zdf_drg from zdfdrg.F90
-                     !  Gives friction velocity sqrt bottom drag/rho_0 i.e. u* = SQRT(rCdU_bot*ub)
-                     ! wet-cell averaging ..
-                     zmsku = 0.5_wp * ( 2.0_wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
-                     zmskv = 0.5_wp * ( 2.0_wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )
-                     zb_coup(ji,jj) = 0.4_wp * SQRT(-1.0_wp * rCdU_bot(ji,jj) *   &
-                        &             SQRT(  ( zmsku*( uu(ji,jj,mbkt(ji,jj),Kbb)+uu(ji-1,jj,mbkt(ji,jj),Kbb) ) )**2   &
-                        &                  + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2  ) )
-
-                     zz_b = -gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)   ! ag 19/03
-                     zc_coup_vis(ji,jj) = -0.5_wp * ( 0.5_wp * zvisml_sc(ji,jj) / zhml(ji,jj) - zb_coup(ji,jj) ) /   &
-                        &                 ( zhml(ji,jj) + zz_b )   ! ag 19/03
-#ifdef key_osm_debug
-                     IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
-                        WRITE(narea+100,'(4(a,g11.3))')' lcoup = T; 1st pz_b= ', zz_b, ' pb_coup ', zb_coup(ji,jj),   &
-                           &                           ' pc_coup_vis ', zc_coup_vis(ji,jj), ' rCdU_bot ',rCdU_bot(ji,jj)
-                        WRITE(narea+100,'(2(a,g11.3))')' zmsku ', zmsku, ' zmskv ', zmskv
-                        FLUSH(narea+100)
-                     END IF
-#endif
-!#ifdef key_osm_debug
-!                     WRITE(narea+400,'(4(a,i7))') ' lcoup = T at ji=',ji,' jj= ',jj,' jig= ', mig(ji), ' jjg= ', mjg(jj)
-!                     WRITE(narea+400,'(3(a,g11.3))') '1st pz_b= ', zz_b, 'pb_coup', zb_coup(ji,jj),   &
-!                        &                            ' pc_coup_vis', zc_coup_vis(ji,jj)
-!                     FLUSH(narea+400)
-!#endif
-                     zz_b = -zhml(ji,jj) + gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)   ! ag 19/03
-                     zbeta_v_sc(ji,jj) = 1.0_wp - 2.0_wp * ( zb_coup(ji,jj) * zz_b + zc_coup_vis(ji,jj) * zz_b**2 ) /   &
-                        &                                  zvisml_sc(ji,jj)   ! ag 19/03
-                     zbeta_d_sc(ji,jj) = 1.0_wp - ( ( zb_coup(ji,jj) * zz_b + zc_coup_vis(ji,jj) * zz_b**2 ) /   &
-                        &                           zdifml_sc(ji,jj) )**p2third
-                     zc_coup_dif(ji,jj) = 0.5_wp * ( -zdifml_sc(ji,jj) / zhml(ji,jj) * ( 1.0_wp - zbeta_d_sc(ji,jj) )**1.5_wp +   &
-                        &                 1.5_wp * ( zdifml_sc(ji,jj) / zhml(ji,jj) ) * zbeta_d_sc(ji,jj) *   &
-                        &                          SQRT( 1.0_wp - zbeta_d_sc(ji,jj) ) - zb_coup(ji,jj) ) / zz_b   ! ag 19/03
-#ifdef key_osm_debug
-                     IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
-                        WRITE(narea+100,'(2(a,g11.3))')' 2nd pz_b= ', zz_b, ' pc_coup_dif', zc_coup_dif(ji,jj)
-                        FLUSH(narea+100)
-                     END IF
-#endif
-!#ifdef key_osm_debug
-!                     WRITE(narea+400,'(3(a,g11.3))') '2nd pz_b= ', pz_b,' pc_coup_dif', zc_coup_dif(ji,jj)
-!                     FLUSH(narea+400)
-!#endif
-                  ELSE   ! ag 19/03
-                    zbeta_d_sc(ji,jj) = 1.0_wp - ( ( zdifpyc_n_sc(ji,jj) + 1.4_wp * zdifpyc_s_sc(ji,jj) ) /   &
-                       &                           ( zdifml_sc(ji,jj) + epsln ) )**p2third   ! ag 19/03
-                    zbeta_v_sc(ji,jj) = 1.0_wp - 2.0_wp * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) /   &
-                       &                         ( zvisml_sc(ji,jj) + epsln )   ! ag 19/03
-                  ENDIF   ! ag 19/03
-               ENDIF      ! ag 19/03
-#ifdef key_osm_debug
-               IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
-                  WRITE(narea+100,'(2(a,g11.3))')'lconv=T: zbeta_d_sc',zbeta_d_sc(ji,jj) ,' zbeta_v_sc=',zbeta_v_sc(ji,jj)
-                  WRITE(narea+100,'(2(a,g11.3))')' Final zdifpyc_n_sc',zdifpyc_n_sc(ji,jj) ,' zvispyc_n_sc=',zvispyc_n_sc(ji,jj)
-                  WRITE(narea+100,'(2(a,g11.3))')' Final zdifpyc_s_sc',zdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',zvispyc_s_sc(ji,jj)
-                  FLUSH(narea+100)
-               END IF
-#endif
-            ELSE
-               zdifml_sc(ji,jj) = svstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( shol(ji,jj) / 0.6_wp )**2 ), 0.2_wp)
-               zvisml_sc(ji,jj) = svstr(ji,jj) * zhbl(ji,jj) * MAX( EXP ( -( shol(ji,jj) / 0.6_wp )**2 ), 0.2_wp)
-#ifdef key_osm_debug
-               IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
-                  WRITE(narea+100,'(a,g11.3)')'End of 1st major loop of osm_diffusivity_viscositys, lconv=F: zdifml_sc=',zdifml_sc(ji,jj),' zvisml_sc=',zvisml_sc(ji,jj)
-                  FLUSH(narea+100)
-               END IF
-#endif
-            END IF
-         END_2D
-         !
-         DO_2D( 0, 0, 0, 0 )
-            IF ( lconv(ji,jj) ) THEN
-               DO jk = 2, imld(ji,jj)   ! mixed layer diffusivity
-                  zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj)
-                  !
-                  zdiffut(ji,jj,jk) = zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_d_sc(ji,jj) * zznd_ml )**1.5
-                  !
-                  zviscos(ji,jj,jk) = zvisml_sc(ji,jj) * zznd_ml * ( 1.0 - zbeta_v_sc(ji,jj) * zznd_ml ) &
-                     &            *                                      ( 1.0 - 0.5 * zznd_ml**2 )
-               END DO
-
-               ! Coupling to bottom
-            
-               IF ( lcoup(ji,jj) ) THEN                                                          ! ag 19/03
-                  DO jk = mbkt(ji,jj), imld(ji,jj), -1                                           ! ag 19/03
-                     zz_b = - ( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) )           ! ag 19/03
-                     zviscos(ji,jj,jk) = zb_coup(ji,jj) * zz_b + zc_coup_vis(ji,jj) * zz_b**2    ! ag 19/03
-                     zdiffut(ji,jj,jk) = zb_coup(ji,jj) * zz_b + zc_coup_dif(ji,jj) * zz_b**2    ! ag 19/03
-                  END DO                                                                         ! ag 19/03
-               ENDIF                                                                             ! ag 19/03
-               ! pycnocline
-               IF ( lpyc(ji,jj) ) THEN 
-                  ! Diffusivity profile in the pycnocline given by cubic polynomial. Note, if lpyc TRUE can't be coupled to seabed.
-                  za_cubic = 0.5
-                  zb_cubic = -1.75 * zdifpyc_s_sc(ji,jj) / zdifpyc_n_sc(ji,jj)
-                  zd_cubic = ( zdh(ji,jj) * zdifml_sc(ji,jj) / zhml(ji,jj) * SQRT( 1.0 - zbeta_d_sc(ji,jj) ) * ( 2.5 * zbeta_d_sc(ji,jj) - 1.0 ) &
-                     & - 0.85 * zdifpyc_s_sc(ji,jj) ) / MAX(zdifpyc_n_sc(ji,jj), 1.e-8)
-                  zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic  - zb_cubic )
-                  zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic
-                  DO jk = imld(ji,jj) , ibld(ji,jj)
-                     zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6)
-                     !
-                     zdiffut(ji,jj,jk) = zdifpyc_n_sc(ji,jj) * ( za_cubic + zb_cubic * zznd_pyc + zc_cubic * zznd_pyc**2 +   zd_cubic * zznd_pyc**3 )
-
-                     zdiffut(ji,jj,jk) = zdiffut(ji,jj,jk) + zdifpyc_s_sc(ji,jj) * ( 1.75 * zznd_pyc - 0.15 * zznd_pyc**2 - 0.2 * zznd_pyc**3 )
-                  END DO
-                  ! viscosity profiles.
-                  za_cubic = 0.5
-                  zb_cubic = -1.75 * zvispyc_s_sc(ji,jj) / zvispyc_n_sc(ji,jj)
-                  zd_cubic = ( 0.5 * zvisml_sc(ji,jj) * zdh(ji,jj) / zhml(ji,jj) - 0.85 * zvispyc_s_sc(ji,jj)  )  / MAX(zvispyc_n_sc(ji,jj), 1.e-8)
-                  zd_cubic = zd_cubic - zb_cubic - 2.0 * ( 1.0 - za_cubic - zb_cubic )
-                  zc_cubic = 1.0 - za_cubic - zb_cubic - zd_cubic
-                  DO jk = imld(ji,jj) , ibld(ji,jj)
-                     zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / MAX(zdh(ji,jj), 1.e-6)
-                     zviscos(ji,jj,jk) = zvispyc_n_sc(ji,jj) * ( za_cubic + zb_cubic * zznd_pyc + zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 )
-                     zviscos(ji,jj,jk) = zviscos(ji,jj,jk) + zvispyc_s_sc(ji,jj) * ( 1.75 * zznd_pyc - 0.15 * zznd_pyc**2 -0.2 * zznd_pyc**3 )
-                  END DO
-!                  IF ( zdhdt(ji,jj) > 0._wp ) THEN
-!                     zdiffut(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w(ji,jj,ibld(ji,jj)+1,Kmm), 1.0e-6 )
-!                     zviscos(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * zdhdt(ji,jj) * e3w(ji,jj,ibld(ji,jj)+1,Kmm), 1.0e-6 )
-!                  ELSE
-!                     zdiffut(ji,jj,ibld(ji,jj)) = 0._wp
-!                     zviscos(ji,jj,ibld(ji,jj)) = 0._wp
-!                  ENDIF
-               ENDIF
-            ELSE
-               ! stable conditions
-               DO jk = 2, ibld(ji,jj)
-                  zznd_ml = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
-                  zdiffut(ji,jj,jk) = 0.75 * zdifml_sc(ji,jj) * zznd_ml * ( 1.0 - zznd_ml )**1.5
-                  zviscos(ji,jj,jk) = 0.375 * zvisml_sc(ji,jj) * zznd_ml * (1.0 - zznd_ml) * ( 1.0 - zznd_ml**2 )
-               END DO
-
-               IF ( zdhdt(ji,jj) > 0._wp ) THEN
-                  zdiffut(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w(ji, jj, ibld(ji,jj), Kmm)
-                  zviscos(ji,jj,ibld(ji,jj)) = MAX(zdhdt(ji,jj), 1.0e-6) * e3w(ji, jj, ibld(ji,jj), Kmm)
-               ENDIF
-            ENDIF   ! end if ( lconv )
-            !
-         END_2D
-         IF( iom_use("pb_coup") ) CALL iom_put( "pb_coup", tmask(:,:,1) * zb_coup(:,:) )   ! BBL-coupling velocity scale
-         IF( ln_timing ) CALL timing_stop('zdf_osm_dv')
-
-      END SUBROUTINE zdf_osm_diffusivity_viscosity
-
-      SUBROUTINE zdf_osm_osbl_state( lconv, lshear, j_ddh, zwb_ent, zwb_min, zshear )
-
-         !!---------------------------------------------------------------------
-         !!                   ***  ROUTINE zdf_osm_osbl_state  ***
-         !!
-         !! ** Purpose : Determines the state of the OSBL, stable/unstable, shear/ noshear. Also determines shear production, entrainment buoyancy flux and interfacial Richardson number
-         !!
-         !! ** Method  :
-         !!
-         !! !!----------------------------------------------------------------------
-
-         INTEGER, DIMENSION(jpi,jpj) :: j_ddh  ! j_ddh = 0, active shear layer; j_ddh=1, shear layer not active; j_ddh=2 shear production low.
-
-         LOGICAL, DIMENSION(jpi,jpj) :: lconv, lshear
-
-         REAL(wp), DIMENSION(jpi,jpj) :: zwb_ent, zwb_min ! Buoyancy fluxes at base of well-mixed layer.
-         REAL(wp), DIMENSION(jpi,jpj) :: zshear  ! production of TKE due to shear across the pycnocline
-
-         ! Local Variables
-
-         INTEGER :: jj, ji
-
-         REAL(wp), DIMENSION(jpi,jpj) :: zekman
-         REAL(wp), DIMENSION(jpi,jpj) :: zri_p, zri_b   ! Richardson numbers
-         REAL(wp) :: zshear_u, zshear_v, zwb_shr
-         REAL(wp) :: zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes
-
-         REAL, PARAMETER :: za_shr = 0.4, zb_shr = 6.5, za_wb_s = 0.8
-         REAL, PARAMETER :: zalpha_c = 0.2, zalpha_lc = 0.03
-         REAL, PARAMETER :: zalpha_ls = 0.06, zalpha_s = 0.15
-         REAL, PARAMETER :: rn_ri_p_thresh = 27.0
-         REAL, PARAMETER :: zri_c = 0.25
-         REAL, PARAMETER :: zek = 4.0
-         REAL, PARAMETER :: zrot=0._wp  ! dummy rotation rate of surface stress.
-
-         IF( ln_timing ) CALL timing_start('zdf_osm_os')
-         ! Determins stability and set flag lconv
-         DO_2D( 0, 0, 0, 0 )
-            IF ( shol(ji,jj) < 0._wp ) THEN
-               lconv(ji,jj) = .TRUE.
-            ELSE
-               lconv(ji,jj) = .FALSE.
-            ENDIF
-         END_2D
-
-         zekman(A2D(0)) = EXP( -1.0_wp * zek * ABS( ff_t(A2D(0)) ) * zhbl(A2D(0)) / MAX(sustar(A2D(0)), 1.e-8 ) )
-
-         zshear(:,:) = zlarge
-         zshear(A2D(0)) = 0._wp
-#ifdef key_osm_debug
-         IF(narea==nn_narea_db) THEN
-            ji=iloc_db; jj=jloc_db
-            WRITE(narea+100,'(a,g11.3)') &
-               & 'zdf_osm_osbl_state start: zekman=', zekman(ji,jj)
-            FLUSH(narea+100)
-         END IF
-#endif
-         j_ddh(A2D(0)) = 1
-
-         DO_2D( 0, 0, 0, 0 )
-            IF ( lconv(ji,jj) ) THEN
-               IF ( zdb_bl(ji,jj) > 0._wp ) THEN
-                  zri_p(ji,jj) = MAX (  SQRT( zdb_bl(ji,jj) * zdh(ji,jj) / MAX( zdu_bl(ji,jj)**2 + zdv_bl(ji,jj)**2, 1.e-8) )  *  ( zhbl(ji,jj) / zdh(ji,jj) ) * ( svstr(ji,jj) / MAX( sustar(ji,jj), 1.e-6 ) )**2 &
-                     & / MAX( zekman(ji,jj), 1.e-6 )  , 5._wp )
-
-                  IF ( ff_t(ji,jj) >= 0.0_wp ) THEN
-                     ! Northern hemisphere
-                     zri_b(ji,jj) = zdb_ml(ji,jj) * zdh(ji,jj) / ( MAX( zdu_ml(ji,jj), 1e-5_wp )**2 + MAX( -1.0_wp * zdv_ml(ji,jj), 1e-5_wp)**2 )
-                  ELSE
-                     ! Southern hemisphere
-                     zri_b(ji,jj) = zdb_ml(ji,jj) * zdh(ji,jj) / ( MAX( zdu_ml(ji,jj), 1e-5_wp )**2 + MAX(           zdv_ml(ji,jj), 1e-5_wp)**2 )
-                  END IF
-                  zshear(ji,jj) = za_shr * zekman(ji,jj) * ( MAX( sustar(ji,jj)**2 * zdu_ml(ji,jj) / zhbl(ji,jj), 0._wp ) + zb_shr * MAX( -ff_t(ji,jj) * sustke(ji,jj) * dstokes(ji,jj) * zdv_ml(ji,jj) / zhbl(ji,jj), 0._wp ) )
-#ifdef key_osm_debug
-                  ! IF(narea==nn_narea_db)THEN
-                  !    WRITE(narea+100,'(2(a,i10.4))')'ji',ji,'jj',jj
-                  !    WRITE(narea+100,'(2(a,i10.4))')'iloc_db',iloc_db,'jloc_db',jloc_db
-                  !    WRITE(narea+100,'(2(a,i10.4))')'iloc_db+',mi0(nn_idb),'jloc_db+',mj0(nn_jdb)
-                  !    FLUSH(narea+100)
-                  ! END IF
-                  IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
-                     WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state 1st zshear: zshear=',zshear(ji,jj)
-                     WRITE(narea+100,'(2(a,g11.3))')'zdf_osm_osbl_state 1st zshear: zri_b=',zri_b(ji,jj),' zri_p=',zri_p(ji,jj)
-                     FLUSH(narea+100)
-                  END IF
-#endif
-                  ! Stability dependence
-                  zshear(ji,jj) = zshear(ji,jj) * EXP( -0.75_wp * MAX( 0.0_wp, ( zri_b(ji,jj) - zri_c ) / zri_c ) )
-#ifdef key_osm_debug
-                  IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
-                     WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state 1st zshear: zshear inc ri part=',zshear(ji,jj)
-                     FLUSH(narea+100)
-                  END IF
-#endif
-
-                  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-                  ! Test ensures j_ddh=0 is not selected. Change to zri_p<27 when  !
-                  ! full code available                                          !
-                  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-                  IF ( zshear(ji,jj) > 1e-10 ) THEN
-                     IF ( zri_p(ji,jj) < rn_ri_p_thresh .AND. MIN( hu(ji,jj,Kmm), hu(ji-1,jj,Kmm), hv(ji,jj,Kmm), hv(ji,jj-1,Kmm) ) > 100.0_wp ) THEN
-                        ! Growing shear layer
-                        j_ddh(ji,jj) = 0
-                        lshear(ji,jj) = .TRUE.
-                     ELSE
-                        j_ddh(ji,jj) = 1
-                        !                 IF ( zri_b <= 1.5 .and. zshear(ji,jj) > 0._wp ) THEN
-                        ! shear production large enough to determine layer charcteristics, but can't maintain a shear layer.
-                        lshear(ji,jj) = .TRUE.
-                        !             ELSE
-                     END IF
-                  ELSE
-                     j_ddh(ji,jj) = 2
-                     lshear(ji,jj) = .FALSE.
-                  END IF
-                  ! Shear production may not be zero, but is small and doesn't determine characteristics of pycnocline.
-                  !               zshear(ji,jj) = 0.5 * zshear(ji,jj)
-                  !               lshear(ji,jj) = .FALSE.
-                  !             ENDIF
-               ELSE                ! zdb_bl test, note zshear set to zero
-                  j_ddh(ji,jj) = 2
-                  lshear(ji,jj) = .FALSE.
-               ENDIF
-            ENDIF
-         END_2D
-
-         ! Calculate entrainment buoyancy flux due to surface fluxes.
-
-         DO_2D( 0, 0, 0, 0 )
-            IF ( lconv(ji,jj) ) THEN
-               zwcor = ABS(ff_t(ji,jj)) * zhbl(ji,jj) + epsln
-               zrf_conv = TANH( ( swstrc(ji,jj) / zwcor )**0.69 )
-               zrf_shear = TANH( ( sustar(ji,jj) / zwcor )**0.69 )
-               zrf_langmuir = TANH( ( swstrl(ji,jj) / zwcor )**0.69 )
-               IF (nn_osm_SD_reduce > 0 ) THEN
-                  ! Effective Stokes drift already reduced from surface value
-                  zr_stokes = 1.0_wp
-               ELSE
-                  ! Effective Stokes drift only reduced by factor rn_zdfodm_adjust_sd,
-                  ! requires further reduction where BL is deep
-                  zr_stokes = 1.0 - EXP( -25.0 * dstokes(ji,jj) / hbl(ji,jj) &
-                     &                  * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) )
-               END IF
-               zwb_ent(ji,jj) = -2.0_wp * zalpha_c * zrf_conv * swbav(ji,jj) &
-                  &                  - zalpha_s * zrf_shear * sustar(ji,jj)**3 / zhml(ji,jj) &
-                  &         + zr_stokes * ( zalpha_s * EXP( -1.5_wp * sla(ji,jj) ) * zrf_shear * sustar(ji,jj)**3 &
-                  &                                         - zrf_langmuir * zalpha_lc * swstrl(ji,jj)**3 ) / zhml(ji,jj)
-               !
-#ifdef key_osm_debug
-               IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
-                  WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state conv+shear0/lang: zwb_ent=',zwb_ent(ji,jj)
-                  FLUSH(narea+100)
-               END IF
-#endif
-
-            ENDIF
-         END_2D
-
-         zwb_min(:,:) = zlarge
-
-         DO_2D( 0, 0, 0, 0 )
-            IF ( lshear(ji,jj) ) THEN
-               IF ( lconv(ji,jj) ) THEN
-                  ! Unstable OSBL
-                  zwb_shr = -1.0_wp * za_wb_s * zri_b(ji,jj) * zshear(ji,jj)
-#ifdef key_osm_debug
-                  IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
-                     WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state 1st zwb_shr: zwb_shr=',zwb_shr
-                     FLUSH(narea+100)
-                  END IF
-#endif
-                  IF ( j_ddh(ji,jj) == 0 ) THEN
-
-                     ! ! Developing shear layer, additional shear production possible.
-
-                     !              zshear_u = MAX( zustar(ji,jj)**2 * MAX( zdu_ml(ji,jj), 0._wp ) /  zhbl(ji,jj), 0._wp )
-                     !              zshear(ji,jj) = zshear(ji,jj) + zshear_u * ( 1.0 - MIN( zri_p(ji,jj) / rn_ri_p_thresh, 1.d0 )**2 )
-                     !              zshear(ji,jj) = MIN( zshear(ji,jj), zshear_u )
-
-                     !              zwb_shr = zwb_shr - 0.25 * MAX ( zshear_u, 0._wp) * ( 1.0 - MIN( zri_p(ji,jj) / rn_ri_p_thresh, 1._wp )**2 )
-                     !              zwb_shr = MAX( zwb_shr, -0.25 * zshear_u )
-#ifdef key_osm_debug
-                     IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
-                        WRITE(narea+100,'(3(a,g11.3))')'zdf_osm_osbl_state j_ddh(ji,jj) == 0:zwb_shr=',zwb_shr, &
-                           & '  zshear=',zshear(ji,jj),'  zshear_u=', zshear_u
-                        FLUSH(narea+100)
-                     END IF
-#endif
-
-                  ENDIF
-                  zwb_ent(ji,jj) = zwb_ent(ji,jj) + zwb_shr
-                  !           zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj)
-               ELSE    ! IF ( lconv ) THEN - ENDIF
-                  ! Stable OSBL  - shear production not coded for first attempt.
-               ENDIF  ! lconv
-            END IF  ! lshear
-            IF ( lconv(ji,jj) ) THEN
-               ! Unstable OSBL
-               zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * 2.0_wp * swbav(ji,jj)
-            END IF  ! lconv
-#ifdef key_osm_debug
-            IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
-               WRITE(narea+100,'(3(a,g11.3))')'end of zdf_osm_osbl_state:zwb_ent=',zwb_ent(ji,jj), &
-                  & '  zwb_min=',zwb_min(ji,jj), '  zwb0tot=', zwb0tot(ji,jj), '  swbav= ', swbav(ji,jj)
-               FLUSH(narea+100)
-            END IF
-#endif
-         END_2D
-         IF( ln_timing ) CALL timing_stop('zdf_osm_os')
-      END SUBROUTINE zdf_osm_osbl_state
-
 
       SUBROUTINE zdf_osm_osbl_state_fk( lpyc, lflux, lmle, zwb_fk )
@@ -2752,4 +2299,497 @@
       !
    END SUBROUTINE zdf_osm_velocity_rotation_3d
+
+   SUBROUTINE zdf_osm_osbl_state( Kmm, ldconv, ldshear, k_ddh, pwb_ent, pwb_min, pshear,   &
+      &                           phbl, phml, pdh, pdb_bl, pdu_bl, pdv_bl, pdb_ml, pdu_ml, pdv_ml )
+      !!---------------------------------------------------------------------
+      !!                 ***  ROUTINE zdf_osm_osbl_state  ***
+      !!
+      !! ** Purpose : Determines the state of the OSBL, stable/unstable,
+      !!              shear/ noshear. Also determines shear production,
+      !!              entrainment buoyancy flux and interfacial Richardson
+      !!              number
+      !!
+      !! ** Method  :
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER,                  INTENT(in   ) ::   Kmm       ! Ocean time-level index
+      LOGICAL,  DIMENSION(:,:), INTENT(  out) ::   ldconv    ! BL stability flags
+      LOGICAL,  DIMENSION(:,:), INTENT(  out) ::   ldshear   ! Shear layers
+      INTEGER,  DIMENSION(:,:), INTENT(  out) ::   k_ddh     ! k_ddh=0: active shear layer
+      !                                                      ! k_ddh=1: shear layer not active
+      !                                                      ! k_ddh=2: shear production low
+      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pwb_ent   ! Buoyancy fluxes at base
+      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pwb_min   !    of well-mixed layer
+      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear    ! Production of TKE due to shear across the pycnocline
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   phbl      ! BL depth
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   phml      ! ML depth
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pdh       ! Pycnocline depth
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pdb_bl    ! Buoyancy diff. between BL average and basal value
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pdu_bl    ! Velocity diff. (u) between BL average and basal value
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pdv_bl    ! Velocity diff. (v) between BL average and basal value
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pdb_ml    ! Buoyancy diff. between mixed-layer average and basal value
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pdu_ml    ! Velocity diff. (u) between mixed-layer average and basal value
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pdv_ml    ! Velocity diff. (v) between mixed-layer average and basal value
+      !
+      ! Local Variables
+      INTEGER :: jj, ji   ! Loop indices
+      !
+      REAL(wp), DIMENSION(A2D(0)) ::   zekman
+      REAL(wp), DIMENSION(A2D(0)) ::   zri_p, zri_b   ! Richardson numbers
+      REAL(wp)                    ::   zshear_u, zshear_v, zwb_shr
+      REAL(wp)                    ::   zwcor, zrf_conv, zrf_shear, zrf_langmuir, zr_stokes
+      !
+      REAL(wp), PARAMETER ::   za_shr         = 0.4_wp,  zb_shr    = 6.5_wp,  za_wb_s = 0.8_wp
+      REAL(wp), PARAMETER ::   zalpha_c       = 0.2_wp,  zalpha_lc = 0.03_wp
+      REAL(wp), PARAMETER ::   zalpha_ls      = 0.06_wp, zalpha_s  = 0.15_wp
+      REAL(wp), PARAMETER ::   rn_ri_p_thresh = 27.0_wp
+      REAL(wp), PARAMETER ::   zri_c          = 0.25_wp
+      REAL(wp), PARAMETER ::   zek            = 4.0_wp
+      REAL(wp), PARAMETER ::   zrot           = 0.0_wp   ! Dummy rotation rate of surface stress
+      REAL(wp), PARAMETER ::   zlarge         = -1e10_wp
+      !
+      IF( ln_timing ) CALL timing_start('zdf_osm_os')
+      !
+      ! Initialise INTENT(  out) arrays
+      ldconv(:,:)  = .FALSE.
+      ldshear(:,:) = .FALSE.
+      k_ddh(:,:)   = 1
+      pwb_ent(:,:) = zlarge
+      pwb_min(:,:) = zlarge
+      pshear(:,:)  = zlarge
+      !
+      ! Determins stability and set flag ldconv
+      DO_2D( 0, 0, 0, 0 )
+         IF ( shol(ji,jj) < 0._wp ) THEN
+            ldconv(ji,jj) = .TRUE.
+         ELSE
+            ldconv(ji,jj) = .FALSE.
+         ENDIF
+      END_2D
+      !
+      pshear(A2D(0)) = 0._wp
+      zekman(:,:) = EXP( -1.0_wp * zek * ABS( ff_t(A2D(0)) ) * phbl(A2D(0)) / MAX(sustar(A2D(0)), 1.e-8 ) )
+      !
+#ifdef key_osm_debug
+      IF(narea==nn_narea_db) THEN
+         ji=iloc_db; jj=jloc_db
+         WRITE(narea+100,'(a,g11.3)') &
+            & 'zdf_osm_osbl_state start: zekman=', zekman(ji,jj)
+         FLUSH(narea+100)
+      END IF
+#endif
+      !
+      DO_2D( 0, 0, 0, 0 )
+         IF ( ldconv(ji,jj) ) THEN
+            IF ( pdb_bl(ji,jj) > 0.0_wp ) THEN
+               zri_p(ji,jj) = MAX (  SQRT( pdb_bl(ji,jj) * pdh(ji,jj) / MAX( pdu_bl(ji,jj)**2 + pdv_bl(ji,jj)**2, 1e-8_wp ) ) *   &
+                  &                  ( phbl(ji,jj) / pdh(ji,jj) ) * ( svstr(ji,jj) / MAX( sustar(ji,jj), 1e-6_wp ) )**2 /         &
+                  &                  MAX( zekman(ji,jj), 1.0e-6_wp )  , 5.0_wp )
+               IF ( ff_t(ji,jj) >= 0.0_wp ) THEN   ! Northern hemisphere
+                  zri_b(ji,jj) = pdb_ml(ji,jj) * pdh(ji,jj) / ( MAX( pdu_ml(ji,jj), 1e-5_wp )**2 +   &
+                     &                                          MAX( -1.0_wp * pdv_ml(ji,jj), 1e-5_wp)**2 )
+               ELSE                                ! Southern hemisphere
+                  zri_b(ji,jj) = pdb_ml(ji,jj) * pdh(ji,jj) / ( MAX( pdu_ml(ji,jj), 1e-5_wp )**2 +   &
+                     &                                          MAX(           pdv_ml(ji,jj), 1e-5_wp)**2 )
+               END IF
+               pshear(ji,jj) = za_shr * zekman(ji,jj) *                                                   &
+                  &            ( MAX( sustar(ji,jj)**2 * pdu_ml(ji,jj) / phbl(ji,jj), 0.0_wp ) +          &
+                  &              zb_shr * MAX( -1.0_wp * ff_t(ji,jj) * sustke(ji,jj) * dstokes(ji,jj) *   &
+                  &                            pdv_ml(ji,jj) / phbl(ji,jj), 0.0_wp ) )
+#ifdef key_osm_debug
+               IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
+                  WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state 1st zshear: zshear=',pshear(ji,jj)
+                  WRITE(narea+100,'(2(a,g11.3))')'zdf_osm_osbl_state 1st zshear: zri_b=',zri_b(ji,jj),' zri_p=',zri_p(ji,jj)
+                  FLUSH(narea+100)
+               END IF
+#endif
+               ! Stability dependence
+               pshear(ji,jj) = pshear(ji,jj) * EXP( -0.75_wp * MAX( 0.0_wp, ( zri_b(ji,jj) - zri_c ) / zri_c ) )
+#ifdef key_osm_debug
+               IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
+                  WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state 1st zshear: zshear inc ri part=',pshear(ji,jj)
+                  FLUSH(narea+100)
+               END IF
+#endif
+               !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+               ! Test ensures k_ddh=0 is not selected. Change to zri_p<27 when  !
+               ! full code available                                          !
+               !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+               IF ( pshear(ji,jj) > 1e-10 ) THEN
+                  IF ( zri_p(ji,jj) < rn_ri_p_thresh .AND. MIN( hu(ji,jj,Kmm), hu(ji-1,jj,Kmm), hv(ji,jj,Kmm), hv(ji,jj-1,Kmm) ) > 100.0_wp ) THEN
+                     ! Growing shear layer
+                     k_ddh(ji,jj) = 0
+                     ldshear(ji,jj) = .TRUE.
+                  ELSE
+                     k_ddh(ji,jj) = 1
+                     !             IF ( zri_b <= 1.5 .and. pshear(ji,jj) > 0._wp ) THEN
+                     ! Shear production large enough to determine layer charcteristics, but can't maintain a shear layer
+                     ldshear(ji,jj) = .TRUE.
+                     !             ELSE
+                  END IF
+               ELSE
+                  k_ddh(ji,jj) = 2
+                  ldshear(ji,jj) = .FALSE.
+               END IF
+               ! Shear production may not be zero, but is small and doesn't determine characteristics of pycnocline
+               !               pshear(ji,jj) = 0.5 * pshear(ji,jj)
+               !               ldshear(ji,jj) = .FALSE.
+               !            ENDIF
+            ELSE   ! pdb_bl test, note pshear set to zero
+               k_ddh(ji,jj) = 2
+               ldshear(ji,jj) = .FALSE.
+            ENDIF
+         ENDIF
+      END_2D
+      !
+      ! Calculate entrainment buoyancy flux due to surface fluxes.
+      DO_2D( 0, 0, 0, 0 )
+         IF ( ldconv(ji,jj) ) THEN
+            zwcor        = ABS( ff_t(ji,jj) ) * phbl(ji,jj) + epsln
+            zrf_conv     = TANH( ( swstrc(ji,jj) / zwcor )**0.69_wp )
+            zrf_shear    = TANH( ( sustar(ji,jj) / zwcor )**0.69_wp )
+            zrf_langmuir = TANH( ( swstrl(ji,jj) / zwcor )**0.69_wp )
+            IF ( nn_osm_SD_reduce > 0 ) THEN
+               ! Effective Stokes drift already reduced from surface value
+               zr_stokes = 1.0_wp
+            ELSE
+               ! Effective Stokes drift only reduced by factor rn_zdfodm_adjust_sd,
+               ! requires further reduction where BL is deep
+               zr_stokes = 1.0 - EXP( -25.0_wp * dstokes(ji,jj) / hbl(ji,jj) * ( 1.0_wp + 4.0_wp * dstokes(ji,jj) / hbl(ji,jj) ) )
+            END IF
+            pwb_ent(ji,jj) = -2.0_wp * zalpha_c * zrf_conv * swbav(ji,jj) -                                          &
+               &             zalpha_s * zrf_shear * sustar(ji,jj)**3 / phml(ji,jj) +                                 &
+               &             zr_stokes * ( zalpha_s * EXP( -1.5_wp * sla(ji,jj) ) * zrf_shear * sustar(ji,jj)**3 -   &
+               &                           zrf_langmuir * zalpha_lc * swstrl(ji,jj)**3 ) / phml(ji,jj)
+#ifdef key_osm_debug
+            IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
+               WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state conv+shear0/lang: zwb_ent=',pwb_ent(ji,jj)
+               FLUSH(narea+100)
+            END IF
+#endif
+         ENDIF
+      END_2D
+      !
+      DO_2D( 0, 0, 0, 0 )
+         IF ( ldshear(ji,jj) ) THEN
+            IF ( ldconv(ji,jj) ) THEN
+               ! Unstable OSBL
+               zwb_shr = -1.0_wp * za_wb_s * zri_b(ji,jj) * pshear(ji,jj)
+#ifdef key_osm_debug
+               IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
+                  WRITE(narea+100,'(a,g11.3)')'zdf_osm_osbl_state 1st zwb_shr: zwb_shr=',zwb_shr
+                  FLUSH(narea+100)
+               END IF
+#endif
+               IF ( k_ddh(ji,jj) == 0 ) THEN
+                  ! Developing shear layer, additional shear production possible.
+
+                  !              pshear_u = MAX( zustar(ji,jj)**2 * MAX( pdu_ml(ji,jj), 0._wp ) /  phbl(ji,jj), 0._wp )
+                  !              pshear(ji,jj) = pshear(ji,jj) + pshear_u * ( 1.0 - MIN( zri_p(ji,jj) / rn_ri_p_thresh, 1.d0 )**2 )
+                  !              pshear(ji,jj) = MIN( pshear(ji,jj), pshear_u )
+
+                  !              zwb_shr = zwb_shr - 0.25 * MAX ( pshear_u, 0._wp) * ( 1.0 - MIN( zri_p(ji,jj) / rn_ri_p_thresh, 1._wp )**2 )
+                  !              zwb_shr = MAX( zwb_shr, -0.25 * pshear_u )
+#ifdef key_osm_debug
+                  IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
+                     WRITE(narea+100,'(3(a,g11.3))')'zdf_osm_osbl_state k_ddh(ji,jj) == 0:zwb_shr=',zwb_shr, &
+                        & '  zshear=',pshear(ji,jj),'  zshear_u=', pshear_u
+                     FLUSH(narea+100)
+                  END IF
+#endif
+               ENDIF
+               pwb_ent(ji,jj) = pwb_ent(ji,jj) + zwb_shr
+               !           pwb_min(ji,jj) = pwb_ent(ji,jj) + pdh(ji,jj) / phbl(ji,jj) * zwb0(ji,jj)
+            ELSE   ! IF ( ldconv ) THEN - ENDIF
+               ! Stable OSBL  - shear production not coded for first attempt.
+            ENDIF   ! ldconv
+         END IF   ! ldshear
+         IF ( ldconv(ji,jj) ) THEN
+            ! Unstable OSBL
+            pwb_min(ji,jj) = pwb_ent(ji,jj) + pdh(ji,jj) / phbl(ji,jj) * 2.0_wp * swbav(ji,jj)
+         END IF  ! ldconv
+#ifdef key_osm_debug
+         IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
+            WRITE(narea+100,'(3(a,g11.3))')'end of zdf_osm_osbl_state:zwb_ent=',pwb_ent(ji,jj), &
+               & '  zwb_min=',pwb_min(ji,jj), '  zwb0tot=', zwb0tot(ji,jj), '  swbav= ', swbav(ji,jj)
+            FLUSH(narea+100)
+         END IF
+#endif
+      END_2D
+      !
+      IF( ln_timing ) CALL timing_stop('zdf_osm_os')
+      !
+   END SUBROUTINE zdf_osm_osbl_state
+
+   SUBROUTINE zdf_osm_diffusivity_viscosity( Kbb, Kmm, kbld, kmld, ldconv, ldshear, ldpyc, ldcoup, k_ddh, pdiffut, pviscos,       &
+      &                                      phbl, phml, pdh, pdhdt, pshear, pb_bl, pu_bl, pv_bl, pb_ml, pu_ml, pv_ml, pwb_ent,   &
+      &                                      pwb_min )
+      !!---------------------------------------------------------------------
+      !!           ***  ROUTINE zdf_osm_diffusivity_viscosity  ***
+      !!
+      !! ** Purpose : Determines the eddy diffusivity and eddy viscosity
+      !!              profiles in the mixed layer and the pycnocline.
+      !!
+      !! ** Method  :
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER,                    INTENT(in   ) ::   Kbb, Kmm       ! Ocean time-level indices
+      INTEGER,  DIMENSION(:,:),   INTENT(in   ) ::   kbld           ! BL base layer
+      INTEGER,  DIMENSION(:,:),   INTENT(in   ) ::   kmld           ! ML base layer
+      LOGICAL,  DIMENSION(:,:),   INTENT(in   ) ::   ldconv         ! BL stability flags
+      LOGICAL,  DIMENSION(:,:),   INTENT(in   ) ::   ldshear        ! Shear layers
+      LOGICAL,  DIMENSION(:,:),   INTENT(in   ) ::   ldpyc          ! Pycnocline flags
+      LOGICAL,  DIMENSION(:,:),   INTENT(in   ) ::   ldcoup         ! Coupling to bottom
+      INTEGER,  DIMENSION(:,:),   INTENT(in   ) ::   k_ddh          ! Type of shear layer
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pdiffut        ! t-diffusivity
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pviscos        ! Viscosity
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   phbl           ! BL depth
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   phml           ! ML depth
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdh            ! Pycnocline depth
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdhdt          ! BL depth tendency
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pshear         ! Shear production
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pb_bl          ! Buoyancy average over the depth of the boundary layer
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pu_bl          ! Velocity average over the depth of the boundary layer
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pv_bl          ! Velocity average over the depth of the boundary layer
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pb_ml          ! Buoyancy average over the depth of the mixed layer
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pu_ml          ! Velocity average over the depth of the mixed layer
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pv_ml          ! Velocity average over the depth of the mixed layer
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pwb_ent        ! Buoyancy entrainment flux
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pwb_min
+      !
+      ! Local variables
+      INTEGER ::   ji, jj, jk   ! Loop indices
+      !
+      ! Scales used to calculate eddy diffusivity and viscosity profiles
+      REAL(wp), DIMENSION(A2D(0)) ::   zdifml_sc,    zvisml_sc
+      REAL(wp), DIMENSION(A2D(0)) ::   zdifpyc_n_sc, zdifpyc_s_sc
+      REAL(wp), DIMENSION(A2D(0)) ::   zvispyc_n_sc, zvispyc_s_sc
+      REAL(wp), DIMENSION(A2D(0)) ::   zbeta_d_sc,   zbeta_v_sc
+      REAL(wp), DIMENSION(A2D(0)) ::   zb_coup,      zc_coup_vis,  zc_coup_dif
+      !
+      REAL(wp) ::   zvel_sc_pyc, zvel_sc_ml, zstab_fac, zz_b
+      REAL(wp) ::   za_cubic, zb_d_cubic, zc_d_cubic, zd_d_cubic,   &   ! Coefficients in cubic polynomial specifying diffusivity
+         &                    zb_v_cubic, zc_v_cubic, zd_v_cubic        ! and viscosity in pycnocline
+      REAL(wp) ::   zznd_ml, zznd_pyc, ztmp
+      REAL(wp) ::   zmsku, zmskv
+      !
+      REAL(wp), PARAMETER :: rn_dif_ml     = 0.8_wp,  rn_vis_ml  = 0.375_wp
+      REAL(wp), PARAMETER :: rn_dif_pyc    = 0.15_wp, rn_vis_pyc = 0.142_wp
+      REAL(wp), PARAMETER :: rn_vispyc_shr = 0.15_wp
+      !
+      IF( ln_timing ) CALL timing_start('zdf_osm_dv')
+      !
+      zb_coup(:,:) = 0.0_wp
+      !
+      DO_2D( 0, 0, 0, 0 )
+         IF ( ldconv(ji,jj) ) THEN
+            !
+            zvel_sc_pyc = ( 0.15_wp * svstr(ji,jj)**3 + swstrc(ji,jj)**3 + 4.25_wp * pshear(ji,jj) * phbl(ji,jj) )**pthird
+            zvel_sc_ml  = ( svstr(ji,jj)**3 + 0.5_wp * swstrc(ji,jj)**3 )**pthird
+            zstab_fac   = ( phml(ji,jj) / zvel_sc_ml *   &
+               &            ( 1.4_wp - 0.4_wp / ( 1.0_wp + EXP(-3.5_wp * LOG10(-shol(ji,jj) ) ) )**1.25_wp ) )**2
+            !
+            zdifml_sc(ji,jj) = rn_dif_ml * phml(ji,jj) * zvel_sc_ml
+            zvisml_sc(ji,jj) = rn_vis_ml * zdifml_sc(ji,jj)
+#ifdef key_osm_debug
+            IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
+               WRITE(narea+100,'(2(a,g11.3))')'Start of 1st major loop of osm_diffusivity_viscositys, ldconv=T: zdifml_sc=',zdifml_sc(ji,jj),' zvisml_sc=',zvisml_sc(ji,jj)
+               WRITE(narea+100,'(3(a,g11.3))')'zvel_sc_pyc=',zvel_sc_pyc,' zvel_sc_ml=',zvel_sc_ml,' zstab_fac=',zstab_fac
+               FLUSH(narea+100)
+            END IF
+#endif
+            !
+            IF ( ldpyc(ji,jj) ) THEN
+               zdifpyc_n_sc(ji,jj) = rn_dif_pyc * zvel_sc_ml * pdh(ji,jj)
+               zvispyc_n_sc(ji,jj) = 0.09_wp * zvel_sc_pyc * ( 1.0_wp - phbl(ji,jj) / pdh(ji,jj) )**2 *                           &
+                  &                  ( 0.005 * ( pu_ml(ji,jj)-pu_bl(ji,jj) )**2 + 0.0075 * ( pv_ml(ji,jj)-pv_bl(ji,jj) )**2 ) /   &
+                  &                  pdh(ji,jj)
+               zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * pdh(ji,jj) + zvispyc_n_sc(ji,jj) * zstab_fac
+#ifdef key_osm_debug
+               IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
+                  WRITE(narea+100,'(2(a,g11.3))')' lpyc=ldconv=T, variables w/o shear contributions: zdifpyc_n_sc',zdifpyc_n_sc(ji,jj) ,' zvispyc_n_sc=',zvispyc_n_sc(ji,jj)
+                  FLUSH(narea+100)
+               END IF
+#endif
+               !
+               IF ( ldshear(ji,jj) .AND. k_ddh(ji,jj) /= 2 ) THEN
+                  ztmp = rn_vispyc_shr * ( pshear(ji,jj) * phbl(ji,jj) )**pthird * phbl(ji,jj)
+                  zdifpyc_n_sc(ji,jj) = zdifpyc_n_sc(ji,jj) + ztmp
+                  zvispyc_n_sc(ji,jj) = zvispyc_n_sc(ji,jj) + ztmp
+               ENDIF
+#ifdef key_osm_debug
+               IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
+                  WRITE(narea+100,'(2(a,g11.3))')' lpyc=ldconv=T, variables w shear contributions: zdifpyc_n_sc',zdifpyc_n_sc(ji,jj) ,' zvispyc_n_sc=',zvispyc_n_sc(ji,jj)
+                  FLUSH(narea+100)
+               END IF
+#endif
+               !
+               zdifpyc_s_sc(ji,jj) = pwb_ent(ji,jj) + 0.0025_wp * zvel_sc_pyc * ( phbl(ji,jj) / pdh(ji,jj) - 1.0_wp ) *   &
+                  &                                   ( pb_ml(ji,jj) - pb_bl(ji,jj) )
+               zvispyc_s_sc(ji,jj) = 0.09_wp * ( pwb_min(ji,jj) + 0.0025_wp * zvel_sc_pyc *                 &
+                  &                                               ( phbl(ji,jj) / pdh(ji,jj) - 1.0_wp ) *   &
+                  &                                               ( pb_ml(ji,jj) - pb_bl(ji,jj) ) )
+#ifdef key_osm_debug
+               IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
+                  WRITE(narea+100,'(2(a,g11.3))')' 1st shot at: zdifpyc_s_sc',zdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',zvispyc_s_sc(ji,jj)
+                  FLUSH(narea+100)
+               END IF
+#endif
+               zdifpyc_s_sc(ji,jj) = 0.09_wp * zdifpyc_s_sc(ji,jj) * zstab_fac
+               zvispyc_s_sc(ji,jj) = zvispyc_s_sc(ji,jj) * zstab_fac
+#ifdef key_osm_debug
+               IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
+                  WRITE(narea+100,'(2(a,g11.3))')' 2nd shot at: zdifpyc_s_sc',zdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',zvispyc_s_sc(ji,jj)
+                  FLUSH(narea+100)
+               END IF
+#endif
+               !
+               zdifpyc_s_sc(ji,jj) = MAX( zdifpyc_s_sc(ji,jj), -0.5_wp * zdifpyc_n_sc(ji,jj) )
+               zvispyc_s_sc(ji,jj) = MAX( zvispyc_s_sc(ji,jj), -0.5_wp * zvispyc_n_sc(ji,jj) )
+               
+               zbeta_d_sc(ji,jj) = 1.0_wp - ( ( zdifpyc_n_sc(ji,jj) + 1.4_wp * zdifpyc_s_sc(ji,jj) ) /   &
+                  &                           ( zdifml_sc(ji,jj) + epsln ) )**p2third
+               zbeta_v_sc(ji,jj) = 1.0_wp - 2.0_wp * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) / ( zvisml_sc(ji,jj) + epsln )
+            ELSE
+               zdifpyc_n_sc(ji,jj) = rn_dif_pyc * zvel_sc_ml * pdh(ji,jj)   ! ag 19/03
+               zdifpyc_s_sc(ji,jj) = 0.0_wp   ! ag 19/03
+               zvispyc_n_sc(ji,jj) = rn_vis_pyc * zvel_sc_ml * pdh(ji,jj)   ! ag 19/03
+               zvispyc_s_sc(ji,jj) = 0.0_wp   ! ag 19/03
+               IF(ldcoup(ji,jj) ) THEN   ! ag 19/03
+                  ! code from SUBROUTINE tke_tke zdftke.F90; uses bottom drag velocity rCdU_bot(ji,jj) = -Cd|ub|
+                  !     already calculated at T-points in SUBROUTINE zdf_drg from zdfdrg.F90
+                  !  Gives friction velocity sqrt bottom drag/rho_0 i.e. u* = SQRT(rCdU_bot*ub)
+                  ! wet-cell averaging ..
+                  zmsku = 0.5_wp * ( 2.0_wp - umask(ji-1,jj,mbkt(ji,jj)) * umask(ji,jj,mbkt(ji,jj)) )
+                  zmskv = 0.5_wp * ( 2.0_wp - vmask(ji,jj-1,mbkt(ji,jj)) * vmask(ji,jj,mbkt(ji,jj)) )
+                  zb_coup(ji,jj) = 0.4_wp * SQRT(-1.0_wp * rCdU_bot(ji,jj) *   &
+                     &             SQRT(  ( zmsku*( uu(ji,jj,mbkt(ji,jj),Kbb)+uu(ji-1,jj,mbkt(ji,jj),Kbb) ) )**2   &
+                     &                  + ( zmskv*( vv(ji,jj,mbkt(ji,jj),Kbb)+vv(ji,jj-1,mbkt(ji,jj),Kbb) ) )**2  ) )
+                  
+                  zz_b = -1.0_wp * gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)   ! ag 19/03
+                  zc_coup_vis(ji,jj) = -0.5_wp * ( 0.5_wp * zvisml_sc(ji,jj) / phml(ji,jj) - zb_coup(ji,jj) ) /   &
+                     &                 ( phml(ji,jj) + zz_b )   ! ag 19/03
+#ifdef key_osm_debug
+                  IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
+                     WRITE(narea+100,'(4(a,g11.3))')' lcoup = T; 1st pz_b= ', zz_b, ' pb_coup ', zb_coup(ji,jj),   &
+                        &                           ' pc_coup_vis ', zc_coup_vis(ji,jj), ' rCdU_bot ',rCdU_bot(ji,jj)
+                     WRITE(narea+100,'(2(a,g11.3))')' zmsku ', zmsku, ' zmskv ', zmskv
+                     FLUSH(narea+100)
+                  END IF
+#endif
+                  zz_b = -1.0_wp * phml(ji,jj) + gdepw(ji,jj,mbkt(ji,jj)+1,Kmm)   ! ag 19/03
+                  zbeta_v_sc(ji,jj) = 1.0_wp - 2.0_wp * ( zb_coup(ji,jj) * zz_b + zc_coup_vis(ji,jj) * zz_b**2 ) /   &
+                     &                                  zvisml_sc(ji,jj)   ! ag 19/03
+                  zbeta_d_sc(ji,jj) = 1.0_wp - ( ( zb_coup(ji,jj) * zz_b + zc_coup_vis(ji,jj) * zz_b**2 ) /   &
+                     &                           zdifml_sc(ji,jj) )**p2third
+                  zc_coup_dif(ji,jj) = 0.5_wp * ( -zdifml_sc(ji,jj) / phml(ji,jj) * ( 1.0_wp - zbeta_d_sc(ji,jj) )**1.5_wp +   &
+                     &                 1.5_wp * ( zdifml_sc(ji,jj) / phml(ji,jj) ) * zbeta_d_sc(ji,jj) *   &
+                     &                          SQRT( 1.0_wp - zbeta_d_sc(ji,jj) ) - zb_coup(ji,jj) ) / zz_b   ! ag 19/03
+#ifdef key_osm_debug
+                  IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
+                     WRITE(narea+100,'(2(a,g11.3))')' 2nd pz_b= ', zz_b, ' pc_coup_dif', zc_coup_dif(ji,jj)
+                     FLUSH(narea+100)
+                  END IF
+#endif
+               ELSE   ! ag 19/03
+                  zbeta_d_sc(ji,jj) = 1.0_wp - ( ( zdifpyc_n_sc(ji,jj) + 1.4_wp * zdifpyc_s_sc(ji,jj) ) /   &
+                     &                           ( zdifml_sc(ji,jj) + epsln ) )**p2third   ! ag 19/03
+                  zbeta_v_sc(ji,jj) = 1.0_wp - 2.0_wp * ( zvispyc_n_sc(ji,jj) + zvispyc_s_sc(ji,jj) ) /   &
+                     &                         ( zvisml_sc(ji,jj) + epsln )   ! ag 19/03
+               ENDIF   ! ag 19/03
+            ENDIF      ! ag 19/03
+#ifdef key_osm_debug
+            IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
+               WRITE(narea+100,'(2(a,g11.3))')'ldconv=T: zbeta_d_sc',zbeta_d_sc(ji,jj) ,' zbeta_v_sc=',zbeta_v_sc(ji,jj)
+               WRITE(narea+100,'(2(a,g11.3))')' Final zdifpyc_n_sc',zdifpyc_n_sc(ji,jj) ,' zvispyc_n_sc=',zvispyc_n_sc(ji,jj)
+               WRITE(narea+100,'(2(a,g11.3))')' Final zdifpyc_s_sc',zdifpyc_s_sc(ji,jj) ,' zvispyc_s_sc=',zvispyc_s_sc(ji,jj)
+               FLUSH(narea+100)
+            END IF
+#endif
+         ELSE
+            zdifml_sc(ji,jj) = svstr(ji,jj) * phbl(ji,jj) * MAX( EXP ( -( shol(ji,jj) / 0.6_wp )**2 ), 0.2_wp)
+            zvisml_sc(ji,jj) = zdifml_sc(ji,jj)
+#ifdef key_osm_debug
+            IF(narea==nn_narea_db.and.ji==iloc_db.and.jj==jloc_db)THEN
+               WRITE(narea+100,'(a,g11.3)')'End of 1st major loop of osm_diffusivity_viscositys, ldconv=F: zdifml_sc=',zdifml_sc(ji,jj),' zvisml_sc=',zvisml_sc(ji,jj)
+               FLUSH(narea+100)
+            END IF
+#endif
+         END IF
+      END_2D
+      !
+      DO_2D( 0, 0, 0, 0 )
+      IF ( ldconv(ji,jj) ) THEN
+         DO jk = 2, kmld(ji,jj)   ! Mixed layer diffusivity
+            zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj)
+            pdiffut(ji,jj,jk) = zdifml_sc(ji,jj) * zznd_ml * ( 1.0_wp - zbeta_d_sc(ji,jj) * zznd_ml )**1.5
+            pviscos(ji,jj,jk) = zvisml_sc(ji,jj) * zznd_ml * ( 1.0_wp - zbeta_v_sc(ji,jj) * zznd_ml ) *   &
+               &                ( 1.0_wp - 0.5_wp * zznd_ml**2 )
+         END DO
+         !
+         ! Coupling to bottom
+         !
+         IF ( ldcoup(ji,jj) ) THEN                                                         ! ag 19/03
+            DO jk = mbkt(ji,jj), kmld(ji,jj), -1                                           ! ag 19/03
+               zz_b = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - gdepw(ji,jj,mbkt(ji,jj)+1,Kmm) )   ! ag 19/03
+               pviscos(ji,jj,jk) = zb_coup(ji,jj) * zz_b + zc_coup_vis(ji,jj) * zz_b**2    ! ag 19/03
+               pdiffut(ji,jj,jk) = zb_coup(ji,jj) * zz_b + zc_coup_dif(ji,jj) * zz_b**2    ! ag 19/03
+            END DO                                                                         ! ag 19/03
+         ENDIF                                                                             ! ag 19/03
+         ! Pycnocline
+         IF ( ldpyc(ji,jj) ) THEN 
+            ! Diffusivity and viscosity profiles in the pycnocline given by
+            ! cubic polynomial. Note, if ldpyc TRUE can't be coupled to seabed.
+            za_cubic = 0.5_wp
+            zb_d_cubic = -1.75_wp * zdifpyc_s_sc(ji,jj) / zdifpyc_n_sc(ji,jj)
+            zd_d_cubic = ( pdh(ji,jj) * zdifml_sc(ji,jj) / phml(ji,jj) * SQRT( 1.0_wp - zbeta_d_sc(ji,jj) ) *   &
+               &           ( 2.5_wp * zbeta_d_sc(ji,jj) - 1.0_wp ) - 0.85_wp * zdifpyc_s_sc(ji,jj) ) /            &
+               &           MAX( zdifpyc_n_sc(ji,jj), 1.0e-8_wp )
+            zd_d_cubic = zd_d_cubic - zb_d_cubic - 2.0_wp * ( 1.0_wp - za_cubic  - zb_d_cubic )
+            zc_d_cubic = 1.0_wp - za_cubic - zb_d_cubic - zd_d_cubic
+            zb_v_cubic = -1.75_wp * zvispyc_s_sc(ji,jj) / zvispyc_n_sc(ji,jj)
+            zd_v_cubic = ( 0.5_wp * zvisml_sc(ji,jj) * pdh(ji,jj) / phml(ji,jj) - 0.85_wp * zvispyc_s_sc(ji,jj) ) /   &
+               &           MAX( zvispyc_n_sc(ji,jj), 1.0e-8_wp )
+            zd_v_cubic = zd_v_cubic - zb_v_cubic - 2.0_wp * ( 1.0_wp - za_cubic - zb_v_cubic )
+            zc_v_cubic = 1.0_wp - za_cubic - zb_v_cubic - zd_v_cubic
+            DO jk = kmld(ji,jj) , kbld(ji,jj)
+               zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / MAX(pdh(ji,jj), 1.0e-6_wp )
+               ztmp = ( 1.75_wp * zznd_pyc - 0.15_wp * zznd_pyc**2 - 0.2_wp * zznd_pyc**3 )
+               !
+               pdiffut(ji,jj,jk) = zdifpyc_n_sc(ji,jj) *   &
+                  &                ( za_cubic + zb_d_cubic * zznd_pyc + zc_d_cubic * zznd_pyc**2 + zd_d_cubic * zznd_pyc**3 )
+               !
+               pdiffut(ji,jj,jk) = pdiffut(ji,jj,jk) + zdifpyc_s_sc(ji,jj) * ztmp
+               pviscos(ji,jj,jk) = zvispyc_n_sc(ji,jj) *   &
+                  &                ( za_cubic + zb_v_cubic * zznd_pyc + zc_v_cubic * zznd_pyc**2 + zd_v_cubic * zznd_pyc**3 )
+               pviscos(ji,jj,jk) = pviscos(ji,jj,jk) + zvispyc_s_sc(ji,jj) * ztmp
+            END DO
+!                  IF ( pdhdt(ji,jj) > 0._wp ) THEN
+!                     zdiffut(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * pdhdt(ji,jj) * e3w(ji,jj,ibld(ji,jj)+1,Kmm), 1.0e-6 )
+!                     zviscos(ji,jj,ibld(ji,jj)+1) = MAX( 0.5 * pdhdt(ji,jj) * e3w(ji,jj,ibld(ji,jj)+1,Kmm), 1.0e-6 )
+!                  ELSE
+!                     zdiffut(ji,jj,ibld(ji,jj)) = 0._wp
+!                     zviscos(ji,jj,ibld(ji,jj)) = 0._wp
+!                  ENDIF
+         ENDIF
+      ELSE
+         ! Stable conditions
+         DO jk = 2, kbld(ji,jj)
+            zznd_ml = gdepw(ji,jj,jk,Kmm) / phbl(ji,jj)
+            pdiffut(ji,jj,jk) = 0.75_wp * zdifml_sc(ji,jj) * zznd_ml * ( 1.0_wp - zznd_ml )**1.5_wp
+            pviscos(ji,jj,jk) = 0.375_wp * zvisml_sc(ji,jj) * zznd_ml * ( 1.0_wp - zznd_ml ) * ( 1.0_wp - zznd_ml**2 )
+         END DO
+         !
+         IF ( pdhdt(ji,jj) > 0.0_wp ) THEN
+            pdiffut(ji,jj,kbld(ji,jj)) = MAX( pdhdt(ji,jj), 1.0e-6_wp) * e3w(ji, jj, kbld(ji,jj), Kmm)
+            pviscos(ji,jj,kbld(ji,jj)) = pdiffut(ji,jj,kbld(ji,jj))
+         ENDIF
+      ENDIF   ! End if ( ldconv )
+      !
+      END_2D
+      IF( iom_use("pb_coup") ) CALL iom_put( "pb_coup", tmask(:,:,1) * zb_coup(:,:) )   ! BBL-coupling velocity scale
+      IF( ln_timing ) CALL timing_stop('zdf_osm_dv')
+      !
+   END SUBROUTINE zdf_osm_diffusivity_viscosity
 
    SUBROUTINE zdf_osm_fgr_terms( Kmm, kbld, kmld, kp_ext, ldconv, ldpyc, k_ddh, phbl, phml, pdh, pdhdt, pshear,             &
