xP: support independent weighting of xP in lengthscale equations

src/turbulence/cmue_a.F90

@@ -85,7 +85,7 @@ subroutine cmue_a(nlev)
 !
 ! !USES:
    use turbulence, only: eps
-   use turbulence, only: P,B,Pb,epsb
+   use turbulence, only: P,B,Px,Pb,epsb
    use turbulence, only: an,as,at,r
    use turbulence, only: cmue1,cmue2,gam
    use turbulence, only: cm0
@@ -158,7 +158,7 @@ subroutine cmue_a(nlev)
 
      do i=1,nlev-1
 
-        Pe   =   ( P(i) + B(i) )/eps(i)
+        Pe   =   ( P(i) + Px(i) + B(i) )/eps(i)
         Pbeb =   Pb(i)/epsb(i)
         r_i  =   1./r(i)
 

src/turbulence/dissipationeq.F90

@@ -17,11 +17,15 @@ subroutine dissipationeq(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS)
 !   =
 !   {\cal D}_\epsilon
 !   + \frac{\epsilon}{k} ( c_{\epsilon 1} P + c_{\epsilon 3} G
+!                        + c_{\epsilon x} P_x
+!                        + c_{\epsilon 4} P_s
 !                        - c_{\epsilon 2} \epsilon )
 !   \comma
 ! \end{equation}
 ! where $\dot{\epsilon}$ denotes the material derivative of $\epsilon$.
-! The production terms $P$ and $G$ follow from \eq{PandG} and
+! The production terms $P$ and $G$ follow from \eq{PandG}.
+! $P_s$ is Stokes shear production defined in \eq{computePs}
+! and $P_x$ accounts for extra turbulence production.
 ! ${\cal D}_\epsilon$ represents the sum of the viscous and turbulent
 ! transport terms.
 !
@@ -62,9 +66,9 @@ subroutine dissipationeq(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS)
 ! by setting {\tt length\_lim = .true.} in {\tt gotm.yaml}.
 !
 ! !USES:
-   use turbulence, only: P,B,PSTK,num
+   use turbulence, only: P,B,Px,PSTK,num
    use turbulence, only: tke,tkeo,k_min,eps,eps_min,L
-   use turbulence, only: ce1,ce2,ce3plus,ce3minus,ce4
+   use turbulence, only: ce1,ce2,ce3plus,ce3minus,cex,ce4
    use turbulence, only: cm0,cde,galp,length_lim
    use turbulence, only: epsilon_bc, psi_ubc, psi_lbc, ubc_type, lbc_type
    use turbulence, only: sig_e,sig_e0,sig_peps
@@ -122,7 +126,7 @@ subroutine dissipationeq(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS)
    if (sig_peps) then          ! With wave breaking
       sig_eff(nlev)=sig_e0
       do i=1,nlev-1
-         peps=(P(i)+B(i))/eps(i)
+         peps=(P(i)+Px(i)+B(i))/eps(i)
          if (peps .gt. 1.) peps=_ONE_
          sig_eff(i)=peps*sig_e+(_ONE_-peps)*sig_e0
       end do
@@ -147,7 +151,7 @@ subroutine dissipationeq(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS)
       end if
 
       EpsOverTke  = eps(i)/tkeo(i)
-      prod        = ce1*EpsOverTke*P(i) + ce4*EpsOverTke*PSTK(i)
+      prod        = EpsOverTke * ( ce1*P(i) + cex*Px(i) + ce4*PSTK(i) )
       buoyan      = ce3*EpsOverTke*B(i)
       diss        = ce2*EpsOverTke*eps(i)
 

src/turbulence/genericeq.F90

@@ -43,12 +43,16 @@ subroutine genericeq(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS)
 !   \label{generic}
 !   \dot{\psi} = {\cal D}_\psi
 !   + \frac{\psi}{k} (  c_{\psi_1} P + c_{\psi_3} G
+!                     + c_{\psi x} P_x
+!                     + c_{\psi 4} P_s
 !    - c_{\psi 2} \epsilon )
 !   \comma
 ! \end{equation}
 ! where $\dot{\psi}$ denotes the material derivative of $\psi$,
 ! see \cite{UmlaufBurchard2003}.
 ! The production terms $P$ and $G$ follow from \eq{PandG}.
