From 3101125e14f0fd0f2c354cec66b705a5ac6d0c03 Mon Sep 17 00:00:00 2001 From: Knut Date: Tue, 29 Nov 2022 13:07:30 +0100 Subject: [PATCH 1/5] xP: support independent weighting of xP in lengthscale equations --- src/turbulence/dissipationeq.F90 | 6 ++-- src/turbulence/genericeq.F90 | 6 ++-- src/turbulence/lengthscaleeq.F90 | 6 ++-- src/turbulence/production.F90 | 22 +++++++------- src/turbulence/q2over2eq.F90 | 4 +-- src/turbulence/tkeeq.F90 | 4 +-- src/turbulence/turbulence.F90 | 49 ++++++++++++++++++++++++-------- 7 files changed, 60 insertions(+), 37 deletions(-) diff --git a/src/turbulence/dissipationeq.F90 b/src/turbulence/dissipationeq.F90 index 0af1189a..97b96876 100644 --- a/src/turbulence/dissipationeq.F90 +++ b/src/turbulence/dissipationeq.F90 @@ -62,9 +62,9 @@ subroutine dissipationeq(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS) ! by setting {\tt length\_lim = .true.} in {\tt gotmturb.nml}. ! ! !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 @@ -147,7 +147,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..da89dc87 100644 --- a/src/turbulence/genericeq.F90 +++ b/src/turbulence/genericeq.F90 @@ -118,9 +118,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 +199,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 4b7f7d3d..41fff980 100644 --- a/src/turbulence/lengthscaleeq.F90 +++ b/src/turbulence/lengthscaleeq.F90 @@ -67,9 +67,9 @@ subroutine lengthscaleeq(nlev,dt,depth,u_taus,u_taub,z0s,z0b,h,NN,SS) ! by setting {\tt length\_lim = .true.} in {\tt gotmturb.nml}. ! ! !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 @@ -159,7 +159,7 @@ subroutine lengthscaleeq(nlev,dt,depth,u_taus,u_taub,z0s,z0b,h,NN,SS) avh(i) = sl*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/production.F90 b/src/turbulence/production.F90 index 44ac37cc..df4a43d9 100644 --- a/src/turbulence/production.F90 +++ b/src/turbulence/production.F90 @@ -51,7 +51,7 @@ subroutine production(nlev,NN,SS,xP, SSCSTK) ! 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 use turbulence, only: alpha,iw_model IMPLICIT NONE @@ -89,19 +89,17 @@ subroutine production(nlev,NN,SS,xP, SSCSTK) 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 if ( PRESENT(SSCSTK) ) then do i=0,nlev diff --git a/src/turbulence/q2over2eq.F90 b/src/turbulence/q2over2eq.F90 index 0880ba2a..60135b0f 100644 --- a/src/turbulence/q2over2eq.F90 +++ b/src/turbulence/q2over2eq.F90 @@ -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 @@ -103,7 +103,7 @@ subroutine q2over2eq(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS) avh(i) = sq*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..1c52412d 100644 --- a/src/turbulence/tkeeq.F90 +++ b/src/turbulence/tkeeq.F90 @@ -58,7 +58,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 +118,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 85f3e067..9573640d 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 @@ -157,7 +160,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 @@ -169,7 +173,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 @@ -178,6 +183,7 @@ module turbulence REALTYPE, public :: e1 REALTYPE, public :: e2 REALTYPE, public :: e3 + REALTYPE, public :: ex REALTYPE, public :: e6 REALTYPE, public :: sq REALTYPE, public :: sl @@ -343,14 +349,14 @@ 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 /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, & @@ -431,6 +437,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 @@ -443,6 +450,7 @@ subroutine init_turbulence_nml(namlst,fn,nlev) ce2=1.92 ce3minus=0.0 ce3plus=1.0 + cex=ce1 ce4=_ZERO_ sig_k=1.0 sig_e=1.3 @@ -452,6 +460,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 @@ -648,6 +657,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', '-', & @@ -670,6 +681,8 @@ subroutine init_turbulence_yaml(branch) default=0._rk) call twig%get(ce3plus, 'ce3plus', 'ce3 for unstable stratification', '-', & default=1.0_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._rk) call twig%get(sig_k, 'sig_k', 'Schmidt number for TKE diffusivity', '-', & @@ -686,6 +699,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)', '-', & @@ -866,6 +881,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)' @@ -2230,7 +2249,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 ' ' @@ -2255,6 +2275,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 @@ -2285,7 +2306,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 ' ' @@ -3721,6 +3743,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) @@ -3781,7 +3804,9 @@ 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',P,B,Pb + LEVEL2 'Px', Px + LEVEL2 'PSTK', PSTK LEVEL2 'num,nuh,nus',num,nuh,nus LEVEL2 'gamu,gamv',gamu,gamv LEVEL2 'gamb,gamh,gams',gamb,gamh,gams @@ -3819,14 +3844,14 @@ subroutine print_state_turbulence() kb_min,epsb_min LEVEL2 'generic namelist', compute_param,gen_m,gen_n,gen_p, & - cpsi1,cpsi2,cpsi3minus,cpsi3plus, & - cpsi4,sig_kpsi,sig_psi, & + cpsi1,cpsi2,cpsi3minus,cpsi3plus,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 '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, & From 1a9b53bc9c6e0b367590322072d7b64b28de4ce7 Mon Sep 17 00:00:00 2001 From: Knut Date: Wed, 9 Apr 2025 17:50:32 +0200 Subject: [PATCH 2/5] xP: some formatting --- src/turbulence/turbulence.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/turbulence/turbulence.F90 b/src/turbulence/turbulence.F90 index 83ea8d10..22e99bf4 100644 --- a/src/turbulence/turbulence.F90 +++ b/src/turbulence/turbulence.F90 @@ -4131,8 +4131,8 @@ subroutine print_state_turbulence() kb_min,epsb_min LEVEL2 'generic namelist', compute_param,gen_m,gen_n,gen_p, & - cpsi1,cpsi2,cpsi3minus,cpsi3plus,cpsix,cpsi4, & - sig_kpsi,sig_psi, & + cpsi1,cpsi2,cpsi3minus,cpsi3plus, & + cpsix,cpsi4,sig_kpsi,sig_psi, & gen_d,gen_alpha,gen_l LEVEL2 'keps namelist', ce1,ce2,ce3minus,ce3plus,cex,ce4, & From bf1c8a2259cb982688f8d1277889372f15e0b27c Mon Sep 17 00:00:00 2001 From: Knut Date: Thu, 10 Apr 2025 14:57:48 +0200 Subject: [PATCH 3/5] xP: added cwx to new omega equation --- src/turbulence/omegaeq.F90 | 6 +++--- src/turbulence/turbulence.F90 | 12 +++++++++--- 2 files changed, 12 insertions(+), 6 deletions(-) diff --git a/src/turbulence/omegaeq.F90 b/src/turbulence/omegaeq.F90 index 0be3fb6e..521a0ff3 100644 --- a/src/turbulence/omegaeq.F90 +++ b/src/turbulence/omegaeq.F90 @@ -57,9 +57,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 +130,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/turbulence.F90 b/src/turbulence/turbulence.F90 index 22e99bf4..3e9a320a 100644 --- a/src/turbulence/turbulence.F90 +++ b/src/turbulence/turbulence.F90 @@ -199,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 @@ -380,7 +381,7 @@ subroutine init_turbulence_nml(namlst,fn,nlev) 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,ex,e6,sq,sl,my_length, new_constr @@ -488,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 @@ -744,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', '-', & @@ -2382,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 ' ' @@ -4138,7 +4144,7 @@ subroutine print_state_turbulence() 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,ex,e6,sq,sl,my_length,new_constr From 14a48018aa2f60f22ca8b6ffcc9a114e4b5c417e Mon Sep 17 00:00:00 2001 From: Knut Date: Mon, 14 Apr 2025 10:24:45 +0200 Subject: [PATCH 4/5] xP: manually re-add Px to other calculations for bwd compatibility --- src/turbulence/cmue_a.F90 | 4 ++-- src/turbulence/dissipationeq.F90 | 2 +- src/turbulence/variances.F90 | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) 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 c8fb082a..f45b2732 100644 --- a/src/turbulence/dissipationeq.F90 +++ b/src/turbulence/dissipationeq.F90 @@ -122,7 +122,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 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 From e8cb186dc004f716e32393221f9e528c3098e72b Mon Sep 17 00:00:00 2001 From: Knut Date: Thu, 24 Apr 2025 02:14:29 +0200 Subject: [PATCH 5/5] xP: some documentation --- src/turbulence/dissipationeq.F90 | 6 +++++- src/turbulence/genericeq.F90 | 4 ++++ src/turbulence/lengthscaleeq.F90 | 7 +++++-- src/turbulence/omegaeq.F90 | 6 +++++- src/turbulence/production.F90 | 20 ++++++++++++++++---- src/turbulence/q2over2eq.F90 | 2 +- src/turbulence/tkeeq.F90 | 7 +++++-- 7 files changed, 41 insertions(+), 11 deletions(-) diff --git a/src/turbulence/dissipationeq.F90 b/src/turbulence/dissipationeq.F90 index f45b2732..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. ! diff --git a/src/turbulence/genericeq.F90 b/src/turbulence/genericeq.F90 index da89dc87..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}. diff --git a/src/turbulence/lengthscaleeq.F90 b/src/turbulence/lengthscaleeq.F90 index f8a518f1..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}. ! diff --git a/src/turbulence/omegaeq.F90 b/src/turbulence/omegaeq.F90 index 521a0ff3..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. ! diff --git a/src/turbulence/production.F90 b/src/turbulence/production.F90 index 4249c471..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. diff --git a/src/turbulence/q2over2eq.F90 b/src/turbulence/q2over2eq.F90 index 6fb87564..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$. diff --git a/src/turbulence/tkeeq.F90 b/src/turbulence/tkeeq.F90 index 1c52412d..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