Index: /NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF/zdfosm.F90
===================================================================
--- /NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF/zdfosm.F90	(revision 14404)
+++ /NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF/zdfosm.F90	(revision 14405)
@@ -288,5 +288,4 @@
       REAL(wp), DIMENSION(jpi,jpj) :: zdh   ! pycnocline depth - grid
       REAL(wp), DIMENSION(jpi,jpj) :: zdhdt ! BL depth tendency
-      REAL(wp), DIMENSION(jpi,jpj) :: zddhdt                                    ! correction to dhdt due to internal structure.
       REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_bl_ext,zdsdz_bl_ext,zdbdz_bl_ext              ! external temperature/salinity and buoyancy gradients
       REAL(wp), DIMENSION(jpi,jpj) :: zdtdz_mle_ext,zdsdz_mle_ext,zdbdz_mle_ext              ! external temperature/salinity and buoyancy gradients
@@ -395,6 +394,5 @@
       ghams(:,:,:)   = 0._wp ; ghamu(:,:,:)   = 0._wp ; ghamv(:,:,:) = 0._wp
 
-      zddhdt(:,:) = 0._wp
-      ! hbl = MAX(hbl,epsln)
+       ! hbl = MAX(hbl,epsln)
       !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
       ! Calculate boundary layer scales
@@ -675,5 +673,5 @@
       CALL zdf_osm_vertical_average(imld-1, ibld-imld+1, zt_ml, zs_ml, zb_ml, zu_ml, zv_ml, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml)
 ! Rate of change of hbl
-      CALL zdf_osm_calculate_dhdt( zdhdt, zddhdt )
+      CALL zdf_osm_calculate_dhdt( zdhdt )
       DO jj = 2, jpjm1
         DO ji = 2, jpim1
@@ -856,19 +854,21 @@
                     ghams(ji,jj,jk) = ghams(ji,jj,jk) - 0.045 * ( ( zws_ent * zdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) * MAX( ( 1.75 * zznd_pyc -0.15 * zznd_pyc**2 - 0.2 * zznd_pyc**3 ), 0.0 )
                   END DO
+                  !
+                  IF ( dh(ji,jj)  <  0.2*hbl(ji,jj) ) THEN
+                     zbuoy_pyc_sc = zalpha_pyc(ji,jj) * zdb_ml(ji,jj) / zdh(ji,jj) + zdbdz_bl_ext(ji,jj)
+                     zdelta_pyc = ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird / SQRT( MAX( zbuoy_pyc_sc, ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**p2third / zdh(ji,jj)**2 ) )
 !
-                  zbuoy_pyc_sc = zalpha_pyc(ji,jj) * zdb_ml(ji,jj) / zdh(ji,jj) + zdbdz_bl_ext(ji,jj)
-                  zdelta_pyc = ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird / SQRT( MAX( zbuoy_pyc_sc, ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**p2third / zdh(ji,jj)**2 ) )
+                     zwt_pyc_sc_1 = 0.325 * ( zalpha_pyc(ji,jj) * zdt_ml(ji,jj) / zdh(ji,jj) + zdtdz_bl_ext(ji,jj) ) * zdelta_pyc**2 / zdh(ji,jj)
 !
-                  zwt_pyc_sc_1 = 0.325 * ( zalpha_pyc(ji,jj) * zdt_ml(ji,jj) / zdh(ji,jj) + zdtdz_bl_ext(ji,jj) ) * zdelta_pyc**2 / zdh(ji,jj)
+                     zws_pyc_sc_1 = 0.325 * ( zalpha_pyc(ji,jj) * zds_ml(ji,jj) / zdh(ji,jj) + zdsdz_bl_ext(ji,jj) ) * zdelta_pyc**2 / zdh(ji,jj)
 !
-                  zws_pyc_sc_1 = 0.325 * ( zalpha_pyc(ji,jj) * zds_ml(ji,jj) / zdh(ji,jj) + zdsdz_bl_ext(ji,jj) ) * zdelta_pyc**2 / zdh(ji,jj)
+                     zzeta_pyc = 0.15 - 0.175 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) ) 
+                     DO jk = 2, ibld(ji,jj)
+                        zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj)
+                        ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.05 * zwt_pyc_sc_1 * EXP( -0.25 * ( zznd_pyc / zzeta_pyc )**2 ) * zdh(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird
 !
-                  zzeta_pyc = 0.15 - 0.175 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) ) 
-                  DO jk = 2, ibld(ji,jj)
-                    zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj)
-                    ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.05 * zwt_pyc_sc_1 * EXP( -0.25 * ( zznd_pyc / zzeta_pyc )**2 ) * zdh(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird
-!
-                    ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.05 * zws_pyc_sc_1 * EXP( -0.25 * ( zznd_pyc / zzeta_pyc )**2 ) * zdh(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird
-                  END DO
+                        ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.05 * zws_pyc_sc_1 * EXP( -0.25 * ( zznd_pyc / zzeta_pyc )**2 ) * zdh(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird
+                     END DO
+                  END IF
                ENDIF ! End of pycnocline                  
              ELSE ! lconv test - stable conditions
@@ -912,6 +912,7 @@
        DO jj = 2, jpjm1
          DO ji = 2, jpim1
-           IF ( lpyc(ji,jj) ) THEN
-             IF ( j_ddh(ji,jj) == 0 ) THEN
+           IF( lconv(ji,jj) ) THEN
+             IF ( lpyc(ji,jj) ) THEN
+               IF ( j_ddh(ji,jj) == 0 ) THEN
 ! Place holding code. Parametrization needs checking for these conditions.
                zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) ))**pthird