+! $P_s$ is Stokes shear production defined in \eq{computePs}
+! and $P_x$ accounts for extra turbulence production.
 ! ${\cal D}_\psi$ represents the sum of the viscous and turbulent
 ! transport terms. The rate of dissipation can computed by solving
 ! \eq{psi_l} for $l$ and inserting the result into \eq{epsilon}.
@@ -118,9 +122,9 @@ subroutine genericeq(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS)
 ! implemented in GOTM and are described in \sect{sec:generate}.
 !
 ! !USES:
-   use turbulence, only: P,B,PSTK,num
+   use turbulence, only: P,B,Px,PSTK,num
    use turbulence, only: tke,tkeo,k_min,eps,eps_min,L
-   use turbulence, only: cpsi1,cpsi2,cpsi3plus,cpsi3minus,cpsi4,sig_psi
+   use turbulence, only: cpsi1,cpsi2,cpsi3plus,cpsi3minus,cpsix,cpsi4,sig_psi
    use turbulence, only: gen_m,gen_n,gen_p
    use turbulence, only: cm0,cde,galp,length_lim
    use turbulence, only: psi_bc, psi_ubc, psi_lbc, ubc_type, lbc_type
@@ -199,7 +203,7 @@ subroutine genericeq(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS)
 
 !     compute production terms in psi-equation
       PsiOverTke  = psi(i)/tkeo(i)
-      prod        = cpsi1*PsiOverTke*P(i) + cpsi4*PsiOverTke*PSTK(i)
+      prod        = PsiOverTke * ( cpsi1*P(i) + cpsix*Px(i) + cpsi4*PSTK(i) )
       buoyan      = cpsi3*PsiOverTke*B(i)
       diss        = cpsi2*PsiOverTke*eps(i)
 

src/turbulence/lengthscaleeq.F90

@@ -14,11 +14,14 @@ subroutine lengthscaleeq(nlev,dt,depth,u_taus,u_taub,z0s,z0b,h,NN,SS)
 ! \begin{equation}
 !   \label{MY}
 !   \dot{\overline{q^2 l}}
-!   = {\cal D}_l + l ( E_1  P + E_3 G - E_2  F \epsilon )
+!   = {\cal D}_l + l ( E_1  P + E_3 G + E_x P_x + E_6 P_s - E_2  F \epsilon )
 !   \comma
 ! \end{equation}
 ! where $\dot{\overline{q^2 l}}$ denotes the material derivative of $q^2 l$.
-! The production terms $P$ and $G$ follow from \eq{PandG}, and $\epsilon$
+! The production terms $P$ and $G$ follow from \eq{PandG}.
+! $P_s$ is Stokes shear production defined in \eq{computePs}
+! and $P_x$ accounts for extra turbulence production.
+! $\epsilon$
 ! can be computed either directly from \eq{epsilonMY}, or from \eq{epsilon}
 ! with the help \eq{B1}.
 !
@@ -67,9 +70,9 @@ subroutine lengthscaleeq(nlev,dt,depth,u_taus,u_taub,z0s,z0b,h,NN,SS)
 ! by setting {\tt length\_lim = .true.} in {\tt gotm.yaml}.
 !
 ! !USES:
-   use turbulence, only: P,B, PSTK
+   use turbulence, only: P,B,Px,PSTK
    use turbulence, only: tke,tkeo,k_min,eps,eps_min,L
-   use turbulence, only: kappa,e1,e2,e3,e6,b1
+   use turbulence, only: kappa,e1,e2,e3,ex,e6,b1
    use turbulence, only: MY_length,cm0,cde,galp,length_lim
    use turbulence, only: q2l_bc, psi_ubc, psi_lbc, ubc_type, lbc_type
    use turbulence, only: sl_var
@@ -159,7 +162,7 @@ subroutine lengthscaleeq(nlev,dt,depth,u_taus,u_taub,z0s,z0b,h,NN,SS)
       avh(i)      =  sl_var(i) * sqrt(2.*tkeo(i))*L(i)
 
 !     compute production terms in q^2 l - equation
-      prod        =  e1*L(i)*P(i) + e6*L(i)*PSTK(i)
+      prod        =  L(i) * ( e1*P(i) + ex*Px(i) + e6*PSTK(i) )
       buoyan      =  e3*L(i)*B(i)
       diss        =  q3(i)/b1*(1.+e2*(L(i)/Lz(i))*(L(i)/Lz(i)))
 

