Implementation of an Eddy-Diffusivity Mass Flux (EDMF) scheme

src/gotm/gotm.F90

@@ -62,6 +62,11 @@ module gotm
    use turbulence,  only: Rig
    use turbulence,  only: kappa
    use turbulence,  only: clean_turbulence
+   use turbulence,  only: do_massflux
+   use turbulence,  only: compute_massflux
+   use turbulence,  only: massflux_on_dynamics 
+   use turbulence,  only: massflux_energy   
+   use turbulence,  only: Fmass,u_p,v_p,T_p,S_p,Pmf,mf_nsub
 
    use mtridiagonal,only: init_tridiagonal,clean_tridiagonal
 
@@ -564,6 +569,7 @@ subroutine initialize_gotm()
                Ti(1:nlev) = gsw_t_from_ct(S(1:nlev),T(1:nlev),-z(1:nlev))
          end select
    end select
+
 !GSW_KB
    u(1:nlev) = uprof_input%data(1:nlev)
    v(1:nlev) = vprof_input%data(1:nlev)
@@ -847,12 +853,16 @@ subroutine integrate_gotm()
          end if
    end select
 !GSW
-!  update shear and stratification
+   if(compute_massflux) then
+      call do_density(nlev,S,T,-z,-zi)
+      call do_massflux     ( nlev,dt,u,v,T,S,rho,h )             ! compute mass flux quantities from the meanflow variables
+      call uvTSequation_mf( nlev,dt,Fmass,u_p,v_p,T_p,S_p,Pmf,mf_nsub,massflux_on_dynamics, massflux_energy ) ! add the mass flux term to the meanflow equation 
+   endif 
+
    call shear(nlev,cnpar)
    call do_density(nlev,S,T,-z,-zi)
    buoy(1:nlev) = -gravity*(rho_p(1:nlev)-rho0)/rho0
    call stratification(nlev)
-
 #ifdef SPM
       if (spm_calc) then
          call set_env_spm(nlev,rho0,depth,u_taub,h,u,v,nuh,tx,ty,Hs,Tz,Phiw)

src/gotm/register_all_variables.F90

@@ -549,10 +549,28 @@ subroutine register_turbulence_variables(nlev)
    call fm%register('nuh', 'm2/s', 'turbulent diffusivity of heat', standard_name='??', dimensions=(/id_dim_zi/), data1d=nuh(0:nlev),category='turbulence', part_of_state=.true.)
    call fm%register('nus', 'm2/s', 'turbulent diffusivity of salt', standard_name='??', dimensions=(/id_dim_zi/), data1d=nus(0:nlev),category='turbulence', part_of_state=.true.)
    call fm%register('nucl', 'm2/s', 'turbulent diffusivity of momentum down Stokes gradient', standard_name='??', dimensions=(/id_dim_zi/), data1d=nucl(0:nlev),category='turbulence', part_of_state=.true.)
