Index: NEMO/branches/2020/dev_r13787_doc_latex_recovery/src/OCE/TRA/tramle.F90
===================================================================
--- NEMO/branches/2020/dev_r13787_doc_latex_recovery/src/OCE/TRA/tramle.F90	(revision 13823)
+++ NEMO/branches/2020/dev_r13787_doc_latex_recovery/src/OCE/TRA/tramle.F90	(revision 14066)
@@ -20,4 +20,8 @@
    USE lib_mpp        ! MPP library
    USE lbclnk         ! lateral boundary condition / mpp link
+
+   ! where OSMOSIS_OBL is used with integrated FK
+   USE zdf_oce, ONLY : ln_zdfosm
+   USE zdfosm, ONLY  : ln_osm_mle, hmle, dbdx_mle, dbdy_mle, mld_prof
 
    IMPLICIT NONE
@@ -79,11 +83,11 @@
       !!             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
-      INTEGER                         , INTENT(in   ) ::   Kmm        ! ocean time level 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                     , INTENT(in   ) ::   kt         ! ocean time-step index
+      INTEGER                     , INTENT(in   ) ::   kit000     ! first time step index
+      INTEGER                     , INTENT(in   ) ::   Kmm        ! ocean time level index
+      CHARACTER(len=3)            , INTENT(in   ) ::   cdtype     ! =TRA or TRC (tracer indicator)
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) ::   pu         ! in : 3 ocean transport components
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) ::   pv         ! out: same 3  transport components
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk), INTENT(inout) ::   pw         !   increased by the MLE induced transport
       !
       INTEGER  ::   ji, jj, jk          ! dummy loop indices
@@ -91,95 +95,147 @@
       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_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 )        ! from the bottom to nlb10 (10m)
-            IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle )   inml_mle(ji,jj) = jk      ! Mixed layer
+      INTEGER , DIMENSION(A2D(nn_hls))     :: inml_mle
+      REAL(wp), DIMENSION(A2D(nn_hls))     :: zpsim_u, zpsim_v, zmld, zbm, zhu, zhv, zn2, zLf_MH
+      REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: zpsi_uw, zpsi_vw
+      ! TEMP: [tiling] These changes not necessary if using XIOS (subdomain support)
+      REAL(wp), DIMENSION(:,:),   ALLOCATABLE, SAVE :: zLf_NH
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, SAVE :: zpsiu_mle, zpsiv_mle
+      !!----------------------------------------------------------------------
+      !
+      !
+      IF(ln_osm_mle.and.ln_zdfosm) 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_2D( 1, 0, 1, 0 )
+               zhu(ji,jj) = MIN( hmle(ji+1,jj), hmle(ji,jj) )
+               zhv(ji,jj) = MIN( hmle(ji,jj+1), hmle(ji,jj) )
+            END_2D
+         CASE ( 1 )                                               != average of the 2 neighbour MLDs
+            DO_2D( 1, 0, 1, 0 )
+               zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) )
+               zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) )
+            END_2D
+         CASE ( 2 )                                               != max of the 2 neighbour MLDs
+            DO_2D( 1, 0, 1, 0 )
+               zhu(ji,jj) = MAX( hmle(ji+1,jj), hmle(ji,jj) )
+               zhv(ji,jj) = MAX( hmle(ji,jj+1), hmle(ji,jj) )
+            END_2D
+         END SELECT
+         IF( nn_mle == 0 ) THEN           ! Fox-Kemper et al. 2010 formulation
+            DO_2D( 1, 0, 1, 0 )
+               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_2D
+            !
+         ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat)
+            DO_2D( 1, 0, 1, 0 )
+               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_2D
+         ENDIF
+
+      ELSE !do not use osn_mle
+         !                                      !==  MLD used for MLE  ==!
+         !                                                ! compute from the 10m density to deal with the diurnal cycle
+         DO_2D( 1, 1, 1, 1 )
+            inml_mle(ji,jj) = mbkt(ji,jj) + 1                    ! init. to number of ocean w-level (T-level + 1)
+         END_2D
+         IF ( nla10 > 0 ) THEN                            ! avoid case where first level is thicker than 10m
+           DO_3DS( 1, 1, 1, 1, jpkm1, nlb10, -1 )        ! from the bottom to nlb10 (10m)
+              IF( rhop(ji,jj,jk) > rhop(ji,jj,nla10) + rn_rho_c_mle )   inml_mle(ji,jj) = jk      ! Mixed layer
+           END_3D
+         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_3D( 1, 1, 1, 1, 1, ikmax )                    ! MLD and mean buoyancy and N2 over the mixed layer
+            zc = e3t(ji,jj,jk,Kmm) * 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 * (rho0 - rhop(ji,jj,jk) ) * r1_rho0
+            zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp
          END_3D
