Index: NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/cfgs/SHARED/field_def_nemo-oce.xml
===================================================================
--- NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/cfgs/SHARED/field_def_nemo-oce.xml	(revision 13824)
+++ NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/cfgs/SHARED/field_def_nemo-oce.xml	(revision 13826)
@@ -234,4 +234,9 @@
         <field id="cfl_cw"       long_name="w-courant number"   unit="#" />
 
+        <!-- variables available with ln_zdfmfc=.true. -->
+         <field id="mf_Tp"       long_name="plume_temperature"      standard_name="plume_temperature"     unit="degC"   grid_ref="grid_T_3D" />
+         <field id="mf_Sp"       long_name="plume_salinity"         standard_name="plume_salinity"        unit="1e-3"   grid_ref="grid_T_3D" />
+         <field id="mf_mf"       long_name="mass flux"              standard_name="mf_mass_flux"          unit="m"      grid_ref="grid_T_3D" />
+
       </field_group> <!-- grid_T -->
 
@@ -649,4 +654,9 @@
 	<field id="avm_evd"      long_name="convective enhancement of vertical viscosity"   standard_name="ocean_vertical_momentum_diffusivity_due_to_convection"   unit="m2/s" />
 
+        <!-- mf_app and mf_wp: available with ln_zdfmfc -->
+         <field id="mf_app"      long_name="convective area"        standard_name="mf_convective_area"    unit="%"      grid_ref="grid_W_3D" />
+         <field id="mf_wp"       long_name="convective velocity"    standard_name="mf_convective_velo"    unit="m/s"    grid_ref="grid_W_3D" />
+
+
         <!-- avt_tide: available with ln_zdfiwm=T -->
         <field id="av_ratio"     long_name="S over T diffusivity ratio"                     standard_name="salinity_over_temperature_diffusivity_ratio"                     unit="1"    />
Index: NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/cfgs/SHARED/namelist_ref
===================================================================
--- NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/cfgs/SHARED/namelist_ref	(revision 13824)
+++ NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/cfgs/SHARED/namelist_ref	(revision 13826)
@@ -1104,4 +1104,5 @@
       nn_npc      =    1         ! frequency of application of npc
       nn_npcp     =  365         ! npc control print frequency
+   ln_zdfmfc   = .false.      !  Mass Flux Convection
    !
    ln_zdfddm   = .false.   ! double diffusive mixing
@@ -1213,4 +1214,15 @@
       !                        !  = 1: Pierson Moskowitz wave spectrum
       !                        !  = 0: Constant La# = 0.3
