From 9b13185330217af99f2538f010947866bb73377c Mon Sep 17 00:00:00 2001 From: lemarie Date: Fri, 13 Mar 2026 18:15:40 +0100 Subject: [PATCH 1/5] Add new subroutines for tke 1st order scheme from NEMO and mass flux convection scheme --- src/meanflow/uvTSequation_mf.F90 | 145 ++++++++++ src/turbulence/cmue_nemo.F90 | 89 ++++++ src/turbulence/lengthscale_gaspar.F90 | 117 ++++++++ src/turbulence/massfluxeq.F90 | 388 ++++++++++++++++++++++++++ 4 files changed, 739 insertions(+) create mode 100644 src/meanflow/uvTSequation_mf.F90 create mode 100644 src/turbulence/cmue_nemo.F90 create mode 100644 src/turbulence/lengthscale_gaspar.F90 create mode 100644 src/turbulence/massfluxeq.F90 diff --git a/src/meanflow/uvTSequation_mf.F90 b/src/meanflow/uvTSequation_mf.F90 new file mode 100644 index 00000000..45ecbd29 --- /dev/null +++ b/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 +!----------------------------------------------------------------------- diff --git a/src/turbulence/cmue_nemo.F90 b/src/turbulence/cmue_nemo.F90 new file mode 100644 index 00000000..6e338c1f --- /dev/null +++ b/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 +!----------------------------------------------------------------------- diff --git a/src/turbulence/lengthscale_gaspar.F90 b/src/turbulence/lengthscale_gaspar.F90 new file mode 100644 index 00000000..8c1bb7c0 --- /dev/null +++ b/src/turbulence/lengthscale_gaspar.F90 @@ -0,0 +1,117 @@ +#include"cppdefs.h" +!------------------------------------------------------------------------- +!BOP +! +! !ROUTINE: Diagnostic length-scale adapted from Gaspar with two master scales \label{sec:lengthscale_gaspar} +! +! !INTERFACE: + subroutine lengthscale_gaspar(nlev,z0s,z0b,h,NN) + +! !DESCRIPTION: +! This subroutine computes the length scale used in the NEMO mode, see Madec et al. (2025). +! Two master length scales $l_u$ and $l_d$ are defined and initialised with the Deardorff +! canonical value $l_{d80}=\sqrt(2k/N^2)$. The Deardorff length scale corresponds to a local +! approximation of the Bougeault and André (1986) length scales. +! Starting from $l_{\rm u}=l_{\rm d}=l_{\rm d80}$, the length scales are then limited not only +! by the distance to the surface and to the top but also by the distance to a strongly stratified portion of the air column. +! This limitation amounts to controlling the vertical gradients of $l_{\rm u}$ and $l_{\rm d}$ so that they do not exceed the depth variations. +! From $l_u$ and $l_d$ two length--scales are defined: $l_k$, +! a characteristic mixing length, and $l_\epsilon$, a characteristic dissipation length. +! They are computed according to Gaspar et al. (1990) +! \begin{equation} +! \begin{array}{l} +! l_k(z)= \text{Min} ( l_d(z),l_u(z)) \comma \\[4mm] +! l_{\epsilon}(z)=\left( l_d(z)l_u(z)\right)^\frac{1}{2} +! \point +! \end{array} +! \end{equation} +! +! $l_k$ is used in {\tt kolpran()} to compute eddy viscosity/difussivity. +! $l_{\epsilon}$ is used to compute the dissipation rate $\epsilon$, +! according to +! \begin{equation} +! \epsilon=C_{\epsilon} k^{3/2} l_{\epsilon}^{-1} +! \comma +! C_{\epsilon}=\frac{\sqrt{2}}{2} +! \point +! \end{equation} +! +! !USES: + use turbulence, only: L,eps,tke,k_min,eps_min + use turbulence, only: kappa + + IMPLICIT NONE +! +! !INPUT PARAMETERS: + +! number of vertical layers + integer, intent(in) :: nlev + +! bottom and surface roughness (m) + REALTYPE, intent(in) :: z0b,z0s + +! layer thickness (m) + REALTYPE, intent(in) :: h(0:nlev) + +! buoyancy frequency (1/s^2) + REALTYPE, intent(in) :: NN(0:nlev) +! +! !REVISION HISTORY: +! Original author(s): Florian Lemarié +! +!EOP +!------------------------------------------------------------------------- +! !LOCAL VARIABLES: + integer :: i,j + REALTYPE :: lup (0:nlev),ldwn(0:nlev), ld80(0:nlev) + REALTYPE :: rn2 + REALTYPE :: L_min, ceps, leps + REALTYPE, parameter :: NNmin =1.e-20 + REALTYPE, parameter :: L_min0=0.04 +!BOC + ceps = 0.5*SQRT(2.) + L_min = ceps*k_min**1.5/eps_min +!------------------------------------------------------------------------- +! Calculation of lu and ld based on Deardorff (1980) and further +! limitations following the NEMO methodology +!------------------------------------------------------------------------- + do i=1,nlev-1 + rn2 = MAX( NN(i),NNmin ) + ld80(i)=SQRT( 2.*tke(i)/rn2 ) + enddo + ! + ld80(0 ) = MAX( L_min , kappa*z0b ) + ld80(nlev) = MAX( L_min0, kappa*z0s ) +!------------------------------------------------------------------------- +! Physical limits for the mixing lengths +!------------------------------------------------------------------------- + ldwn(0 ) = ld80(0) ! bottom boundary condition + do i = 1, nlev + ldwn(i) = MIN( ldwn(i-1) + h(i ) , ld80(i) ) + enddo + ! + lup(nlev) = ld80(nlev) ! surface boundary condition + do i = nlev-1,0,-1 + lup (i) = MIN( lup (i+1) + h(i+1) , ld80(i) ) + enddo +!------------------------------------------------------------------------- +! Calculation of L and eps +!------------------------------------------------------------------------- + L(0 )=kappa*z0b + L(nlev)=kappa*z0s + ! + do i=1,nlev-1 + L(i) = MIN( lup(i), ldwn(i) ) ! used in kolpran + end do + ! + do i=0,nlev + leps = SQRT( lup(i)*ldwn(i) ) + eps(i) = (ceps/leps)*sqrt(tke(i)*tke(i)*tke(i)) ! used in tkeeq + enddo + + return + end subroutine lengthscale_gaspar +!EOC +!----------------------------------------------------------------------- +! Copyright by the GOTM-team under the GNU Public License - www.gnu.org +!----------------------------------------------------------------------- diff --git a/src/turbulence/massfluxeq.F90 b/src/turbulence/massfluxeq.F90 new file mode 100644 index 00000000..616c3b22 --- /dev/null +++ b/src/turbulence/massfluxeq.F90 @@ -0,0 +1,388 @@ +#include"cppdefs.h" +!----------------------------------------------------------------------- +!BOP +! +! !ROUTINE: The convective plume nonlinear ODEs \label{sec:massfluxeq} +! +! !INTERFACE: + subroutine massfluxeq(nlev,dt,u,v,T,S,rho,h) + +! !DESCRIPTION: +! This subroutine solves a set of nonlinear ODEs to determine +! the plume fractional area $a_p$, the plume temperature $\Theta_p$ and salinity $S_p$, +! the plume horizontal momentum $U_p$ and $V_p$, and the plume vertical velocity $w_p$. +! The plume equations are (Perrot and Lemarié, 2025): +! \begin{align*} +! \partial_z(a_ w_p) &= E-D \\ +! a_p w_p \partial_z \Theta_p &= E( \Theta - \Theta_p ) \\ +! a_p w_p \partial_z S_p &= E( S - S_p ) \\ +! a_p w_p \partial_z U_p &= E( U - U_p ) + a_p w_p C_u \partial_z U \\ +! a_p w_p \partial_z V_p &= E( V - V_p ) + a_p w_p C_u \partial_z V \\ +! a_p w_p \partial_z w_p &= -b E w_p + a_p \left( a B_p + \frac{b'}{h} w_p^2 \right) \\ +! B_p &= b_{\rm eos}( \Theta_p, S_p ) - b_{\rm eos}( \Theta, S ) \\ +! E &= a_p C_{ent} \max(0,\partial_z w_p) \\ +! D &= -a_p C_{det} \min(0,\partial_z w_p) - a_p w_p \frac{\delta_0}{h} +! \end{align*} +! where $a,b,C_u,b',C_{ent},C_{det},\delta_0$ are user provided parameters. +! +! !USES: + use turbulence, only: Fmass,a_p,w_p,u_p,v_p,T_p,S_p,EmD + use turbulence, only: mf_ap0,mf_wp0,mf_Cent,mf_Cdet, mf_d0 + use turbulence, only: mf_aa,mf_bb,mf_bp,mf_uv,mf_dbkg + use turbulence, only: mf_zinv, Bmf, massflux_energy, mf_nsub + use meanflow, only: gravity, z, zi + use density, only: rho0, get_rho +! + IMPLICIT NONE +! +! !INPUT PARAMETERS: + +! number of vertical layers + integer, intent(in) :: nlev + +! time step (s) + REALTYPE, intent(in) :: dt + +! horizontal velocity (m/s) + REALTYPE, intent(in) :: u(0:nlev),v(0:nlev) + +! temperature and salinity + REALTYPE, intent(in) :: T(0:nlev),S(0:nlev) + +! + REALTYPE, intent(in) :: rho(0:nlev) + +! layer thickness (m) + REALTYPE, intent(in) :: h(0:nlev) + +! !REVISION HISTORY: +! Original author(s): Florian Lemarié +!EOP +!------------------------------------------------------------------------ +! +! !LOCAL VARIABLES: + REALTYPE, parameter :: wp_min = 0.5E-08 ! min value for plume velocity + REALTYPE :: mf_one_p_bCent, cff, cup, cdwn, cff1, Crt_number(0:nlev) + REALTYPE :: bp,d0,izinv,zinvMin, dpth, B_p, rhs_wp2, Crt_max + REALTYPE :: app,wpp,wpp2,wpm2,wpm,temp_p,salt_p,rho_p, dz, nrj + REALTYPE :: mf_EmD, hEnt, hDet, dwp, apm, ap, rhs, zinv, relax + REALTYPE :: totalDepth, penal, del + integer :: i,im1 + logical :: found +! +!------------------------------------------------------------------------ +!BOC +! + mf_one_p_bCent = 1. + mf_bb*mf_Cent + zinvMin = -1.0*h(nlev) + mf_zinv = MIN( mf_zinv,zinvMin ) + izinv = -1./mf_zinv + bp = izinv*mf_bp + mf_d0 = izinv*mf_dbkg + relax = 0.6 + nrj = 0. + if(massflux_energy) nrj = 1. + totalDepth = -zi(0) + del = 1./(0.02*totalDepth) +! + a_p(0:nlev) = 0. + w_p(0:nlev) = -wp_min + Fmass(0:nlev) = 0. + +! surface BC for plume quantities + a_p(nlev) = mf_ap0 + w_p(nlev) = mf_wp0 + Fmass(nlev) = 0. + cff = 1./( h(nlev-1)+h(nlev) ) + cup = cff*(2.*h(nlev)+h(nlev-1)) + cdwn = cff*h(nlev) + t_p(nlev) = cup*T(nlev)-cdwn*T(nlev-1) ! extrapolation toward the surface + s_p(nlev) = cup*S(nlev)-cdwn*S(nlev-1) + u_p(nlev) = cup*u(nlev)-cdwn*u(nlev-1) + v_p(nlev) = cup*v(nlev)-cdwn*v(nlev-1) + +! +! Main loop for the computation of plume quantities (from top to bottom) +! +! found = .false. + + do i=nlev,1,-1 + +! 1- Initialize values at the top of the current grid cell + + app = a_p(i) + wpp = w_p(i) + wpp2 = wpp*wpp + penal = 0.5*( 1.+TANH( del*(zi(i-1)+0.9*totalDepth) ) ) + +! 2- Compute the buoyancy forcing term B_p + + dpth = - z(i) + B_p = - gravity*( get_rho(S_p(i),T_p(i),p=dpth) - rho(i) )/rho0 + +! 3- w_p equation + +! >>> iteration 1 + rhs_wp2 = bp*(wpp2+wpp2)+2.*mf_aa*B_p + if( rhs_wp2 < 0. ) then + cff = mf_one_p_bCent + else + cff = 1. + endif + + cff1 = 1./(cff+h(i)*bp) + wpm2 = cff1*((cff-h(i)*bp)*wpp2-mf_aa*h(i)*2.*B_p) + + +! >>> iteration 2 + rhs_wp2 = bp*(wpp2+wpm2)+2.*mf_aa*B_p + if( rhs_wp2 < 0. ) then + cff = mf_one_p_bCent + else + cff = 1. + dz = (wpp*wpp-wp_min*wp_min) / (2.*mf_aa*B_p+bp*(wpp*wpp+wp_min*wp_min)) + if( dz>0. .and. dz= 0.) then + hEnt = mf_Cent*dwp + hDet = -d0 + else + hEnt = 0. + hDet = -( mf_Cdet*dwp + d0 ) + endif + + mf_EmD = hEnt-hDet + cff = 1. / (2.*wpm+mf_EmD) + cff1 = app *(2.*wpp-mf_EmD) + apm = MAX(cff*cff1,0.) + +! 5- Plume tracer equations & Plume horizontal velocity equations + + ap = 0.5*(app+apm) + if( apm > 0. .and. wpm <-wp_min ) then + + rhs = ap*(hEnt*T(i)-hDet*T_p(i)) + t_p(i-1) = (app*wpp*T_p(i)-rhs)/(apm*wpm) + rhs = ap*(hEnt*S(i)-hDet*S_p(i)) + s_p(i-1) = (app*wpp*S_p(i)-rhs)/(apm*wpm) + im1 = MAX(i-1,1) + rhs = ap*(hEnt*u(i)-hDet*u_p(i)) + u_p(i-1) = (app*wpp*u_p(i)-rhs)/(apm*wpm) - mf_uv*(u(i)-u(im1)) + rhs = ap*(hEnt*v(i)-hDet*v_p(i)) + v_p(i-1) = (app*wpp*v_p(i)-rhs)/(apm*wpm) - mf_uv*(v(i)-v(im1)) + EmD(i ) = - 2.*(app*wpp-apm*wpm)/(h(i)*(app*wpp+apm*wpm)) + else + t_p(i-1) = t_p(i) + s_p(i-1) = s_p(i) + u_p(i-1) = u_p(i) + v_p(i-1) = v_p(i) + EmD(i ) = 0. + endif + + w_p(i-1) = wpm + a_p(i-1) = apm + Fmass(i-1) = -apm*wpm + Bmf (i ) = -nrj*Fmass(i)*B_p + +! if(found) exit + + end do + +! Compute Courant number to use substepping if necessary + + do i = 1,nlev-1 + Crt_number(i) = Fmass(i)*dt/h(i) + enddo + Crt_max = MAXVAL(Crt_number(:)) + mf_nsub = MAX( ceiling(Crt_max), 1 ) + + return + end subroutine massfluxeq +!EOC +!----------------------------------------------------------------------- +!BOP +! +! !ROUTINE: Computation of mass flux contribution to the turbulent transport of TKE \label{sec:massfluxeq_Witek11} +! +! !INTERFACE: + subroutine tke_transport_mf(nlev,u,v,h) +! !DESCRIPTION: +! This subroutine computes the term \eqref{TurbTransportTKE} to be added in the right hand side of the TKE equation. +! +! !USES: + use turbulence, only: Fmass,w_p,u_p,v_p + use turbulence, only: normVp,wk_mf + use turbulence, only: massflux_on_dynamics +! + IMPLICIT NONE +! +! !INPUT PARAMETERS: + +! number of vertical layers + integer, intent(in) :: nlev + +! horizontal velocity (m/s) + REALTYPE, intent(in) :: u(0:nlev),v(0:nlev) + +! layer thickness (m) + REALTYPE, intent(in) :: h(0:nlev) + +! !REVISION HISTORY: +! Original author(s): Florian Lemarié +!EOP +!------------------------------------------------------------------------ +! +! !LOCAL VARIABLES: + integer :: i + REALTYPE :: cff, FC(0:nlev) +! +!------------------------------------------------------------------------ +!BOC +! + FC(nlev)= 0. + cff = 0. + if( massflux_on_dynamics ) cff = 1. + do i = 1,nlev-1 + normVp(i) = 0.5*( w_p(i)*w_p(i) + cff*( u_p(i)-u(i) )**2 & + + cff*( v_p(i)-v(i) )**2 ) + FC(i) = Fmass(i )*normVp(i ) + enddo + ! + do i=1,nlev + wk_mf(i) = 2.*( FC(i+1) - FC(i ))/(h(i+1)+h(i)) + enddo + ! + return + end subroutine tke_transport_mf +!EOC + + + +!----------------------------------------------------------------------- +!BOP +! + +! !ROUTINE: Computation of a mass flux term on TKE \label{sec:massfluxeq_mfTKE} +! +! !INTERFACE + subroutine tkeeq_mf(nlev,dt,h) +! !DESCRIPTION: +! This subroutine computes the subplume TKE k_p from the equation +! \[ +! a_p w_p \partial_z k_p = E\left( (k-k_p) + \frac{1}{2} \| \mathbf{u}_p - \mathbf{u} \|^2 \right) - a_p \epsilon_p +! \] +! with $\epsilon_p = (c_\epsilon/l_\epsilon) k_p^{3/2}$. Then a mass flux term is added to the TKE equation +! \[ +! \partial_t k = - \partial_z \left( a_p w_p ( k_p - k ) \right) +! \] +! This term is added as a correction to the already computed $k$ from tkeeq subroutine \ref{sec:tkeeq}. +! +! +! !USES: + use turbulence, only: normVp,Fmass,a_p,w_p + use turbulence, only: tke,tkeo,tke_p,k_min,eps + use turbulence, only: mf_Cent,mf_Cdet, mf_d0, mf_nsub +! +! IMPLICIT NONE +! +! !INPUT PARAMETERS: +! +! number of vertical layers + integer, intent(in) :: nlev +! +! time step (s) + REALTYPE, intent(in) :: dt +!! +! layer thickness (m) + REALTYPE, intent(in) :: h(0:nlev) +! +! !REVISION HISTORY: +! Original author(s): Florian Lemarié +!EOP +!------------------------------------------------------------------------ +! +! !LOCAL VARIABLES: + REALTYPE, parameter :: wp_min = 0.5E-08 ! min value for plume velocity + REALTYPE :: dt_sub, tke_star(0:nlev) + REALTYPE :: ap,dissp,kpp,kpm,dwp,k_env,d0,rhs,FC(0:nlev) + integer :: i,isub +! +!------------------------------------------------------------------------ +!BOC +! +! 1- Compute subplume TKE tke_p + + tke_p(nlev-1) = tke(nlev-1) ! upper BC + + do i = nlev-1,1,-1 + kpp = a_p(i)*w_p(i)*tke_p(i) + ap = 0.5*( a_p(i ) + a_p(i-1) ) + dissp = eps(i)*tke_p(i)*SQRT(tke_p(i)) & + / ( tkeo (i)*SQRT(tkeo (i)) ) ! (ceps/leps) tke_p**(3/2) + dwp = w_p(i)-w_p(i-1) + ! Compute entrainment and detrainment rates multiplied by h + d0 = MIN( 0.5*h(i)*mf_d0*(w_p(i)+w_p(i-1)), -2.*wp_min ) + if( dwp >= 0.) then + hEnt = ap*mf_Cent*dwp + hDet = -ap*d0 + else + hEnt = 0. + hDet = -ap*( mf_Cdet*dwp + d0 ) + endif + ! + k_env = tke(i-1) + rhs = hEnt*k_env-hDet*tke_p(i) + rhs = rhs + hEnt*normVp(i)-ap*dissp + kpm = kpp - rhs + ! + if( a_p(i-1)>0. .and. w_p(i-1) < -wpmin ) then + tke_p(i-1) = MAX( kpm/(a_p(i-1)*w_p(i-1)), k_min ) + else + tke_p(i-1) = k_min + endif + enddo + + tke_p(nlev) = tke_p(nlev-1) + +! 2- Apply mass flux term to the TKE equation + dt_sub = dt / real(mf_nsub) + + do isub=1,mf_nsub + ! save old value + tke_star = tke + FC(1:nlev) = 0. + + do i = 1,nlev-1 + FC(i) = Fmass(i)*( tke_p(i)-tke_star(i-1) ) + enddo + + do i = 1,nlev-1 + tke(i) = tke_star(i) + dt_sub*(FC(i+1)-FC(i))*2./(h(i+1)+h(i)) + enddo + ! + enddo + +! clip at k_min + do i=0,nlev + tke(i) = max(tke(i),k_min) + enddo + + return + end subroutine tkeeq_mf +!EOC + +!----------------------------------------------------------------------- +! Copyright by the GOTM-team under the GNU Public License - www.gnu.org +!----------------------------------------------------------------------- From 403e30736db36133c8b655a7225e909b4e940cf0 Mon Sep 17 00:00:00 2001 From: lemarie Date: Fri, 13 Mar 2026 18:22:18 +0100 Subject: [PATCH 2/5] Mass flux scheme with energetic consistency --- src/gotm/gotm.F90 | 14 ++- src/meanflow/CMakeLists.txt | 1 + src/turbulence/CMakeLists.txt | 3 + src/turbulence/tkeeq.F90 | 17 ++- src/turbulence/turbulence.F90 | 222 ++++++++++++++++++++++++++++++++-- 5 files changed, 241 insertions(+), 16 deletions(-) diff --git a/src/gotm/gotm.F90 b/src/gotm/gotm.F90 index a10572d2..e75a00bc 100644 --- a/src/gotm/gotm.F90 +++ b/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) diff --git a/src/meanflow/CMakeLists.txt b/src/meanflow/CMakeLists.txt index f115df55..b307bece 100644 --- a/src/meanflow/CMakeLists.txt +++ b/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) diff --git a/src/turbulence/CMakeLists.txt b/src/turbulence/CMakeLists.txt index 7a317702..7e5a1288 100644 --- a/src/turbulence/CMakeLists.txt +++ b/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 diff --git a/src/turbulence/tkeeq.F90 b/src/turbulence/tkeeq.F90 index 82a95981..32d7c63c 100644 --- a/src/turbulence/tkeeq.F90 +++ b/src/turbulence/tkeeq.F90 @@ -16,16 +16,17 @@ subroutine tkeeq(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS) ! \label{tkeA} ! \dot{k} ! = -! {\cal D}_k + P + G + P_x + P_s - \epsilon +! {\cal D}_k + P + G + P_x + P_s + P_mf + G_mf - \epsilon ! \comma ! \end{equation} ! where $\dot{k}$ denotes the material derivative of $k$. $P$ and $G$ are ! the production of $k$ by mean shear and buoyancy, respectively, and ! $\epsilon$ the rate of dissipation. -! $P_s$ is Stokes shear production defined in \eq{computePs} -! and $P_x$ accounts for extra turbulence production. +! $P_s$ is Stokes shear production defined in \eq{computePs}, +! $P_x$ accounts for extra turbulence production, $P_{mf}$ and G_{mf} are energy +! transfers associated to convective plumes when a mass flux scheme is used. ! ${\cal D}_k$ represents the sum of -! the viscous and turbulent transport terms. +! the viscous and turbulent transport terms. ! For horizontally homogeneous flows, the transport term ${\cal D}_k$ ! appearing in \eq{tkeA} is presently expressed by a simple ! gradient formulation, @@ -35,7 +36,12 @@ subroutine tkeeq(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS) ! \comma ! \end{equation} ! where $\sigma_k$ is the constant Schmidt-number for $k$. -! +! When a mass flux scheme is used, additional terms are present in the +! turbulent fluxes of TKE (Perrot and Lemarié, 2025), in particular +! \begin{equation} +! \label{TurbTransportTKE} +! {\cal D}_k^{mf} = \frstder{z} \left( \frac{a_p w_p}{2} \| \mathbf{u}_{h,p} - \left< \mathbf{u}_h \right> \|^2 \right). +! \end{equation} ! In horizontally homogeneous flows, the shear and the buoyancy ! production, $P$ and $G$, can be written as ! \begin{equation} @@ -65,6 +71,7 @@ subroutine tkeeq(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS) use turbulence, only: tke,tkeo,k_min,eps use turbulence, only: k_bc, k_ubc, k_lbc, ubc_type, lbc_type use turbulence, only: sig_k + use turbulence, only: wk_mf,Pmf,Bmf use util, only: Dirichlet,Neumann IMPLICIT NONE diff --git a/src/turbulence/turbulence.F90 b/src/turbulence/turbulence.F90 index db32367f..b2e96817 100644 --- a/src/turbulence/turbulence.F90 +++ b/src/turbulence/turbulence.F90 @@ -36,7 +36,7 @@ module turbulence private ! ! !PUBLIC MEMBER FUNCTIONS: - public init_turbulence,post_init_turbulence,do_turbulence + public init_turbulence,post_init_turbulence,do_turbulence, do_massflux public k_bc,q2over2_bc,epsilon_bc,omega_bc,psi_bc,q2l_bc public clean_turbulence #ifdef _PRINTSTATE_ @@ -68,6 +68,9 @@ module turbulence ! Stokes production REALTYPE, public, dimension(:), allocatable :: PSTK +! Extra TKE production terms by convective plumes with mass flux scheme + REALTYPE, public, dimension(:), allocatable :: Pmf,Bmf + ! turbulent diffusivities ! of momentum, temperature, salinity REALTYPE, public, dimension(:), allocatable, target :: num @@ -79,6 +82,14 @@ module turbulence ! algebraic Reynolds stress and flux models. REALTYPE, public, dimension(:), allocatable :: nucl +! plumes properties for convective mass flux scheme + REALTYPE, public, dimension(:), allocatable :: a_p,w_p,Fmass + REALTYPE, public, dimension(:), allocatable :: T_p,S_p,wk_mf + REALTYPE, public, dimension(:), allocatable :: normVp,tke_p + REALTYPE, public, dimension(:), allocatable :: u_p,v_p,EmD + REALTYPE, public :: mf_zinv, mf_d0 + integer , public :: mf_nsub + ! non-local fluxes of momentum REALTYPE, public, dimension(:), allocatable :: gamu,gamv @@ -238,6 +249,20 @@ module turbulence REALTYPE, public :: numiw REALTYPE, public :: nuhiw REALTYPE, public :: numshear + +! the 'mfconvec' namelist + logical , public :: compute_massflux + logical , public :: massflux_on_dynamics + logical , public :: massflux_energy + REALTYPE, public :: mf_ap0 + REALTYPE, public :: mf_wp0 + REALTYPE, public :: mf_Cent + REALTYPE, public :: mf_Cdet + REALTYPE, public :: mf_aa + REALTYPE, public :: mf_bb + REALTYPE, public :: mf_bp + REALTYPE, public :: mf_uv + REALTYPE, public :: mf_dbkg ! ! !DEFINED PARAMETERS: @@ -256,6 +281,7 @@ module turbulence integer, parameter, public :: Constant=1 integer, parameter, public :: Munk_Anderson=2 integer, parameter, public :: Schumann_Gerz=3 + integer, parameter, public :: Prandtl_nemo=4 ! method to update length scale integer, parameter, public :: Parabolic=1 @@ -264,6 +290,7 @@ module turbulence integer, parameter, public :: Robert_Ouellet=4 integer, parameter, public :: Blackadar=5 integer, parameter, public :: Bougeault_Andre=6 + integer, parameter, public :: Gaspar_nemo=7 integer, parameter, public :: diss_eq=8 integer, parameter, public :: omega_eq=11 integer, parameter, public :: length_eq=9 @@ -393,6 +420,10 @@ subroutine init_turbulence_nml(namlst,fn,nlev) namelist /iw/ iw_model,alpha,klimiw,rich_cr, & numiw,nuhiw,numshear + + namelist /mfconvec/ compute_massflux,mf_ap0,mf_wp0,mf_Cent, & + mf_Cdet,mf_aa,mf_bb,mf_bp,mf_uv,mf_dbkg + !EOP !----------------------------------------------------------------------- !BOC @@ -531,6 +562,20 @@ subroutine init_turbulence_nml(namlst,fn,nlev) nuhiw=5.e-5 numshear=5.e-3 +! the 'mfconvec' namelist + compute_massflux=.false. + massflux_on_dynamics=.false. + massflux_energy=.false. + mf_ap0 = 0.2 + mf_wp0 = -0.5e-8 + mf_Cent = 0.99 + mf_Cdet = 1.99 + mf_aa = 1.0 + mf_bb = 1.0 + mf_bp = 0.75 + mf_uv = 0.5 + mf_dbkg = 1.125 + ! read the variables from the namelist file open(namlst,file=fn,status='old',action='read',err=80) read(namlst,nml=turbulence,err=81) @@ -550,6 +595,7 @@ subroutine init_turbulence_nml(namlst,fn,nlev) read(namlst,nml=my,err=87) read(namlst,nml=scnd,err=88) read(namlst,nml=iw,err=89) + read(namlst,nml=mfconvec,err=90) close (namlst) endif LEVEL2 'done.' @@ -578,6 +624,8 @@ subroutine init_turbulence_nml(namlst,fn,nlev) stop 'init_turbulence' 89 FATAL 'I could not read "iw" namelist' stop 'init_turbulence' +90 FATAL 'I could not read "mfconvec" namelist' + stop 'init_turbulence' end subroutine init_turbulence_nml !EOC @@ -632,11 +680,17 @@ subroutine init_turbulence_yaml(branch) call branch%get(turb_method, 'turb_method', 'turbulence closure', & options=(/option(no_model, 'constant nuh and num', 'no_model'), option(first_order, 'first-order', 'first_order'), option(second_order, 'second-order', 'second_order'), option(100, 'cvmix', 'cvmix')/),default=second_order) call branch%get(tke_method, 'tke_method', 'turbulent kinetic energy equation', & - options=(/option(tke_local_eq, 'algebraic length scale equation', 'local_eq'), option(tke_keps, 'differential equation for tke (k-eps or k-w style)', 'tke'), option(tke_MY, 'differential equation for q^2/2 (Mellor-Yamada style)', 'Mellor_Yamada')/),default=tke_keps) + options=(/option(tke_local_eq, 'algebraic length scale equation', 'local_eq'), option(tke_keps, 'differential equation for tke (k-eps, k-w, or 1eq k style)', 'tke'), & + option(tke_MY, 'differential equation for q^2/2 (Mellor-Yamada style)', 'Mellor_Yamada')/),default=tke_keps) call branch%get(len_scale_method, 'len_scale_method', 'dissipative length scale', & - options=(/option(parabolic, 'parabolic', 'parabolic'), option(triangular, 'triangular', 'triangular'), option(Xing_Davies, 'Xing and Davies (1995)', 'Xing_Davies'), option(Robert_Ouellet, 'Robert and Ouellet (1987)', 'Robert_Ouellet'), option(Blackadar, 'Blackadar (two boundaries) (1962)', 'Blackadar'), option(Bougeault_Andre, 'Bougeault and Andre (1986)', 'Bougeault_Andre'), option(diss_eq, 'dynamic dissipation rate equation', 'dissipation'), option(omega_eq, 'dynamic omega equation', 'omega'), option(length_eq, 'dynamic Mellor-Yamada q^2 l-equation', 'Mellor_Yamada'), option(generic_eq, 'generic length scale (GLS)', 'gls')/),default=diss_eq) + options=(/option(parabolic, 'parabolic', 'parabolic'), option(triangular, 'triangular', 'triangular'), & + option(Xing_Davies, 'Xing and Davies (1995)', 'Xing_Davies'), option(Robert_Ouellet, 'Robert and Ouellet (1987)', 'Robert_Ouellet'), & + option(Blackadar, 'Blackadar (two boundaries) (1962)', 'Blackadar'), option(Bougeault_Andre, 'Bougeault and Andre (1986)', 'Bougeault_Andre'), & + option(Gaspar_nemo, 'Gaspar (1990)', 'Gaspar_nemo'), & + option(diss_eq, 'dynamic dissipation rate equation', 'dissipation'), option(omega_eq, 'dynamic omega equation', 'omega'), option(length_eq, 'dynamic Mellor-Yamada q^2 l-equation', 'Mellor_Yamada'), option(generic_eq, 'generic length scale (GLS)', 'gls')/),default=diss_eq) call branch%get(stab_method, 'stab_method', 'stability functions', & - options=(/option(1, 'constant', 'constant'), option(Munk_Anderson, 'Munk and Anderson (1954)', 'Munk_Anderson'), option(Schumann_Gerz, 'Schumann and Gerz (1995)', 'Schumann_Gerz')/),default=1) + options=(/option(1, 'constant', 'constant'), option(Munk_Anderson, 'Munk and Anderson (1954)', 'Munk_Anderson'), & + option(Schumann_Gerz, 'Schumann and Gerz (1995)', 'Schumann_Gerz'), option(Prandtl_nemo, 'Madec et al. (2025)', 'Prandtl_nemo')/),default=1) twig => branch%get_child('bc', 'boundary conditions') call twig%get(k_ubc, 'k_ubc', 'upper boundary condition for k-equation', & @@ -823,6 +877,28 @@ subroutine init_turbulence_yaml(branch) default=1.e-4_rk) call twig%get(nuhiw, 'nuh', 'background diffusivity for internal wave breaking', 'm^2/s', & default=1.e-5_rk) + + twig => branch%get_child('MFconvec', 'mass flux convection scheme', display=display_advanced) + call twig%get(compute_massflux, 'compute_massflux', 'Compute convective mass flux term', default=.false.) + call twig%get(massflux_on_dynamics, 'massflux_on_dynamics', 'Apply mass flux to horizontal velocities', default=.false.) + call twig%get(massflux_energy, 'massflux_energy', 'Add contribution of mass flux terms to TKE equation', default=.false.) + call twig%get(mf_ap0, 'mf_ap0', 'surface bdy condition for plume area', '-', default=0.2_rk) + call twig%get(mf_wp0, 'mf_wp0', 'surface bdy condition for plume vertical velocity', 'm/s', & + default=-0.5e-8_rk) + call twig%get(mf_Cent, 'mf_Cent', 'coefficient for entrainment param', '-', & + default=0.99_rk) + call twig%get(mf_Cdet, 'mf_Cdet', 'coefficient for detrainment param', '-', & + default=1.99_rk) + call twig%get(mf_aa, 'mf_aa', 'a coefficient in wp eqn', '-', & + default=1.0_rk) + call twig%get(mf_bb, 'mf_bb', 'b coefficient in wp eqn', '-', & + default=1.0_rk) + call twig%get(mf_bp, 'mf_bp', 'b prime coeficient in wp eqn', '-', & + default=0.75_rk) + call twig%get(mf_uv, 'mf_uv', 'Cu coefficient in the up eqn', '-', & + default=0.5_rk) + call twig%get(mf_dbkg, 'mf_dbkg', 'background detrainment', '-', & + default=1.125_rk) LEVEL2 'done.' return @@ -953,6 +1029,16 @@ subroutine post_init_turbulence(nlev) allocate(PSTK(0:nlev),stat=rc) if (rc /= 0) stop 'init_turbulence: Error allocating (PSTK)' PSTK = _ZERO_ + + ! Mass flux production + allocate(Pmf(0:nlev),stat=rc) + if (rc /= 0) stop 'init_turbulence: Error allocating (Pmf)' + Pmf = _ZERO_ + + ! Mass flux production + allocate(Bmf(0:nlev),stat=rc) + if (rc /= 0) stop 'init_turbulence: Error allocating (Bmf)' + Bmf = _ZERO_ allocate(num(0:nlev),stat=rc) if (rc /= 0) stop 'init_turbulence: Error allocating (num)' @@ -970,6 +1056,54 @@ subroutine post_init_turbulence(nlev) if (rc /= 0) stop 'init_turbulence: Error allocating (nucl)' nucl = _ZERO_ + allocate(a_p(0:nlev),stat=rc) + if (rc /= 0) stop 'init_turbulence: Error allocating (a_p)' + a_p = _ZERO_ + + allocate(w_p(0:nlev),stat=rc) + if (rc /= 0) stop 'init_turbulence: Error allocating (w_p)' + w_p = _ZERO_ + + allocate(Fmass(0:nlev),stat=rc) + if (rc /= 0) stop 'init_turbulence: Error allocating (Fmass)' + Fmass = _ZERO_ + + allocate(u_p(0:nlev),stat=rc) + if (rc /= 0) stop 'init_turbulence: Error allocating (u_p)' + u_p = _ZERO_ + + allocate(v_p(0:nlev),stat=rc) + if (rc /= 0) stop 'init_turbulence: Error allocating (v_p)' + v_p = _ZERO_ + + allocate(T_p(0:nlev),stat=rc) + if (rc /= 0) stop 'init_turbulence: Error allocating (T_p)' + T_p = _ZERO_ + + allocate(S_p(0:nlev),stat=rc) + if (rc /= 0) stop 'init_turbulence: Error allocating (S_p)' + S_p = _ZERO_ + + allocate(normVp(0:nlev),stat=rc) + if (rc /= 0) stop 'init_turbulence: Error allocating (normVp)' + normVp = _ZERO_ + + allocate(wk_mf(0:nlev),stat=rc) + if (rc /= 0) stop 'init_turbulence: Error allocating (wk_mf)' + wk_mf = _ZERO_ + + allocate(tke_p(0:nlev),stat=rc) + if (rc /= 0) stop 'init_turbulence: Error allocating (tke_p)' + tke_p = _ZERO_ + + allocate(EmD(0:nlev),stat=rc) + if (rc /= 0) stop 'init_turbulence: Error allocating (EmD)' + EmD = _ZERO_ + + mf_zinv = _ZERO_ + mf_d0 = _ZERO_ + mf_nsub = 1 + allocate(gamu(0:nlev),stat=rc) if (rc /= 0) stop 'init_turbulence: Error allocating (gamu)' gamu = _ZERO_ @@ -2351,6 +2485,17 @@ subroutine report_model LEVEL3 'Schmidt-number for k, sig_k =', sig_k LEVEL3 'von Karman constant, kappa =', kappa LEVEL2 ' ' + case(Gaspar_nemo) ! Florian + LEVEL2 ' ' + LEVEL2 '--------------------------------------------------------' + LEVEL2 'You are using a one-equation model' + LEVEL2 'with the length-scale of Gaspar (1990)' + LEVEL2 'as implemented in the NEMO model (Madec et al.)' + LEVEL2 'The properties of this model are:' + LEVEL2 ' ' + LEVEL3 'cm =', cm0_fix + LEVEL3 'ceps =', 0.5*SQRT(2.) + LEVEL2 ' ' case(diss_eq) LEVEL2 ' ' LEVEL2 '--------------------------------------------------------' @@ -2614,16 +2759,16 @@ end subroutine alpha_mnb ! empirical stability function call production(nlev,NN,SS, xP=xP, SSCSTK=SSCSTK, SSSTK=SSSTK) - call alpha_mnb(nlev,NN,SS) call stabilityfunctions(nlev) call do_tke(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS) + if(compute_massflux .and. massflux_energy) then + call tkeeq_mf(nlev,dt,h) + endif call do_lengthscale(nlev,dt,depth,u_taus,u_taub,z0s,z0b,h,NN,SS) call kolpran(nlev) - call internal_wave(nlev,NN,SS) - case (second_order) ! solve a model for the second moments @@ -2740,9 +2885,9 @@ subroutine do_tke(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS) select case (tke_method) case(tke_local_eq) ! use algebraic length scale equation - call tkealgebraic(nlev,u_taus,u_taub,NN,SS) + call tkealgebraic(nlev,u_taus,u_taub,NN,SS) case(tke_keps) - ! use differential equation for tke (k-epsilon style) + ! use differential equation for tke (k-epsilon style) - works also for 1 equation case (with Gaspar_nemo mixing lengths) call tkeeq(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS) case(tke_MY) ! use differential equation for q^2/2 (Mellor-Yamada style) @@ -2861,6 +3006,8 @@ subroutine do_lengthscale(nlev,dt,depth,u_taus,u_taub, z0s,z0b,h,NN,SS) call lengthscaleeq(nlev,dt,depth,u_taus,u_taub,z0s,z0b,h,NN,SS) case(Bougeault_Andre) call potentialml(nlev,z0b,z0s,h,depth,NN) + case(Gaspar_nemo) + call lengthscale_gaspar(nlev,z0s,z0b,h,NN) case default call algebraiclength(len_scale_method,nlev,z0b,z0s,depth,h,NN) end select @@ -2871,6 +3018,55 @@ end subroutine do_lengthscale !EOC + +!----------------------------------------------------------------------- +!BOP +! +! !IROUTINE: Handle convection for first-order TKE schemes using a mass flux formalism\label{sec:convec} +! +! !INTERFACE: + subroutine do_massflux(nlev,dt,umean,vmean,Tmean,Smean,Rmean,h) +! +! !DESCRIPTION: +! Based on user input, this routine calls the appropriate routines for +! calculating the plume properties ($U_p$, $V_p$, $\Theta_p$, $S_p$, $w_p$, $a_p$) +! from a set of nonlinear ordinary differential equations. +! This is the basis for the representation of the effect of +! coherent structures with a mass flux formalism (Perrot and Lemarié, 2025). +! When the energetic consistency is required, an extra term ${\cal D}_k^{mf}$ +! is computed following \eqref{TurbTransportTKE}. +! +! !USES: + IMPLICIT NONE +! +! !INPUT PARAMETERS: + integer, intent(in) :: nlev + REALTYPE, intent(in) :: dt + REALTYPE, intent(in) :: umean(0:nlev) + REALTYPE, intent(in) :: vmean(0:nlev) + REALTYPE, intent(in) :: Tmean(0:nlev) + REALTYPE, intent(in) :: Smean(0:nlev) + REALTYPE, intent(in) :: Rmean(0:nlev) + REALTYPE, intent(in) :: h(0:nlev) +! +! !REVISION HISTORY: +! Original author(s): Florian Lemarié +! +!EOP +!----------------------------------------------------------------------- +!BOC + call massfluxeq(nlev,dt,umean,vmean,tmean,smean,Rmean,h) + if( massflux_energy ) call tke_transport_mf(nlev,umean,vmean,h) + return + end subroutine do_massflux +!----------------------------------------------------------------------- +!EOC + + + + + + !----------------------------------------------------------------------- !BOP ! @@ -3031,6 +3227,8 @@ subroutine stabilityfunctions(nlev) call cmue_ma(nlev) case(Schumann_Gerz) call cmue_sg(nlev) + case(Prandtl_nemo) + call cmue_nemo(nlev) case default STDERR '... not a valid stability function' @@ -3146,6 +3344,9 @@ subroutine compute_cm0(turb_method,stab_method,scnd_method) case(Schumann_Gerz) cm0 = cm0_fix cmsf = cm0_fix + case(Prandtl_nemo) + cm0 = 0.1 + cmsf = 0.1 case default STDERR '... not a valid stability function to compute cm0' @@ -4162,6 +4363,9 @@ subroutine print_state_turbulence() LEVEL2 'iw namelist', iw_model,alpha,klimiw,rich_cr, & numiw,nuhiw,numshear + LEVEL2 'mfconvec namelist', compute_massflux,massflux_on_dynamics,massflux_energy, & + mf_ap0,mf_wp0,mf_Cent,mf_Cdet,mf_aa,mf_bb,mf_bp,mf_uv, mf_dbkg + end subroutine print_state_turbulence !EOC #endif From 4a12459e8ed34e086ae0a27cc26118f152bae558 Mon Sep 17 00:00:00 2001 From: lemarie Date: Tue, 17 Mar 2026 14:16:23 +0100 Subject: [PATCH 3/5] Add extra terms for energetic consistency in TKE equation (EDMF case) --- src/turbulence/tkeeq.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/turbulence/tkeeq.F90 b/src/turbulence/tkeeq.F90 index 32d7c63c..ca48d992 100644 --- a/src/turbulence/tkeeq.F90 +++ b/src/turbulence/tkeeq.F90 @@ -129,7 +129,7 @@ subroutine tkeeq(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS) ! compute production terms in k-equation prod = P(i) + Px(i) + PSTK(i) - buoyan = B(i) + buoyan = B(i) + Bmf(i) + Pmf(i) + wk_mf(i) ! the contributions Bmf, Pmf, and wk_mf are not positive-definite that's the reason why they are added to buoyan diss = eps(i) From 48c596611ae19b059fda7a6b71791beafa34ffec Mon Sep 17 00:00:00 2001 From: lemarie Date: Tue, 17 Mar 2026 14:16:59 +0100 Subject: [PATCH 4/5] Add extra variables in output files for EDMF case --- src/gotm/register_all_variables.F90 | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/src/gotm/register_all_variables.F90 b/src/gotm/register_all_variables.F90 index 58aa79c5..9bc06d00 100644 --- a/src/gotm/register_all_variables.F90 +++ b/src/gotm/register_all_variables.F90 @@ -549,10 +549,24 @@ 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_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_ From d57059c39890042bf8387c9481dac31cb6c58848 Mon Sep 17 00:00:00 2001 From: lemarie Date: Tue, 17 Mar 2026 15:03:45 +0100 Subject: [PATCH 5/5] Add extra variables in output files for EDMF case --- src/gotm/register_all_variables.F90 | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/gotm/register_all_variables.F90 b/src/gotm/register_all_variables.F90 index 9bc06d00..0b56dde1 100644 --- a/src/gotm/register_all_variables.F90 +++ b/src/gotm/register_all_variables.F90 @@ -556,6 +556,10 @@ subroutine register_turbulence_variables(nlev) 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')