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 14315)
+++ NEMO/branches/2020/dev_r14122_ENHANCE-14_smueller_OSMOSIS_streamlining/src/OCE/ZDF/zdfosm.F90	(revision 14316)
@@ -277,7 +277,4 @@
       INTEGER, DIMENSION(jpi, jpj) :: j_ddh ! Type of shear layer
 
-      REAL(wp) :: ztgrad,zsgrad,zbgrad ! Temporary variables used to calculate pycnocline gradients
-      REAL(wp) :: zugrad,zvgrad        ! temporary variables for calculating pycnocline shear
-
       REAL(wp), DIMENSION(jpi,jpj) :: zhbl  ! bl depth - grid
       REAL(wp), DIMENSION(jpi,jpj) :: zhml  ! ml depth - grid
@@ -299,21 +296,9 @@
       REAL(wp), DIMENSION(jpi,jpj) :: zdt_ml,zds_ml,zdu_ml,zdv_ml,zdb_ml ! difference between mixed layer average and parameter at base of blayer
 !      REAL(wp), DIMENSION(jpi,jpj) :: zwth_ent,zws_ent ! heat and salinity fluxes at the top of the pycnocline
-      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) ::   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) ::   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
       REAL(wp), DIMENSION(jpi, jpj) :: zshear, zri_i ! Shear production and interfacial richardon number.
 
-      REAL(wp) :: zl_c,zl_l,zl_eps  ! Used to calculate turbulence length scale.
-
-      REAL(wp) :: za_cubic, zb_cubic, zc_cubic, zd_cubic ! coefficients in cubic polynomial specifying diffusivity in pycnocline.
-      REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_1,zsc_ws_1 ! Temporary scales used to calculate scalar non-gradient terms.
-      REAL(wp), DIMENSION(jpi,jpj) :: zsc_wth_pyc, zsc_ws_pyc ! Scales for pycnocline transport term/
-      REAL(wp), DIMENSION(jpi,jpj) :: zsc_uw_1,zsc_uw_2,zsc_vw_1,zsc_vw_2 ! Temporary scales for non-gradient momentum flux terms.
       REAL(wp), DIMENSION(jpi,jpj) :: zhbl_t ! holds boundary layer depth updated by full timestep
 
@@ -334,7 +319,4 @@
       REAL(wp), DIMENSION(jpi,jpj,jpk) :: zdiffut ! t-diffusivity
       REAL(wp), DIMENSION(jpi,jpj) :: zalpha_pyc
-      REAL(wp), DIMENSION(jpi,jpj) :: ztau_sc_u ! dissipation timescale at baes of WML.
-      REAL(wp) :: zdelta_pyc, zwt_pyc_sc_1, zws_pyc_sc_1, zzeta_pyc
-      REAL(wp) :: zbuoy_pyc_sc, zomega, zvw_max
       INTEGER :: ibld_ext=0                          ! does not have to be zero for modified scheme
       REAL(wp) :: zgamma_b_nd, zgamma_b, zdhoh, ztau
@@ -344,6 +326,4 @@
       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
@@ -354,5 +334,5 @@
       ibld(:,:)   = 0     ; imld(:,:)  = 0
       zustar(:,:) = 0.0_wp
-      zwstrl(:,:) = 0._wp ; zvstr(:,:) = 0._wp ; zwstrc(:,:)    = 0._wp ; zuw0(:,:)      = 0._wp
+      zwstrl(:,:) = 0._wp ; zvstr(:,:) = 0._wp ; zwstrc(:,:)    = 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
@@ -372,5 +352,4 @@
       zdt_ml(:,:)  = 0._wp ; zds_ml(:,:)  = 0._wp ; zdu_ml(:,:)   = 0._wp ; zdv_ml(:,:)  = 0._wp
       zdb_ml(:,:)  = 0._wp
-      zwth_ent = 0._wp ; zws_ent = 0._wp
       !
       zdbdz_pyc(:,:,:) = 0.0_wp
@@ -386,8 +365,5 @@
       zwb_fk_b(:,:) = 0._wp   ! must be initialised even with ln_osm_mle=F as used in zdf_osm_calculate_dhdt
 
-      ! Flux-Gradient arrays.
-      zsc_wth_1(:,:)  = 0._wp ; zsc_ws_1(:,:)   = 0._wp ; zsc_uw_1(:,:)   = 0._wp
-      zsc_uw_2(:,:)   = 0._wp ; zsc_vw_1(:,:)   = 0._wp ; zsc_vw_2(:,:)   = 0._wp
-      zhbl_t(:,:)     = 0._wp ; zdhdt(:,:)      = 0._wp
+      zhbl_t(:,:)     = 0._wp
 
       zdiffut(:,:,:) = 0._wp ; zviscos(:,:,:) = 0._wp ; ghamt(:,:,:) = 0._wp
@@ -706,451 +682,11 @@
        CALL zdf_osm_diffusivity_viscosity( zdiffut, zviscos )
 