@@ -940,5 +941,6 @@
              END DO
            ENDIF  ! lpyc
-         END DO ! ji loop
+           ENDIF ! lconv
+        END DO ! ji loop
        END DO  ! jj loop
 
@@ -1029,4 +1031,5 @@
                        & + 0.3 * 0.1 * ( EXP( -zznd_d ) + EXP( -5.0 * ( 1.0 - zznd_ml ) ) ) * zsc_vw_1(ji,jj)
                END DO
+
              ELSE
                DO jk = 2, ibld(ji,jj)
@@ -1052,5 +1055,26 @@
        END DO
 
-       IF(ln_dia_osm) THEN
+!        DO jj = 1, jpjm1
+!          DO ji = 1, jpim1
+!            IF ( lconv(ji,jj) ) THEN
+!              IF ( lpyc(ji,jj) ) THEN
+!                zd_cubic = ( 0.948 - 2.13 * zdh(ji,jj) / zhml(ji,jj) ) * zustar(ji,jj)**2
+!                zc_cubic = -0.474 * zustar(ji,jj)**2 - zd_cubic
+!                DO jk = imld(ji,jj), ibld(ji,jj)
+!                  zznd_pyc = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zdh(ji,jj)
+!                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 )
+!                END DO
+!                zc_cubic= 3.0 * ff_t(ji,jj) * zustar(ji,jj) * zhml(ji,jj)
+!                zd_cubic = -2.0 * ff_t(ji,jj) * zustar(ji,jj) * zhml(ji,jj)
+!                DO jk = imld(ji,jj), ibld(ji,jj)
+!                  zznd_pyc = -( gdepw_n(ji,jj,jk)-zhbl(ji,jj) ) / zdh(ji,jj)
+!                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0.3 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 )
+!                END DO
+!              ENDIF
+!            ENDIF
+!          END DO
+!        END DO
+
+       IF(ln_dSia_osm) THEN
           IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask*ghamu )
           IF ( iom_use("ghamv_f") ) CALL iom_put( "ghamv_f", wmask*ghamv )
@@ -1388,5 +1412,5 @@
                  zdifpyc_n_sc(ji,jj) =  rn_dif_pyc * zvel_sc_ml * zdh(ji,jj)
 
-                 IF ( lshear(ji,jj) .and. j_ddh(ji,jj) == 1 ) THEN
+                 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)
                  ENDIF
@@ -1398,5 +1422,5 @@
                  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
-                 IF ( lshear(ji,jj) .and. j_ddh(ji,jj) == 1 ) THEN
+                 IF ( lshear(ji,jj) .AND. j_ddh(ji,jj) /= 2 ) THEN
                    zvispyc_n_sc(ji,jj) = zvispyc_n_sc(ji,jj) + rn_vispyc_shr * ( zshear(ji,jj) * zhbl(ji,jj ) )**pthird * zhbl(ji,jj)
                  ENDIF
@@ -1508,13 +1532,14 @@
      
      REAL(wp), DIMENSION(jpi,jpj) :: zekman
-     REAL(wp) :: zri_p, zri_b   ! Richardson numbers
+     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.1
-     REAL, PARAMETER :: rn_ri_thres_a = 0.5, rn_ri_thresh_b = 0.59
-     REAL, PARAMETER :: zalpha_c = 0.2, zalpha_lc = 0.04     
+     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.
      
@@ -1530,5 +1555,5 @@
      END DO       
  
-     zekman(:,:) = EXP( - 4.0 * ABS( ff_t(:,:) ) * zhbl(:,:) / MAX(zustar(:,:), 1.e-8 ) )
+     zekman(:,:) = EXP( - zek * ABS( ff_t(:,:) ) * zhbl(:,:) / MAX(zustar(:,:), 1.e-8 ) )
      
      WHERE ( lconv )