+/
+!-----------------------------------------------------------------------
+&namzdf_mfc     !   Mass Flux Convection
+!-----------------------------------------------------------------------
+   ln_edmfuv  = .false.        ! Activate on velocity fields (Not available yet)
+   rn_cemf    =  1.            ! entrain/detrain coef. (<0 => cte; >0 % depending on dW/dz
+   rn_cwmf    = -0.            ! entrain/detrain coef. (<0 => cte; >0 % depending on dW/dz
+   rn_cent    = 2.e-5          ! entrain of convective area
+   rn_cdet    = 3.e-5          ! detrain of convective area 
+   rn_cap     = 0.9            ! Coef. for CAP estimation
+   App_max    = 0.1            ! Maximum convection area (% of the cell)
 /
 !-----------------------------------------------------------------------
Index: NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/C1D/step_c1d.F90
===================================================================
--- NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/C1D/step_c1d.F90	(revision 13824)
+++ NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/C1D/step_c1d.F90	(revision 13826)
@@ -104,4 +104,5 @@
       IF( ln_tradmp )   CALL tra_dmp( kstp, Nbb, Nnn, ts, Nrhs  )  ! internal damping trends- tracers
       IF(.NOT.ln_linssh)CALL tra_adv( kstp, Nbb, Nnn, ts, Nrhs  )  ! horizontal & vertical advection
+      IF( ln_zdfmfc  )  CALL tra_mfc( kstp, Nbb     , ts, Nrhs  )  ! Mass Flux Convection
       IF( ln_zdfosm  )  CALL tra_osm( kstp, Nnn     , ts, Nrhs  )  ! OSMOSIS non-local tracer fluxes
                         CALL tra_zdf( kstp, Nbb, Nnn, Nrhs, ts, Naa   )         ! vertical mixing
@@ -122,7 +123,4 @@
                         CALL dyn_atf    ( kstp, Nbb, Nnn, Naa , uu, vv, e3t, e3u, e3v )  ! time filtering of "now" fields
       IF(.NOT.ln_linssh)CALL ssh_atf    ( kstp, Nbb, Nnn, Naa , ssh )                    ! time filtering of "now" sea surface height
-      IF( kstp == nit000 .AND. ln_linssh) THEN
-         ssh(:,:,Naa) = ssh(:,:,Nnn)  ! init ssh after in ln_linssh case
-      ENDIF
       !
       ! Swap time levels
Index: NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/TRA/eosbn2.F90
===================================================================
--- NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/TRA/eosbn2.F90	(revision 13824)
+++ NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/TRA/eosbn2.F90	(revision 13826)
@@ -55,5 +55,5 @@
    !                  !! * Interface
    INTERFACE eos
-      MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d
+      MODULE PROCEDURE eos_insitu, eos_insitu_pot, eos_insitu_2d, eos_insitu_pot_2d
    END INTERFACE
    !
@@ -538,4 +538,75 @@
 
 
+   SUBROUTINE eos_insitu_pot_2d( pts, prhop )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE eos_insitu_pot  ***
+      !!
+      !! ** Purpose :   Compute the in situ density (ratio rho/rho0) and the
+      !!      potential volumic mass (Kg/m3) from potential temperature and
+      !!      salinity fields using an equation of state selected in the
+      !!     namelist.
+      !!
+      !! ** Action  :
+      !!              - prhop, the potential volumic mass (Kg/m3)
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj,jpts), INTENT(in   ) ::   pts    ! 1 : potential temperature  [Celsius]
+      !                                                                ! 2 : salinity               [psu]
+      REAL(wp), DIMENSION(jpi,jpj     ), INTENT(  out) ::   prhop  ! potential density (surface referenced)
+      !
+      INTEGER  ::   ji, jj, jk, jsmp             ! dummy loop indices
+      INTEGER  ::   jdof
+      REAL(wp) ::   zt , zh , zstemp, zs , ztm   ! local scalars
+      REAL(wp) ::   zn , zn0, zn1, zn2, zn3      !   -      -
+      REAL(wp), DIMENSION(:), ALLOCATABLE :: zn0_sto, zn_sto, zsign    ! local vectors
+      !!----------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('eos-pot')
+      !
+      SELECT CASE ( neos )
+      !
+      CASE( np_teos10, np_eos80 )                !==  polynomial TEOS-10 / EOS-80 ==!
+         !
+            DO_2D( 1, 1, 1, 1 )
+               !
+               zt  = pts (ji,jj,jp_tem) * r1_T0                           ! temperature
+               zs  = SQRT( ABS( pts(ji,jj,jp_sal) + rdeltaS ) * r1_S0 )   ! square root salinity
+               ztm = tmask(ji,jj,1)                                         ! tmask
+               !
+               zn0 = (((((EOS060*zt   &
+                  &   + EOS150*zs+EOS050)*zt   &
+                  &   + (EOS240*zs+EOS140)*zs+EOS040)*zt   &
+                  &   + ((EOS330*zs+EOS230)*zs+EOS130)*zs+EOS030)*zt   &
+                  &   + (((EOS420*zs+EOS320)*zs+EOS220)*zs+EOS120)*zs+EOS020)*zt   &
+                  &   + ((((EOS510*zs+EOS410)*zs+EOS310)*zs+EOS210)*zs+EOS110)*zs+EOS010)*zt   &
+                  &   + (((((EOS600*zs+EOS500)*zs+EOS400)*zs+EOS300)*zs+EOS200)*zs+EOS100)*zs+EOS000
+                  !
+               !
+               prhop(ji,jj) = zn0 * ztm                           ! potential density referenced at the surface
+               !
+            END_2D
+
+      CASE( np_seos )                !==  simplified EOS  ==!
+         !
+         DO_2D( 1, 1, 1, 1 )
+            zt  = pts  (ji,jj,jp_tem) - 10._wp
+            zs  = pts  (ji,jj,jp_sal) - 35._wp
+            ztm = tmask(ji,jj,1)
+            !                                                     ! potential density referenced at the surface
+            zn =  - rn_a0 * ( 1._wp + 0.5_wp*rn_lambda1*zt ) * zt   &
+               &  + rn_b0 * ( 1._wp - 0.5_wp*rn_lambda2*zs ) * zs   &
+               &  - rn_nu * zt * zs
+            prhop(ji,jj) = ( rho0 + zn ) * ztm
+            !
+         END_2D
+         !
+      END SELECT
+      !RBB IF(sn_cfctl%l_prtctl)   CALL prt_ctl( tab2d_1=prd, clinfo1=' eos-pot: ', tab2d_2=prhop, clinfo2=' pot : ', kdim=1 )
+      !
+      IF( ln_timing )   CALL timing_stop('eos-pot')
+      !
+   END SUBROUTINE eos_insitu_pot_2d
+
+
    SUBROUTINE rab_3d( pts, pab, Kmm )
       !!----------------------------------------------------------------------
Index: NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/TRA/trazdf.F90
===================================================================
--- NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/TRA/trazdf.F90	(revision 13824)
+++ NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/TRA/trazdf.F90	(revision 13826)
@@ -17,4 +17,5 @@
    USE phycst         ! physical constant
    USE zdf_oce        ! ocean vertical physics variables
+   USE zdfmfc         ! Mass FLux Convection 
    USE sbc_oce        ! surface boundary condition: ocean
    USE ldftra         ! lateral diffusion: eddy diffusivity
@@ -189,4 +190,9 @@
             ENDIF
             !
+            ! Modification of diagonal to add MF scheme
+            IF ( ln_zdfmfc ) THEN
+               CALL diag_mfc( zwi, zwd, zws, p2dt, Kaa )
+            END IF
+            !
             !! Matrix inversion from the first level
             !!----------------------------------------------------------------------
@@ -217,7 +223,12 @@
          ENDIF 
          !         
+         ! Modification of rhs to add MF scheme
+         IF ( ln_zdfmfc ) THEN
+            CALL rhs_mfc( pt(:,:,:,jn,Krhs), jn )
+         END IF
+         !
          DO_2D( 0, 0, 0, 0 )         !* 2nd recurrence:    Zk = Yk - Ik / Tk-1  Zk-1
             pt(ji,jj,1,jn,Kaa) =        e3t(ji,jj,1,Kbb) * pt(ji,jj,1,jn,Kbb)    &
-               &               + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs)
+               &               + p2dt * e3t(ji,jj,1,Kmm) * pt(ji,jj,1,jn,Krhs) 
          END_2D
          DO_3D( 0, 0, 0, 0, 2, jpkm1 )