src/turbulence/omegaeq.F90

@@ -17,11 +17,15 @@ subroutine omegaeq(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS)
 !   =
 !   {\cal D}_\omega
 !   + \frac{\omega}{k} ( c_{\omega 1} P + c_{\omega 3} G
+!                        + c_{\omega x} P_x
+!                        + c_{\omega 4} P_s
 !                        - c_{\omega 2} \varepsilon )
 !   \comma
 ! \end{equation}
 ! where $\dot{\omega}$ denotes the material derivative of $\omega$.
-! The production terms $P$ and $G$ follow from \eq{PandG} and
+! The production terms $P$ and $G$ follow from \eq{PandG}.
+! $P_s$ is Stokes shear production defined in \eq{computePs}
+! and $P_x$ accounts for extra turbulence production.
 ! ${\cal D}_\omega$ represents the sum of the viscous and turbulent
 ! transport terms.
 !
@@ -57,9 +61,9 @@ subroutine omegaeq(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS)
 ! by setting {\tt length\_lim = .true.} in {\tt gotm.yaml}.     
 !
 ! !USES:
-   use turbulence, only: P,B,PSTK,num
+   use turbulence, only: P,B,Px,PSTK,num
    use turbulence, only: tke,tkeo,k_min,eps,eps_min,L
-   use turbulence, only: cw1,cw2,cw3plus,cw3minus,cw4
+   use turbulence, only: cw1,cw2,cw3plus,cw3minus,cwx,cw4
    use turbulence, only: cm0,cde,galp,length_lim
    use turbulence, only: omega_bc, psi_ubc, psi_lbc, ubc_type, lbc_type
    use turbulence, only: sig_w
@@ -130,7 +134,7 @@ subroutine omegaeq(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS)
       end if
 
       OmgOverTke  = omega(i)/tkeo(i)
-      prod        = cw1*OmgOverTke*P(i) + cw4*OmgOverTke*PSTK(i)
+      prod        = OmgOverTke * ( cw1*P(i) + cwx*Px(i) + cw4*PSTK(i) )
       buoyan      = cw3*OmgOverTke*B(i)
       diss        = cw2*OmgOverTke*eps(i)
 

src/turbulence/production.F90

@@ -9,20 +9,32 @@ subroutine production(nlev,NN,SS,xP, SSCSTK, SSSTK)
 !
 ! !DESCRIPTION:
 !  This subroutine calculates the production terms of turbulent kinetic
-!  energy as defined in \eq{PandG} and the production of buoayancy
+!  energy as defined in \eq{PandG} and the production of buoyancy
 !  variance as defined in \eq{Pbvertical}.
-!  The shear-production is computed according to
+!  The Eulerian shear-production is computed according to
 !  \begin{equation}
 !    \label{computeP}
-!     P = \nu_t (M^2 + \alpha_w N^2) + X_P
+!     P = \nu_t (M^2 + \alpha_w N^2) + \nu^S_t S_c^2
 !    \comma
 !  \end{equation}
 !  with the turbulent diffusivity of momentum, $\nu_t$, defined in
 !  \eq{nu}. The shear-frequency, $M$, is discretised as described
 !  in \sect{sec:shear}.
 !   The term multiplied by $\alpha_w$ traces back to
 !  a parameterisation of breaking internal waves suggested by
-!  \cite{Mellor89}. $X_P$ is an extra production term, connected for
+!  \cite{Mellor89}.
+!  The turbulent momentum fluxes due to Stokes velocities induce the
+!  Stokes-Eulerian cross-shear term
+!  $S_c^2 = \frac{\partial u}{\partial z}\frac{\partial u_s}{\partial z} + \frac{\partial v}{\partial z}\frac{\partial v_s}{\partial z}$
+!  with corresponding diffusivity $\nu^S_t$, and the additional
+!  Stokes shear-production
+!  \begin{equation}
+!    \label{computePs}
+!     P_s = \nu_t S_c^2 + \nu^S_t S_s^2
+!  \end{equation}
+!  with squared Stokes shear
+!  $S_s^2 = \frac{\partial u_s}{\partial z}^2 + \frac{\partial v_s}{\partial z}^2$.
+!  $X_P$ is an extra production term, connected for
 !  example with turbulence production caused by sea-grass, see
 !  \eq{sgProduction} in  \sect{sec:seagrass}. {\tt xP} is an {\tt optional}
 !  argument in the FORTRAN code.