-      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_3D( 1, 1, 1, 1, 1, ikmax )                    ! MLD and mean buoyancy and N2 over the mixed layer
-         zc = e3t(ji,jj,jk,Kmm) * 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 * (rho0 - rhop(ji,jj,jk) ) * r1_rho0
-         zn2 (ji,jj) = zn2 (ji,jj) + zc * (rn2(ji,jj,jk)+rn2(ji,jj,jk+1))*0.5_wp
-      END_3D
-
-      SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts
-      CASE ( 0 )                                               != min of the 2 neighbour MLDs
-         DO_2D( 1, 0, 1, 0 )
-            zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) )
-            zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) )
+   
+         SELECT CASE( nn_mld_uv )                         ! MLD at u- & v-pts
+         CASE ( 0 )                                               != min of the 2 neighbour MLDs
+            DO_2D( 1, 0, 1, 0 )
+               zhu(ji,jj) = MIN( zmld(ji+1,jj), zmld(ji,jj) )
+               zhv(ji,jj) = MIN( zmld(ji,jj+1), zmld(ji,jj) )
+            END_2D
+         CASE ( 1 )                                               != average of the 2 neighbour MLDs
+            DO_2D( 1, 0, 1, 0 )
+               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_2D
+         CASE ( 2 )                                               != max of the 2 neighbour MLDs
+            DO_2D( 1, 0, 1, 0 )
+               zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) )
+               zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) )
+            END_2D
+         END SELECT
+         !                                                ! convert density into buoyancy
+         DO_2D( 1, 1, 1, 1 )
+            zbm(ji,jj) = + grav * zbm(ji,jj) / MAX( e3t(ji,jj,1,Kmm), zmld(ji,jj) )
          END_2D
