Index: NEMO/branches/2020/dev_r14122_ENHANCE-14_smueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90
===================================================================
--- NEMO/branches/2020/dev_r14122_ENHANCE-14_smueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90	(revision 14287)
+++ NEMO/branches/2020/dev_r14122_ENHANCE-14_smueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90	(revision 14305)
@@ -113,4 +113,6 @@
    INTEGER  ::   nn_osm_SD_reduce   ! = 0/1/2 flag for getting effective stokes drift from surface value
    LOGICAL  ::   ln_dia_osm         ! Use namelist  rn_osm_la
+   LOGICAL  ::   ln_dia_pyc_scl = .FALSE.   ! Output of pycnocline scalar-gradient profiles
+   LOGICAL  ::   ln_dia_pyc_shr = .FALSE.   ! Output of pycnocline velocity-shear  profiles
 
 
@@ -215,9 +217,9 @@
       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::  p_avm, p_avt   ! momentum and tracer Kz (w-points)
       !!
-      INTEGER ::   ji, jj, jk                   ! dummy loop indices
+      INTEGER ::   ji, jj, jk, jkflt                               ! dummy loop indices
 
       INTEGER ::   jl                   ! dummy loop indices
 
-      INTEGER ::   ikbot, jkmax, jkm1, jkp2     !
+      INTEGER ::   ikbot, jkm1, jkp2     !
 
       REAL(wp) ::   ztx, zty, zflageos, zstabl, zbuofdep,zucube      !
@@ -233,7 +235,7 @@
       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)   :: zrad0                       ! Surface solar temperature flux (deg m/s)
+      REAL(wp)   :: zradh                       ! Radiative flux at bl base (Buoyancy units)
+      REAL(wp)   :: zradav                      ! Radiative flux, bl average (Buoyancy Units)
       REAL(wp), DIMENSION(jpi,jpj) :: zustar    ! friction velocity
       REAL(wp), DIMENSION(jpi,jpj) :: zwstrl    ! Langmuir velocity scale
@@ -241,5 +243,5 @@
       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)   :: zvw0                        ! Surface v-momentum flux
       REAL(wp), DIMENSION(jpi,jpj) :: zwth0     ! Surface heat flux (Kinematic)
       REAL(wp), DIMENSION(jpi,jpj) :: zws0      ! Surface freshwater flux
@@ -299,9 +301,9 @@
       REAL(wp) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline
       REAL(wp) :: zuw_bse,zvw_bse  ! momentum fluxes at the top of the pycnocline
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz_pyc    ! parametrized gradient of temperature in pycnocline
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdsdz_pyc    ! parametrised gradient of salinity in pycnocline
+      REAL(wp) ::   zdtdz_pyc    ! Parametrized gradient of temperature in pycnocline
+      REAL(wp) ::   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
+      REAL(wp) ::   zdudz_pyc    ! u-shear across the pycnocline
+      REAL(wp) ::   zdvdz_pyc    ! v-shear across the pycnocline
       REAL(wp), DIMENSION(jpi,jpj) :: zdbds_mle    ! Magnitude of horizontal buoyancy gradient.
       ! Flux-gradient relationship variables
@@ -317,7 +319,7 @@
 
       ! 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
+      REAL(wp), DIMENSION(jpi,jpj) ::   z2du     ! u-shear^2
+      REAL(wp), DIMENSION(jpi,jpj) ::   z2dv     ! v-shear^2
+      REAL(wp)                     ::   zrimix   ! Spatial form of ri#-induced diffusion
 
       ! Temporary variables
@@ -342,4 +344,6 @@
       REAL(wp) :: zsqrtpi, z_two_thirds, zproportion, ztransp, zthickness
       REAL(wp) :: z2k_times_thickness, zsqrt_depth, zexp_depth, zdstokes0, zf, zexperfc
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   z3ddz_pyc_1, z3ddz_pyc_2   ! Pycnocline gradient/shear profiles
+      INTEGER  ::   istat                                                     ! Memory allocation status
 
       ! For debugging
@@ -349,7 +353,7 @@
       IF( ln_timing ) CALL timing_start('zdf_osm')
       ibld(:,:)   = 0     ; imld(:,:)  = 0
-      zrad0(:,:)  = 0._wp ; zradh(:,:) = 0._wp ; zradav(:,:)    = 0._wp ; zustar(:,:)    = 0._wp
+      zustar(:,:) = 0.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
+      zwth0(:,:)  = 0.0_wp ; zws0(:,:)  = 0.0_wp ; zwb0(:,:)      = 0.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
@@ -370,6 +374,5 @@
       zwth_ent = 0._wp ; zws_ent = 0._wp
       !
