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 14529)
+++ NEMO/branches/2021/dev_r14122_HPC-08_Mueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90	(revision 14550)
@@ -988,5 +988,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)
@@ -2154,5 +2154,5 @@
       REAL(wp), DIMENSION(A2D(0)) ::   zzeta_pyc                               !
       REAL(wp)                    ::   zomega, zvw_max                         !
-      REAL(wp)                    ::   zuw_bse,zvw_bse                         ! Momentum, heat, and salinity fluxes
+      REAL(wp), DIMENSION(A2D(0)) ::   zuw_bse,zvw_bse                         ! Momentum, heat, and salinity fluxes
       REAL(wp), DIMENSION(A2D(0)) ::   zwth_ent,zws_ent                        ! at the top of the pycnocline
       REAL(wp), DIMENSION(A2D(0)) ::   zsc_wth_pyc, zsc_ws_pyc                 ! Scales for pycnocline transport term
@@ -2351,17 +2351,17 @@
             IF ( k_ddh(ji,jj) == 0 ) THEN
                ! Place holding code. Parametrization needs checking for these conditions.
-               zomega = ( 0.15_wp * pwstrl(ji,jj)**3 + pwstrc(ji,jj)**3 + 4.75_wp * ( pshear(ji,jj)* phbl(ji,jj) )**pthird )**pthird
-               zuw_bse = -0.0035_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdu_ml(ji,jj)
-               zvw_bse = -0.0075_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdv_ml(ji,jj)
+               zomega         = ( 0.15_wp * pwstrl(ji,jj)**3 + pwstrc(ji,jj)**3 + 4.75_wp * ( pshear(ji,jj) * phbl(ji,jj) ) )**pthird
+               zuw_bse(ji,jj) = -0.0035_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdu_ml(ji,jj)
+               zvw_bse(ji,jj) = -0.0075_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdv_ml(ji,jj)
             ELSE
-               zomega = ( 0.15_wp * pwstrl(ji,jj)**3 + pwstrc(ji,jj)**3 + 4.75_wp * ( pshear(ji,jj)* phbl(ji,jj) )**pthird )**pthird
-               zuw_bse = -0.0035_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdu_ml(ji,jj)
-               zvw_bse = -0.0075_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdv_ml(ji,jj)
+               zomega         = ( 0.15_wp * pwstrl(ji,jj)**3 + pwstrc(ji,jj)**3 + 4.75_wp * ( pshear(ji,jj) * phbl(ji,jj) ) )**pthird
+               zuw_bse(ji,jj) = -0.0035_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdu_ml(ji,jj)
+               zvw_bse(ji,jj) = -0.0075_wp * zomega * ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdv_ml(ji,jj)
             ENDIF
-            zb_cubic(ji,jj) = pdh(ji,jj) / phbl(ji,jj) * puw0(ji,jj) - ( 2.0 + pdh(ji,jj) / phml(ji,jj) ) * zuw_bse
-            za_cubic(ji,jj) = zuw_bse - zb_cubic(ji,jj)
+            zb_cubic(ji,jj) = pdh(ji,jj) / phbl(ji,jj) * puw0(ji,jj) - ( 2.0 + pdh(ji,jj) / phml(ji,jj) ) * zuw_bse(ji,jj)
+            za_cubic(ji,jj) = zuw_bse(ji,jj) - zb_cubic(ji,jj)
             zvw_max = 0.7_wp * ff_t(ji,jj) * ( pustke(ji,jj) * dstokes(ji,jj) + 0.7_wp * pustar(ji,jj) * phml(ji,jj) )
-            zd_cubic(ji,jj) = zvw_max * pdh(ji,jj) / phml(ji,jj) - ( 2.0_wp + pdh(ji,jj) / phml(ji,jj) ) * zvw_bse
-            zc_cubic(ji,jj) = zvw_bse - zd_cubic(ji,jj)
+            zd_cubic(ji,jj) = zvw_max * pdh(ji,jj) / phml(ji,jj) - ( 2.0_wp + pdh(ji,jj) / phml(ji,jj) ) * zvw_bse(ji,jj)
+            zc_cubic(ji,jj) = zvw_bse(ji,jj) - zd_cubic(ji,jj)
          END IF
       END_2D
@@ -2369,8 +2369,8 @@
          IF ( ldpyc(ji,jj) .AND. ( jk >= kmld(ji,jj) ) .AND. ( jk <= kbld(ji,jj) ) ) THEN
             zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj)
-            ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045_wp * ztau_sc_u(ji,jj)**2 *                                      &
+            ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045_wp * ( ztau_sc_u(ji,jj)**2 ) * zuw_bse(ji,jj) *                 &
                &                                ( za_cubic(ji,jj) * zznd_pyc**2 + zb_cubic(ji,jj) * zznd_pyc**3 ) *   &
                &                                ( 0.75_wp + 0.25_wp * zznd_pyc )**2 * pdbdz_pyc(ji,jj,jk)
-            ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.045_wp * ztau_sc_u(ji,jj)**2 *                                      &
+            ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.045_wp * ( ztau_sc_u(ji,jj)**2 ) * zvw_bse(ji,jj) *                 &
                &                                ( zc_cubic(ji,jj) * zznd_pyc**2 + zd_cubic(ji,jj) * zznd_pyc**3 ) *   &
                &                                ( 0.75_wp + 0.25_wp * zznd_pyc )**2 * pdbdz_pyc(ji,jj,jk)
@@ -2483,5 +2483,5 @@
       DO_3D( 0, 0, 0, 0, 2, jkm_bld )
          IF ( ( .NOT. ldconv(ji,jj) ) .AND. ( jk <= kbld(ji,jj) ) ) THEN
-            znd = ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / phbl(ji,jj)   ! ALMG to think about
+            znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / phbl(ji,jj)   ! ALMG to think about
             IF ( znd >= 0.0_wp ) THEN
                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0_wp - EXP( -10.0_wp * znd**2 ) )
