Index: /NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/TRA/tramle.F90
===================================================================
--- /NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/TRA/tramle.F90	(revision 12216)
+++ /NEMO/branches/2019/dev_r11078_OSMOSIS_IMMERSE_Nurser/src/OCE/TRA/tramle.F90	(revision 12217)
@@ -21,4 +21,6 @@
    USE lbclnk         ! lateral boundary condition / mpp link
 
+   USE zdfosm, ONLY  : ln_osm_mle, hmle, dbdx_mle, dbdy_mle, mld_prof
+
    IMPLICIT NONE
    PRIVATE
@@ -56,197 +58,257 @@
 CONTAINS
 
-   SUBROUTINE tra_mle_trp( kt, kit000, pu, pv, pw, cdtype )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE tra_mle_trp  ***
-      !!
-      !! ** Purpose :   Add to the transport the Mixed Layer Eddy induced transport
-      !!
-      !! ** Method  :   The 3 components of the Mixed Layer Eddy (MLE) induced
-      !!              transport are computed as follows :
-      !!                zu_mle = dk[ zpsi_uw ]
-      !!                zv_mle = dk[ zpsi_vw ]
-      !!                zw_mle = - di[ zpsi_uw ] - dj[ zpsi_vw ]
-      !!                where zpsi is the MLE streamfunction at uw and vw points (see the doc)
-      !!              and added to the input velocity :
-      !!                p.n = p.n + z._mle
-      !!
-      !! ** Action  : - (pun,pvn,pwn) increased by the mle transport
-      !!                CAUTION, the transport is not updated at the last line/raw
-      !!                         this may be a problem for some advection schemes
-      !!
-      !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008
-      !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008
-      !!----------------------------------------------------------------------
-      INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
-      INTEGER                         , INTENT(in   ) ::   kit000     ! first time step index
-      CHARACTER(len=3)                , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu         ! in : 3 ocean transport components
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pv         ! out: same 3  transport components
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pw         !   increased by the MLE induced transport
-      !
-      INTEGER  ::   ji, jj, jk          ! dummy loop indices
-      INTEGER  ::   ii, ij, ik, ikmax   ! local integers
-      REAL(wp) ::   zcuw, zmuw, zc      ! local scalar
-      REAL(wp) ::   zcvw, zmvw          !   -      -
-      INTEGER , DIMENSION(jpi,jpj)     :: inml_mle
-      REAL(wp), DIMENSION(jpi,jpj)     :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw
-      !!----------------------------------------------------------------------
-      !
-      !                                      !==  MLD used for MLE  ==!
-      !                                                ! compute from the 10m density to deal with the diurnal cycle
-      inml_mle(:,:) = mbkt(:,:) + 1                    ! init. to number of ocean w-level (T-level + 1)
-      IF ( nla10 > 0 ) THEN                            ! avoid case where first level is thicker than 10m
-         DO jk = jpkm1, nlb10, -1                      ! from the bottom to nlb10 (10m)
-            DO jj = 1, jpj
-               DO ji = 1, jpi                          ! index of the w-level at the ML based
-                  IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle )   inml_mle(ji,jj) = jk      ! Mixed layer
-               END DO
-            END DO
-         END DO
-      ENDIF
-      ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 )                  ! max level of the computation
-      !
-      !
-      zmld(:,:) = 0._wp                      !==   Horizontal shape of the MLE  ==!
-      zbm (:,:) = 0._wp
-      zn2 (:,:) = 0._wp
-      DO jk = 1, ikmax                                 ! MLD and mean buoyancy and N2 over the mixed layer
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points
-               zmld(ji,jj) = zmld(ji,jj) + zc
-               zbm (ji,jj) = zbm (ji,jj) + zc * (rau0 - rhop(ji,jj,jk) ) * r1_rau0
-               zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp
-            END DO
-         END DO
-      END DO
-
-      SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts
-      CASE ( 0 )                                               != min of the 2 neighbour MLDs
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) )
-               zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) )
-            END DO
-         END DO
-      CASE ( 1 )                                               != average of the 2 neighbour MLDs
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp
-               zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp
-            END DO
-         END DO
-      CASE ( 2 )                                               != max of the 2 neighbour MLDs
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) )
-               zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) )
-            END DO
-         END DO
-      END SELECT
-      !                                                ! convert density into buoyancy
-      zbm(:,:) = + grav * zbm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) )
-      !
-      !
-      !                                      !==  Magnitude of the MLE stream function  ==!
-      !
-      !                 di[bm]  Ds
-      ! Psi = Ce  H^2 ---------------- e2u  mu(z)   where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 )
-      !                  e1u   Lf fu                      and the e2u for the "transport"
-      !                                                      (not *e3u as divided by e3u at the end)
-      !
-      IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2_e1u(ji,jj)                                            &
-                  &           * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )   &
-                  &           / (  MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) )   )
-                  !
-               zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj)  * e1_e2v(ji,jj)                                            &
-                  &           * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) )   &
-                  &           / (  MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) )   )
-            END DO
-         END DO
-         !
-      ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat)
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2_e1u(ji,jj)               &
-                  &                  * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )
-                  !
-               zpsim_v(ji,jj) = rc_f *   zhv(ji,jj)   * zhv(ji,jj)   * e1_e2v(ji,jj)               &
-                  &                  * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) )
-            END DO
-         END DO
-      ENDIF
-      !
-      IF( nn_conv == 1 ) THEN              ! No MLE in case of convection
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp )   zpsim_u(ji,jj) = 0._wp
-               IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp )   zpsim_v(ji,jj) = 0._wp
-            END DO
-         END DO
-      ENDIF
-      !
-      !                                      !==  structure function value at uw- and vw-points  ==!
-      DO jj = 1, jpjm1
-         DO ji = 1, fs_jpim1   ! vector opt.
-            zhu(ji,jj) = 1._wp / zhu(ji,jj)                   ! hu --> 1/hu
-            zhv(ji,jj) = 1._wp / zhv(ji,jj)
-         END DO
-      END DO
-      !
-      zpsi_uw(:,:,:) = 0._wp
-      zpsi_vw(:,:,:) = 0._wp
-      !
-      DO jk = 2, ikmax                                ! start from 2 : surface value = 0
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               zcuw = 1._wp - ( gdepw_n(ji+1,jj,jk) + gdepw_n(ji,jj,jk) ) * zhu(ji,jj)
-               zcvw = 1._wp - ( gdepw_n(ji,jj+1,jk) + gdepw_n(ji,jj,jk) ) * zhv(ji,jj)
-               zcuw = zcuw * zcuw
-               zcvw = zcvw * zcvw
-               zmuw = MAX(  0._wp , ( 1._wp - zcuw ) * ( 1._wp + r5_21 * zcuw )  )
-               zmvw = MAX(  0._wp , ( 1._wp - zcvw ) * ( 1._wp + r5_21 * zcvw )  )
-               !
-               zpsi_uw(ji,jj,jk) = zpsim_u(ji,jj) * zmuw * umask(ji,jj,jk)
-               zpsi_vw(ji,jj,jk) = zpsim_v(ji,jj) * zmvw * vmask(ji,jj,jk)
-            END DO
-         END DO
-      END DO
-      !
-      !                                      !==  transport increased by the MLE induced transport ==!
-      DO jk = 1, ikmax
-         DO jj = 1, jpjm1                          ! CAUTION pu,pv must be defined at row/column i=1 / j=1
-            DO ji = 1, fs_jpim1   ! vector opt.
-               pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) )
-               pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) )
-            END DO
-         END DO
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk)   &
+  SUBROUTINE tra_mle_trp( kt, kit000, pu, pv, pw, cdtype )
+    !!----------------------------------------------------------------------
+    !!                  ***  ROUTINE tra_mle_trp  ***
+    !!
+    !! ** Purpose :   Add to the transport the Mixed Layer Eddy induced transport
+    !!
+    !! ** Method  :   The 3 components of the Mixed Layer Eddy (MLE) induced
+    !!              transport are computed as follows :
+    !!                zu_mle = dk[ zpsi_uw ]
+    !!                zv_mle = dk[ zpsi_vw ]
+    !!                zw_mle = - di[ zpsi_uw ] - dj[ zpsi_vw ]
+    !!                where zpsi is the MLE streamfunction at uw and vw points (see the doc)
+    !!              and added to the input velocity :
+    !!                p.n = p.n + z._mle
+    !!
+    !! ** Action  : - (pun,pvn,pwn) increased by the mle transport
+    !!                CAUTION, the transport is not updated at the last line/raw
+    !!                         this may be a problem for some advection schemes
+    !!
+    !! References: Fox-Kemper et al., JPO, 38, 1145-1165, 2008
+    !!             Fox-Kemper and Ferrari, JPO, 38, 1166-1179, 2008
+    !!----------------------------------------------------------------------
+    INTEGER                         , INTENT(in   ) ::   kt         ! ocean time-step index
+    INTEGER                         , INTENT(in   ) ::   kit000     ! first time step index
+    CHARACTER(len=3)                , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
+    REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pu         ! in : 3 ocean transport components
+    REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pv         ! out: same 3  transport components
+    REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pw         !   increased by the MLE induced transport
+    !
+    INTEGER  ::   ji, jj, jk          ! dummy loop indices
+    INTEGER  ::   ii, ij, ik, ikmax   ! local integers
+    REAL(wp) ::   zcuw, zmuw, zc      ! local scalar
+    REAL(wp) ::   zcvw, zmvw          !   -      -
+    INTEGER , DIMENSION(jpi,jpj)     :: inml_mle
+    REAL(wp), DIMENSION(jpi,jpj)     :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_NH, zLf_MH
+    REAL(wp), DIMENSION(jpi,jpj,jpk) :: zpsi_uw, zpsi_vw
+    !!----------------------------------------------------------------------
+    !
+    !
+    IF(ln_osm_mle) THEN
+       ikmax = MIN( MAXVAL( mld_prof(:,:) ), jpkm1 )                  ! max level of the computation
+       !
+       !
+       SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts
+       CASE ( 0 )                                               != min of the 2 neighbour MLDs
+          DO jj = 1, jpjm1
+             DO ji = 1, fs_jpim1   ! vector opt.
+                zhu(ji,jj) = MIN( hmle(ji+1,jj), hmle(ji,jj) )
+                zhv(ji,jj) = MIN( hmle(ji,jj+1), hmle(ji,jj) )
+             END DO
+          END DO
+       CASE ( 1 )                                               != average of the 2 neighbour MLDs
+          DO jj = 1, jpjm1
+             DO ji = 1, fs_jpim1   ! vector opt.
+                zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) )
+                zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) )
+             END DO
+          END DO
+       CASE ( 2 )                                               != max of the 2 neighbour MLDs
+          DO jj = 1, jpjm1
+             DO ji = 1, fs_jpim1   ! vector opt.
+                zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) )
+                zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) )
+             END DO
+          END DO
+       END SELECT
+       IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation
+          DO jj = 1, jpjm1
+             DO ji = 1, fs_jpim1   ! vector opt.
+                zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2u(ji,jj)                                            &
+                     &           * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) )   &
+                     &           / (  MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) )   )
+                !
+                zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj)  * e1v(ji,jj)                                            &
+                     &           * dbdy_mle(ji,jj)  * MIN( 111.e3_wp , e2v(ji,jj) )   &
+                     &           / (  MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) )   )
+             END DO
+          END DO
+          !
+       ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat)
+          DO jj = 1, jpjm1
+             DO ji = 1, fs_jpim1   ! vector opt.
+                zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2u(ji,jj)               &
+                     &                  * dbdx_mle(ji,jj) * MIN( 111.e3_wp , e1u(ji,jj) )
+                !
+                zpsim_v(ji,jj) = rc_f *   zhv(ji,jj)   * zhv(ji,jj)   * e1v(ji,jj)               &
+                     &                  * dbdy_mle(ji,jj) * MIN( 111.e3_wp , e2v(ji,jj) )
+             END DO
+          END DO
+       ENDIF
+
+    ELSE !do not use osn_mle
+       !                                      !==  MLD used for MLE  ==!
+       !                                                ! convert density into buoyancy
+       inml_mle(:,:) = mbkt(:,:) + 1                    ! init. to number of ocean w-level (T-level + 1)
+       IF ( nla10 > 0 ) THEN                            ! avoid case where first level is thicker than 10m
+          DO jk = jpkm1, nlb10, -1                      ! from the bottom to nlb10 (10m)
+             DO jj = 1, jpj
+                DO ji = 1, jpi                          ! index of the w-level at the ML based
+                   IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle )   inml_mle(ji,jj) = jk      ! Mixed layer
+                END DO
+             END DO
+          END DO
+       ENDIF
+
+       !
+       !
+       zmld(:,:) = 0._wp                      !==   Horizontal shape of the MLE  ==!
+       zbm (:,:) = 0._wp
+       zn2 (:,:) = 0._wp
+       DO jk = 1, ikmax                                 ! MLD and mean buoyancy and N2 over the mixed layer
+          DO jj = 1, jpj
+             DO ji = 1, jpi
+                zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, inml_mle(ji,jj)-jk ) , 1  )  )    ! zc being 0 outside the ML t-points
+                zmld(ji,jj) = zmld(ji,jj) + zc
+                zbm (ji,jj) = zbm (ji,jj) + zc * (rau0 - rhop(ji,jj,jk) ) * r1_rau0
+                zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp
+             END DO
+          END DO
+       END DO
+
+       SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts
+       CASE ( 0 )                                               != min of the 2 neighbour MLDs
+          DO jj = 1, jpjm1
+             DO ji = 1, fs_jpim1   ! vector opt.
+                zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) )
+                zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) )
+             END DO
+          END DO
+       CASE ( 1 )                                               != average of the 2 neighbour MLDs
+          DO jj = 1, jpjm1
+             DO ji = 1, fs_jpim1   ! vector opt.
+                zhu(ji,jj) = ( zmld(ji+1,jj) + zmld(ji,jj) ) * 0.5_wp
+                zhv(ji,jj) = ( zmld(ji,jj+1) + zmld(ji,jj) ) * 0.5_wp
+             END DO
+          END DO
+       CASE ( 2 )                                               != max of the 2 neighbour MLDs
+          DO jj = 1, jpjm1
+             DO ji = 1, fs_jpim1   ! vector opt.
+                zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) )
+                zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) )
+             END DO
+          END DO
+       END SELECT
+       !                                                ! convert density into buoyancy
+       ikmax = MIN( MAXVAL( inml_mle(:,:) ), jpkm1 )                  ! max level of the computation
+       zbm(:,:) = + grav * zbm(:,:) / MAX( e3t_n(:,:,1), zmld(:,:) )
+       !
+       !
+       !                                      !==  Magnitude of the MLE stream function  ==!
+       !
+       !                 di[bm]  Ds
+       ! Psi = Ce  H^2 ---------------- e2u  mu(z)   where fu Lf = MAX( fu*rn_fl , (Db H)^1/2 )
+       !                  e1u   Lf fu                      and the e2u for the "transport"
+       !                                                      (not *e3u as divided by e3u at the end)
+       !
+       IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation
+          DO jj = 1, jpjm1
+             DO ji = 1, fs_jpim1   ! vector opt.
+                zpsim_u(ji,jj) = rn_ce * zhu(ji,jj) * zhu(ji,jj)  * e2_e1u(ji,jj)                                            &
+                     &           * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )   &
+                     &           / (  MAX( rn_lf * rfu(ji,jj) , SQRT( rb_c * zhu(ji,jj) ) )   )
+                !
+                zpsim_v(ji,jj) = rn_ce * zhv(ji,jj) * zhv(ji,jj)  * e1_e2v(ji,jj)                                            &
+                     &           * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) )   &
+                     &           / (  MAX( rn_lf * rfv(ji,jj) , SQRT( rb_c * zhv(ji,jj) ) )   )
+             END DO
+          END DO
+          !
+       ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat)
+          DO jj = 1, jpjm1
+             DO ji = 1, fs_jpim1   ! vector opt.
+                zpsim_u(ji,jj) = rc_f *   zhu(ji,jj)   * zhu(ji,jj)   * e2_e1u(ji,jj)               &
+                     &                  * ( zbm(ji+1,jj) - zbm(ji,jj) ) * MIN( 111.e3_wp , e1u(ji,jj) )
+                !
+                zpsim_v(ji,jj) = rc_f *   zhv(ji,jj)   * zhv(ji,jj)   * e1_e2v(ji,jj)               &
+                     &                  * ( zbm(ji,jj+1) - zbm(ji,jj) ) * MIN( 111.e3_wp , e2v(ji,jj) )
+             END DO
+          END DO
+       ENDIF
+       !
+       IF( nn_conv == 1 ) THEN              ! No MLE in case of convection
+          DO jj = 1, jpjm1
+             DO ji = 1, fs_jpim1   ! vector opt.
+                IF( MIN( zn2(ji,jj) , zn2(ji+1,jj) ) < 0._wp )   zpsim_u(ji,jj) = 0._wp
+                IF( MIN( zn2(ji,jj) , zn2(ji,jj+1) ) < 0._wp )   zpsim_v(ji,jj) = 0._wp
+             END DO
+          END DO
+       ENDIF
+       !
+    ENDIF  ! end of lm_osm_mle loop
+    !                                      !==  structure function value at uw- and vw-points  ==!
+    DO jj = 1, jpjm1
+       DO ji = 1, fs_jpim1   ! vector opt.
+          zhu(ji,jj) = 1._wp / zhu(ji,jj)                   ! hu --> 1/hu
+          zhv(ji,jj) = 1._wp / zhv(ji,jj)
+       END DO
+    END DO
+    !
+    zpsi_uw(:,:,:) = 0._wp
+    zpsi_vw(:,:,:) = 0._wp
+    !
+    DO jk = 2, ikmax                                ! start from 2 : surface value = 0
+       DO jj = 1, jpjm1
+          DO ji = 1, fs_jpim1   ! vector opt.
+             zcuw = 1._wp - ( gdepw_n(ji+1,jj,jk) + gdepw_n(ji,jj,jk) ) * zhu(ji,jj)
+             zcvw = 1._wp - ( gdepw_n(ji,jj+1,jk) + gdepw_n(ji,jj,jk) ) * zhv(ji,jj)
+             zcuw = zcuw * zcuw
+             zcvw = zcvw * zcvw
+             zmuw = MAX(  0._wp , ( 1._wp - zcuw ) * ( 1._wp + r5_21 * zcuw )  )
+             zmvw = MAX(  0._wp , ( 1._wp - zcvw ) * ( 1._wp + r5_21 * zcvw )  )
+             !
+             zpsi_uw(ji,jj,jk) = zpsim_u(ji,jj) * zmuw * umask(ji,jj,jk)
+             zpsi_vw(ji,jj,jk) = zpsim_v(ji,jj) * zmvw * vmask(ji,jj,jk)
+          END DO
+       END DO
+    END DO
+    !
+    !                                      !==  transport increased by the MLE induced transport ==!
+    DO jk = 1, ikmax
+       DO jj = 1, jpjm1                          ! CAUTION pu,pv must be defined at row/column i=1 / j=1
+          DO ji = 1, fs_jpim1   ! vector opt.
+             pu(ji,jj,jk) = pu(ji,jj,jk) + ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji,jj,jk+1) )
+             pv(ji,jj,jk) = pv(ji,jj,jk) + ( zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj,jk+1) )
+          END DO
+       END DO
+       DO jj = 2, jpjm1
+          DO ji = fs_2, fs_jpim1   ! vector opt.
+             pw(ji,jj,jk) = pw(ji,jj,jk) - ( zpsi_uw(ji,jj,jk) - zpsi_uw(ji-1,jj,jk)   &
                   &                          + zpsi_vw(ji,jj,jk) - zpsi_vw(ji,jj-1,jk) )