-   call fm%register('gamu', 'm2/s2', 'non-local flux of u-momentum', standard_name='??', dimensions=(/id_dim_zi/), data1d=gamu(0:nlev),category='turbulence')
-   call fm%register('gamv', 'm2/s2', 'non-local flux of v-momentum', standard_name='??', dimensions=(/id_dim_zi/), data1d=gamv(0:nlev),category='turbulence')
-   call fm%register('gamh', 'K m/s', 'non-local heat flux', standard_name='??', dimensions=(/id_dim_zi/), data1d=gamh(0:nlev),category='turbulence')
-   call fm%register('gams', 'g/kg m/s', 'non-local salinity flux', standard_name='??', dimensions=(/id_dim_zi/), data1d=gams(0:nlev),category='turbulence')
+   if ( compute_massflux ) then 
+      call fm%register('a_p', 'm2/s2', 'convective plume fractional area', standard_name='??', dimensions=(/id_dim_zi/), data1d=a_p(0:nlev),category='turbulence')
+      call fm%register('w_p', 'm/s', 'convective plume vertical velocity', standard_name='??', dimensions=(/id_dim_zi/), data1d=w_p(0:nlev),category='turbulence')
+      call fm%register('t_p', 'Celsius', 'convective plume temperature', standard_name='??', dimensions=(/id_dim_zi/), data1d=T_p(0:nlev),category='turbulence')
+      call fm%register('fmass', 'm/s', 'convective mass flux', standard_name='??', dimensions=(/id_dim_zi/), data1d=Fmass(0:nlev),category='turbulence')
+      call fm%register('EmD', '', 'Plume entrainment minus detrainment', standard_name='??', dimensions=(/id_dim_zi/), data1d=EmD(0:nlev),category='turbulence') 
+      call fm%register('zinv', 'm', 'Penetration depth for convective plumes', standard_name='??', data0d=mf_zinv,category='surface') 
+      if(massflux_on_dynamics) then
+         call fm%register('u_p', 'm/s', 'convective plume x-velocity', standard_name='??', dimensions=(/id_dim_zi/), data1d=u_p(0:nlev),category='turbulence')
+         call fm%register('v_p', 'm/s', 'convective plume y-velocity', standard_name='??', dimensions=(/id_dim_zi/), data1d=v_p(0:nlev),category='turbulence')
+      endif
+      if(massflux_energy) then
+         call fm%register('Pmf', 'm2/s3', 'shear production by convective plumes', standard_name='??', dimensions=(/id_dim_zi/), data1d=Pmf(0:nlev),category='turbulence/shear')
+         call fm%register('Gmf', 'm2/s3', 'buoyancy production by convective plumes', standard_name='??', dimensions=(/id_dim_zi/), data1d=Bmf(0:nlev),category='turbulence/buoyancy')
+         call fm%register('tke_p', 'm2/s2', 'Subplume turbulent kinetic energy', standard_name='??', dimensions=(/id_dim_zi/), data1d=tke_p(0:nlev),category='turbulence')
+      endif 
+   else 
+      call fm%register('gamu', 'm2/s2', 'non-local flux of u-momentum', standard_name='??', dimensions=(/id_dim_zi/), data1d=gamu(0:nlev),category='turbulence')
+      call fm%register('gamv', 'm2/s2', 'non-local flux of v-momentum', standard_name='??', dimensions=(/id_dim_zi/), data1d=gamv(0:nlev),category='turbulence')
+      call fm%register('gamh', 'K m/s', 'non-local heat flux', standard_name='??', dimensions=(/id_dim_zi/), data1d=gamh(0:nlev),category='turbulence')
+      call fm%register('gams', 'g/kg m/s', 'non-local salinity flux', standard_name='??', dimensions=(/id_dim_zi/), data1d=gams(0:nlev),category='turbulence')
+   endif
    call fm%register('Rig', '', 'gradient Richardson number', standard_name='??', dimensions=(/id_dim_zi/), data1d=Rig(0:nlev),category='turbulence')
 
 #ifdef _CVMIX_

src/meanflow/CMakeLists.txt

@@ -12,6 +12,7 @@ add_library(meanflow
    updategrid.F90
    vequation.F90
    wequation.F90
+   uvTSequation_mf.F90
 )
 target_link_libraries(meanflow PRIVATE gsw_static stokes_drift observations airsea_driver config)
 set_property(TARGET meanflow PROPERTY FOLDER gotm)

src/meanflow/uvTSequation_mf.F90