-      zdtdz_pyc(:,:,:) = 0._wp ; zdsdz_pyc(:,:,:) = 0._wp ; zdbdz_pyc(:,:,:) = 0._wp
-      zdudz_pyc(:,:,:) = 0._wp ; zdvdz_pyc(:,:,:) = 0._wp
+      zdbdz_pyc(:,:,:) = 0.0_wp
       !
       zdtdz_bl_ext(:,:) = 0._wp ; zdsdz_bl_ext(:,:) = 0._wp ; zdbdz_bl_ext(:,:) = 0._wp
@@ -396,40 +399,38 @@
       ! 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
+      !
+      ! Turbulent surface fluxes and fluxes averaged over depth of the OSBL
+      zz0 =           rn_abs   ! Assume two-band radiation model for depth of OSBL - surface equi-partition in 2-bands
+      zz1 =  1.0_wp - rn_abs
+      DO_2D( 0, 0, 0, 0 )
+         zrad0         = qsr(ji,jj) * r1_rho0_rcp                          ! Surface downward irradiance (so always +ve)
+         zradh         = zrad0 *                                       &   ! Downwards irradiance at base of boundary layer
+            &            ( zz0 * EXP( -1.0_wp * hbl(ji,jj) / rn_si0 ) + zz1 * EXP( -1.0_wp * hbl(ji,jj) / rn_si1 ) )
+         zradav        = zrad0 *                                       &   ! Downwards irradiance averaged over depth of the OSBL
+            &     ( zz0 * ( 1.0_wp - EXP( -hbl(ji,jj)/rn_si0 ) ) * rn_si0 +   &
+            &       zz1 * ( 1.0_wp - EXP( -hbl(ji,jj)/rn_si1 ) ) * rn_si1 ) / hbl(ji,jj)
+         zwth0(ji,jj)  = - qns(ji,jj) * r1_rho0_rcp * tmask(ji,jj,1)       ! Upwards surface Temperature flux for non-local term
+         zwthav(ji,jj) = 0.5_wp * zwth0(ji,jj) -                       &   ! Turbulent heat flux averaged over depth of OSBL
+            &            ( 0.5_wp * ( zrad0 + zradh ) - zradav )
+     END_2D
      DO_2D( 0, 0, 0, 0 )
-        ! Surface downward irradiance (so always +ve)
-        zrad0(ji,jj) = qsr(ji,jj) * r1_rho0_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)
+        zws0(ji,jj)    = -1.0_wp *                                     &   ! Upwards surface salinity flux for non-local term
+           &          ( ( emp(ji,jj) - rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm) + sfx(ji,jj) ) * r1_rho0 * tmask(ji,jj,1)
+        zthermal       = rab_n(ji,jj,1,jp_tem)
+        zbeta          = rab_n(ji,jj,1,jp_sal)
+        zwb0(ji,jj)    = grav * zthermal * zwth0(ji,jj) -              &   ! Non radiative upwards surface buoyancy flux
+           &             grav * zbeta * zws0(ji,jj)
+        zwsav(ji,jj)   = 0.5 * zws0(ji,jj)                                 ! Turbulent salinity flux averaged over depth of the OBSL
+        zwbav(ji,jj)   = grav  * zthermal * zwthav(ji,jj) -            &   ! Turbulent buoyancy flux averaged over the depth of the
+           &             grav  * zbeta * zwsav(ji,jj)                      ! OBSBL
      END_2D
-     ! Turbulent surface fluxes and fluxes averaged over depth of the OSBL
      DO_2D( 0, 0, 0, 0 )
-        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_rho0_rcp * tmask(ji,jj,1)
-        ! Upwards surface salinity flux for non-local term
-        zws0(ji,jj) = - ( ( emp(ji,jj)-rnf(ji,jj) ) * ts(ji,jj,1,jp_sal,Kmm)  + sfx(ji,jj) ) * r1_rho0 * 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) = - 0.5 * (utau(ji-1,jj) + utau(ji,jj)) * r1_rho0 * tmask(ji,jj,1)
-        zvw0(ji,jj) = - 0.5 * (vtau(ji,jj-1) + vtau(ji,jj)) * r1_rho0 * tmask(ji,jj,1)
-        ! Friction velocity (zustar), at T-point : LMD94 eq. 2
-        zustar(ji,jj) = MAX( SQRT( SQRT( zuw0(ji,jj) * zuw0(ji,jj) + zvw0(ji,jj) * zvw0(ji,jj) ) ), 1.0e-8 )
+        zuw0(ji,jj)    = - 0.5 * (utau(ji-1,jj) + utau(ji,jj)) *       &   ! Surface upward velocity fluxes
+           &             r1_rho0 * tmask(ji,jj,1)
+        zvw0           = - 0.5 * (vtau(ji,jj-1) + vtau(ji,jj)) * r1_rho0 * tmask(ji,jj,1)
+        zustar(ji,jj)  = MAX( SQRT( SQRT( zuw0(ji,jj) *                &   ! Friction velocity (zustar), at T-point : LMD94 eq. 2
+           &                              zuw0(ji,jj) + zvw0 * zvw0 ) ), 1.0e-8_wp )
         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) )
