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 13858)
+++ NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser_4.0/src/OCE/ZDF/zdfosm.F90	(revision 13859)
@@ -1049,4 +1049,24 @@
        END DO
 
+       DO jj = 2, jpjm1
+          DO ji = 2, jpim1
+            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= 0.2 * ff_t(ji,jj) * zustke(ji,jj) * zhml(ji,jj)
+              zd_cubic = -0.1 * ff_t(ji,jj) * zustke(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
+          END DO
+       END DO
+
+
        IF(ln_dia_osm) THEN
           IF ( iom_use("ghamu_f") ) CALL iom_put( "ghamu_f", wmask*ghamu )
@@ -1067,5 +1087,5 @@
             IF ( .not. lconv(ji,jj) ) THEN
                DO jk = 2, ibld(ji,jj)
-                  znd = ( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zhbl(ji,jj) !ALMG to think about
+                  znd = -( gdepw_n(ji,jj,jk) - zhbl(ji,jj) ) / zhbl(ji,jj) !ALMG to think about
                   IF ( znd >= 0.0 ) THEN
                      ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) )
@@ -1448,5 +1468,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 - zd_cubic )
+                   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)
@@ -2305,19 +2325,43 @@
          IF ( lshear(ji,jj) ) THEN
             IF ( lconv(ji,jj) ) THEN
-              IF ( j_ddh(ji,jj) == 0 ) THEN
+               IF( ln_osm_mle ) THEN
+                  IF ( ( zwb_ent(ji,jj) + 2.0 * zwb_fk_b(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
-              ELSE
+                      dh(ji,jj) = dh(ji,jj) + zddhdt(ji,jj) * rn_rdt
+                    ELSE
 ! Temporary (probably) Recalculate dh_ref to ensure dh doesn't go negative. Can't do this using zddhdt from calculate_dhdt 
-                IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN
-                  zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )
+                      IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN
+                        zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )
+                      ELSE
+                        zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )
+                      ENDIF
+                      ztau = MAX( zdb_bl(ji,jj) * ( zari * hbl(ji,jj) ) / ( a_ddh_2 * MAX(-zwb_ent(ji,jj), 1.e-12) ), 2.0 * rn_rdt )
+                      dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zari * zhbl(ji,jj) * ( 1.0 - EXP( -rn_rdt / ztau ) )
+ !                    IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zari * zhbl(ji,jj)!
+                    ENDIF ! j_ddh
+                    IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zari * zhbl(ji,jj)
+			      ELSE  ! dhdt inc FK <0                      
+                     ztau = 0.2 * hbl(ji,jj) / MAX(epsln, (zvstr(ji,jj)**3 + 0.5 *zwstrc(ji,jj)**3)**pthird)
+                     dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + 0.2 * zhbl(ji,jj) * ( 1.0 - EXP( -rn_rdt / ztau ) )
+ !                    IF ( zdhdt(ji,jj) < 0._wp .and. dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zari * zhbl(ji,jj)!
+                     IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = 0.2 * zhbl(ji,jj)
+				  ENDIF ! dhdt<0
+               ELSE ! ln_osm_mle
+                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
                 ELSE
-                  zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )
+! Temporary (probably) Recalculate dh_ref to ensure dh doesn't go negative. Can't do this using zddhdt from calculate_dhdt 
+                  IF ( ( zwstrc(ji,jj) / zvstr(ji,jj) )**3 <= 0.5 ) THEN
+                    zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zvstr(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )
+                  ELSE
+                    zari = MIN( 1.5 * zdb_bl(ji,jj) / ( zhbl(ji,jj) * ( MAX(zdbdz_bl_ext(ji,jj),0._wp) + zdb_bl(ji,jj)**2 / MAX(4.5 * zwstrc(ji,jj)**2 , 1.e-12 )) ), 0.2d0 )
+                  ENDIF
+                  ztau = MAX( zdb_bl(ji,jj) * ( zari * hbl(ji,jj) ) / ( a_ddh_2 * MAX(-zwb_ent(ji,jj), 1.e-12) ), 2.0 * rn_rdt )
+                  dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zari * zhbl(ji,jj) * ( 1.0 - EXP( -rn_rdt / ztau ) )
                 ENDIF
-                ztau = MAX( zdb_bl(ji,jj) * ( zari * hbl(ji,jj) ) / ( a_ddh_2 * MAX(-zwb_ent(ji,jj), 1.e-12) ), 2.0 * rn_rdt )
-                dh(ji,jj) = dh(ji,jj) * EXP( -rn_rdt / ztau ) + zari * zhbl(ji,jj) * ( 1.0 - EXP( -rn_rdt / ztau ) )
                 IF ( dh(ji,jj) >= hbl(ji,jj) ) dh(ji,jj) = zari * zhbl(ji,jj)
-              ENDIF
-               
+	       ENDIF   ! ln_osm_mle
             ELSE ! lconv
 ! Initially shear only for entraining OSBL. Stable code will be needed if extended to stable OSBL 
@@ -2415,5 +2459,5 @@
 
     END SUBROUTINE zdf_osm_pycnocline_thickness
-
+ 
 
    SUBROUTINE zdf_osm_zmld_horizontal_gradients( zmld, zdtdx, zdtdy, zdsdx, zdsdy, dbdx_mle, dbdy_mle, zdbds_mle )
@@ -2578,5 +2622,5 @@
            ENDIF
            hmle(ji,jj) = MIN(hmle(ji,jj), ht_n(ji,jj))
-          IF(ln_osm_hmle_limit) hmle(ji,jj) = MIN(hmle(ji,jj), MAX(rn_osm_hmle_limit,1.2*hbl(ji,jj)) )
+           IF(ln_osm_hmle_limit) hmle(ji,jj) = MIN(hmle(ji,jj), MAX(rn_osm_hmle_limit,1.2*hbl(ji,jj)) )
         END DO
       END DO