-            END DO
-         END DO
-      END DO
-
-      IF( cdtype == 'TRA') THEN              !==  outputs  ==!
-         !
-         zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:)      ! Lf = N H / f
-         CALL iom_put( "Lf_NHpf" , zLf_NH  )    ! Lf = N H / f
-         !
-         ! divide by cross distance to give streamfunction with dimensions m^2/s
-         DO jk = 1, ikmax+1
-            zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk) * r1_e2u(:,:)
-            zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk) * r1_e1v(:,:)
-         END DO
-         CALL iom_put( "psiu_mle", zpsi_uw )    ! i-mle streamfunction
-         CALL iom_put( "psiv_mle", zpsi_vw )    ! j-mle streamfunction
-      ENDIF
-      !
-   END SUBROUTINE tra_mle_trp
+          END DO
+       END DO
+    END DO
+
+    IF( cdtype == 'TRA') THEN              !==  outputs  ==!
+       !
+       IF (ln_osm_mle) THEN
+          zLf_NH(:,:) = SQRT( rb_c * hmle(:,:) ) * r1_ft(:,:)      ! Lf = N H / f
+       ELSE
+          zLf_NH(:,:) = SQRT( rb_c * zmld(:,:) ) * r1_ft(:,:)      ! Lf = N H / f
+       END IF
+       CALL iom_put( "Lf_NHpf" , zLf_NH  )    ! Lf = N H / f
+       !
+       ! divide by cross distance to give streamfunction with dimensions m^2/s
+       DO jk = 1, ikmax+1
+          zpsi_uw(:,:,jk) = zpsi_uw(:,:,jk) * r1_e2u(:,:)
+          zpsi_vw(:,:,jk) = zpsi_vw(:,:,jk) * r1_e1v(:,:)
+       END DO
+       CALL iom_put( "psiu_mle", zpsi_uw )    ! i-mle streamfunction
+       CALL iom_put( "psiv_mle", zpsi_vw )    ! j-mle streamfunction
+    ENDIF
+    !
+  END SUBROUTINE tra_mle_trp
 
 
@@ -290,6 +352,12 @@
          WRITE(numout,*) '         Density difference used to define ML for FK              rn_rho_c_mle  = ', rn_rho_c_mle
       ENDIF
-      !
-      IF(lwp) THEN
+
+
+      IF( ln_osm_mle .AND. ln_mle ) THEN
+         WRITE(numout,*) 'WARNING: You are running with both OSM-FK and default FK'
+         ! CALL ctl_stop('STOP in zdf_osm_init: Cannot run with both OSM-FK and default FK')
+     END IF
+
+     IF(lwp) THEN
          WRITE(numout,*)
          IF( ln_mle ) THEN