-       !
-       ! calculate non-gradient components of the flux-gradient relationships
-       !
-! Stokes term in scalar flux, flux-gradient relationship
-       WHERE ( lconv )
-          zsc_wth_1 = zwstrl**3 * zwth0 / ( zvstr**3 + 0.5 * zwstrc**3 + epsln)
-          !
-          zsc_ws_1 = zwstrl**3 * zws0 / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
-       ELSEWHERE
-          zsc_wth_1 = 2.0 * zwthav
-          !
-          zsc_ws_1 = 2.0 * zwsav
-       ENDWHERE
-
-
-       DO_2D( 0, 0, 0, 0 )
-         IF ( lconv(ji,jj) ) THEN
-           DO jk = 2, imld(ji,jj)
-              zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
-              ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_wth_1(ji,jj)
-              !
-              ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.35 * EXP ( -zznd_d ) * ( 1.0 - EXP ( -2.0 * zznd_d ) ) *  zsc_ws_1(ji,jj)
-           END DO ! end jk loop
-         ELSE     ! else for if (lconv)
- ! Stable conditions
-            DO jk = 2, ibld(ji,jj)
-               zznd_d=gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
-               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) &
-                    &          *                 ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_wth_1(ji,jj)
-               !
-               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5 * EXP ( -0.9 * zznd_d ) &
-                    &          *                 ( 1.0 - EXP ( -4.0 * zznd_d ) ) *  zsc_ws_1(ji,jj)
-            END DO
-         ENDIF               ! endif for check on lconv
-
-       END_2D
-
-! Stokes term in flux-gradient relationship (note in zsc_uw_n don't use zvstr since term needs to go to zero as zwstrl goes to zero)
-       WHERE ( lconv )
-          zsc_uw_1 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MAX( ( 1.0 - 1.0 * 6.5 * zla**(8.0/3.0) ), 0.2 )
-          zsc_uw_2 = ( zwstrl**3 + 0.5 * zwstrc**3 )**pthird * zustke / MIN( zla**(8.0/3.0) + epsln, 0.12 )
-          zsc_vw_1 = ff_t * zhml * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / ( ( zvstr**3 + 0.5 * zwstrc**3 )**(2.0/3.0) + epsln )
-       ELSEWHERE
-          zsc_uw_1 = zustar**2
-          zsc_vw_1 = ff_t * zhbl * zustke**3 * MIN( zla**(8.0/3.0), 0.12 ) / (zvstr**2 + epsln)
-       ENDWHERE
-       IF(ln_dia_osm) THEN
-          IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu )
-          IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv )
-       END IF
-       DO_2D( 0, 0, 0, 0 )
-          IF ( lconv(ji,jj) ) THEN
-             DO jk = 2, imld(ji,jj)
-                zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
-                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) +      ( -0.05 * EXP ( -0.4 * zznd_d )   * zsc_uw_1(ji,jj)   &
-                     &          +                        0.00125 * EXP (      - zznd_d )   * zsc_uw_2(ji,jj) ) &
-                     &          *                          ( 1.0 - EXP ( -2.0 * zznd_d ) )
-!
-                ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.65 *  0.15 * EXP (      - zznd_d )                       &
-                     &          *                          ( 1.0 - EXP ( -2.0 * zznd_d ) ) * zsc_vw_1(ji,jj)
-             END DO   ! end jk loop
-          ELSE
-! Stable conditions
-             DO jk = 2, ibld(ji,jj) ! corrected to ibld
-                zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
-                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75 *   1.3 * EXP ( -0.5 * zznd_d )                       &
-                     &                                   * ( 1.0 - EXP ( -4.0 * zznd_d ) ) * zsc_uw_1(ji,jj)
-                ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0._wp
-             END DO   ! end jk loop
-          ENDIF
-       END_2D
-
-! Buoyancy term in flux-gradient relationship [note : includes ROI ratio (X0.3) and pressure (X0.5)]
-
-       WHERE ( lconv )
-          zsc_wth_1 = zwbav * zwth0 * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
-          zsc_ws_1  = zwbav * zws0  * ( 1.0 + EXP ( 0.2 * zhol ) ) / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
-       ELSEWHERE
-          zsc_wth_1 = 0._wp
-          zsc_ws_1 = 0._wp
-       ENDWHERE
-
-       DO_2D( 0, 0, 0, 0 )
-          IF (lconv(ji,jj) ) THEN
-             DO jk = 2, imld(ji,jj)
-                zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj)
-                ! calculate turbulent length scale
-                zl_c = 0.9 * ( 1.0 - EXP ( - 7.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           &
-                     &     * ( 1.0 - EXP ( -15.0 * (     1.1 - zznd_ml          ) ) )
-                zl_l = 2.0 * ( 1.0 - EXP ( - 2.0 * ( zznd_ml - zznd_ml**3 / 3.0 ) ) )                                           &
-                     &     * ( 1.0 - EXP ( - 5.0 * (     1.0 - zznd_ml          ) ) ) * ( 1.0 + dstokes(ji,jj) / zhml (ji,jj) )
-                zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0 + EXP ( -3.0 * LOG10 ( - zhol(ji,jj) ) ) ) ** (3.0 / 2.0)
-                ! non-gradient buoyancy terms
-                ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * 0.5 * zsc_wth_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml )
-                ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * 0.5 *  zsc_ws_1(ji,jj) * zl_eps * zhml(ji,jj) / ( 0.15 + zznd_ml )
-             END DO
-
-             IF ( lpyc(ji,jj) ) THEN
-               ztau_sc_u(ji,jj) = zhml(ji,jj) / ( zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird
-               ztau_sc_u(ji,jj) = ztau_sc_u(ji,jj) * ( 1.4 -0.4 / ( 1.0 + EXP( -3.5 * LOG10( -zhol(ji,jj) ) ) )**1.5 )
-               zwth_ent =  -0.003 * ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird * ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zdt_ml(ji,jj)
-               zws_ent =  -0.003 * ( 0.15 * zvstr(ji,jj)**3 + zwstrc(ji,jj)**3 )**pthird * ( 1.0 - zdh(ji,jj) /zhbl(ji,jj) ) * zds_ml(ji,jj)
-! Cubic profile used for buoyancy term
-               za_cubic = 0.755 * ztau_sc_u(ji,jj)
-               zb_cubic = 0.25 * ztau_sc_u(ji,jj)
-               DO jk = 2, ibld(ji,jj)
-                 zznd_pyc = -( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj)
-                 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) - 0.045 * ( ( zwth_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 )
-
-                 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
-!
-               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)
-!
-               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(ji,jj,jk,Kmm) - 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
-            ENDIF ! End of pycnocline
-          ELSE ! lconv test - stable conditions
-             DO jk = 2, ibld(ji,jj)
-                ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj)
-                ghams(ji,jj,jk) = ghams(ji,jj,jk) +  zsc_ws_1(ji,jj)
-             END DO
-          ENDIF
-       END_2D
-
-       WHERE ( lconv )
-          zsc_uw_1 = -zwb0 * zustar**2 * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )
-          zsc_uw_2 =  zwb0 * zustke    * zhml / ( zvstr**3 + 0.5 * zwstrc**3 + epsln )**(2.0/3.0)
-          zsc_vw_1 = 0._wp
-       ELSEWHERE
-         zsc_uw_1 = 0._wp
-         zsc_vw_1 = 0._wp
-       ENDWHERE
-
-       DO_2D( 0, 0, 0, 0 )
-          IF ( lconv(ji,jj) ) THEN
-             DO jk = 2 , imld(ji,jj)
-                zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
-                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3 * 0.5 * ( zsc_uw_1(ji,jj) +   0.125 * EXP( -0.5 * zznd_d )     &
-                     &                                                            * (   1.0 - EXP( -0.5 * zznd_d ) )   &
-                     &                                          * zsc_uw_2(ji,jj)                                    )
-                ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj)
-             END DO  ! jk loop
-          ELSE
-          ! stable conditions
-             DO jk = 2, ibld(ji,jj)
-                ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj)
-                ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj)
-             END DO
-          ENDIF
-       END_2D
-
-       DO_2D( 0, 0, 0, 0 )
-        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 )**pthird
-            zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj)
-            zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj)
-          ELSE
-            zomega = ( 0.15 * zwstrl(ji,jj)**3 + zwstrc(ji,jj)**3 + 4.75 * ( zshear(ji,jj)* zhbl(ji,jj) )**pthird )**pthird
-            zuw_bse = -0.0035 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdu_ml(ji,jj)
-            zvw_bse = -0.0075 * zomega * ( 1.0 - zdh(ji,jj) / zhbl(ji,jj) ) * zdv_ml(ji,jj)
-          ENDIF
-          zd_cubic = zdh(ji,jj) / zhbl(ji,jj) * zuw0(ji,jj) - ( 2.0 + zdh(ji,jj) /zhml(ji,jj) ) * zuw_bse
-          zc_cubic = zuw_bse - zd_cubic
-! need ztau_sc_u to be available. Change to array.
-          DO jk = imld(ji,jj), ibld(ji,jj)
-             zznd_pyc = - ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj)
-             ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.045 * ztau_sc_u(ji,jj)**2 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk)
-          END DO
-          zvw_max = 0.7 * ff_t(ji,jj) * ( zustke(ji,jj) * dstokes(ji,jj) + 0.75 * zustar(ji,jj) * zhml(ji,jj) )
-          zd_cubic = zvw_max * zdh(ji,jj) / zhml(ji,jj) - ( 2.0 + zdh(ji,jj) /zhml(ji,jj) ) * zvw_bse
-          zc_cubic = zvw_bse - zd_cubic
-          DO jk = imld(ji,jj), ibld(ji,jj)
-            zznd_pyc = -( gdepw(ji,jj,jk,Kmm) -zhbl(ji,jj) ) / zdh(ji,jj)
-            ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.045 * ztau_sc_u(ji,jj)**2 * ( zc_cubic * zznd_pyc**2 + zd_cubic * zznd_pyc**3 ) * ( 0.75 + 0.25 * zznd_pyc )**2 * zdbdz_pyc(ji,jj,jk)
-          END DO
-        ENDIF  ! lpyc
-       END_2D
-
-       IF(ln_dia_osm) THEN
-          IF ( iom_use("ghamu_0") ) CALL iom_put( "ghamu_0", wmask*ghamu )
-          IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 )
-       END IF
-! Transport term in flux-gradient relationship [note : includes ROI ratio (X0.3) ]
-
-       DO_2D( 1, 0, 1, 0 )
-
-         IF ( lconv(ji,jj) ) THEN
-           zsc_wth_1(ji,jj) = zwth0(ji,jj) / ( 1.0 - 0.56 * EXP( zhol(ji,jj) ) )
-           zsc_ws_1(ji,jj) = zws0(ji,jj) / (1.0 - 0.56 *EXP( zhol(ji,jj) ) )
-           IF ( lpyc(ji,jj) ) THEN
-! Pycnocline scales
-              zsc_wth_pyc(ji,jj) = -0.2 * zwb0(ji,jj) * zdt_bl(ji,jj) / zdb_bl(ji,jj)
-              zsc_ws_pyc(ji,jj) = -0.2 * zwb0(ji,jj) * zds_bl(ji,jj) / zdb_bl(ji,jj)
-            ENDIF
-         ELSE
-           zsc_wth_1(ji,jj) = 2.0 * zwthav(ji,jj)
-           zsc_ws_1(ji,jj) = zws0(ji,jj)
-         ENDIF
-       END_2D
-
-       DO_2D( 0, 0, 0, 0 )
-         IF ( lconv(ji,jj) ) THEN
-            DO jk = 2, imld(ji,jj)
-               zznd_ml=gdepw(ji,jj,jk,Kmm) / zhml(ji,jj)
-               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * zsc_wth_1(ji,jj)                &
-                    &          * ( -2.0 + 2.75 * (       ( 1.0 + 0.6 * zznd_ml**4 )      &
-                    &                               - EXP(     - 6.0 * zznd_ml    ) ) )  &
-                    &          * ( 1.0 - EXP( - 15.0 * (         1.0 - zznd_ml    ) ) )
-               !
-               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * zsc_ws_1(ji,jj)  &
-                    &          * ( -2.0 + 2.75 * (       ( 1.0 + 0.6 * zznd_ml**4 )      &
-                    &                               - EXP(     - 6.0 * zznd_ml    ) ) )  &
-                    &          * ( 1.0 - EXP ( -15.0 * (         1.0 - zznd_ml    ) ) )
-            END DO
-!
-            IF ( lpyc(ji,jj) ) THEN
-! pycnocline
-              DO jk = imld(ji,jj), ibld(ji,jj)
-                zznd_pyc = - ( gdepw(ji,jj,jk,Kmm) - zhbl(ji,jj) ) / zdh(ji,jj)
-                ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 4.0 * zsc_wth_pyc(ji,jj) * ( 0.48 - EXP( -1.5 * ( zznd_pyc -0.3)**2 ) )
-                ghams(ji,jj,jk) = ghams(ji,jj,jk) + 4.0 * zsc_ws_pyc(ji,jj) * ( 0.48 - EXP( -1.5 * ( zznd_pyc -0.3)**2 ) )
-              END DO
-           ENDIF
-         ELSE
-            IF( zdhdt(ji,jj) > 0. ) THEN
-              DO jk = 2, ibld(ji,jj)
-                 zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
-                 znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
-                 ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + &
-              &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_wth_1(ji,jj)
-                 ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3 * ( -4.06 * EXP( -2.0 * zznd_d ) * (1.0 - EXP( -4.0 * zznd_d ) ) + &
-               &  7.5 * EXP ( -10.0 * ( 0.95 - znd )**2 ) * ( 1.0 - znd ) ) * zsc_ws_1(ji,jj)
-              END DO
-            ENDIF
-         ENDIF
-       END_2D
-
-       WHERE ( lconv )
-          zsc_uw_1 = zustar**2
-          zsc_vw_1 = ff_t * zustke * zhml
-       ELSEWHERE
-          zsc_uw_1 = zustar**2
-          zsc_uw_2 = (2.25 - 3.0 * ( 1.0 - EXP( -1.25 * 2.0 ) ) ) * ( 1.0 - EXP( -4.0 * 2.0 ) ) * zsc_uw_1
-          zsc_vw_1 = ff_t * zustke * zhbl
-          zsc_vw_2 = -0.11 * SIN( 3.14159 * ( 2.0 + 0.4 ) ) * EXP(-( 1.5 + 2.0 )**2 ) * zsc_vw_1
-       ENDWHERE
-
-       DO_2D( 0, 0, 0, 0 )
-          IF ( lconv(ji,jj) ) THEN
-            DO jk = 2, imld(ji,jj)
-               zznd_ml = gdepw(ji,jj,jk,Kmm) / zhml(ji,jj)
-               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
-               ghamu(ji,jj,jk) = ghamu(ji,jj,jk)&
-                    & + 0.3 * ( -2.0 + 2.5 * ( 1.0 + 0.1 * zznd_ml**4 ) - EXP ( -8.0 * zznd_ml ) ) * zsc_uw_1(ji,jj)
-               !
-               ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
-                    & + 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)
-               znd = gdepw(ji,jj,jk,Kmm) / zhbl(ji,jj)
-               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
-               IF ( zznd_d <= 2.0 ) THEN
-                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5 * 0.3 &
-                       &*  ( 2.25 - 3.0  * ( 1.0 - EXP( - 1.25 * zznd_d ) ) * ( 1.0 - EXP( -2.0 * zznd_d ) ) ) * zsc_uw_1(ji,jj)
-                  !
-               ELSE
-                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk)&
-                       & + 0.5 * 0.3 * ( 1.0 - EXP( -5.0 * ( 1.0 - znd ) ) ) * zsc_uw_2(ji,jj)
-                  !
-               ENDIF
-
-               ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
-                    & + 0.3 * 0.15 * SIN( 3.14159 * ( 0.65 * zznd_d ) ) * EXP( -0.25 * zznd_d**2 ) * zsc_vw_1(ji,jj)
-               ghamv(ji,jj,jk) = ghamv(ji,jj,jk)&
-                    & + 0.3 * 0.15 * EXP( -5.0 * ( 1.0 - znd ) ) * ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj)
-            END DO
-          ENDIF
-       END_2D
-
-       IF(ln_dia_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 )
-          IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 )
-          IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 )
-          IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 )
-          IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 )
-       END IF
-!
-! Make surface forced velocity non-gradient terms go to zero at the base of the mixed layer.
-
-
- ! Make surface forced velocity non-gradient terms go to zero at the base of the boundary layer.
-
-      DO_2D( 0, 0, 0, 0 )
-         IF ( .not. lconv(ji,jj) ) THEN
-            DO jk = 2, ibld(ji,jj)
-               znd = ( gdepw(ji,jj,jk,Kmm) - 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 ) )
-                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0 - EXP( -10.0 * znd**2 ) )
-               ELSE
-                  ghamu(ji,jj,jk) = 0._wp
-                  ghamv(ji,jj,jk) = 0._wp
-               ENDIF
-            END DO
-         ENDIF
-      END_2D
-
-      ! 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
-      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 ( 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 )
-          ghamt(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
-          ghams(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
-          ghamu(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
-          ghamv(ji,jj,ibld(ji,jj)+ibld_ext) = 0._wp
-       END_2D
-
-       IF(ln_dia_osm) THEN
-          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("zviscos") ) CALL iom_put( "zviscos", wmask*zviscos )
-       END IF
+      !
+      ! Calculate non-gradient components of the flux-gradient relationships
+      ! --------------------------------------------------------------------
+      CALL zdf_osm_fgr_terms( Kmm, ibld, imld, jp_ext, ibld_ext, lconv, lpyc, j_ddh, zhbl, zhml, zdh, zdhdt, zhol, zshear,   &
+         &                    zustar, zwstrl, zvstr, zwstrc, zuw0, zwth0, zws0, zwb0, zwthav, zwsav, zwbav, zustke, zla,     &
+         &                    zdt_bl, zds_bl, zdb_bl, zdu_bl, zdv_bl, zdt_ml, zds_ml, zdb_ml, zdu_ml, zdv_ml,                &
+         &                    zdtdz_bl_ext, zdsdz_bl_ext, zdbdz_bl_ext, zdbdz_pyc, zalpha_pyc, zdiffut, zviscos )
 
        !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
@@ -1329,7 +865,5 @@
          IF ( iom_use("zhol") ) CALL iom_put( "zhol", tmask(:,:,1)*zhol )               ! ML depth internal to zdf_osm routine
          IF ( iom_use("zwthav") ) CALL iom_put( "zwthav", tmask(:,:,1)*zwthav )         ! upward BL-avged turb temp flux
-         IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent )   ! upward turb temp entrainment flux
          IF ( iom_use("zwb_ent") ) CALL iom_put( "zwb_ent", tmask(:,:,1)*zwb_ent )      ! upward turb buoyancy entrainment flux
-         IF ( iom_use("zws_ent") ) CALL iom_put( "zws_ent", tmask(:,:,1)*zws_ent )      ! upward turb salinity entrainment flux
          IF ( iom_use("zt_ml") ) CALL iom_put( "zt_ml", tmask(:,:,1)*zt_ml )            ! average T in ML
 
@@ -1374,4 +908,5 @@
 !
       REAL(wp) :: zvel_sc_pyc, zvel_sc_ml, zstab_fac
+      REAL(wp) ::   za_cubic, zb_cubic, zc_cubic, zd_cubic   ! Coefficients in cubic polynomial specifying diffusivity in pycnocline
 
       REAL(wp), PARAMETER :: rn_dif_ml = 0.8, rn_vis_ml = 0.375
@@ -1883,4 +1418,8 @@
       END_2D
       !
+      IF ( ln_dia_pyc_scl ) THEN   ! Output of pycnocline gradient profiles
+         IF ( iom_use("zdbdz_pyc") ) CALL iom_put( "zdbdz_pyc", wmask(:,:,:) * pdbdz(:,:,:) )
+      END IF
+      !
       IF( ln_timing ) CALL timing_stop('zdf_osm_pscp')
       !
@@ -2542,4 +2081,545 @@
    END SUBROUTINE zdf_osm_vertical_average
 
+   SUBROUTINE zdf_osm_fgr_terms( Kmm, kbld, kmld, kp_ext, kbld_ext, ldconv, ldpyc, k_ddh, phbl, phml, pdh, pdhdt, phol, pshear,   &
+      &                          pustar, pwstrl, pvstr, pwstrc, puw0, pwth0, pws0, pwb0, pwthav, pwsav, pwbav, pustke, pla,       &
+      &                          pdt_bl, pds_bl, pdb_bl, pdu_bl, pdv_bl, pdt_ml, pds_ml, pdb_ml, pdu_ml, pdv_ml,                  &
+      &                          pdtdz_bl_ext, pdsdz_bl_ext, pdbdz_bl_ext, pdbdz_pyc, palpha_pyc, pdiffut, pviscos )
+      !!---------------------------------------------------------------------
+      !!                 ***  ROUTINE zdf_osm_fgr_terms ***
+      !!
+      !! ** Purpose : Compute non-gradient terms in flux-gradient relationship
+      !!
+      !! ** Method  :
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER,                    INTENT(in   ) ::   Kmm            ! Time-level index
+      INTEGER,  DIMENSION(:,:),   INTENT(in   ) ::   kbld           ! BL base layer
+      INTEGER,  DIMENSION(:,:),   INTENT(in   ) ::   kmld           ! ML base layer
+      INTEGER,  DIMENSION(:,:),   INTENT(in   ) ::   kp_ext         ! Offset for external level
+      INTEGER,                    INTENT(in   ) ::   kbld_ext       ! Offset for external level
+      LOGICAL,  DIMENSION(:,:),   INTENT(in   ) ::   ldconv         ! BL stability flags
+      LOGICAL,  DIMENSION(:,:),   INTENT(in   ) ::   ldpyc          ! Pycnocline flags
+      INTEGER,  DIMENSION(:,:),   INTENT(in   ) ::   k_ddh          ! Type of shear layer
+      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   ) ::   phol           ! Stability parameter for boundary layer
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pshear         ! Shear production
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pustar         ! Friction velocity
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pwstrl         ! Langmuir velocity scale
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pvstr          ! Velocity scale (approaches zustar for large Langmuir number)
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pwstrc         ! Convective velocity scale
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   puw0           ! Surface u-momentum flux
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pwth0          ! Surface heat flux
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pws0           ! Surface freshwater flux
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pwb0           ! Surface buoyancy flux
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pwthav         ! BL average heat flux
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pwsav          ! BL average freshwater flux
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pwbav          ! BL average buoyancy flux
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pustke         ! Surface Stokes drift
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pla            ! Langmuir number
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdt_bl         ! Temperature diff. between BL average and basal value
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pds_bl         ! Salinity diff. between BL average and basal value
+      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. (u) between BL average and basal value
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdt_ml         ! Temperature diff. between mixed-layer average and basal value
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pds_ml         ! Salinity diff. between mixed-layer 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
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdtdz_bl_ext   ! External temperature gradients
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdsdz_bl_ext   ! External salinity gradients
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   pdbdz_bl_ext   ! External buoyancy gradients
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdbdz_pyc      ! Pycnocline buoyancy gradients
+      REAL(wp), DIMENSION(:,:),   INTENT(in   ) ::   palpha_pyc     !
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pdiffut        ! t-diffusivity
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pviscos        ! Viscosity
+      !
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   z3ddz_pyc_1, z3ddz_pyc_2   ! Pycnocline gradient/shear profiles
+      !
+      INTEGER                     ::   ji, jj, jk, jkm_bld, jkf_mld, jkm_mld   ! Loop indices
+      INTEGER                     ::   istat                                   ! Memory allocation status
+      REAL(wp)                    ::   zznd_d, zznd_ml, zznd_pyc, znd          ! Temporary non-dimensional depths
+      REAL(wp), DIMENSION(A2D(0)) ::   zsc_wth_1,zsc_ws_1                      ! Temporary scales
+      REAL(wp), DIMENSION(A2D(0)) ::   zsc_uw_1, zsc_uw_2                      ! Temporary scales
+      REAL(wp), DIMENSION(A2D(0)) ::   zsc_vw_1, zsc_vw_2                      ! Temporary scales
+      REAL(wp), DIMENSION(A2D(0)) ::   ztau_sc_u                               ! Dissipation timescale at base of WML
+      REAL(wp)                    ::   zbuoy_pyc_sc, zdelta_pyc                !
+      REAL(wp)                    ::   zl_c,zl_l,zl_eps                        ! Used to calculate turbulence length scale
+      REAL(wp), DIMENSION(A2D(0)) ::   za_cubic, zb_cubic                      ! Coefficients in cubic polynomial specifying
+      REAL(wp), DIMENSION(A2D(0)) ::   zc_cubic, zd_cubic                      ! diffusivity in pycnocline
+      REAL(wp), DIMENSION(A2D(0)) ::   zwt_pyc_sc_1, zws_pyc_sc_1              !
+      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)) ::   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
+      REAL(wp)                    ::   ztmp                                    !
+      REAL(wp)                    ::   ztgrad, zsgrad, zbgrad                  ! Variables used to calculate pycnocline gradients
+      REAL(wp)                    ::   zugrad, zvgrad                          ! Variables for calculating pycnocline shear
+      REAL(wp)                    ::   zdtdz_pyc                               ! Parametrized gradient of temperature in pycnocline
+      REAL(wp)                    ::   zdsdz_pyc                               ! Parametrised gradient of salinity in pycnocline
+      REAL(wp)                    ::   zdudz_pyc                               ! u-shear across the pycnocline
+      REAL(wp)                    ::   zdvdz_pyc                               ! v-shear across the pycnocline
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing ) CALL timing_start('zdf_osm_ft')
+      !
+      ! Auxiliary indices
+      ! -----------------
+      jkm_bld = 0
+      jkf_mld = jpk
+      jkm_mld = 0
+      DO_2D( 0, 0, 0, 0 )
+         IF ( kbld(ji,jj) > jkm_bld ) jkm_bld = kbld(ji,jj)
+         IF ( kbld(ji,jj) < jkf_mld ) jkf_mld = kbld(ji,jj)
+         IF ( kmld(ji,jj) > jkm_mld ) jkm_mld = kmld(ji,jj)
+      END_2D
+      !
+      ! Stokes term in scalar flux, flux-gradient relationship
+      ! ------------------------------------------------------
+      WHERE ( ldconv(A2D(0)) )
+         zsc_wth_1(:,:) = pwstrl(A2D(0))**3 * pwth0(A2D(0)) / ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln )
+         zsc_ws_1(:,:)  = pwstrl(A2D(0))**3 * pws0(A2D(0))  / ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln )
+      ELSEWHERE
+         zsc_wth_1(:,:) = 2.0_wp * pwthav(A2D(0))
+         zsc_ws_1(:,:)  = 2.0_wp * pwsav(A2D(0))
+      ENDWHERE
+      DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) )
+         IF ( ldconv(ji,jj) ) THEN
+            IF ( jk <= kmld(ji,jj) ) THEN
+               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
+               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.35_wp * EXP( -1.0_wp * zznd_d ) *   &
+                  &                                ( 1.0_wp - EXP( -2.0_wp * zznd_d ) ) * zsc_wth_1(ji,jj)
+               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.35_wp * EXP( -1.0_wp * zznd_d ) *   &
+                  &                                ( 1.0_wp - EXP( -2.0_wp * zznd_d ) ) * zsc_ws_1(ji,jj)
+            END IF
+         ELSE   ! Stable conditions
+            IF ( jk <= kbld(ji,jj) ) THEN
+               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
+               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 1.5_wp * EXP( -0.9_wp * zznd_d ) *   &
+                  &                                ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) * zsc_wth_1(ji,jj)
+               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 1.5_wp * EXP( -0.9_wp * zznd_d ) *   &
+                  &                                ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) * zsc_ws_1(ji,jj)
+            END IF
+         END IF   ! Check on ldconv
+      END_3D
+      !
+      IF ( ln_dia_osm ) THEN
+         IF ( iom_use("ghamu_00") ) CALL iom_put( "ghamu_00", wmask*ghamu )
+         IF ( iom_use("ghamv_00") ) CALL iom_put( "ghamv_00", wmask*ghamv )
+      END IF
+      !
+      ! Stokes term in flux-gradient relationship (note in zsc_uw_n don't use
+      ! zvstr since term needs to go to zero as zwstrl goes to zero)
+      ! ---------------------------------------------------------------------
+      WHERE ( ldconv(A2D(0)) )
+         zsc_uw_1(:,:) = ( pwstrl(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 )**pthird * pustke(A2D(0)) /   &
+            &                                  MAX( ( 1.0_wp - 1.0_wp * 6.5_wp * pla(A2D(0))**( 8.0_wp / 3.0_wp ) ), 0.2_wp )
+         zsc_uw_2(:,:) = ( pwstrl(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 )**pthird * pustke(A2D(0)) /   &
+            &                                  MIN( pla(A2D(0))**( 8.0_wp / 3.0_wp ) + epsln, 0.12_wp )
+         zsc_vw_1(:,:) = ff_t(A2D(0)) * phml(A2D(0)) * pustke(A2D(0))**3 * MIN( pla(A2D(0))**( 8.0_wp / 3.0_wp ), 0.12_wp ) /   &
+            &            ( ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 )**( 2.0_wp / 3.0_wp ) + epsln )
+      ELSEWHERE
+         zsc_uw_1(:,:) = pustar(A2D(0))**2
+         zsc_vw_1(:,:) = ff_t(A2D(0)) * phbl(A2D(0)) * pustke(A2D(0))**3 * MIN( pla(A2D(0))**( 8.0_wp / 3.0_wp ), 0.12_wp ) /   &
+            &            ( pvstr(A2D(0))**2 + epsln )
+      ENDWHERE
+      DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) )
+         IF ( ldconv(ji,jj) ) THEN
+            IF ( jk <= kmld(ji,jj) ) THEN
+               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
+               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + ( -0.05_wp   * EXP( -0.4_wp * zznd_d ) * zsc_uw_1(ji,jj) +     &
+                  &                                  0.00125_wp * EXP( -1.0_wp * zznd_d ) * zsc_uw_2(ji,jj) ) *   &
+                  &                                ( 1.0_wp - EXP( -2.0_wp * zznd_d ) )
+               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) - 0.65_wp *  0.15_wp * EXP( -1.0_wp * zznd_d ) *                 &
+                  &                                ( 1.0_wp - EXP( -2.0_wp * zznd_d ) ) * zsc_vw_1(ji,jj)
+            END IF
+         ELSE   ! Stable conditions
+            IF ( jk <= kbld(ji,jj) ) THEN   ! Corrected to ibld
+               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
+               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) - 0.75_wp * 1.3_wp * EXP( -0.5_wp * zznd_d ) *             &
+                  &                                ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) * zsc_uw_1(ji,jj)
+            END IF
+         END IF
+      END_3D
+      !
+      ! Buoyancy term in flux-gradient relationship [note : includes ROI ratio
+      ! (X0.3) and pressure (X0.5)]
+      ! ----------------------------------------------------------------------
+      WHERE ( ldconv(A2D(0)) )
+         zsc_wth_1(:,:) = pwbav(A2D(0)) * pwth0(A2D(0)) * ( 1.0_wp + EXP( 0.2_wp * phol(A2D(0)) ) ) /   &
+            &             ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln )
+         zsc_ws_1(:,:)  = pwbav(A2D(0)) * pws0(A2D(0))  * ( 1.0_wp + EXP( 0.2_wp * phol(A2D(0)) ) ) /   &
+            &             ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln )
+      ELSEWHERE
+         zsc_wth_1(:,:) = 0.0_wp
+         zsc_ws_1(:,:)  = 0.0_wp
+      ENDWHERE
+      DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) )
+         IF ( ldconv(ji,jj) ) THEN
+            IF ( jk <= kmld(ji,jj) ) THEN
+               zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj)
+               ! Calculate turbulent length scale
+               zl_c   = 0.9_wp * ( 1.0_wp - EXP( -7.0_wp * ( zznd_ml - zznd_ml**3 / 3.0_wp ) ) ) *                         &
+                  &     ( 1.0_wp - EXP( -15.0_wp * ( 1.1_wp - zznd_ml ) ) )
+               zl_l   = 2.0_wp * ( 1.0_wp - EXP( -2.0_wp * ( zznd_ml - zznd_ml**3 / 3.0_wp ) ) ) *                         &
+                  &     ( 1.0_wp - EXP( -5.0_wp  * ( 1.0_wp - zznd_ml ) ) ) * ( 1.0_wp + dstokes(ji,jj) / phml (ji,jj) )
+               zl_eps = zl_l + ( zl_c - zl_l ) / ( 1.0_wp + EXP( -3.0_wp * LOG10( -1.0_wp * phol(ji,jj) ) ) )**( 3.0_wp / 2.0_wp )
+               ! Non-gradient buoyancy terms
+               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3_wp * 0.5_wp * zsc_wth_1(ji,jj) * zl_eps * phml(ji,jj) / ( 0.15_wp + zznd_ml )
+               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3_wp * 0.5_wp *  zsc_ws_1(ji,jj) * zl_eps * phml(ji,jj) / ( 0.15_wp + zznd_ml )
+            END IF
+         ELSE   ! Stable conditions
+            IF ( jk <= kbld(ji,jj) ) THEN
+               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + zsc_wth_1(ji,jj)
+               ghams(ji,jj,jk) = ghams(ji,jj,jk) +  zsc_ws_1(ji,jj)
+            END IF
+         END IF
+      END_3D
+      DO_2D( 0, 0, 0, 0 )
+         IF ( ldconv(ji,jj) .AND. ldpyc(ji,jj) ) THEN
+            ztau_sc_u(ji,jj)    = phml(ji,jj) / ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird *                             &
+               &                ( 1.4_wp - 0.4_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * phol(ji,jj) ) ) )**1.5_wp )
+            zwth_ent(ji,jj)     = -0.003_wp * ( 0.15_wp * pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird *   &
+               &                ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pdt_ml(ji,jj)
+            zws_ent(ji,jj)      = -0.003_wp * ( 0.15_wp * pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird *   &
+               &                ( 1.0_wp - pdh(ji,jj) / phbl(ji,jj) ) * pds_ml(ji,jj)
+            zbuoy_pyc_sc        = palpha_pyc(ji,jj) * pdb_ml(ji,jj) / pdh(ji,jj) + pdbdz_bl_ext(ji,jj)
+            zdelta_pyc          = ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird /   &
+               &                  SQRT( MAX( zbuoy_pyc_sc, ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**p2third / pdh(ji,jj)**2 ) )
+            zwt_pyc_sc_1(ji,jj) = 0.325_wp * ( palpha_pyc(ji,jj) * pdt_ml(ji,jj) / pdh(ji,jj) + pdtdz_bl_ext(ji,jj) ) *   &
+               &                  zdelta_pyc**2 / pdh(ji,jj)
+            zws_pyc_sc_1(ji,jj) = 0.325_wp * ( palpha_pyc(ji,jj) * pds_ml(ji,jj) / pdh(ji,jj) + pdsdz_bl_ext(ji,jj) ) *   &
+               &                  zdelta_pyc**2 / pdh(ji,jj)
+            zzeta_pyc(ji,jj)    = 0.15_wp - 0.175_wp / ( 1.0_wp + EXP( -3.5_wp * LOG10( -1.0_wp * phol(ji,jj) ) ) )
+         END IF
+      END_2D
+      DO_3D( 0, 0, 0, 0, 2, jkm_bld )
+         IF ( ldconv(ji,jj) .AND. ldpyc(ji,jj) .AND. ( jk <= kbld(ji,jj) ) ) THEN
+            zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj)
+            ghamt(ji,jj,jk) = ghamt(ji,jj,jk) -                                                                                &
+               &              0.045_wp * ( ( zwth_ent(ji,jj) * pdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) *                 &
+               &                         MAX( ( 1.75_wp * zznd_pyc -0.15_wp * zznd_pyc**2 - 0.2_wp * zznd_pyc**3 ), 0.0_wp )
+            ghams(ji,jj,jk) = ghams(ji,jj,jk) -                                                                                &
+               &              0.045_wp * ( ( zws_ent(ji,jj)  * pdbdz_pyc(ji,jj,jk) ) * ztau_sc_u(ji,jj)**2 ) *                 &
+               &                         MAX( ( 1.75_wp * zznd_pyc -0.15_wp * zznd_pyc**2 - 0.2_wp * zznd_pyc**3 ), 0.0_wp )
+            ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.05_wp  * zwt_pyc_sc_1(ji,jj) *                              &
+               &                                EXP( -0.25_wp * ( zznd_pyc / zzeta_pyc(ji,jj) )**2 ) *        &
+               &                                pdh(ji,jj) / ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird
+            ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.05_wp  * zws_pyc_sc_1(ji,jj) *                              &
+               &                                EXP( -0.25_wp * ( zznd_pyc / zzeta_pyc(ji,jj) )**2 ) *        &
+               &                                pdh(ji,jj) / ( pvstr(ji,jj)**3 + pwstrc(ji,jj)**3 )**pthird
+         END IF   ! End of pycnocline
+      END_3D
+      !
+      IF(ln_dia_osm) THEN
+         IF ( iom_use("zwth_ent") ) CALL iom_put( "zwth_ent", tmask(:,:,1)*zwth_ent )   ! Upward turb. temperature entrainment flux
+         IF ( iom_use("zws_ent")  ) CALL iom_put( "zws_ent",  tmask(:,:,1)*zws_ent  )   ! Upward turb. salinity    entrainment flux
+      END IF
+      !
+      zsc_vw_1(:,:) = 0.0_wp
+      WHERE ( ldconv(A2D(0)) )
+         zsc_uw_1(:,:) = -1.0_wp * pwb0(A2D(0)) * pustar(A2D(0))**2 * phml(A2D(0)) /   &
+            &            ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln )
+         zsc_uw_2(:,:) =           pwb0(A2D(0)) * pustke(A2D(0))    * phml(A2D(0)) /   &
+            &            ( pvstr(A2D(0))**3 + 0.5_wp * pwstrc(A2D(0))**3 + epsln )**( 2.0_wp / 3.0_wp )
+      ELSEWHERE
+         zsc_uw_1(:,:) = 0.0_wp
+      ENDWHERE
+      DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) )
+         IF ( ldconv(ji,jj) ) THEN
+            IF ( jk <= kmld(ji,jj) ) THEN
+               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
+               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.3_wp * 0.5_wp *   &
+                  &                                ( zsc_uw_1(ji,jj) + 0.125_wp * EXP( -0.5_wp * zznd_d ) *       &
+                  &                                  (   1.0_wp - EXP( -0.5_wp * zznd_d ) ) * zsc_uw_2(ji,jj) )
+               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj)
+            END IF
+         ELSE   ! Stable conditions
+            IF ( jk <= kbld(ji,jj) ) THEN
+               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + zsc_uw_1(ji,jj)
+               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + zsc_vw_1(ji,jj)
+            END IF
+         ENDIF
+      END_3D
+      !
+      DO_2D( 0, 0, 0, 0 )
+         IF ( ldpyc(ji,jj) ) THEN
+            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)
+            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)
+            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)
+            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)
+         END IF
+      END_2D
+      DO_3D( 0, 0, 0, 0, jkf_mld, jkm_bld )   ! Need ztau_sc_u to be available. Change to array.
+         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 *                                      &
+               &                                ( 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 *                                      &
+               &                                ( 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)
+         END IF   ! ldpyc
+      END_3D
+      !
+      IF(ln_dia_osm) THEN
+         IF ( iom_use("ghamu_0") )    CALL iom_put( "ghamu_0",    wmask*ghamu           )
+         IF ( iom_use("zsc_uw_1_0") ) CALL iom_put( "zsc_uw_1_0", tmask(:,:,1)*zsc_uw_1 )
+      END IF
+      !
+      ! Transport term in flux-gradient relationship [note : includes ROI ratio
+      ! (X0.3) ]
+      ! -----------------------------------------------------------------------
+      WHERE ( ldconv(A2D(0)) )
+         zsc_wth_1(:,:) = pwth0(A2D(0)) / ( 1.0_wp - 0.56_wp * EXP( phol(A2D(0)) ) )
+         zsc_ws_1(:,:)  = pws0(A2D(0))  / ( 1.0_wp - 0.56_wp * EXP( phol(A2D(0)) ) )
+         WHERE ( ldpyc(A2D(0)) )   ! Pycnocline scales
+            zsc_wth_pyc(:,:) = -0.2_wp * pwb0(A2D(0)) * pdt_bl(A2D(0)) / pdb_bl(A2D(0))
+            zsc_ws_pyc(:,:)  = -0.2_wp * pwb0(A2D(0)) * pds_bl(A2D(0)) / pdb_bl(A2D(0))
+         END WHERE
+      ELSEWHERE
+         zsc_wth_1(:,:) = 2.0 * pwthav(A2D(0))
+         zsc_ws_1(:,:)  =       pws0(A2D(0))
+      END WHERE
+      DO_3D( 0, 0, 0, 0, 1, MAX( jkm_mld, jkm_bld ) )
+         IF ( ldconv(ji,jj) ) THEN
+            IF ( ( jk > 1 ) .AND. ( jk <= kmld(ji,jj) ) ) THEN
+               zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj)
+               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3_wp * zsc_wth_1(ji,jj) * ( -2.0_wp + 2.75_wp * ( ( 1.0_wp + 0.6_wp * zznd_ml**4 ) - EXP( -6.0_wp * zznd_ml ) ) ) *   &
+                  &      ( 1.0_wp - EXP( -15.0_wp * ( 1.0_wp - zznd_ml ) ) )
+               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3_wp * zsc_ws_1(ji,jj)  * ( -2.0_wp + 2.75_wp * ( ( 1.0_wp + 0.6_wp * zznd_ml**4 ) - EXP( -6.0_wp * zznd_ml ) ) ) *   &
+                  &      ( 1.0_wp - EXP( -15.0_wp * ( 1.0_wp - zznd_ml ) ) )
+            END IF
+            IF ( ldpyc(ji,jj) .AND. ( jk >= kmld(ji,jj) ) .AND. ( jk <= kbld(ji,jj) ) ) THEN   ! Pycnocline
+               zznd_pyc = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phbl(ji,jj) ) / pdh(ji,jj)
+               ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 4.0_wp * zsc_wth_pyc(ji,jj) * ( 0.48_wp - EXP( -1.5_wp * ( zznd_pyc - 0.3_wp )**2 ) )
+               ghams(ji,jj,jk) = ghams(ji,jj,jk) + 4.0_wp * zsc_ws_pyc(ji,jj)  * ( 0.48_wp - EXP( -1.5_wp * ( zznd_pyc - 0.3_wp )**2 ) )
+            END IF
+         ELSE
+            IF( pdhdt(ji,jj) > 0. ) THEN
+               IF ( ( jk > 1 ) .AND. ( jk <= kbld(ji,jj) ) ) THEN
+                  zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
+                  znd    = gdepw(ji,jj,jk,Kmm) / phbl(ji,jj)
+                  ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + 0.3_wp * ( -4.06_wp * EXP( -2.0_wp * zznd_d ) * ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) +   &
+                                      7.5_wp * EXP ( -10.0_wp * ( 0.95_wp - znd )**2 ) * ( 1.0_wp - znd ) ) * zsc_wth_1(ji,jj)
+                  ghams(ji,jj,jk) = ghams(ji,jj,jk) + 0.3_wp * ( -4.06_wp * EXP( -2.0_wp * zznd_d ) * ( 1.0_wp - EXP( -4.0_wp * zznd_d ) ) +   &
+                                      7.5_wp * EXP ( -10.0_wp * ( 0.95_wp - znd )**2 ) * ( 1.0_wp - znd ) ) * zsc_ws_1(ji,jj)
+               END IF
+            ENDIF
+         ENDIF
+      END_3D
+      !
+      WHERE ( ldconv(A2D(0)) )
+         zsc_uw_1(:,:) = pustar(A2D(0))**2
+         zsc_vw_1(:,:) = ff_t(A2D(0)) * pustke(A2D(0)) * phml(A2D(0))
+      ELSEWHERE
+         zsc_uw_1(:,:) = pustar(A2D(0))**2
+         zsc_uw_2(:,:) = ( 2.25_wp - 3.0_wp * ( 1.0_wp - EXP( -1.25_wp * 2.0_wp ) ) ) * ( 1.0_wp - EXP( -4.0_wp * 2.0_wp ) ) *   &
+            &            zsc_uw_1(:,:)
+         zsc_vw_1(:,:) = ff_t * pustke(A2D(0)) * phbl(A2D(0))
+         zsc_vw_2(:,:) = -0.11_wp * SIN( 3.14159_wp * ( 2.0_wp + 0.4_wp ) ) * EXP( -1.0_wp * ( 1.5_wp + 2.0_wp )**2 ) *   &
+            &            zsc_vw_1(:,:)
+      ENDWHERE
+      DO_3D( 0, 0, 0, 0, 2, MAX( jkm_mld, jkm_bld ) )
+         IF ( ldconv(ji,jj) ) THEN
+            IF ( jk <= kmld(ji,jj) ) THEN
+               zznd_ml = gdepw(ji,jj,jk,Kmm) / phml(ji,jj)
+               zznd_d  = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
+               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) +   &
+                  &              0.3_wp * ( -2.0_wp + 2.5_wp * ( 1.0_wp + 0.1_wp * zznd_ml**4 ) - EXP( -8.0_wp * zznd_ml ) ) *   &
+                  &              zsc_uw_1(ji,jj)
+               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) +   &
+                  &              0.3_wp * 0.1_wp * ( EXP( -1.0_wp * zznd_d ) + EXP( -5.0_wp * ( 1.0_wp - zznd_ml ) ) ) *   &
+                  &              zsc_vw_1(ji,jj)
+            END IF
+         ELSE
+            IF ( jk <= kbld(ji,jj) ) THEN
+               znd    = gdepw(ji,jj,jk,Kmm) / phbl(ji,jj)
+               zznd_d = gdepw(ji,jj,jk,Kmm) / dstokes(ji,jj)
+               IF ( zznd_d <= 2.0 ) THEN
+                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5_wp * 0.3_wp *                                              &
+                     &                                ( 2.25_wp - 3.0_wp * ( 1.0_wp - EXP( -1.25_wp * zznd_d ) ) *   &
+                     &                                  ( 1.0_wp - EXP( -2.0_wp * zznd_d ) ) ) * zsc_uw_1(ji,jj)
+               ELSE
+                  ghamu(ji,jj,jk) = ghamu(ji,jj,jk) + 0.5_wp * 0.3_wp *   &
+                     &                                ( 1.0_wp - EXP( -5.0_wp * ( 1.0_wp - znd ) ) ) * zsc_uw_2(ji,jj)
+               ENDIF
+               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0.3_wp * 0.15_wp * SIN( 3.14159_wp * ( 0.65_wp * zznd_d ) ) *   &
+                  &                                EXP( -0.25_wp * zznd_d**2 ) * zsc_vw_1(ji,jj)
+               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + 0.3_wp * 0.15_wp * EXP( -5.0 * ( 1.0 - znd ) ) * ( 1.0 - EXP( -20.0 * ( 1.0 - znd ) ) ) * zsc_vw_2(ji,jj)
+            END IF
+         END IF
+      END_3D
+      !
+      IF(ln_dia_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    )
+         IF ( iom_use("zsc_uw_1_f") ) CALL iom_put( "zsc_uw_1_f", tmask(:,:,1)*zsc_uw_1 )
+         IF ( iom_use("zsc_vw_1_f") ) CALL iom_put( "zsc_vw_1_f", tmask(:,:,1)*zsc_vw_1 )
+         IF ( iom_use("zsc_uw_2_f") ) CALL iom_put( "zsc_uw_2_f", tmask(:,:,1)*zsc_uw_2 )
+         IF ( iom_use("zsc_vw_2_f") ) CALL iom_put( "zsc_vw_2_f", tmask(:,:,1)*zsc_vw_2 )
+      END IF
+      !
+      ! Make surface forced velocity non-gradient terms go to zero at the base
+      ! of the mixed layer.
+      !
+      ! Make surface forced velocity non-gradient terms go to zero at the base
+      ! of the boundary layer.
+      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
+            IF ( znd >= 0.0_wp ) THEN
+               ghamu(ji,jj,jk) = ghamu(ji,jj,jk) * ( 1.0_wp - EXP( -10.0_wp * znd**2 ) )
+               ghamv(ji,jj,jk) = ghamv(ji,jj,jk) * ( 1.0_wp - EXP( -10.0_wp * znd**2 ) )
+            ELSE
+               ghamu(ji,jj,jk) = 0.0_wp
+               ghamv(ji,jj,jk) = 0.0_wp
+            ENDIF
+         END IF
+      END_3D
+      !
+      ! Pynocline contributions
+      !
+      IF ( ln_dia_pyc_scl .OR. ln_dia_pyc_shr ) THEN   ! Allocate arrays for output of pycnocline gradient/shear profiles
+         ALLOCATE( z3ddz_pyc_1(jpi,jpj,jpk), z3ddz_pyc_2(jpi,jpj,jpk), STAT=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_3D( 0, 0, 0, 0, 2, jkm_bld )
+         IF ( ldconv (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 ( kbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN
+               ! Pycnocline profile only defined when depth steady of increasing.
+               IF ( pdhdt(ji,jj) > 0.0_wp ) THEN   ! Depth increasing, or steady.
+                  IF ( pdb_bl(ji,jj) > 0.0_wp ) THEN
+                     IF ( phol(ji,jj) >= 0.5_wp ) THEN   ! Very stable - 'thick' pycnocline
+                        ztmp = 1.0_wp / MAX( phbl(ji,jj), epsln )
+                        ztgrad = pdt_bl(ji,jj) * ztmp
+                        zsgrad = pds_bl(ji,jj) * ztmp
+                        zbgrad = pdb_bl(ji,jj) * ztmp
+                        IF ( jk <= kbld(ji,jj) ) THEN
+                           znd = gdepw(ji,jj,jk,Kmm) * ztmp
+                           zdtdz_pyc =  ztgrad * EXP( -15.0_wp * ( znd - 0.9_wp )**2 )
+                           zdsdz_pyc =  zsgrad * EXP( -15.0_wp * ( znd - 0.9_wp )**2 )
+                           ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + pdiffut(ji,jj,jk) * zdtdz_pyc
+                           ghams(ji,jj,jk) = ghams(ji,jj,jk) + pdiffut(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 IF
+                     ELSE   ! Slightly stable - 'thin' pycnoline - needed when stable layer begins to form.
+                        ztmp = 1.0_wp / MAX( pdh(ji,jj), epsln )
+                        ztgrad = pdt_bl(ji,jj) * ztmp
+                        zsgrad = pds_bl(ji,jj) * ztmp
+                        zbgrad = pdb_bl(ji,jj) * ztmp
+                        IF ( jk <= kbld(ji,jj) ) THEN
+                           znd = -1.0_wp * ( gdepw(ji,jj,jk,Kmm) - phml(ji,jj) ) * ztmp
+                           zdtdz_pyc =  ztgrad * EXP( -1.75_wp * ( znd + 0.75_wp )**2 )
+                           zdsdz_pyc =  zsgrad * EXP( -1.75_wp * ( znd + 0.75_wp )**2 )
+                           ghamt(ji,jj,jk) = ghamt(ji,jj,jk) + pdiffut(ji,jj,jk) * zdtdz_pyc
+                           ghams(ji,jj,jk) = ghams(ji,jj,jk) + pdiffut(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 IF
+                     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_3D
+      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(:,:,:) )
+      END IF
+      DO_3D( 0, 0, 0, 0, 2, jkm_bld )
+         IF ( .NOT. ldconv (ji,jj) ) THEN
+            IF ( kbld(ji,jj) + kp_ext(ji,jj) < mbkt(ji,jj) ) THEN
+               zugrad = 3.25_wp * pdu_bl(ji,jj) / phbl(ji,jj)
+               zvgrad = 2.75_wp * pdv_bl(ji,jj) / phbl(ji,jj)
+               IF ( jk <= kbld(ji,jj) ) THEN
+                  znd = gdepw(ji,jj,jk,Kmm) / phbl(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) + pviscos(ji,jj,jk) * zdudz_pyc
+                  ghamv(ji,jj,jk) = ghamv(ji,jj,jk) + pviscos(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 IF
+            END IF
+         END IF
+      END_3D
+      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 ( 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 )
+         ghamt(ji,jj,kbld(ji,jj)+kbld_ext) = 0.0_wp
+         ghams(ji,jj,kbld(ji,jj)+kbld_ext) = 0.0_wp
+         ghamu(ji,jj,kbld(ji,jj)+kbld_ext) = 0.0_wp
+         ghamv(ji,jj,kbld(ji,jj)+kbld_ext) = 0.0_wp
+      END_2D
+      !
+      IF(ln_dia_osm) THEN
+         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("zviscos") ) CALL iom_put( "zviscos", wmask*pviscos )
+      END IF
+      !
+      IF( ln_timing ) CALL timing_stop('zdf_osm_ft')
+      !
+   END SUBROUTINE zdf_osm_fgr_terms
+
    SUBROUTINE zdf_osm_init( Kmm )
      !!----------------------------------------------------------------------