@@ -1543,15 +1568,18 @@
          IF ( lconv(ji,jj) ) THEN
             IF ( zdb_bl(ji,jj) > 0._wp ) THEN
-              zri_p = 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) ) * ( zvstr(ji,jj) / MAX( zustar(ji,jj), 1.e-6 ) )**2 &
+              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) ) * ( zvstr(ji,jj) / MAX( zustar(ji,jj), 1.e-6 ) )**2 &
                    & / MAX( zekman(ji,jj), 1.e-6 )  , 5._wp )
             
-              zri_b = zdb_ml(ji,jj) * zdh(ji,jj) / MAX( zdu_ml(ji,jj)**2 + zdv_ml(ji,jj)**2, 1.e-8 )
+              zri_b(ji,jj) = zdb_ml(ji,jj) * zdh(ji,jj) / ( MAX( zdu_ml(ji,jj), 1.e-5 )**2 + MAX( zdv_ml(ji,jj), 1.e-5)**2 )
                         
               zshear(ji,jj) = za_shr * zekman(ji,jj) * ( MAX( zustar(ji,jj)**2 * zdu_ml(ji,jj) / zhbl(ji,jj), 0._wp ) + zb_shr * MAX( -ff_t(ji,jj) * zustke(ji,jj) * dstokes(ji,jj) * zdv_ml(ji,jj) / zhbl(ji,jj), 0._wp ) )
+! Stability Dependence
+              zshear(ji,jj) = zshear(ji,jj) * EXP( -0.75 * MAX(0._wp,( zri_b(ji,jj) - zri_c ) / zri_c ) )
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 ! Test ensures j_ddh=0 is not selected. Change to zri_p<27 when  !
 ! full code available                                          !
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
-              IF ( zri_p < -rn_ri_p_thresh .and. zshear(ji,jj) > 0._wp ) THEN
+            IF ( zshear(ji,jj) > 1.e-10 ) THEN
+              IF ( zri_p(ji,jj) < rn_ri_p_thresh ) THEN
 ! Growing shear layer
                 j_ddh(ji,jj) = 0
@@ -1559,13 +1587,17 @@
               ELSE
                 j_ddh(ji,jj) = 1
-                IF ( zri_b <= 1.5 .and. zshear(ji,jj) > 0._wp ) THEN
+!                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
+                lshear(ji,jj) = .TRUE.
+!                ELSE
+              ENDIF
+            ELSE
+                j_ddh(ji,jj) = 2
+                lshear(ji,jj) = .FALSE.
+            ENDIF
 ! 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 
-              ENDIF                
+!                  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
@@ -1594,8 +1626,8 @@
             &                  * ( 1.0 + 4.0 * dstokes(ji,jj) / hbl(ji,jj) ) )
            END IF
-           zwb_ent(ji,jj) = - 2.0 * 0.2 * zrf_conv * zwbav(ji,jj) &
-                  &                  - 0.15 * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) &
-                  &         + zr_stokes * ( 0.15 * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 &
-                  &                                         - zrf_langmuir * 0.03 * zwstrl(ji,jj)**3 ) / zhml(ji,jj)
+           zwb_ent(ji,jj) = - 2.0 * zalpha_c * zrf_conv * zwbav(ji,jj) &
+                  &                  - zalpha_s * zrf_shear * zustar(ji,jj)**3 /zhml(ji,jj) &
+                  &         + zr_stokes * ( zalpha_s * EXP( -1.5 * zla(ji,jj) ) * zrf_shear * zustar(ji,jj)**3 &
+                  &                                         - zrf_langmuir * zalpha_lc * zwstrl(ji,jj)**3 ) / zhml(ji,jj)
              !
          ENDIF
@@ -1610,29 +1642,27 @@
            IF ( lconv(ji,jj) ) THEN
 ! Unstable OSBL
-              zwb_shr = -za_wb_s * zshear(ji,jj)
+              zwb_shr = -za_wb_s * zri_b(ji,jj) * zshear(ji,jj)
               IF ( j_ddh(ji,jj) == 0 ) THEN
 
-! Developing shear layer, additional shear production possible.
-
-                zshear_u = MAX( zustar(ji,jj)**2 * zdu_ml(ji,jj) /  zhbl(ji,jj), 0._wp )
-                zshear(ji,jj) = zshear(ji,jj) + zshear_u * ( 1.0 - MIN( zri_p / rn_ri_p_thresh, 1.d0 ) )
-                zshear(ji,jj) = MIN( zshear(ji,jj), zshear_u )
+! ! 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 = -za_wb_s * zshear(ji,jj)
+!                 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 )
                 
               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)