+        zsin_wind(ji,jj) = -zvw0        / ( zustar(ji,jj) * zustar(ji,jj) )
      END_2D
      ! Calculate Stokes drift in direction of wind (zustke) and Stokes penetration depth (dstokes)
@@ -699,6 +700,5 @@
 
       CALL zdf_osm_external_gradients( ibld+2, zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext )
-      CALL zdf_osm_pycnocline_scalar_profiles( zdtdz_pyc, zdsdz_pyc, zdbdz_pyc, zalpha_pyc )
-      CALL zdf_osm_pycnocline_shear_profiles( zdudz_pyc, zdvdz_pyc )
+      CALL zdf_osm_pycnocline_buoyancy_profiles( zdbdz_pyc, zalpha_pyc )
        !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
        ! Eddy viscosity/diffusivity and non-gradient terms in the flux-gradient relationship
@@ -1028,22 +1028,116 @@
       END_2D
 
-      ! pynocline contributions
-       DO_2D( 0, 0, 0, 0 )
-         IF ( .not. lconv(ji,jj) ) THEN
-          IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
-             DO jk= 2, ibld(ji,jj)
-                znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
-                ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc(ji,jj,jk)
-                ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc(ji,jj,jk)
-                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc(ji,jj,jk)
-                ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc(ji,jj,jk)
-             END DO
-          END IF
+      ! Pynocline contributions
+      !
+      IF ( ln_dia_pyc_scl .OR. ln_dia_pyc_shr ) THEN   ! Allocate arrays used for output of pycnocline gradient/shear profiles
+         ALLOCATE( z3ddz_pyc_1(jpi,jpj,jpk), z3ddz_pyc_2(jpi,jpj,jpk), STAT=istat )
+         CALL mpp_sum( 'zdfosm', istat )
+         IF ( istat /= 0 ) CALL ctl_stop( 'zdf_osm: failed to allocate temporary arrays' )
+         z3ddz_pyc_1(:,:,:) = 0.0_wp
+         z3ddz_pyc_2(:,:,:) = 0.0_wp
+      END IF
+      DO_2D( 0, 0, 0, 0 )
+         IF ( lconv (ji,jj) ) THEN
+            ! Unstable conditions. Shouldn;t be needed with no pycnocline code.
+            !                  zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / &
+            !                       &      ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * &
+            !                      &      MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 ))
+            !Alan is this right?
+            !                  zvgrad = ( 0.7 * zdv_ml(ji,jj) + &
+            !                       &    2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / &
+            !                       &          ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird  + epsln ) &
+            !                       &      )/ (zdh(ji,jj)  + epsln )
+            !                  DO jk = 2, ibld(ji,jj) - 1 + ibld_ext
+            !                     znd = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v
+            !                     IF ( znd <= 0.0 ) THEN
+            !                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd )
+            !                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd )
+            !                     ELSE
+            !                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd )
+            !                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd )
+            !                     ENDIF
+            !                  END DO
+         ELSE   ! Stable conditions
+            IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
+               ! If pycnocline profile only defined when depth steady of increasing.
+               IF ( zdhdt(ji,jj) > 0.0_wp ) THEN   ! Depth increasing, or steady.
+                  IF ( zdb_bl(ji,jj) > 0.0_wp ) THEN
+                     IF ( zhol(ji,jj) >= 0.5_wp ) THEN   ! Very stable - 'thick' pycnocline
+                        ztmp = 1.0_wp / MAX( zhbl(ji,jj), epsln )
+                        ztgrad = zdt_bl(ji,jj) * ztmp
+                        zsgrad = zds_bl(ji,jj) * ztmp
+                        zbgrad = zdb_bl(ji,jj) * ztmp
+                        DO jk = 2, ibld(ji,jj)
+                           znd = gdepw(ji,jj,jk,Kmm) * ztmp
+                           zdtdz_pyc =  ztgrad * EXP( -15.0 * ( znd - 0.9_wp )**2 )
+                           zdsdz_pyc =  zsgrad * EXP( -15.0 * ( znd - 0.9_wp )**2 )
+                           ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc
+                           ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc
+                           IF ( ln_dia_pyc_scl ) THEN
+                              z3ddz_pyc_1(ji,jj,jk) = zdtdz_pyc
+                              z3ddz_pyc_2(ji,jj,jk) = zdsdz_pyc
+                           END IF
+                        END DO
+                     ELSE   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form.
+                        ztmp = 1.0_wp / MAX( zdh(ji,jj), epsln )
+                        ztgrad = zdt_bl(ji,jj) * ztmp
+                        zsgrad = zds_bl(ji,jj) * ztmp
+                        zbgrad = zdb_bl(ji,jj) * ztmp
+                        DO jk = 2, ibld(ji,jj)
+                           znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) * ztmp
+                           zdtdz_pyc =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
+                           zdsdz_pyc =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
+                           ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zdiffut(ji,jj,jk) * zdtdz_pyc
+                           ghams(ji,jj,jk) = ghams(ji,jj,jk) + zdiffut(ji,jj,jk) * zdsdz_pyc
+                           IF ( ln_dia_pyc_scl ) THEN
+                              z3ddz_pyc_1(ji,jj,jk) = zdtdz_pyc
+                              z3ddz_pyc_2(ji,jj,jk) = zdsdz_pyc
+                           END IF
+                        END DO
+                     ENDIF   ! IF (zhol >=0.5)
+                  ENDIF      ! IF (zdb_bl> 0.)
+               ENDIF         ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero
+            END IF
          END IF