@@ -0,0 +1,145 @@
+#include"cppdefs.h"
+!-----------------------------------------------------------------------
+!BOP
+!
+! !ROUTINE: Add mass flux term to the meanflow equations\label{sec:meanflow_massflux}
+!
+! !INTERFACE:
+   subroutine uvTSequation_mf( nlev,dt,Fmass,u_p,v_p,T_p,S_p,Pmf,nsub,mf_dyn,mf_energy )
+!
+! !DESCRIPTION:
+!  This subroutine computes the mass flux contribution to the tracer and momentum meanflow equations 
+!  \begin{align*}
+!    \partial_t\Theta &= {\cal F}_{\Theta} + \partial_z\left( F_M (\Theta_p - \Theta)  \right) \\
+!    \partial_t S     &= {\cal F}_{S} + \partial_z\left( F_M (S_p - S)  \right) \\
+!    \partial_t S     &= {\cal F}_{U} + \partial_z\left( F_M (U_p - U)  \right) \\
+!    \partial_t S     &= {\cal F}_{V} + \partial_z\left( F_M (V_p - V)  \right) 
+!    \comma
+!  \end{align*}
+!  where $F_M = -a_p w_p$, and ${\cal F}_{X}$ groups all the right-hand-side terms already added in  
+!  the temperature, salinity, uequation, and vequation subroutines. 
+!
+! !USES:
+   use meanflow,     only: h,u,v,T,S,uo,vo 
+
+   IMPLICIT NONE
+!
+! !INPUT PARAMETERS:
+
+!  number of vertical layers
+   integer, intent(in)                 :: nlev
+
+!  time step (s)
+   REALTYPE, intent(in)                :: dt
+
+!  mass flux (m/s)
+   REALTYPE, intent(in)                :: Fmass(0:nlev)
+
+!  convective plumes u-velocity (m/s)
+   REALTYPE, intent(in)                :: u_p(0:nlev)
+
+!  convective plumes v-velocity (m/s)
+   REALTYPE, intent(in)                :: v_p(0:nlev)
+
+!  convective plumes temperature (Celsius)
+   REALTYPE, intent(in)                :: T_p(0:nlev)
+
+!  convective plumes salinity (psu)
+   REALTYPE, intent(in)                :: S_p(0:nlev)
+
+!  TKE production term associated with convective plumes  
+   REALTYPE, intent(out)               :: Pmf(0:nlev)
+
+!  
+   integer, intent(in)                :: nsub
+
+!  apply mass flux to the dynamics 
+   logical, intent(in)                :: mf_dyn, mf_energy
+
+! !REVISION HISTORY:
+!  Original author(s): Florian Lemarié 
+!EOP
+!
+! !LOCAL VARIABLES:
+   integer                   :: i,isub
+   REALTYPE                  :: VFlux(0:nlev,2),dt_sub 
+   REALTYPE                  :: u_ED(0:nlev), v_ED(0:nlev)
+   REALTYPE                  :: u_star(0:nlev), v_star(0:nlev)
+   REALTYPE                  :: T_star(0:nlev), S_star(0:nlev)
+   REALTYPE                  :: idz,SSUnp12,SSVnp12,cff,normVel
+!
+!-----------------------------------------------------------------------
+!BOC
+!  
+   if( mf_dyn .and. mf_energy) then 
+      u_ED(:) = u(:)   ! Save the current value of u,v for which only ED has been applied 
+      v_ED(:) = v(:)
+   endif 
+
+   dt_sub  = dt / real(nsub)
+
+!  Substepping 
+
+   do isub = 1, nsub
+      !  save old value
+      T_star = T 
+      S_star = S
+      !
+      VFlux(0,:) = 0. !; VFlux(nlev) = 0.
+      ! Compute fluxes associated to mass flux
+      do i = 1,nlev 
+         VFlux(i,1) = Fmass(i)*(t_p(i)-T_star(i)) 
+         VFlux(i,2) = Fmass(i)*(s_p(i)-s_star(i)) 
+      enddo 
+      ! Apply flux divergence
+      do i = 1,nlev
+         T(i)=T_star(i)+dt_sub*(VFlux(i,1)-VFlux(i-1,1))/h(i) 
+         S(i)=S_star(i)+dt_sub*(VFlux(i,2)-VFlux(i-1,2))/h(i) 
+      enddo
+      !
+   enddo
+   !
+   if( mf_dyn ) then 
+      do isub = 1, nsub
+         !  save old value
+         u_star = u
+         v_star = v
+         !
+         VFlux(0,:) = 0. !; VFlux(nlev) = 0.
+         ! Compute fluxes associated to mass flux
+         do i = 1,nlev 
+            VFlux(i,1) = Fmass(i)*(u_p(i)-u_star(i)) 
+            VFlux(i,2) = Fmass(i)*(v_p(i)-v_star(i)) 
+         enddo 
+         ! Apply flux divergence
+         do i = 1,nlev
+            u(i)=u_star(i)+dt_sub*(VFlux(i,1)-VFlux(i-1,1))/h(i) 
+            v(i)=v_star(i)+dt_sub*(VFlux(i,2)-VFlux(i-1,2))/h(i) 
+         enddo
+      !
+      enddo
+   endif
+
+! Compute additional TKE production term for energetic consistency 
+
+   if( mf_dyn .and. mf_energy ) then 
+      do i = 1,nlev-1
+         idz     = 2./(h(i+1)+h(i))
+         SSUnp12 = 0.5*( u(i+1)-u(i) + uo(i+1)-uo(i) )
+         SSVnp12 = 0.5*( v(i+1)-v(i) + vo(i+1)-vo(i) ) 
+         Pmf(i)  =   idz*SSUnp12*Fmass(i)*( u_p(i)-u_ED(i) )     & 
+                   + idz*SSVnp12*Fmass(i)*( v_p(i)-v_ED(i) )
+      enddo  
+   else 
+      do i = 1,nlev-1
+         Pmf(i)    =  0. 
+      enddo 
+   endif 
+   !
+   return
+   end subroutine uvTSequation_mf
+!EOC
+
+!-----------------------------------------------------------------------
+! Copyright by the GOTM-team under the GNU Public License - www.gnu.org
+!-----------------------------------------------------------------------

src/turbulence/CMakeLists.txt