+!              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
-         ELSE  ! lshear
-           IF ( lconv(ji,jj) ) THEN
+         ENDIF  ! lshear
+         IF ( lconv(ji,jj) ) THEN
 ! Unstable OSBL
-              zwb_shr = -za_wb_s * zshear(ji,jj)
-              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)
-           ENDIF  ! lconv
-         ENDIF    ! lshear
+            zwb_min(ji,jj) = zwb_ent(ji,jj) + zdh(ji,jj) / zhbl(ji,jj) * zwb0(ji,jj)
+         ENDIF  ! lconv
        END DO   ! ji
      END DO     ! jj
@@ -2033,5 +2063,5 @@
     END SUBROUTINE zdf_osm_pycnocline_shear_profiles
 
-   SUBROUTINE zdf_osm_calculate_dhdt( zdhdt, zddhdt )
+   SUBROUTINE zdf_osm_calculate_dhdt( zdhdt )
      !!---------------------------------------------------------------------
      !!                   ***  ROUTINE zdf_osm_calculate_dhdt  ***
@@ -2043,13 +2073,13 @@
      !!----------------------------------------------------------------------
 
-    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt, zddhdt        ! Rate of change of hbl
+    REAL(wp), DIMENSION(jpi,jpj) :: zdhdt        ! Rate of change of hbl
 
     INTEGER :: jj, ji
     REAL(wp) :: zgamma_b_nd, zgamma_dh_nd, zpert, zpsi
-    REAL(wp) :: zvel_max!, zwb_min
-    REAL(wp) :: zzeta_m = 0.3
-    REAL(wp) :: zgamma_c = 2.0
-    REAL(wp) :: zdhoh = 0.1
-    REAL(wp) :: alpha_bc = 0.5
+    REAL(wp) :: zvel_max,zddhdt
+    REAL(wp), PARAMETER :: zzeta_m = 0.3
+    REAL(wp), PARAMETER :: zgamma_c = 2.0
+    REAL(wp), PARAMETER :: zdhoh = 0.1
+    REAL(wp), PARAMETER :: zalpha_b = 0.3
     REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth
  
@@ -2073,14 +2103,9 @@
                    ! OSBL is deepening, entrainment > restratification
                    IF ( zdb_bl(ji,jj) > 0.0 .and. zdbdz_bl_ext(ji,jj) > 0.0 ) THEN
-! *** Used for shear Needs to be changed to work stabily
-!                zgamma_b_nd = zdbdz_bl_ext * dh / zdb_ml
-!                zalpha_b = 6.7 * zgamma_b_nd / ( 1.0 + zgamma_b_nd )
-!                zgamma_b = zgamma_b_nd / ( 0.12 * ( 1.25 + zgamma_b_nd ) )
-!                za_1 = 1.0 / zgamma_b**2 - 0.017
-!                za_2 = 1.0 / zgamma_b**3 - 0.0025
-!                zpsi = zalpha_b * ( 1.0 + zgamma_b_nd ) * ( za_1 - 2.0 * za_2 * dh / hbl )
-                      zpsi = 0._wp
-                      zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) )
-                      zdhdt(ji,jj) = zdhdt(ji,jj)! - zpsi * ( -1.0 / zhml(ji,jj) + 2.4 * zdbdz_bl_ext(ji,jj) / zdb_ml(ji,jj) ) * zwb_min(ji,jj) * zdh(ji,jj) / zdb_bl(ji,jj)
+                      zgamma_b_nd = zdbdz_bl_ext(ji,jj) + zdh(ji,jj) / zdb_ml(ji,jj)
+                      zpsi = -zalpha_pyc(ji,jj) * ( zwb0(ji,jj) - ( zwb_min(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) ) * zdh(ji,jj) / zhbl(ji,jj)
+                      zpsi = zpsi - 4.0 * ( 1.0 + zdh(ji,jj) /zhbl(ji,jj) ) * zgamma_b_nd * ( zwb_min(ji,jj) + 2.0 * zwb_fk_b(ji,jj) )
+                      zpsi = zalpha_b * MAX ( zpsi, 0._wp )
+                      zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) / ( zvel_max + MAX(zdb_bl(ji,jj), 1.0e-15) ) + zpsi / MAX( zdb_bl(ji,jj), 1.e-15 )
                       IF ( j_ddh(ji,jj) == 1 ) THEN
                         IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN
@@ -2090,13 +2115,16 @@
                         ENDIF
 ! Relaxation to dh_ref = zari * hbl