-       END_2D
+      END_2D
+      IF ( ln_dia_pyc_scl ) THEN   ! Output of pycnocline gradient profiles
+         IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask(:,:,:) * z3ddz_pyc_1(:,:,:) )
+         IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask(:,:,:) * z3ddz_pyc_2(:,:,:) )
+         IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask(:,:,:) * zdbdz_pyc(:,:,:)   )
+      END IF
+      DO_2D( 0, 0, 0, 0 )
+         IF ( .NOT. lconv (ji,jj) ) THEN
+            IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
+               zugrad = 3.25_wp * zdu_bl(ji,jj) / zhbl(ji,jj)
+               zvgrad = 2.75_wp * zdv_bl(ji,jj) / zhbl(ji,jj)
+               DO jk = 2, ibld(ji,jj)
+                  znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
+                  IF ( znd < 1.0 ) THEN
+                     zdudz_pyc = zugrad * EXP( -40.0_wp * ( znd - 1.0_wp )**2 )
+                  ELSE
+                     zdudz_pyc = zugrad * EXP( -20.0_wp * ( znd - 1.0_wp )**2 )
+                  ENDIF
+                  zdvdz_pyc = zvgrad * EXP( -20.0_wp * ( znd - 0.85_wp )**2 )
+                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zviscos(ji,jj,jk) * zdudz_pyc
+                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zviscos(ji,jj,jk) * zdvdz_pyc
+                  IF ( ln_dia_pyc_shr ) THEN
+                     z3ddz_pyc_1(ji,jj,jk) = zdudz_pyc
+                     z3ddz_pyc_2(ji,jj,jk) = zdvdz_pyc
+                  END IF
+               END DO
+            END IF
+         END IF
+      END_2D
+      IF ( ln_dia_pyc_shr ) THEN   ! Output of pycnocline shear profiles
+         IF ( iom_use("dudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask(:,:,:) * z3ddz_pyc_1(:,:,:) )
+         IF ( iom_use("dvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask(:,:,:) * z3ddz_pyc_2(:,:,:) )
+      END IF
       IF(ln_dia_osm) THEN
-          IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu )
-          IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv )
-       END IF
+         IF ( iom_use("ghamu_b") ) CALL iom_put( "ghamu_b", wmask*ghamu )
+         IF ( iom_use("ghamv_b") ) CALL iom_put( "ghamv_b", wmask*ghamv )
+      END IF
+      IF ( ln_dia_pyc_scl .OR. ln_dia_pyc_shr ) THEN   ! Deallocate arrays used for output of pycnocline gradient/shear profiles
+         DEALLOCATE( z3ddz_pyc_1, z3ddz_pyc_2 )
+      END IF
 
        DO_2D( 0, 0, 0, 0 )
@@ -1057,8 +1151,7 @@
           IF ( iom_use("ghamu_1") ) CALL iom_put( "ghamu_1", wmask*ghamu )
           IF ( iom_use("ghamv_1") ) CALL iom_put( "ghamv_1", wmask*ghamv )
-          IF ( iom_use("zdudz_pyc") ) CALL iom_put( "zdudz_pyc", wmask*zdudz_pyc )
-          IF ( iom_use("zdvdz_pyc") ) CALL iom_put( "zdvdz_pyc", wmask*zdvdz_pyc )
           IF ( iom_use("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos )
        END IF
+
        !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
        ! Need to put in code for contributions that are applied explicitly to
@@ -1080,39 +1173,39 @@
        END_2D
 
-       IF(ln_dia_osm) THEN
-          IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc )
-          IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc )
-          IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc )
-       END IF
-
-! KPP-style Ri# mixing
-       IF( ln_kpprimix) THEN
-          DO_3D( 1, 0, 1, 0, 2, jpkm1 )      !* Shear production at uw- and vw-points (energy conserving form)
-             z3du(ji,jj,jk) = 0.5 * (  uu(ji,jj,jk-1,Kmm) -  uu(ji  ,jj,jk,Kmm) )   &
-                  &                 * (  uu(ji,jj,jk-1,Kbb) -  uu(ji  ,jj,jk,Kbb) ) * wumask(ji,jj,jk) &
-                  &                 / (  e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) )
-             z3dv(ji,jj,jk) = 0.5 * (  vv(ji,jj,jk-1,Kmm) -  vv(ji,jj  ,jk,Kmm) )   &
-                  &                 * (  vv(ji,jj,jk-1,Kbb) -  vv(ji,jj  ,jk,Kbb) ) * wvmask(ji,jj,jk) &
-                  &                 / (  e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) )
-          END_3D
-      !
-         DO_3D( 0, 0, 0, 0, 2, jpkm1 )
-            !                                          ! 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_3D
-
+      ! KPP-style Ri# mixing
+      IF ( ln_kpprimix ) THEN
+          jkflt = jpk
           DO_2D( 0, 0, 0, 0 )