@@ -51,7 +63,7 @@ subroutine production(nlev,NN,SS,xP, SSCSTK, SSSTK)
 !  kinetic and potential energy in \eq{tkeA} and \eq{kbeq}, respectively.
 !
 ! !USES:
-   use turbulence, only: P,B,Pb, PSTK
+   use turbulence, only: P,B,Pb,Px,PSTK
    use turbulence, only: num,nuh, nucl
    use turbulence, only: alpha,iw_model
    IMPLICIT NONE
@@ -92,19 +104,17 @@ subroutine production(nlev,NN,SS,xP, SSCSTK, SSSTK)
       alpha_eff=alpha
    end if
 
+   do i=0,nlev
+      P(i)    =  num(i)*( SS(i)+alpha_eff*NN(i) )
+      B(i)    = -nuh(i)*NN(i)
+      Pb(i)   = -  B(i)*NN(i)
+   end do
+
    if ( PRESENT(xP) ) then
       do i=0,nlev
-         P(i)    =  num(i)*( SS(i)+alpha_eff*NN(i) ) + xP(i)
-         B(i)    = -nuh(i)*NN(i)
-         Pb(i)   = -  B(i)*NN(i)
-      enddo
-   else
-      do i=0,nlev
-         P(i)    =  num(i)*( SS(i)+alpha_eff*NN(i) )
-         B(i)    = -nuh(i)*NN(i)
-         Pb(i)   = -  B(i)*NN(i)
-      enddo
-   endif
+         Px(i) = xP(i)
+      end do
+   end if
 
 !  P is -<u'w'>du/dz for q2l production with e1,
 !  PSTK is -<u'w'>dus/dz for q2l prodcution with e6,

src/turbulence/q2over2eq.F90

@@ -13,7 +13,7 @@ subroutine q2over2eq(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS)
 !   \label{tkeB}
 !   \dot{\overline{q^2/2}}
 !   =
-!   {\cal D}_q +  P + G  - \epsilon
+!   {\cal D}_q +  P + G + P_x + P_s - \epsilon
 !   \comma
 ! \end{equation}
 ! where $\dot{\overline{q^2/2}}$ denotes the material derivative of $q^2/2$.
@@ -44,7 +44,7 @@ subroutine q2over2eq(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS)
 ! \end{equation}
 !
 ! !USES:
-   use turbulence,   only: P,B, PSTK
+   use turbulence,   only: P,B,Px,PSTK
    use turbulence,   only: tke,tkeo,k_min,eps,L
    use turbulence,   only: q2over2_bc, k_ubc, k_lbc, ubc_type, lbc_type
    use turbulence,   only: sq_var
@@ -103,7 +103,7 @@ subroutine q2over2eq(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS)
       avh(i) = sq_var(i) * sqrt( 2.*tke(i) )*L(i)
 
 !     compute production terms in q^2/2-equation
-      prod     = P(i) + PSTK(i)
+      prod     = P(i) + Px(i) + PSTK(i)
       buoyan   = B(i)
       diss     = eps(i)
 

src/turbulence/tkeeq.F90

@@ -16,12 +16,15 @@ subroutine tkeeq(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS)
 !   \label{tkeA}
 !   \dot{k}
 !   =
-!   {\cal D}_k +  P + G  - \epsilon
+!   {\cal D}_k +  P + G + P_x + P_s - \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. ${\cal D}_k$ represents the sum of
+! $\epsilon$ the rate of dissipation.
+! $P_s$ is Stokes shear production defined in \eq{computePs}
+! and $P_x$ accounts for extra turbulence production.
+! ${\cal D}_k$ represents the sum of
 ! 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
@@ -58,7 +61,7 @@ subroutine tkeeq(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS)
 ! where $c_\mu^0$ is a constant of the model.
 !
 ! !USES:
-   use turbulence,   only: P,B,PSTK,num
+   use turbulence,   only: P,B,Px,PSTK,num
    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
@@ -118,7 +121,7 @@ subroutine tkeeq(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS)
       avh(i) = num(i)/sig_k
 
 !     compute production terms in k-equation
-      prod     = P(i) + PSTK(i)
+      prod     = P(i) + Px(i) + PSTK(i)
       buoyan   = B(i)
       diss     = eps(i)