-                        zddhdt(ji,jj) = -a_ddh_2 * ( 1.0 - zdh(ji,jj) / ( zari * zhbl(ji,jj) ) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj)
+                        zddhdt = -a_ddh_2 * ( 1.0 - zdh(ji,jj) / ( zari * zhbl(ji,jj) ) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj)
                         
-                      ELSE  ! j_ddh == 0
+                      ELSE IF ( j_ddh(ji,jj) == 0 ) THEN
 ! Growing shear layer
-                        zddhdt(ji,jj) = -a_ddh * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj)
+                        zddhdt = -a_ddh * ( 1.0 - 1.6 * zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj)
+                        zddhdt = EXP( - 4.0 * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX(zustar(ji,jj), 1.e-8 ) ) * zddhdt
+                      ELSE
+                        zddhdt = 0._wp
                       ENDIF ! j_ddh
-                        zdhdt(ji,jj) = zdhdt(ji,jj) ! + zpsi * zddhdt(ji,jj)
+                      zdhdt(ji,jj) = zdhdt(ji,jj) + zalpha_b * ( zalpha_pyc(ji,jj) * zdb_ml(ji,jj) + 2.0 * zdbdz_bl_ext(ji,jj) * zdh(ji,jj) ) * zddhdt / zdb_bl(ji,jj)
                    ELSE    ! zdb_bl >0
-                      zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) /  MAX( zvel_max, 1.0e-15)
+                       zdhdt(ji,jj) = -( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(ji,jj) ) /  MAX( zvel_max, 1.0e-15)
                    ENDIF
                 ELSE   ! zwb_min + 2*zwb_fk_b < 0
@@ -2302,6 +2330,6 @@
       INTEGER :: jj, ji
       INTEGER :: inhml
-      REAL(wp) :: zari, ztau, zdh_ref
-      REAL, PARAMETER :: a_ddh_2 = 3.5 ! also in pycnocline_depth
+      REAL(wp) :: zari, ztau, zdh_ref, zddhdt
+      REAL, PARAMETER :: a_ddh = 2.5, a_ddh_2 = 3.5 ! also in pycnocline_depth
 
     DO jj = 2, jpjm1
@@ -2310,7 +2338,14 @@
          IF ( lshear(ji,jj) ) THEN
             IF ( lconv(ji,jj) ) THEN
+             IF ( zdb_bl(ji,jj) > 0._wp ) THEN
               IF ( j_ddh(ji,jj) == 0 ) THEN
 ! ddhdt for pycnocline determined in osm_calculate_dhdt
-                dh(ji,jj) = dh(ji,jj) + zddhdt(ji,jj) * rn_rdt
+!                zddhdt = -a_ddh * ( 1.0 - 1.6 * zdh(ji,jj) / zhbl(ji,jj) ) * zwb_ent(ji,jj) / zdb_bl(ji,jj)
+!                zddhdt = EXP( - 4.0 * ABS( ff_t(ji,jj) ) * zhbl(ji,jj) / MAX(zustar(ji,jj), 1.e-8 ) ) * zddhdt
+! maximum limit for how thick the shear layer can grow relative to the thickness of the boundary kayer
+!                dh(ji,jj) = MIN( dh(ji,jj) + zddhdt * rn_rdt, 0.625 * hbl(ji,jj) )
+                ztau = MAX( zdb_bl(ji,jj) * ( 0.625 * hbl(ji,jj) ) / ( a_ddh * MAX(-zwb_ent(ji,jj), 1.e-12) ), 2.0 * rn_rdt )
+                dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + 0.625 * zhbl(ji,jj) * ( 1.0 - EXP( -rn_rdt / ztau ) )
+                IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = 0.8 * zhbl(ji,jj)
               ELSE
 ! Temporary (probably) Recalculate dh_ref to ensure dh doesn't go negative. Can't do this using zddhdt from calculate_dhdt 
@@ -2324,5 +2359,9 @@
                 IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zari * zhbl(ji,jj)
               ENDIF
-               
+             ELSE
+              ztau = MAX( MAX( hbl(ji,jj) / ( zvstr(ji,jj)**3 + 0.5 * zwstrc(ji,jj)**3 )**pthird, epsln), 2.0 * rn_rdt )
+              dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + 0.2 * zhbl(ji,jj) * ( 1.0 - EXP( -rn_rdt / ztau ) )
+              IF ( dh(ji,jj) > hbl(ji,jj) ) dh(ji,jj) = 0.2 * hbl(ji,jj)
+             ENDIF
             ELSE ! lconv
 ! Initially shear only for entraining OSBL. Stable code will be needed if extended to stable OSBL 