-             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
+             IF ( ibld(ji,jj) < jkflt ) jkflt = ibld(ji,jj)
           END_2D
-
+          DO jk = jkflt+1, jpkm1
+             ! Shear production at uw- and vw-points (energy conserving form)
+             DO_2D( 1, 0, 1, 0 )
+                IF ( jk > MIN( ibld(ji,jj), ibld(ji+1,jj) ) ) THEN
+                   z2du(ji,jj) = 0.5_wp * ( uu(ji,jj,jk-1,Kmm) - uu(ji,jj,jk,Kmm) ) *                      &
+                      &                   ( uu(ji,jj,jk-1,Kbb) - uu(ji,jj,jk,Kbb) ) * wumask(ji,jj,jk) /   &
+                      &                   ( e3uw(ji,jj,jk,Kmm) * e3uw(ji,jj,jk,Kbb) )
+                END IF
+                IF ( jk > MIN( ibld(ji,jj), ibld(ji,jj+1) ) ) THEN
+                   z2dv(ji,jj) = 0.5_wp * ( vv(ji,jj,jk-1,Kmm) - vv(ji,jj,jk,Kmm) ) *                      &
+                      &                   ( vv(ji,jj,jk-1,Kbb) - vv(ji,jj,jk,Kbb) ) * wvmask(ji,jj,jk) /   &
+                      &                   ( e3vw(ji,jj,jk,Kmm) * e3vw(ji,jj,jk,Kbb) )
+                END IF
+             END_2D
+             DO_2D( 0, 0, 0, 0 )
+                IF ( jk > ibld(ji,jj) ) THEN
+                   ! Shear prod. at w-point weightened by mask
+                   zesh2  =  ( z2du(ji-1,jj) + z2du(ji,jj) ) / MAX( 1.0_wp , umask(ji-1,jj,jk) + umask(ji,jj,jk) )   &
+                      &    + ( z2dv(ji,jj-1) + z2dv(ji,jj) ) / MAX( 1.0_wp , vmask(ji,jj-1,jk) + vmask(ji,jj,jk) )
+                   ! Local Richardson number
+                   zri   = MAX( rn2b(ji,jj,jk), 0.0_wp ) / MAX(zesh2, epsln)
+                   zfri =  MIN( zri / rn_riinfty , 1.0_wp )
+                   zfri  = ( 1.0_wp - zfri * zfri )
+                   zrimix  =  zfri * zfri  * zfri * wmask(ji, jj, jk)
+                   zdiffut(ji,jj,jk) = zrimix*rn_difri
+                   zviscos(ji,jj,jk) = zrimix*rn_difri
+                END IF
+             END_2D
+          END DO
        END IF ! ln_kpprimix = .true.
 
@@ -1157,11 +1250,4 @@
           END_2D
        ENDIF
-
-       IF(ln_dia_osm) THEN
-          IF ( iom_use("zdtdz_pyc") ) CALL iom_put( "zdtdz_pyc", wmask*zdtdz_pyc )
-          IF ( iom_use("zdsdz_pyc") ) CALL iom_put( "zdsdz_pyc", wmask*zdsdz_pyc )
-          IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask*zdbdz_pyc )
-       END IF
-
 
        ! Lateral boundary conditions on zvicos (sign unchanged), needed to caclulate viscosities on u and v grids
@@ -1727,145 +1813,77 @@
     END SUBROUTINE zdf_osm_external_gradients
 