@@ -8,6 +8,7 @@ add_library(turbulence
    cmue_d_h15.F90
    cmue_ma.F90
    cmue_sg.F90
+   cmue_nemo.F90
    compute_cpsi3.F90
    compute_rist.F90
    dissipationeq.F90
@@ -20,6 +21,8 @@ add_library(turbulence
    kbalgebraic.F90
    kbeq.F90
    lengthscaleeq.F90
+   lengthscale_gaspar.F90
+   massfluxeq.F90
    potentialml.F90
    production.F90
    q2over2eq.F90

src/turbulence/cmue_nemo.F90

@@ -0,0 +1,89 @@
+#include"cppdefs.h"
+!-----------------------------------------------------------------------
+!BOP
+!
+! !ROUTINE: The NEMO stability function\label{sec:stabNEMO}
+!
+! !INTERFACE:
+   subroutine cmue_nemo(nlev)
+!
+! !DESCRIPTION:
+!  This subroutine computes stability functions according to
+! \begin{equation}
+! c_{\mu}=c_{\mu}^{0}=0.1,\qquad c'_{\mu}=c_{\mu}^0 Pr_t^{-1}
+! \end{equation}
+! with constant $c_{\mu}^0 = 0.1$ and $Pr_t$ the turbulent Prandtl number. 
+! The empirical relation for the inverse turbulent Prandtl--number is 
+! \begin{equation}
+!   Pr_t^{-1} = \max\left(  0.1,  \mathrm{Ri}_c / \max( \mathrm{Ri}, \mathrm{Ri}_c ) \right)
+!   \comma
+! \end{equation}
+! where $Ri$ is the gradient Richardson--number and \mathrm{Ri}_c its critical 
+! value. This value is estimated from the local equilibrium $P+G = \epsilon$, 
+! which translates into (considering $l_\epsilon=\sqrt{2k/N^2}$):
+!\[
+!\mathrm{num} S^2 - \mathrm{nuh} N^2 = c_\epsilon \frac{k \sqrt{N^2}}{\sqrt{2}}
+!\]
+!which amounts to
+!\[
+!c_{\mu}^0 l_k \sqrt{k} \left( S^2 - Pr_t^{-1} N^2 \right) = c_\epsilon \frac{k \sqrt{N^2}}{\sqrt{2}}
+!\]
+!Considering that $l_k=\sqrt{2k/N^2}$ and $\mathrm{Ri} = N^2 / S^2$ we get 
+!\[
+!1 - Pr_t^{-1} \mathrm{Ri} = \frac{c_\epsilon}{2 c_{\mu}^0} \mathrm{Ri} \qquad \longrightarrow \qquad \mathrm{Ri}_c = \frac{1}{1+\frac{1}{2}\frac{c_\epsilon}{c_{\mu}^0}} \approx 0.22
+!\]
+!
+! !USES:
+
+   use turbulence, only: cmue1,cmue2
+   use turbulence, only: P,B  
+   use turbulence, only: num,nuh
+   IMPLICIT NONE
+!
+! !INPUT PARAMETERS:
+   integer, intent(in)                 :: nlev
+!
+! !REVISION HISTORY:
+!  Original author(s): Florian Lemarié 
+!
+!EOP
+!
+! !LOCAL VARIABLES:
+   integer                   :: i
+   REALTYPE                  :: Ri(0:nlev)
+   REALTYPE                  :: Ri_loc,Ri_cri,Denom
+   REALTYPE,parameter        :: limit=0.1
+   REALTYPE,parameter        :: bshear=1.e-16
+!
+!-----------------------------------------------------------------------
+!BOC
+   cmue1(0:nlev) = 0.1
+   Ri_cri = 2.*cmue1(1)/( 2.*cmue1(1) + 0.5*SQRT(2.) )
+
+   do i=1,nlev-1
+      if( B(i) >= 0.0 ) then
+         Ri(i) = 0.
+      else
+         Denom = nuh(i)*P(i)
+         if (Denom == 0.0) then
+            Ri(i) = -B(i)*num(i) / bshear
+         else
+            Ri(i) = -B(i)*num(i) / Denom
+         endif
+      endif
+   end do
+   !
+   Ri(0)=Ri(1); Ri(nlev)=Ri(nlev-1)
+   !
+   do i=1,nlev-1
+      Ri_loc   = 0.25*(2.*Ri(i)+Ri(i+1)+Ri(i-1))
+      cmue2(i) = cmue1(i)*MAX( limit, Ri_cri / MAX( Ri_cri , Ri_loc ) )
+   enddo   
+   !
+   return
+   end subroutine cmue_nemo
+!EOC
+
+!-----------------------------------------------------------------------
+! Copyright by the GOTM-team under the GNU Public License - www.gnu.org
+!-----------------------------------------------------------------------