-      CASE ( 1 )                                               != average of the 2 neighbour MLDs
-         DO_2D( 1, 0, 1, 0 )
-            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_2D
-      CASE ( 2 )                                               != max of the 2 neighbour MLDs
-         DO_2D( 1, 0, 1, 0 )
-            zhu(ji,jj) = MAX( zmld(ji+1,jj), zmld(ji,jj) )
-            zhv(ji,jj) = MAX( zmld(ji,jj+1), zmld(ji,jj) )
-         END_2D
-      END SELECT
-      !                                                ! convert density into buoyancy
-      zbm(:,:) = + grav * zbm(:,:) / MAX( e3t(:,:,1,Kmm), 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_2D( 1, 0, 1, 0 )
-            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) ) )   )
+         !
+         !
+         !                                      !==  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_2D( 1, 0, 1, 0 )
+               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_2D
-         !
-      ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat)
-         DO_2D( 1, 0, 1, 0 )
-            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) = 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_2D
+            !
+         ELSEIF( nn_mle == 1 ) THEN       ! New formulation (Lf = 5km fo/ff with fo=Coriolis parameter at latitude rn_lat)
+            DO_2D( 1, 0, 1, 0 )
+               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_2D
-      ENDIF
-      !
-      IF( nn_conv == 1 ) THEN              ! No MLE in case of convection
-         DO_2D( 1, 0, 1, 0 )
-            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_2D
-      ENDIF
-      !
-      !                                      !==  structure function value at uw- and vw-points  ==!
-      DO_2D( 1, 0, 1, 0 )
-         zhu(ji,jj) = 1._wp / zhu(ji,jj)                   ! hu --> 1/hu
-         zhv(ji,jj) = 1._wp / zhv(ji,jj)
-      END_2D
-      !
-      zpsi_uw(:,:,:) = 0._wp
-      zpsi_vw(:,:,:) = 0._wp
-      !
+               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_2D
+         ENDIF
+         !
+         IF( nn_conv == 1 ) THEN              ! No MLE in case of convection
+            DO_2D( 1, 0, 1, 0 )
+               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_2D
+         ENDIF
+         !
+      ENDIF  ! end of ln_osm_mle conditional
+    !                                      !==  structure function value at uw- and vw-points  ==!
+    DO_2D( 1, 0, 1, 0 )
+       zhu(ji,jj) = 1._wp / MAX(zhu(ji,jj), rsmall)                   ! hu --> 1/hu
+       zhv(ji,jj) = 1._wp / MAX(zhv(ji,jj), rsmall) 
+    END_2D
+    !
+    zpsi_uw(:,:,:) = 0._wp
+    zpsi_vw(:,:,:) = 0._wp
+    !
       DO_3D( 1, 0, 1, 0, 2, ikmax )                ! start from 2 : surface value = 0
          zcuw = 1._wp - ( gdepw(ji+1,jj,jk,Kmm) + gdepw(ji,jj,jk,Kmm) ) * zhu(ji,jj)
@@ -206,20 +262,36 @@
       END DO
 
+      ! TEMP: [tiling] These changes not necessary if using XIOS (subdomain support)
       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
+         IF( ntile == 0 .OR. ntile == 1 ) THEN                             ! Do only on the first tile
+            ALLOCATE( zLf_NH(jpi,jpj), zpsiu_mle(jpi,jpj,jpk), zpsiv_mle(jpi,jpj,jpk) )
+            zpsiu_mle(:,:,:) = 0._wp ; zpsiv_mle(:,:,:) = 0._wp
+         ENDIF
+         !
+         IF (ln_osm_mle.and.ln_zdfosm) THEN
+            DO_2D( 0, 0, 0, 0 )
+               zLf_NH(ji,jj) = SQRT( rb_c * hmle(ji,jj) ) * r1_ft(ji,jj)      ! Lf = N H / f
+            END_2D
+         ELSE
+            DO_2D( 0, 0, 0, 0 )
+               zLf_NH(ji,jj) = SQRT( rb_c * zmld(ji,jj) ) * r1_ft(ji,jj)      ! Lf = N H / f
+            END_2D
+         ENDIF
          !
          ! 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
+         DO_3D( 0, 0, 0, 0, 1, ikmax+1 )
+            zpsiu_mle(ji,jj,jk) = zpsi_uw(ji,jj,jk) * r1_e2u(ji,jj)
+            zpsiv_mle(ji,jj,jk) = zpsi_vw(ji,jj,jk) * r1_e1v(ji,jj)
+         END_3D
+
+         IF( ntile == 0 .OR. ntile == nijtile ) THEN                       ! Do only on the last tile
+            CALL iom_put( "Lf_NHpf" , zLf_NH  )    ! Lf = N H / f
+            CALL iom_put( "psiu_mle", zpsiu_mle )    ! i-mle streamfunction
+            CALL iom_put( "psiv_mle", zpsiv_mle )    ! j-mle streamfunction
+            DEALLOCATE( zLf_NH, zpsiu_mle, zpsiv_mle )
+         ENDIF
       ENDIF
       !
    END SUBROUTINE tra_mle_trp
-
 
    SUBROUTINE tra_mle_init