-    SUBROUTINE zdf_osm_pycnocline_scalar_profiles( zdtdz, zdsdz, zdbdz, zalpha )
-
-     REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdtdz, zdsdz, zdbdz      ! gradients in the pycnocline
-     REAL(wp), DIMENSION(jpi,jpj) :: zalpha
-
-     INTEGER :: jk, jj, ji
-     REAL(wp) :: ztgrad, zsgrad, zbgrad
-     REAL(wp) :: zgamma_b_nd, znd
-     REAL(wp) :: zzeta_m, zzeta_en, zbuoy_pyc_sc
-     REAL(wp), PARAMETER :: zgamma_b = 2.25, zzeta_sh = 0.15
-
-     IF( ln_timing ) CALL timing_start('zdf_osm_pscp')
-     DO_2D( 0, 0, 0, 0 )
-        IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
-           IF ( lconv(ji,jj) ) THEN  ! convective conditions
-             IF ( lpyc(ji,jj) ) THEN
-                zzeta_m = 0.1 + 0.3 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) )
-                zalpha(ji,jj) = 2.0 * ( 1.0 - ( 0.80 * zzeta_m + 0.5 * SQRT( 3.14159 / zgamma_b ) ) * zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) ) / ( 0.723 + SQRT( 3.14159 / zgamma_b ) )
-                zalpha(ji,jj) = MAX( zalpha(ji,jj), 0._wp )
-
-                ztmp = 1._wp/MAX(zdh(ji,jj), epsln)
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-! Commented lines in this section are not needed in new code, once tested !
-! can be removed                                                          !
-!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-!                   ztgrad = zalpha * zdt_ml(ji,jj) * ztmp + zdtdz_bl_ext(ji,jj)
-!                   zsgrad = zalpha * zds_ml(ji,jj) * ztmp + zdsdz_bl_ext(ji,jj)
-                zbgrad = zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp + zdbdz_bl_ext(ji,jj)
-                zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln)
-                DO jk = 2, ibld(ji,jj)+ibld_ext
-                  znd = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) * ztmp
-                  IF ( znd <= zzeta_m ) THEN
-!                        zdtdz(ji,jj,jk) = zdtdz_bl_ext(ji,jj) + zalpha * zdt_ml(ji,jj) * ztmp * &
-!                &        EXP( -6.0 * ( znd -zzeta_m )**2 )
-!                        zdsdz(ji,jj,jk) = zdsdz_bl_ext(ji,jj) + zalpha * zds_ml(ji,jj) * ztmp * &
-!                                  & EXP( -6.0 * ( znd -zzeta_m )**2 )
-                     zdbdz(ji,jj,jk) = zdbdz_bl_ext(ji,jj) + zalpha(ji,jj) * zdb_ml(ji,jj) * ztmp * &
-                               & EXP( -6.0 * ( znd -zzeta_m )**2 )
-                  ELSE
-!                         zdtdz(ji,jj,jk) =  ztgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 )
-!                         zdsdz(ji,jj,jk) =  zsgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 )
-                      zdbdz(ji,jj,jk) =  zbgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 )
-                  ENDIF
-               END DO
-            ENDIF ! if no pycnocline pycnocline gradients set to zero
-           ELSE
-              ! stable conditions
-              ! if pycnocline profile only defined when depth steady of increasing.
-              IF ( zdhdt(ji,jj) > 0.0 ) THEN        ! Depth increasing, or steady.
-                 IF ( zdb_bl(ji,jj) > 0._wp ) THEN
-                    IF ( zhol(ji,jj) >= 0.5 ) THEN      ! Very stable - 'thick' pycnocline
-                       ztmp = 1._wp/MAX(zhbl(ji,jj), epsln)
-                       ztgrad = zdt_bl(ji,jj) * ztmp
-                       zsgrad = zds_bl(ji,jj) * ztmp
-                       zbgrad = zdb_bl(ji,jj) * ztmp
-                       DO jk = 2, ibld(ji,jj)
-                          znd = gdepw(ji,jj,jk,Kmm) * ztmp
-                          zdtdz(ji,jj,jk) =  ztgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
-                          zdbdz(ji,jj,jk) =  zbgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
-                          zdsdz(ji,jj,jk) =  zsgrad * EXP( -15.0 * ( znd - 0.9 )**2 )
-                       END DO
-                    ELSE                                   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form.
-                       ztmp = 1._wp/MAX(zdh(ji,jj), epsln)
-                       ztgrad = zdt_bl(ji,jj) * ztmp
-                       zsgrad = zds_bl(ji,jj) * ztmp
-                       zbgrad = zdb_bl(ji,jj) * ztmp
-                       DO jk = 2, ibld(ji,jj)
-                          znd = -( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) * ztmp
-                          zdtdz(ji,jj,jk) =  ztgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
-                          zdbdz(ji,jj,jk) =  zbgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
-                          zdsdz(ji,jj,jk) =  zsgrad * EXP( -1.75 * ( znd + 0.75 )**2 )
-                       END DO
-                    ENDIF ! IF (zhol >=0.5)
-                 ENDIF    ! IF (zdb_bl> 0.)
-              ENDIF       ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero
-           ENDIF          ! IF (lconv)
-        ENDIF      ! IF ( ibld(ji,jj) < mbkt(ji,jj) )
-     END_2D
-     IF( ln_timing ) CALL timing_stop('zdf_osm_pscp')
-
-    END SUBROUTINE zdf_osm_pycnocline_scalar_profiles
-
-    SUBROUTINE zdf_osm_pycnocline_shear_profiles( zdudz, zdvdz )
-      !!---------------------------------------------------------------------
-      !!                   ***  ROUTINE zdf_osm_pycnocline_shear_profiles  ***
-      !!
-      !! ** Purpose : Calculates velocity shear in the pycnocline
-      !!
-      !! ** Method  :
-      !!
-      !!----------------------------------------------------------------------
-
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdudz, zdvdz
-
-      INTEGER :: jk, jj, ji
-      REAL(wp) :: zugrad, zvgrad, znd
-      REAL(wp) :: zzeta_v = 0.45
+   SUBROUTINE zdf_osm_pycnocline_buoyancy_profiles( pdbdz, palpha )
+      REAL(wp), DIMENSION(:,:,:), INTENT( inout ) ::   pdbdz   ! Gradients in the pycnocline
+      REAL(wp), DIMENSION(:,:),   INTENT( inout ) ::   palpha
+      INTEGER                                     ::   jk, jj, ji
+      REAL(wp)                                    ::   zbgrad
+      REAL(wp)                                    ::   zgamma_b_nd, znd
+      REAL(wp)                                    ::   zzeta_m
+      REAL(wp), PARAMETER                         ::   ppgamma_b = 2.25_wp
       !