Index: NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/ZDF/zdf_oce.F90
===================================================================
--- NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/ZDF/zdf_oce.F90	(revision 13824)
+++ NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/ZDF/zdf_oce.F90	(revision 13826)
@@ -40,4 +40,5 @@
    LOGICAL , PUBLIC ::   ln_zdfswm   !: surface  wave-induced mixing flag
    LOGICAL , PUBLIC ::   ln_zdfiwm   !: internal wave-induced mixing flag
+   LOGICAL , PUBLIC ::   ln_zdfmfc   !: convection: eddy diffusivity Mass Flux Convection
    !                             ! coefficients 
    REAL(wp), PUBLIC ::   rn_avm0     !: vertical eddy viscosity (m2/s)
Index: NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/ZDF/zdfmfc.F90
===================================================================
--- NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/ZDF/zdfmfc.F90	(revision 13826)
+++ NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/ZDF/zdfmfc.F90	(revision 13826)
@@ -0,0 +1,472 @@
+ MODULE zdfmfc
+   !!======================================================================
+   !!                       ***  MODULE  zdfmfc  ***
+   !! Ocean physics: Mass-Flux scheme parameterization of Convection:
+   !!                Non-local transport for the convective ocean boundary
+   !!                layer. Subgrid-scale large eddies are represented by a
+   !!                mass-flux contribution (ln_zdfmfc = .TRUE.)
+   !!======================================================================
+   !! History : NEMO  !
+   !!            3.6  !  2016-06  (H. Giordani, R. Bourdallé-Badie)  Original code
+   !!            4.2  !  2020-12  (H. Giordani, R. Bourdallé-Badie)  adapt to NEM04.2
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   tra_mfc       : Compute the Mass Flux and trends of T/S
+   !!   diag_mfc      : Modify diagonal of trazdf Matrix
+   !!   rhs_mfc       : Modify RHS of trazdf Matrix
+   !!   zdf_mfc_init  : initialization, namelist read, and parameters control
+   !!----------------------------------------------------------------------
+   !
+   USE oce            ! ocean dynamics and active tracers
+   USE dom_oce        ! ocean space and time domain
+   USE domvvl         ! ocean space and time domain : variable volume layer
+   USE domzgr
+   USE zdf_oce        ! ocean vertical physics
+   USE sbc_oce        ! surface boundary condition: ocean
+   USE phycst         ! physical constants
+   USE eosbn2         ! equation of state (eos routine)
+   USE zdfmxl         ! mixed layer
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE lib_mpp        ! MPP manager
+   USE prtctl         ! Print control
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O manager library
+   USE timing         ! Timing
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   tra_mfc         ! routine called in step module
+   PUBLIC   diag_mfc        ! routine called in trazdf module
+   PUBLIC   rhs_mfc         ! routine called in trazdf module
+   PUBLIC   zdf_mfc_init    ! routine called in nemo module
+   !
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::  edmfa, edmfb, edmfc   !: diagonal term of the matrix.
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::  edmftra               !: y term for matrix inversion
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::  edmfm               !: y term for matrix inversion
+   !
+   !! ** Namelist  namzdf_edmf  **
+   REAL(wp) ::   rn_cemf           ! entrain of T/S
+   REAL(wp) ::   rn_cwmf           ! detrain of T/S
+   REAL(wp) ::   rn_cent           ! entrain of the convective mass flux
+   REAL(wp) ::   rn_cdet           ! detrain of the convective mass flux
+   REAL(wp) ::   rn_cap            ! Factor of computation for convective area (negative => area constant)
+   REAL(wp) ::   App_max           ! Maximum of the convective area
+   LOGICAL, PUBLIC, SAVE  ::   ln_edmfuv         !: EDMF flag for velocity  !
+   !
+   !! * Substitutions
+#  include "do_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OCE 4.2 , NEMO Consortium (2018)
+   !! $Id: zdfmfc.F90 13783 2020-20-02 15:30:22Z rbourdal $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   INTEGER FUNCTION zdf_mfc_alloc()
+      !!----------------------------------------------------------------------
+      !!                ***  FUNCTION zdf_edmf_alloc  ***
+      !!----------------------------------------------------------------------
+      ALLOCATE( edmfa(jpi,jpj,jpk) , edmfb(jpi,jpj,jpk) , edmfc(jpi,jpj,jpk)      &
+         &      , edmftra(jpi,jpj,jpk,2), edmfm(jpi,jpj,jpk) ,  STAT= zdf_mfc_alloc )
+         !
+      IF( lk_mpp             )   CALL mpp_sum ( 'zdfmfc', zdf_mfc_alloc )
+      IF( zdf_mfc_alloc /= 0 )   CALL ctl_warn('zdf_mfc_alloc: failed to allocate arrays')
+   END FUNCTION zdf_mfc_alloc
+
+
+   SUBROUTINE tra_mfc( kt, Kmm, pts, Krhs )
+      !!----------------------------------------------------------------------
+      !!                   ***  ROUTINE zdf_mfc  ***
+      !!
+      !! ** Purpose :      Compute a mass flux, depending on surface flux, over
+      !!            the instable part of the water column.
+      !!
+      !! ** Method  :     Compute surface instability and mix tracers until stable level
+      !!           
+      !!
+      !! ** Action  :      Compute convection plume and (ta,sa)-trends for trazdf (EDMF scheme)
+      !!
+      !! References :      
+      !!                   Giordani, Bourdallé-Badie and Madec JAMES 2020
+      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      INTEGER                                  , INTENT(in)    :: Kmm, Krhs ! time level indices
+      REAL(wp), DIMENSION(jpi,jpj,jpk,jpts,jpt), INTENT(inout) :: pts       ! active tracers and RHS of tracer equation
+      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   ztsp         ! T/S of the plume
+      REAL(wp), DIMENSION(jpi,jpj,jpk,2) ::   ztse         ! T/S at W point
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrwp          !
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zrwp2         !
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zapp          !
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zedmf         !
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zepsT, zepsW  !
+      !
+      REAL(wp), DIMENSION(jpi,jpj) :: zustar, zustar2   !
+      REAL(wp), DIMENSION(jpi,jpj) :: zuws, zvws, zsws, zfnet          !
+      REAL(wp), DIMENSION(jpi,jpj) :: zfbuo, zrautbm1, zrautb, zraupl
+      REAL(wp), DIMENSION(jpi,jpj) :: zwpsurf            !
+      REAL(wp), DIMENSION(jpi,jpj) :: zop0 , zsp0 !
+      REAL(wp), DIMENSION(jpi,jpj) :: zrwp_0, zrwp2_0  !
+      REAL(wp), DIMENSION(jpi,jpj) :: zapp0           !
+      REAL(wp), DIMENSION(jpi,jpj) :: zphp, zph, zphpm1, zphm1, zNHydro
+      REAL(wp), DIMENSION(jpi,jpj) :: zhcmo          !
+      !
+      REAL(wp), DIMENSION(jpi,jpj,jpk)   ::   zn2    ! N^2
+      REAL(wp), DIMENSION(jpi,jpj,2  ) ::   zab, zabm1, zabp ! alpha and beta
+     
+      REAL(wp) :: zrho, zrhop
+      REAL(wp) :: zcnh, znum, zden, zcoef1, zcoef2
+      REAL(wp) :: zca, zcb, zcd, zrw, zxl, zcdet, zctre
+      REAL(wp) :: zaw, zbw, zxw
+      REAL(wp) :: alpha
+     !
+      INTEGER, INTENT(in   )    ::   kt   ! ocean time-step index      !
+      !
+      INTEGER  ::   ji, jj, jk  ! dummy  loop arguments   
+      !
+      !------------------------------------------------------------------
+      ! Initialisation of coefficients
+      !------------------------------------------------------------------
+      zca          = 1._wp
+      zcb          = 1._wp
+      zcd          = 1._wp
+
+      !------------------------------------------------------------------
+      ! Surface boundary condition
+      !------------------------------------------------------------------
+      ! surface Stress
+      !--------------------
+      zuws(:,:) = utau(:,:) * r1_rho0 
+      zvws(:,:) = vtau(:,:) * r1_rho0 
+      zustar2(:,:) = SQRT(zuws(:,:)*zuws(:,:)+zvws(:,:)*zvws(:,:))
+      zustar(:,:)  = SQRT(zustar2(:,:))
+
+      ! Heat Flux
+      !--------------------
+      zfnet(:,:) = qns(:,:) + qsr(:,:)
+      zfnet(:,:) = zfnet(:,:) / (rho0 * rcp)
+
+      ! Water Flux
+      !---------------------
+      zsws(:,:) = emp(:,:)
+
+      !-------------------------------------------
+      ! Initialisation of prognostic variables
+      !-------------------------------------------
+      zrwp (:,:,:) =  0. ; zrwp2(:,:,:) =  0. ; zedmf(:,:,:) =  0.
+      zph  (:,:)   =  0. ; zphm1(:,:)   =  0. ; zphpm1(:,:)  =  0.
+      ztsp(:,:,:,:)=  0.
+
+      ! Tracers inside plume (ztsp) and environment (ztse)
+      ztsp(:,:,1,jp_tem) = pts(:,:,1,jp_tem,Kmm) * tmask(:,:,1)
+      ztsp(:,:,1,jp_sal) = pts(:,:,1,jp_sal,Kmm) * tmask(:,:,1)
+      ztse(:,:,1,jp_tem) = pts(:,:,1,jp_tem,Kmm) * tmask(:,:,1)
+      ztse(:,:,1,jp_sal) = pts(:,:,1,jp_sal,Kmm) * tmask(:,:,1)
+
+      CALL eos( ztse(:,:,1,:) ,  zrautb(:,:) )
+      CALL eos( ztsp(:,:,1,:) ,  zraupl(:,:) )
+
+      !-------------------------------------------
+      ! Boundary Condition of Mass Flux (plume velo.; convective area, entrain/detrain)
+      !-------------------------------------------
+      zhcmo(:,:) = e3t(:,:,1,Kmm)
+      zfbuo(:,:)   = 0.
+      WHERE ( ABS(zrautb(:,:)) > 1.e-20 ) zfbuo(:,:)   =   &
+         &      grav * ( 2.e-4_wp *zfnet(:,:) - 7.6E-4_wp*pts(:,:,1,jp_sal,Kmm)*zsws(:,:)/zrautb(:,:)) * zhcmo(:,:)
+
+      zedmf(:,:,1) = -0.065_wp*(ABS(zfbuo(:,:)))**(1._wp/3._wp)*SIGN(1.,zfbuo(:,:))
+      zedmf(:,:,1) = MAX(0., zedmf(:,:,1))
+
+      zwpsurf(:,:) = 2._wp/3._wp*zustar(:,:) + 2._wp/3._wp*ABS(zfbuo(:,:))**(1._wp/3._wp)
+      zwpsurf(:,:) = MAX(1.e-5_wp,zwpsurf(:,:))
+      zwpsurf(:,:) = MIN(1.,zwpsurf(:,:))
+
+      zapp(:,:,:)  = App_max
+      WHERE(zwpsurf .NE. 0.) zapp(:,:,1)   = MIN(MAX(0.,zedmf(:,:,1)/zwpsurf(:,:)), App_max)
+
+      zedmf(:,:,1) = -zedmf(:,:,1)
+      zrwp(:,:,1)  = zwpsurf(:,:)
+      zrwp2(:,:,1) = zwpsurf(:,:)**2
+      zepsT(:,:,:) = 0.001
+      zepsW(:,:,:) = 0.001
+
+
+      !--------------------------------------------------------------
+      ! Compute plume properties 
+      ! In the same loop on vert. levels computation of:
+      !    - Vertical velocity: zWp
+      !    - Convective Area: zAp
+      !    - Tracers properties inside the plume (if necessary): ztp
+      !---------------------------------------------------------------
+
+      DO jk= 2, jpk
+
+         ! Compute the buoyancy acceleration on T-points at jk-1
+         zrautbm1(:,:) = zrautb(:,:)
+         !CALL eos( pts (:,:,jk-1,:,Kmm) ,  zrautbm1(:,:) )
+         CALL eos( pts (:,:,jk  ,:,Kmm) ,  zrautb(:,:)   )
+         CALL eos( ztsp(:,:,jk-1,:    ) ,  zraupl(:,:)   )
+
+         zphm1(:,:)  = zphm1(:,:)  + grav * zrautbm1(:,:) * e3t(:,:,jk-1, Kmm)
+         zphpm1(:,:) = zphpm1(:,:) + grav * zraupl(:,:)   * e3t(:,:,jk-1, Kmm)
+         zph(:,:)    = zphm1(:,:)  + grav * zrautb(:,:)   * e3t(:,:,jk  , Kmm)
+
+         WHERE(zrautbm1 .NE. 0.) zfbuo(:,:)  =  grav * (zraupl(:,:) - zrautbm1(:,:)) / zrautbm1(:,:)
+
+         DO_2D( 0, 0, 0, 0 )
+
+            ! Compute Environment of Plume. Interpolation T/S (before time step) on W-points
+            zrw              =  (gdept(ji,jj,jk,Kmm) - gdepw(ji,jj,jk,Kmm)) &
+               &              / (gdept(ji,jj,jk,Kmm) - gdept(ji,jj,jk-1,Kmm))
+            ztse(ji,jj,jk,:) = (pts(ji,jj,jk,:,Kmm) * zrw + pts(ji,jj,jk-1,:,Kmm)*(1._wp - zrw) )*tmask(ji,jj,jk)
+
+            !---------------------------------------------------------------
+            ! Compute the vertical velocity on W-points
+            !---------------------------------------------------------------
+
+            ! Non-hydrostatic pressure terms in the wp2 equation
+            zcnh = 0.2 
+            znum = 0.5_wp  + zcnh - &
+                   (zcnh*grav*zraupl(ji,jj)/zph(ji,jj)+zcb*zepsW(ji,jj,jk-1)) &
+                   *e3t(ji,jj,jk-1,Kmm)*0.5_wp   
+            zden = 0.5_wp + zcnh + &
+                   (zcnh*grav*zraupl(ji,jj)/zph(ji,jj)+zcb*zepsW(ji,jj,jk-1)) &
+                   *e3t(ji,jj,jk-1,Kmm)*0.5_wp   
+
+            zcoef1 = zca*e3t(ji,jj,jk-1,Kmm) / zden
+            zcoef2 = znum/zden
+
+            ! compute wp2 
+            zrwp2(ji,jj,jk) = zcoef1*zfbuo(ji,jj) &
+                            + zcoef2*zrwp2(ji,jj,jk-1)
+            zrwp2(ji,jj,jk) = MAX ( zrwp2(ji,jj,jk)*wmask(ji,jj,jk) , 0.)
+            zrwp (ji,jj,jk) = SQRT( zrwp2(ji,jj,jk) )
+
+            !----------------------------------------------------------------------------------
+            ! Compute convective area on W-point
+            ! Compute vertical profil of the convective area with mass conservation hypothesis
+            ! If rn_cap negative => constant value on the water column.
+            !----------------------------------------------------------------------------------
+            IF( rn_cap .GT. 0. ) THEN
+
+               zxw = MAX(zrwp(ji,jj,jk-1), zrwp(ji,jj,jk) )
+               IF( zxw > 0. ) THEN
+
+                  zxl = (zrwp(ji,jj,jk-1)-zrwp(ji,jj,jk))/(e3t(ji,jj,jk-1,Kmm)*zxw)
+                  IF (zxl .LT. 0._wp) THEN
+                     zctre  = -1.*rn_cap*zxl 
+                     zcdet  =  0.
+                  ELSE
+                     zctre  =  0.
+                     zcdet  =  rn_cap*zxl 
+                  END IF
+                     zapp(ji,jj,jk) = zapp(ji,jj,jk-1)*     &
+                     &                (1._wp + (zxl + zctre - zcdet )*e3t(ji,jj,jk-1,Kmm))
+                  ELSE
+                     zapp(ji,jj,jk) = App_max
+                  END IF
+                  zapp(ji,jj,jk) = MIN( MAX(zapp(ji,jj,jk),0.), App_max)
+               ELSE
+                  zapp(ji,jj,jk) = -1. * rn_cap
+               END IF
+
+            ! Compute Mass Flux on W-point
+            zedmf(ji,jj,jk)   = -zapp(ji,jj,jk) * zrwp(ji,jj,jk)* wmask(ji,jj,jk)
+
+            ! Compute Entrainment coefficient
+            IF(rn_cemf .GT. 0.) THEN
+               zxw = 0.5*(zrwp(ji,jj,jk-1)+ zrwp(ji,jj,jk) )
+               zepsT(ji,jj,jk)  =  0.01
+               IF( zxw > 0.  ) THEN
+                  zepsT(ji,jj,jk)  =  zepsT(ji,jj,jk) +                       &
+                                   &  1.*ABS( zrwp(ji,jj,jk-1)-zrwp(ji,jj,jk) )  &
+                                   &  / ( e3t(ji,jj,jk-1,Kmm) * zxw )
+                  zepsT(ji,jj,jk)  = zepsT(ji,jj,jk) * rn_cemf * wmask(ji,jj,jk)
+               ENDIF
+            ELSE
+               zepsT(ji,jj,jk)  = -rn_cemf
+            ENDIF
+
+            ! Compute the detrend coef for velocity (on W-point and not T-points, bug ???)
+            IF(rn_cwmf .GT. 0.) THEN
+               zepsW(ji,jj,jk)  =  rn_cwmf * zepsT(ji,jj,jk)
+            ELSE
+               zepsW(ji,jj,jk)  = -rn_cwmf
+            ENDIF
+
+            !---------------------------------------------------------------
+            ! Compute the plume properties on T-points
+            !---------------------------------------------------------------
+            IF(zrwp (ji,jj,jk) .LT. 1.e-12 .AND. zrwp (ji,jj,jk-1) .LT. 1.e-12) THEN
+               ztsp(ji,jj,jk-1,jp_tem) = pts(ji,jj,jk-1,jp_tem,Kmm)
+               ztsp(ji,jj,jk-1,jp_sal) = pts(ji,jj,jk-1,jp_sal,Kmm)
+            ENDIF
+
+            zcoef1 =  (1._wp-zepsT(ji,jj,jk)*(1._wp-zrw)*e3w(ji,jj,jk,Kmm)*wmask(ji,jj,jk ) ) &
+            &       / (1._wp+zepsT(ji,jj,jk)*zrw*e3w(ji,jj,jk,Kmm)*wmask(ji,jj,jk) )
+            !
+            zcoef2 =  zepsT(ji,jj,jk)*e3w(ji,jj,jk,Kmm)*wmask(ji,jj,jk)                       &
+            &       / (1._wp+zepsT(ji,jj,jk)*zrw*e3w(ji,jj,jk,Kmm)*wmask(ji,jj,jk))
+            !
+            ztsp(ji,jj,jk,jp_tem) = (zcoef1 * ztsp(ji,jj,jk-1,jp_tem) +  &
+            &                        zcoef2 * ztse(ji,jj,jk  ,jp_tem) )*tmask(ji,jj,jk)
+            ztsp(ji,jj,jk,jp_sal) = (zcoef1 * ztsp(ji,jj,jk-1,jp_sal) +  &
+            &                        zcoef2 * ztse(ji,jj,jk  ,jp_sal) )*tmask(ji,jj,jk)
+
+         END_2D 
+      END DO ! end of loop on jpk
+
+      ! Compute Mass Flux on T-point
+      DO jk=1,jpk-1
+         edmfm(:,:,jk) = (zedmf(:,:,jk+1)  + zedmf(:,:,jk) )*0.5
+      END DO
+
+      ! Save variable (on T point)
+      CALL iom_put( "mf_Tp" , ztsp(:,:,:,jp_tem) )  ! Save plume temperature
+      CALL iom_put( "mf_Sp" , ztsp(:,:,:,jp_sal) )  ! Save plume salinity
+      CALL iom_put( "mf_mf" , edmfm(:,:,:)       )  ! Save Mass Flux
+      ! Save variable (on W point)
+      CALL iom_put( "mf_wp" , zrwp (:,:,:)       )  ! Save convective velocity in the plume
+      CALL iom_put( "mf_app", zapp (:,:,:)       )  ! Save convective area
+
+      !=================================================================================
+      !  Computation of a tridiagonal matrix and right hand side terms of the linear system
+      !=================================================================================
+      edmfa(:,:,:)     = 0.
+      edmfb(:,:,:)     = 0.
+      edmfc(:,:,:)     = 0.  
+      edmftra(:,:,:,:) = 0.
+
+      !---------------------------------------------------------------
+      ! Diagonal terms 
+      !---------------------------------------------------------------
+      DO jk=1,jpk-1
+         edmfa(:,:,jk) =  0.
+         edmfb(:,:,jk) = -edmfm(:,:,jk  ) / e3w(:,:,jk+1,Kmm)
+         edmfc(:,:,jk) =  edmfm(:,:,jk+1) / e3w(:,:,jk+1,Kmm)
+      END DO
+      edmfa(:,:,jpk) = -edmfm(:,:,jpk-1) / e3w(:,:,jpk,Kmm)
+      edmfb(:,:,jpk) =  edmfm(:,:,jpk  ) / e3w(:,:,jpk,Kmm)
+      edmfc(:,:,jpk) =  0.
+
+      !---------------------------------------------------------------
+      ! right hand side term for Temperature
+      !---------------------------------------------------------------
+      DO jk=1,jpk-1
+        edmftra(:,:,jk,1) = - edmfm(:,:,jk  ) * ztsp(:,:,jk  ,jp_tem) / e3w(:,:,jk+1,Kmm) &
+                          & + edmfm(:,:,jk+1) * ztsp(:,:,jk+1,jp_tem) / e3w(:,:,jk+1,Kmm)
+      END DO
+      edmftra(:,:,jpk,1) = - edmfm(:,:,jpk-1) * ztsp(:,:,jpk-1,jp_tem) / e3w(:,:,jpk,Kmm) &
+                         & + edmfm(:,:,jpk  ) * ztsp(:,:,jpk  ,jp_tem) / e3w(:,:,jpk,Kmm)
+                        
+      !---------------------------------------------------------------
+      ! Right hand side term for Salinity
+      !---------------------------------------------------------------
+      DO jk=1,jpk-1
+         edmftra(:,:,jk,2) =  - edmfm(:,:,jk  ) * ztsp(:,:,jk  ,jp_sal) / e3w(:,:,jk+1,Kmm) &
+                           &  + edmfm(:,:,jk+1) * ztsp(:,:,jk+1,jp_sal) / e3w(:,:,jk+1,Kmm)
+      END DO
+      edmftra(:,:,jpk,2) = - edmfm(:,:,jpk-1) * ztsp(:,:,jpk-1,jp_sal) / e3w(:,:,jpk,Kmm) &
+                         & + edmfm(:,:,jpk  ) * ztsp(:,:,jpk  ,jp_sal) / e3w(:,:,jpk,Kmm)
+      !
+      !
+      CALL lbc_lnk_multi( 'zdfmfc', edmfm,'T',1., edmfa,'T',1., edmfb,'T',1., edmfc,'T',1., edmftra(:,:,:,jp_tem),'T',1., edmftra(:,:,:,jp_sal),'T',1.)
+      !
+   END SUBROUTINE tra_mfc
+
+   
+   SUBROUTINE diag_mfc( zdiagi, zdiagd, zdiags, p2dt, Kaa )
+
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::  zdiagi, zdiagd, zdiags  ! inout: tridaig. terms 
+      REAL(wp)                        , INTENT(in   ) ::   p2dt                   ! tracer time-step
+      INTEGER                         , INTENT(in   ) ::   Kaa                    ! ocean time level indices
+
+      INTEGER  ::   ji, jj, jk  ! dummy  loop arguments   
+
+         DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+            zdiagi(ji,jj,jk) = zdiagi(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfa(ji,jj,jk)
+            zdiags(ji,jj,jk) = zdiags(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfc(ji,jj,jk)
+            zdiagd(ji,jj,jk) = zdiagd(ji,jj,jk) + e3t(ji,jj,jk,Kaa) * p2dt *edmfb(ji,jj,jk)
+         END_3D
+
+   END SUBROUTINE diag_mfc
+
+   SUBROUTINE rhs_mfc( zrhs, jjn )
+
+      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   zrhs                   ! inout: rhs trend 
+      INTEGER                         , INTENT(in   ) ::   jjn                    ! tracer indices
+
+      INTEGER  ::   ji, jj, jk  ! dummy  loop arguments   
+
+      DO_3D( 0, 0, 0, 0, 1, jpkm1 )
+         zrhs(ji,jj,jk) = zrhs(ji,jj,jk) + edmftra(ji,jj,jk,jjn)
+      END_3D
+
+   END SUBROUTINE rhs_mfc
+
+
+
+   SUBROUTINE zdf_mfc_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_mfc_init  ***
+      !!                    
+      !! ** Purpose :   Initialization of the vertical eddy diffivity and
+      !!      mass flux
+      !!
+      !! ** Method  :   Read the namzdf_mfc namelist and check the parameters
+      !!      called at the first timestep (nit000)
+      !!
+      !! ** input   :   Namlist namzdf_mfc
+      !!
+      !! ** Action  :   Increase by 1 the nstop flag is setting problem encounter
+      !!
+      !!----------------------------------------------------------------------
+      !
+      INTEGER ::   jk    ! dummy loop indices
+      INTEGER ::   ios   ! Local integer output status for namelist read
+      REAL(wp)::   zcr   ! local scalar
+      !!
+      NAMELIST/namzdf_mfc/ ln_edmfuv, rn_cemf, rn_cwmf, rn_cent, rn_cdet, rn_cap, App_max
+      !!----------------------------------------------------------
+      !
+      !
+!      REWIND( numnam_ref )              ! Namelist namzdf_mfc in reference namelist : Vertical eddy diffivity mass flux
+      READ  ( numnam_ref, namzdf_mfc, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_edmf in reference namelist' )
+
+!      REWIND( numnam_cfg )              ! Namelist namzdf_mfc in configuration namelist : Vertical eddy diffivity mass flux
+      READ  ( numnam_cfg, namzdf_mfc, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_edmf in configuration namelist' )
+      IF(lwm) WRITE ( numond, namzdf_mfc )
+
+      IF(lwp) THEN                     !* Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'zdf_mfc_init'
+         WRITE(numout,*) '~~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namzdf_mfc : set eddy diffusivity Mass Flux Convection'
+         WRITE(numout,*) '   Apply mass flux on velocities (Not yet avail.)     ln_edmfuv = ', ln_edmfuv
+         WRITE(numout,*) '   Coeff for entrain/detrain T/S of plume (Neg => cte) rn_cemf  = ', rn_cemf
+         WRITE(numout,*) '   Coeff for entrain/detrain Wp of plume  (Neg => cte) rn_cwmf  = ', rn_cwmf
+         WRITE(numout,*) '   Coeff for entrain/detrain area of plume             rn_cap   = ', rn_cap
+         WRITE(numout,*) '   Coeff for entrain area of plume                     rn_cent  = ', rn_cent
+         WRITE(numout,*) '   Coeff for detrain area of plume                     rn_cdet  = ', rn_cdet
+         WRITE(numout,*) '   Max convective area                                 App_max  = ', App_max
+       ENDIF
+                                     !* allocate edmf arrays
+      IF( zdf_mfc_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_edmf_init : unable to allocate arrays' )
+      edmfa(:,:,:)     = 0.
+      edmfb(:,:,:)     = 0.
+      edmfc(:,:,:)     = 0.
+      edmftra(:,:,:,:) = 0.
+      !
+   END SUBROUTINE zdf_mfc_init
+
+   !!======================================================================
+ 
+   !!======================================================================
+END MODULE zdfmfc
+
+
+
Index: NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/ZDF/zdfphy.F90
===================================================================
--- NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/ZDF/zdfphy.F90	(revision 13824)
+++ NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/ZDF/zdfphy.F90	(revision 13826)
@@ -21,4 +21,5 @@
    USE zdfddm         ! vertical physics: double diffusion mixing      
    USE zdfevd         ! vertical physics: convection via enhanced vertical diffusion  
+   USE zdfmfc         ! vertical physics: Mass Flux Convection 
    USE zdfiwm         ! vertical physics: internal wave-induced mixing  
    USE zdfswm         ! vertical physics: surface  wave-induced mixing
@@ -78,4 +79,5 @@
       NAMELIST/namzdf/ ln_zdfcst, ln_zdfric, ln_zdftke, ln_zdfgls,   &     ! type of closure scheme
          &             ln_zdfosm,                                    &     ! type of closure scheme
+         &             ln_zdfmfc,                                    &     ! convection : mass flux
          &             ln_zdfevd, nn_evdm, rn_evd ,                  &     ! convection : evd
          &             ln_zdfnpc, nn_npc , nn_npcp,                  &     ! convection : npc
@@ -112,4 +114,5 @@
          WRITE(numout,*) '         OSMOSIS-OBL closure (OSM)               ln_zdfosm = ', ln_zdfosm
          WRITE(numout,*) '      convection: '
+         WRITE(numout,*) '         convection mass flux (mfc)              ln_zdfmfc = ', ln_zdfmfc
          WRITE(numout,*) '         enhanced vertical diffusion             ln_zdfevd = ', ln_zdfevd
          WRITE(numout,*) '            applied on momentum (=1/0)             nn_evdm = ', nn_evdm
@@ -172,4 +175,7 @@
       IF( ln_zdfnpc .AND. ln_zdfevd )   CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfnpc and ln_zdfevd' )
       IF( ln_zdfosm .AND. ln_zdfevd )   CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfosm and ln_zdfevd' )
+      IF( ln_zdfmfc .AND. ln_zdfevd )   CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfmfc and ln_zdfevd' )
+      IF( ln_zdfmfc .AND. ln_zdfnpc )   CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfmfc and ln_zdfnpc' )
+      IF( ln_zdfmfc .AND. ln_zdfosm )   CALL ctl_stop( 'zdf_phy_init: chose between ln_zdfmfc and ln_zdfosm' )
       IF( lk_top    .AND. ln_zdfnpc )   CALL ctl_stop( 'zdf_phy_init: npc scheme is not working with key_top' )
       IF( lk_top    .AND. ln_zdfosm )   CALL ctl_stop( 'zdf_phy_init: osmosis scheme is not working with key_top' )
@@ -178,4 +184,5 @@
          IF    ( ln_zdfnpc ) THEN  ;   WRITE(numout,*) '   ==>>>   convection: use non penetrative convective scheme'
          ELSEIF( ln_zdfevd ) THEN  ;   WRITE(numout,*) '   ==>>>   convection: use enhanced vertical diffusion scheme'
+         ELSEIF( ln_zdfmfc ) THEN  ;   WRITE(numout,*) '   ==>>>   convection: use Mass Flux scheme'
          ELSE                      ;   WRITE(numout,*) '   ==>>>   convection: no specific scheme used'
          ENDIF
@@ -205,5 +212,7 @@
       ELSE                   ;   l_zdfsh2 = .TRUE.
       ENDIF
-
+      !                          !== Mass Flux Convectiive algorithm  ==!
+      IF( ln_zdfmfc )   CALL zdf_mfc_init       ! Convection computed with eddy diffusivity mass flux
+      !
       !                          !== gravity wave-driven mixing  ==!
       IF( ln_zdfiwm )   CALL zdf_iwm_init       ! internal wave-driven mixing
Index: NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/step.F90
===================================================================
--- NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/step.F90	(revision 13824)
+++ NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/step.F90	(revision 13826)
@@ -262,4 +262,5 @@
 #endif
                          CALL tra_adv    ( kstp, Nbb, Nnn, ts, Nrhs )  ! hor. + vert. advection	==> RHS
+      IF( ln_zdfmfc  )   CALL tra_mfc    ( kstp, Nbb,      ts, Nrhs )  ! Mass Flux Convection 
       IF( ln_zdfosm  )   CALL tra_osm    ( kstp,      Nnn, ts, Nrhs )  ! OSMOSIS non-local tracer fluxes ==> RHS
       IF( lrst_oce .AND. ln_zdfosm ) &
Index: NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/step_oce.F90
===================================================================
--- NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/step_oce.F90	(revision 13824)
+++ NEMO/branches/2020/dev_r2052_ENHANCE-09_rbourdal_massfluxconvection/src/OCE/step_oce.F90	(revision 13826)
@@ -69,4 +69,5 @@
    USE zdfphy          ! vertical physics manager      (zdf_phy_init routine)
    USE zdfosm  , ONLY : osm_rst, dyn_osm, tra_osm      ! OSMOSIS routines used in step.F90
+   USE zdfmfc          ! Mass FLux Convection routine used in step.F90
 
    USE diu_layers      ! diurnal SST bulk and coolskin routines
