diff --git a/src/turbulence/cmue_a.F90 b/src/turbulence/cmue_a.F90 index 79471980..32ba0fb5 100644 --- a/src/turbulence/cmue_a.F90 +++ b/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) diff --git a/src/turbulence/dissipationeq.F90 b/src/turbulence/dissipationeq.F90 index ae61a168..9de73ea1 100644 --- a/src/turbulence/dissipationeq.F90 +++ b/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) diff --git a/src/turbulence/genericeq.F90 b/src/turbulence/genericeq.F90 index 3ebb04c8..e75a829c 100644 --- a/src/turbulence/genericeq.F90 +++ b/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) diff --git a/src/turbulence/lengthscaleeq.F90 b/src/turbulence/lengthscaleeq.F90 index 0ac79236..2ef19488 100644 --- a/src/turbulence/lengthscaleeq.F90 +++ b/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))) diff --git a/src/turbulence/omegaeq.F90 b/src/turbulence/omegaeq.F90 index 0be3fb6e..05a49553 100644 --- a/src/turbulence/omegaeq.F90 +++ b/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) diff --git a/src/turbulence/production.F90 b/src/turbulence/production.F90 index 0a0cfc17..68935468 100644 --- a/src/turbulence/production.F90 +++ b/src/turbulence/production.F90 @@ -9,12 +9,12 @@ 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 @@ -22,7 +22,19 @@ subroutine production(nlev,NN,SS,xP, SSCSTK, SSSTK) ! 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 -du/dz for q2l production with e1, ! PSTK is -dus/dz for q2l prodcution with e6, diff --git a/src/turbulence/q2over2eq.F90 b/src/turbulence/q2over2eq.F90 index 40fa2b3a..4244b481 100644 --- a/src/turbulence/q2over2eq.F90 +++ b/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) diff --git a/src/turbulence/tkeeq.F90 b/src/turbulence/tkeeq.F90 index 917f52ed..82a95981 100644 --- a/src/turbulence/tkeeq.F90 +++ b/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) diff --git a/src/turbulence/turbulence.F90 b/src/turbulence/turbulence.F90 index c29ce1f0..3e9a320a 100644 --- a/src/turbulence/turbulence.F90 +++ b/src/turbulence/turbulence.F90 @@ -62,6 +62,9 @@ module turbulence ! of tke and buoyancy variance REALTYPE, public, dimension(:), allocatable :: P,B,Pb +! extra production term + REALTYPE, public, dimension(:), allocatable :: Px + ! Stokes production REALTYPE, public, dimension(:), allocatable :: PSTK @@ -172,7 +175,8 @@ module turbulence REALTYPE, public :: cpsi2 REALTYPE, public :: cpsi3minus REALTYPE, public :: cpsi3plus - REALTYPE, public :: cpsi4 + REALTYPE, public :: cpsix + REALTYPE, public :: cpsi4 REALTYPE :: sig_kpsi REALTYPE, public :: sig_psi REALTYPE :: gen_d @@ -184,7 +188,8 @@ module turbulence REALTYPE, public :: ce2 REALTYPE, public :: ce3minus REALTYPE, public :: ce3plus - REALTYPE, public :: ce4 + REALTYPE, public :: cex + REALTYPE, public :: ce4 REALTYPE, public :: sig_k REALTYPE, public :: sig_e logical, public :: sig_peps @@ -194,6 +199,7 @@ module turbulence REALTYPE, public :: cw2 REALTYPE, public :: cw3minus REALTYPE, public :: cw3plus + REALTYPE, public :: cwx REALTYPE, public :: cw4 REALTYPE, public :: sig_kw REALTYPE, public :: sig_w @@ -202,6 +208,7 @@ module turbulence REALTYPE, public :: e1 REALTYPE, public :: e2 REALTYPE, public :: e3 + REALTYPE, public :: ex REALTYPE, public :: e6 REALTYPE, public :: sq REALTYPE, public :: sl @@ -367,17 +374,17 @@ subroutine init_turbulence_nml(namlst,fn,nlev) kb_min,epsb_min namelist /generic/ compute_param,gen_m,gen_n,gen_p, & - cpsi1,cpsi2,cpsi3minus,cpsi3plus,cpsi4, & + cpsi1,cpsi2,cpsi3minus,cpsi3plus,cpsix,cpsi4, & sig_kpsi,sig_psi, & gen_d,gen_alpha,gen_l - namelist /keps/ ce1,ce2,ce3minus,ce3plus,ce4, & + namelist /keps/ ce1,ce2,ce3minus,ce3plus,cex,ce4, & sig_k,sig_e,sig_peps - namelist /kw/ cw1,cw2,cw3minus,cw3plus,cw4, & + namelist /kw/ cw1,cw2,cw3minus,cw3plus,cwx,cw4, & sig_kw,sig_w - namelist /my/ e1,e2,e3,e6,sq,sl,my_length, new_constr + namelist /my/ e1,e2,e3,ex,e6,sq,sl,my_length, new_constr namelist /scnd/ scnd_method,kb_method,epsb_method, & scnd_coeff, & @@ -458,6 +465,7 @@ subroutine init_turbulence_nml(namlst,fn,nlev) cpsi2=1.92 cpsi3minus=0.0 cpsi3plus=1.0 + cpsix=cpsi1 cpsi4=_ZERO_ sig_kpsi=1.0 sig_psi=1.3 @@ -470,6 +478,7 @@ subroutine init_turbulence_nml(namlst,fn,nlev) ce2=1.92 ce3minus=0.0 ce3plus=1.5 + cex=ce1 ce4=_ZERO_ sig_k=1.0 sig_e=1.3 @@ -480,6 +489,7 @@ subroutine init_turbulence_nml(namlst,fn,nlev) cw2=0.833 cw3minus=0.0 cw3plus=0.5 + cwx=cw1 cw4=0.15 sig_kw=2.0 sig_w=2.0 @@ -488,6 +498,7 @@ subroutine init_turbulence_nml(namlst,fn,nlev) e1=1.8 e2=1.33 e3=1.8 + ex=e1 e6=4.0 sq=0.2 sl=0.2 @@ -691,6 +702,8 @@ subroutine init_turbulence_yaml(branch) default=0.0_rk) call twig%get(cpsi3plus, 'cpsi3plus', 'cpsi3 for unstable stratification', '-', & default=1._rk) + call twig%get(cpsix, 'cpsix', 'empirical coefficient cpsix in psi equation', '-', & + default=cpsi1) call twig%get(cpsi4, 'cpsi4', 'empirical coefficient cpsi4 in psi equation', '-', & default=0.00_rk) call twig%get(sig_kpsi, 'sig_kpsi', 'Schmidt number for TKE diffusivity', '-', & @@ -713,6 +726,8 @@ subroutine init_turbulence_yaml(branch) default=0._rk) call twig%get(ce3plus, 'ce3plus', 'ce3 for unstable stratification', '-', & default=1.5_rk) + call twig%get(cex, 'cex', 'empirical coefficient cex in dissipation equation', '-', & + default=ce1) call twig%get(ce4, 'ce4', 'empirical coefficient ce4 in dissipation equation', '-', & default=0.0_rk) call twig%get(sig_k, 'sig_k', 'Schmidt number for TKE diffusivity', '-', & @@ -731,6 +746,9 @@ subroutine init_turbulence_yaml(branch) default=0.0_rk) call twig%get(cw3plus, 'cw3plus', 'cw3 for unstable stratification', '-', & default=0.5_rk) + call twig%get(cwx, 'cwx', 'empirical coefficient cwx in omega equation', '-', & + default=cw1) + call twig%get(cw4, 'cw4', 'empirical coefficient cw4 in omega equation', '-', & default=0.15_rk) call twig%get(sig_kw, 'sig_kw', 'Schmidt number for TKE diffusivity', '-', & @@ -745,6 +763,8 @@ subroutine init_turbulence_yaml(branch) default=1.33_rk) call twig%get(e3, 'e3', 'coefficient e3 in q^2 l equation', '-', & default=1.8_rk) + call twig%get(ex, 'ex', 'coefficient ex in q^2 l equation', '-', & + default=e1) call twig%get(e6, 'e6', 'coefficient e6 in q^2 l equation', '-', & default=4.0_rk) call twig%get(sq, 'sq', 'turbulent diffusivities of q^2 (= 2k)', '-', & @@ -925,6 +945,10 @@ subroutine post_init_turbulence(nlev) if (rc /= 0) stop 'init_turbulence: Error allocating (Pb)' Pb = _ZERO_ + allocate(Px(0:nlev),stat=rc) + if (rc /= 0) stop 'init_turbulence: Error allocating (Px)' + Px = _ZERO_ + ! Stokes production allocate(PSTK(0:nlev),stat=rc) if (rc /= 0) stop 'init_turbulence: Error allocating (PSTK)' @@ -2337,7 +2361,8 @@ subroutine report_model LEVEL3 'ce2 =', ce2 LEVEL3 'ce3minus =', ce3minus LEVEL3 'ce3plus =', ce3plus - LEVEL3 'ce4 =', ce4 + LEVEL3 'cex =', cex + LEVEL3 'ce4 =', ce4 LEVEL3 'sig_k =', sig_k LEVEL3 'sig_e =', sig_e LEVEL2 ' ' @@ -2362,7 +2387,8 @@ subroutine report_model LEVEL3 'cw2 =', cw2 LEVEL3 'cw3minus =', cw3minus LEVEL3 'cw3plus =', cw3plus - LEVEL3 'cw4 =', cw4 + LEVEL3 'cwx =', cwx + LEVEL3 'cw4 =', cw4 LEVEL3 'sig_k =', sig_k LEVEL3 'sig_w =', sig_w LEVEL2 ' ' @@ -2376,7 +2402,7 @@ subroutine report_model LEVEL3 'length-scale slope (no shear), L =', gen_l LEVEL3 'steady-state Richardson-number, Ri_st=', ri_st LEVEL2 '--------------------------------------------------------' - LEVEL2 ' ' + LEVEL2 ' ' case(length_eq) LEVEL2 ' ' LEVEL2 '--------------------------------------------------------' @@ -2387,6 +2413,7 @@ subroutine report_model LEVEL3 'E1 =', e1 LEVEL3 'E2 =', e2 LEVEL3 'E3 =', e3 + LEVEL3 'Ex =', ex LEVEL3 'E6 =', e6 LEVEL3 'Sq =', sq LEVEL3 'Sl =', sl @@ -2417,7 +2444,8 @@ subroutine report_model LEVEL3 'cpsi2 =', cpsi2 LEVEL3 'cpsi3minus =', cpsi3minus LEVEL3 'cpsi3plus =', cpsi3plus - LEVEL3 'cpsi4 =', cpsi4 + LEVEL3 'cpsix =', cpsix + LEVEL3 'cpsi4 =', cpsi4 LEVEL3 'sig_k =', sig_kpsi LEVEL3 'sig_psi =', sig_psi LEVEL2 ' ' @@ -2585,13 +2613,7 @@ end subroutine alpha_mnb ! solve a model for tke, length scale ! empirical stability function - if ( PRESENT(xP) ) then -! with seagrass turbulence - call production(nlev,NN,SS,xP, SSCSTK=SSCSTK, SSSTK=SSSTK) - else -! without - call production(nlev,NN,SS, SSCSTK=SSCSTK, SSSTK=SSSTK) - end if + call production(nlev,NN,SS, xP=xP, SSCSTK=SSCSTK, SSSTK=SSSTK) call alpha_mnb(nlev,NN,SS) call stabilityfunctions(nlev) @@ -2606,13 +2628,7 @@ end subroutine alpha_mnb ! solve a model for the second moments - if ( PRESENT(xP) ) then -! with seagrass turbulence - call production(nlev,NN,SS,xP, SSCSTK=SSCSTK, SSSTK=SSSTK) - else -! without - call production(nlev,NN,SS, SSCSTK=SSCSTK, SSSTK=SSSTK) - end if + call production(nlev,NN,SS, xP=xP, SSCSTK=SSCSTK, SSSTK=SSSTK) select case(scnd_method) case (quasi_Eq) @@ -3843,7 +3859,7 @@ REALTYPE function psi_bc(bc,type,zi,ki,z0,u_tau) psi_bc = - (gen_m*gen_alpha+gen_n)*cmsf*cm0**gen_p/sig_psi & *K**(gen_m+0.5)*gen_l**(gen_n+1.) & *(zi+z0)**((gen_m+0.5)*gen_alpha+gen_n) - endif + endif case default end select end function psi_bc @@ -4015,6 +4031,7 @@ subroutine clean_turbulence() if (allocated(P)) deallocate(P) if (allocated(B)) deallocate(B) if (allocated(Pb)) deallocate(Pb) + if (allocated(Px)) deallocate(Px) if (allocated(PSTK)) deallocate(PSTK) if (allocated(num)) deallocate(num) if (allocated(nuh)) deallocate(nuh) @@ -4080,7 +4097,7 @@ subroutine print_state_turbulence() LEVEL2 'tke,eps,L',tke,eps,L LEVEL2 'tkeo',tkeo LEVEL2 'kb,epsb',kb,epsb - LEVEL2 'P,B,Pb,PSTK',P,B,Pb, PSTK + LEVEL2 'P,B,Pb,Px,PSTK',P,B,Pb,Px, PSTK LEVEL2 'num,nuh,nus,nucl',num,nuh,nus, nucl LEVEL2 'gamu,gamv',gamu,gamv LEVEL2 'gamb,gamh,gams',gamb,gamh,gams @@ -4121,16 +4138,16 @@ subroutine print_state_turbulence() LEVEL2 'generic namelist', compute_param,gen_m,gen_n,gen_p, & cpsi1,cpsi2,cpsi3minus,cpsi3plus, & - cpsi4,sig_kpsi,sig_psi, & + cpsix,cpsi4,sig_kpsi,sig_psi, & gen_d,gen_alpha,gen_l - LEVEL2 'keps namelist', ce1,ce2,ce3minus,ce3plus,ce4, & + LEVEL2 'keps namelist', ce1,ce2,ce3minus,ce3plus,cex,ce4, & sig_k,sig_e,sig_peps - LEVEL2 'kw namelist', cw1,cw2,cw3minus,cw3plus,cw4, & + LEVEL2 'kw namelist', cw1,cw2,cw3minus,cw3plus,cwx,cw4, & sig_kw,sig_w - LEVEL2 'my namelist', e1,e2,e3, e6, sq,sl,my_length,new_constr + LEVEL2 'my namelist', e1,e2,e3,ex,e6,sq,sl,my_length,new_constr LEVEL2 'scnd namelist', scnd_method,kb_method,epsb_method, & scnd_coeff, & diff --git a/src/turbulence/variances.F90 b/src/turbulence/variances.F90 index 27910b0d..3feb7898 100644 --- a/src/turbulence/variances.F90 +++ b/src/turbulence/variances.F90 @@ -45,7 +45,7 @@ subroutine variances(nlev,SSU,SSV) ! ! !USES: use turbulence, only: uu,vv,ww - use turbulence, only: tke,eps,P,B,num + use turbulence, only: tke,eps,P,B,Px,num use turbulence, only: cc1,ct1,a2,a3,a5 IMPLICIT NONE ! @@ -90,7 +90,7 @@ subroutine variances(nlev,SSU,SSV) vv(i) = tke(i)*( fac1 + fac2*( fac3*num(i)*SSV(i) & - fac5*num(i)*SSU(i) - 4./3.*a5*B(i) ) ) - ww(i) = tke(i)*( fac1 + fac2*( fac4*P(i) + 8./3.*a5*B(i) ) ) + ww(i) = tke(i)*( fac1 + fac2*( fac4*(P(i)+Px(i)) + 8./3.*a5*B(i) ) ) enddo