-      IF( ln_timing ) CALL timing_start('zdf_osm_pshp')
+      IF( ln_timing ) CALL timing_start('zdf_osm_pscp')
+      !
       DO_2D( 0, 0, 0, 0 )
-         !
          IF ( ibld(ji,jj) + jp_ext(ji,jj) < mbkt(ji,jj) ) THEN
-            IF ( lconv (ji,jj) ) THEN
-               ! Unstable conditions. Shouldn;t be needed with no pycnocline code.
-!                  zugrad = 0.7 * zdu_ml(ji,jj) / zdh(ji,jj) + 0.3 * zustar(ji,jj)*zustar(ji,jj) / &
-!                       &      ( ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird * zhml(ji,jj) ) * &
-!                      &      MIN(zla(ji,jj)**(8.0/3.0) + epsln, 0.12 ))
-               !Alan is this right?
-!                  zvgrad = ( 0.7 * zdv_ml(ji,jj) + &
-!                       &    2.0 * ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) / &
-!                       &          ( ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird  + epsln ) &
-!                       &      )/ (zdh(ji,jj)  + epsln )
-!                  DO jk = 2, ibld(ji,jj) - 1 + ibld_ext
-!                     znd = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / (zdh(ji,jj) + epsln ) - zzeta_v
-!                     IF ( znd <= 0.0 ) THEN
-!                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( 3.0 * znd )
-!                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( 3.0 * znd )
-!                     ELSE
-!                        zdudz(ji,jj,jk) = 1.25 * zugrad * EXP( -2.0 * znd )
-!                        zdvdz(ji,jj,jk) = 1.25 * zvgrad * EXP( -2.0 * znd )
-!                     ENDIF
-!                  END DO
-            ELSE
-               ! stable conditions
-               zugrad = 3.25 * zdu_bl(ji,jj) / zhbl(ji,jj)
-               zvgrad = 2.75 * zdv_bl(ji,jj) / zhbl(ji,jj)
-               DO jk = 2, ibld(ji,jj)
-                  znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
-                  IF ( znd < 1.0 ) THEN
-                     zdudz(ji,jj,jk) = zugrad * EXP( -40.0 * ( znd - 1.0 )**2 )
-                  ELSE
-                     zdudz(ji,jj,jk) = zugrad * EXP( -20.0 * ( znd - 1.0 )**2 )
-                  ENDIF
-                  zdvdz(ji,jj,jk) = zvgrad * EXP( -20.0 * ( znd - 0.85 )**2 )
-               END DO
-            ENDIF
-            !
-         END IF      ! IF ( ibld(ji,jj) + ibld_ext < mbkt(ji,jj) )
+            IF ( lconv(ji,jj) ) THEN   ! convective conditions
+               IF ( lpyc(ji,jj) ) THEN
+                  zzeta_m = 0.1_wp + 0.3_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * zhol(ji,jj) ) ) )
+                  palpha(ji,jj) = 2.0_wp * ( 1.0_wp - ( 0.80_wp * zzeta_m + 0.5_wp * SQRT( 3.14159_wp / ppgamma_b ) ) *   &
+                     &                                zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / zdb_ml(ji,jj) ) /                &
+                     &            ( 0.723_wp + SQRT( 3.14159_wp / ppgamma_b ) )
+                  palpha(ji,jj) = MAX( palpha(ji,jj), 0.0_wp )
+                  ztmp = 1.0_wp / MAX( zdh(ji,jj), epsln )
+                  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+                  ! Commented lines in this section are not needed in new code, once tested !
+                  ! can be removed                                                          !
+                  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+                  ! ztgrad = zalpha * zdt_ml(ji,jj) * ztmp + zdtdz_bl_ext(ji,jj)
+                  ! zsgrad = zalpha * zds_ml(ji,jj) * ztmp + zdsdz_bl_ext(ji,jj)
+                  zbgrad = palpha(ji,jj) * zdb_ml(ji,jj) * ztmp + zdbdz_bl_ext(ji,jj)
+                  zgamma_b_nd = zdbdz_bl_ext(ji,jj) * zdh(ji,jj) / MAX(zdb_ml(ji,jj), epsln)
+                  DO jk = 2, ibld(ji,jj)+ibld_ext
+                     znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) * ztmp
+                     IF ( znd <= zzeta_m ) THEN
+                        ! zdtdz(ji,jj,jk) = zdtdz_bl_ext(ji,jj) + zalpha * zdt_ml(ji,jj) * ztmp * &
+                        ! &        EXP( -6.0 * ( znd -zzeta_m )**2 )
+                        ! zdsdz(ji,jj,jk) = zdsdz_bl_ext(ji,jj) + zalpha * zds_ml(ji,jj) * ztmp * &
+                        ! & EXP( -6.0 * ( znd -zzeta_m )**2 )
+                        pdbdz(ji,jj,jk) = zdbdz_bl_ext(ji,jj) + palpha(ji,jj) * zdb_ml(ji,jj) * ztmp * &
+                           & EXP( -6.0_wp * ( znd -zzeta_m )**2 )
+                     ELSE
+                        ! zdtdz(ji,jj,jk) =  ztgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 )
+                        ! zdsdz(ji,jj,jk) =  zsgrad * EXP( -zgamma_b * ( znd - zzeta_m )**2 )
+                        pdbdz(ji,jj,jk) =  zbgrad * EXP( -1.0_wp * ppgamma_b * ( znd - zzeta_m )**2 )
+                     ENDIF
+                  END DO
+               ENDIF   ! If no pycnocline pycnocline gradients set to zero
+            ELSE   ! Stable conditions
+               ! If pycnocline profile only defined when depth steady of increasing.
+               IF ( zdhdt(ji,jj) > 0.0_wp ) THEN   ! Depth increasing, or steady.
+                  IF ( zdb_bl(ji,jj) > 0.0_wp ) THEN
+                     IF ( zhol(ji,jj) >= 0.5_wp ) THEN   ! Very stable - 'thick' pycnocline
+                        ztmp = 1.0_wp / MAX( zhbl(ji,jj), epsln )
+                        zbgrad = zdb_bl(ji,jj) * ztmp
+                        DO jk = 2, ibld(ji,jj)
+                           znd = gdepw(ji,jj,jk,Kmm) * ztmp
+                           pdbdz(ji,jj,jk) =  zbgrad * EXP( -15.0_wp * ( znd - 0.9_wp )**2 )
+                        END DO
+                     ELSE   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form.
+                        ztmp = 1.0_wp / MAX( zdh(ji,jj), epsln )
+                        zbgrad = zdb_bl(ji,jj) * ztmp
+                        DO jk = 2, ibld(ji,jj)
+                           znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - zhml(ji,jj) ) * ztmp
+                           pdbdz(ji,jj,jk) =  zbgrad * EXP( -1.75_wp * ( znd + 0.75_wp )**2 )
+                        END DO
+                     ENDIF   ! IF (zhol >=0.5)
+                  ENDIF      ! IF (zdb_bl> 0.)
+               ENDIF         ! IF (zdhdt >= 0) zdhdt < 0 not considered since pycnocline profile is zero and profile arrays are intialized to zero
+            ENDIF            ! IF (lconv)
+         ENDIF   ! IF ( ibld(ji,jj) < mbkt(ji,jj) )
       END_2D
-      IF( ln_timing ) CALL timing_stop('zdf_osm_pshp')
-    END SUBROUTINE zdf_osm_pycnocline_shear_profiles
+      !
+      IF( ln_timing ) CALL timing_stop('zdf_osm_pscp')
+      !
+   END SUBROUTINE zdf_osm_pycnocline_buoyancy_profiles
 
    SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zddhdt )
@@ -2609,4 +2627,8 @@
      END IF
 
+     ! Flags associated with diagnostic output
+     IF ( ln_dia_osm .AND. ( iom_use("zdudz_pyc") .OR. iom_use("zdvdz_pyc") ) )                            ln_dia_pyc_shr = .TRUE.
+     IF ( ln_dia_osm .AND. ( iom_use("zdtdz_pyc") .OR. iom_use("zdsdz_pyc") .OR. iom_use("zdbdz_pyc" ) ) ) ln_dia_pyc_scl = .TRUE.
+     
      !                              ! allocate zdfosm arrays
      IF( zdf_osm_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_osm_init : unable to allocate arrays' )
