diff --git a/src/gotm/diagnostics.F90 b/src/gotm/diagnostics.F90 index 1acfb380..a412b3ea 100644 --- a/src/gotm/diagnostics.F90 +++ b/src/gotm/diagnostics.F90 @@ -92,9 +92,10 @@ subroutine do_diagnostics(nlev) use density, only: rho0,cp use meanflow, only: gravity,drag use meanflow, only: h,u,v,s,t,NN,SS,buoy,rad + use stokes_drift, only: dusdz, dvsdz use turbulence, only: turb_method use turbulence, only: kappa - use turbulence, only: num,nuh + use turbulence, only: num,nuh, nucl use turbulence, only: tke,P,B use turbulence, only: Rig use observations, only: tprof_input,b_obs_sbf @@ -151,8 +152,10 @@ subroutine do_diagnostics(nlev) taux(nlev) = -tx tauy(nlev) = -ty do i=nlev-1,1,-1 - taux(i)=-num(i)*(u(i+1)-u(i))/(0.5*(h(i+1)+h(i))) - tauy(i)=-num(i)*(v(i+1)-v(i))/(0.5*(h(i+1)+h(i))) + taux(i)=-num(i)*(u(i+1)-u(i))/(0.5*(h(i+1)+h(i))) & + -nucl(i)*dusdz%data(i) + tauy(i)=-num(i)*(v(i+1)-v(i))/(0.5*(h(i+1)+h(i))) & + -nucl(i)*dvsdz%data(i) end do taux(0)=-drag(1)*u(1)*sqrt(u(1)**2+v(1)**2) tauy(0)=-drag(1)*v(1)*sqrt(u(1)**2+v(1)**2) diff --git a/src/gotm/gotm.F90 b/src/gotm/gotm.F90 index c05aa480..57b86f9d 100644 --- a/src/gotm/gotm.F90 +++ b/src/gotm/gotm.F90 @@ -56,7 +56,7 @@ module gotm use turbulence, only: turb_method use turbulence, only: init_turbulence,post_init_turbulence,do_turbulence - use turbulence, only: num,nuh,nus + use turbulence, only: num,nuh,nus, nucl use turbulence, only: const_num,const_nuh use turbulence, only: gamu,gamv,gamh,gams use turbulence, only: Rig @@ -792,8 +792,8 @@ subroutine integrate_gotm() call coriolis(nlev,dt) ! update velocity - call uequation(nlev,dt,cnpar,tx,num,gamu,ext_press_mode) - call vequation(nlev,dt,cnpar,ty,num,gamv,ext_press_mode) + call uequation(nlev,dt,cnpar,tx,num, nucl, gamu,ext_press_mode) + call vequation(nlev,dt,cnpar,ty,num, nucl, gamv,ext_press_mode) call external_pressure(ext_press_mode,nlev) call internal_pressure(nlev) call friction(nlev,kappa,avmolu,tx,ty,plume_type) @@ -889,10 +889,10 @@ subroutine integrate_gotm() ! update one-point models # ifdef SEAGRASS call do_turbulence(nlev,dt,depth,u_taus,u_taub,z0s,z0b,h, & - NN,SS,xP, SSCSTK=SSCSTK) + NN,SS,xP, SSCSTK=SSCSTK, SSSTK=SSSTK) # else call do_turbulence(nlev,dt,depth,u_taus,u_taub,z0s,z0b,h, & - NN,SS, SSCSTK=SSCSTK) + NN,SS, SSCSTK=SSCSTK, SSSTK=SSSTK) # endif end select diff --git a/src/gotm/register_all_variables.F90 b/src/gotm/register_all_variables.F90 index 437460cb..aa0ecf25 100644 --- a/src/gotm/register_all_variables.F90 +++ b/src/gotm/register_all_variables.F90 @@ -513,6 +513,7 @@ subroutine register_turbulence_variables(nlev) call fm%register('num', 'm2/s', 'turbulent diffusivity of momentum', standard_name='??', dimensions=(/id_dim_zi/), data1d=num(0:nlev),category='turbulence', part_of_state=.true.) 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') diff --git a/src/meanflow/meanflow.F90 b/src/meanflow/meanflow.F90 index 9014f30d..3c64cac6 100644 --- a/src/meanflow/meanflow.F90 +++ b/src/meanflow/meanflow.F90 @@ -60,7 +60,7 @@ module meanflow REALTYPE, public, dimension(:), allocatable :: SS,SSU,SSV ! Stokes-Eulerian cross-shear and Stokes shear squared - REALTYPE, public, dimension(:), allocatable :: SSCSTK + REALTYPE, public, dimension(:), allocatable :: SSCSTK, SSSTK ! buoyancy, short-wave radiation, ! extra production of tke by see-grass etc @@ -364,6 +364,10 @@ subroutine post_init_meanflow(nlev,latitude) if (rc /= 0) STOP 'init_meanflow: Error allocating (SSCSTK)' SSCSTK = _ZERO_ + allocate(SSSTK(0:nlev),stat=rc) + if (rc /= 0) STOP 'init_meanflow: Error allocating (SSSTK)' + SSSTK = _ZERO_ + allocate(xP(0:nlev),stat=rc) if (rc /= 0) STOP 'init_meanflow: Error allocating (xP)' xP = _ZERO_ @@ -464,6 +468,7 @@ subroutine clean_meanflow() if (allocated(SSU)) deallocate(SSU) if (allocated(SSV)) deallocate(SSV) if (allocated(SSCSTK)) deallocate(SSCSTK) + if (allocated(SSSTK)) deallocate(SSSTK) if (allocated(xP)) deallocate(xP) if (allocated(buoy)) deallocate(buoy) if (allocated(rad)) deallocate(rad) @@ -523,6 +528,7 @@ subroutine print_state_meanflow() if (allocated(SSU)) LEVEL2 'SSU',SSU if (allocated(SSV)) LEVEL2 'SSV',SSV if (allocated(SSCSTK)) LEVEL2 'SSCSTK',SSCSTK + if (allocated(SSSTK)) LEVEL2 'SSSTK',SSSTK if (allocated(buoy)) LEVEL2 'buoy',buoy if (allocated(rad)) LEVEL2 'rad',rad if (allocated(xp)) LEVEL2 'xP',xP diff --git a/src/meanflow/shear.F90 b/src/meanflow/shear.F90 index e3a2cf81..cb38870f 100644 --- a/src/meanflow/shear.F90 +++ b/src/meanflow/shear.F90 @@ -111,7 +111,7 @@ subroutine shear(nlev,cnpar) ! !USES: use meanflow, only: h,u,v,uo,vo use meanflow, only: SS,SSU,SSV - use meanflow, only: SSCSTK + use meanflow, only: SSCSTK, SSSTK use stokes_drift, only: dusdz, dvsdz IMPLICIT NONE @@ -164,6 +164,9 @@ subroutine shear(nlev,cnpar) SSCSTK(i) = dusdz%data(i) * (u(i+1)-u(i)) / (0.5*(h(i+1)+h(i))) & + dvsdz%data(i) * (v(i+1)-v(i)) / (0.5*(h(i+1)+h(i))) + ! Stokes shear squared + SSSTK(i) = dusdz%data(i)**2 + dvsdz%data(i)**2 + end do SSU(0 ) = SSU(1 ) @@ -178,6 +181,9 @@ subroutine shear(nlev,cnpar) SSCSTK(0 ) = SSCSTK(1 ) SSCSTK(nlev) = SSCSTK(nlev-1) + SSSTK (0 ) = SSSTK (1 ) + SSSTK (nlev) = SSSTK (nlev-1) + return end subroutine shear !EOC diff --git a/src/meanflow/uequation.F90 b/src/meanflow/uequation.F90 index e4bfe0fe..c59c6789 100644 --- a/src/meanflow/uequation.F90 +++ b/src/meanflow/uequation.F90 @@ -5,7 +5,7 @@ ! !ROUTINE: The U-momentum equation\label{sec:uequation} ! ! !INTERFACE: - subroutine uequation(nlev,dt,cnpar,tx,num,gamu,ext_method) + subroutine uequation(nlev,dt,cnpar,tx,num, nucl, gamu,ext_method) ! ! !DESCRIPTION: ! This subroutine computes the transport of momentum in @@ -66,6 +66,7 @@ subroutine uequation(nlev,dt,cnpar,tx,num,gamu,ext_method) use meanflow, only: gravity,avmolu use meanflow, only: h,u,uo,v,w,avh use meanflow, only: drag,SS,runtimeu + use stokes_drift, only: dusdz use observations, only: w_adv_input,w_adv_discr use observations, only: uprof_input,vel_relax_tau,vel_relax_ramp use observations, only: int_press_type @@ -94,6 +95,9 @@ subroutine uequation(nlev,dt,cnpar,tx,num,gamu,ext_method) ! diffusivity of momentum (m^2/s) REALTYPE, intent(in) :: num(0:nlev) +! eddy coefficient of momentum flux down the Stokes gradient (m^2/s) + REALTYPE, intent(in) :: nucl(0:nlev) + ! non-local flux of momentum (m^2/s^2) REALTYPE, intent(in) :: gamu(0:nlev) @@ -175,6 +179,10 @@ subroutine uequation(nlev,dt,cnpar,tx,num,gamu,ext_method) ! add non-local fluxes ! Qsour(i) = Qsour(i) - ( gamu(i) - gamu(i-1) )/h(i) +! add down Stokes gradient fluxes + Qsour(i) = Qsour(i) + ( nucl(i )*dusdz%data(i ) & + -nucl(i-1)*dusdz%data(i-1) )/h(i) + end do ! implement bottom friction as source term diff --git a/src/meanflow/vequation.F90 b/src/meanflow/vequation.F90 index 6923d8f7..c019c3ee 100644 --- a/src/meanflow/vequation.F90 +++ b/src/meanflow/vequation.F90 @@ -5,7 +5,7 @@ ! !ROUTINE: The V-momentum equation\label{sec:vequation} ! ! !INTERFACE: - subroutine vequation(nlev,dt,cnpar,ty,num,gamv,ext_method) + subroutine vequation(nlev,dt,cnpar,ty,num, nucl, gamv,ext_method) ! ! !DESCRIPTION: ! This subroutine computes the transport of momentum in @@ -45,6 +45,7 @@ subroutine vequation(nlev,dt,cnpar,ty,num,gamv,ext_method) use meanflow, only: gravity,avmolu use meanflow, only: h,v,vo,u,w,avh use meanflow, only: drag,SS,runtimev + use stokes_drift, only: dvsdz use observations, only: w_adv_input,w_adv_discr use observations, only: vprof_input,vel_relax_tau,vel_relax_ramp use observations, only: int_press_type @@ -73,6 +74,9 @@ subroutine vequation(nlev,dt,cnpar,ty,num,gamv,ext_method) ! diffusivity of momentum (m^2/s) REALTYPE, intent(in) :: num(0:nlev) +! eddy coefficient of momentum flux down the Stokes gradient (m^2/s) + REALTYPE, intent(in) :: nucl(0:nlev) + ! non-local flux of momentum (m^2/s^2) REALTYPE, intent(in) :: gamv(0:nlev) @@ -154,6 +158,10 @@ subroutine vequation(nlev,dt,cnpar,ty,num,gamv,ext_method) ! add non-local fluxes ! Qsour(i) = Qsour(i) - ( gamv(i) - gamv(i-1) )/h(i) +! add down Stokes gradient fluxes + Qsour(i) = Qsour(i) + ( nucl(i )*dvsdz%data(i ) & + -nucl(i-1)*dvsdz%data(i-1) )/h(i) + end do ! implement bottom friction as source term diff --git a/src/stokes_drift/stokes_drift_theory.F90 b/src/stokes_drift/stokes_drift_theory.F90 index c928b44f..080a7bae 100644 --- a/src/stokes_drift/stokes_drift_theory.F90 +++ b/src/stokes_drift/stokes_drift_theory.F90 @@ -45,8 +45,11 @@ subroutine stokes_drift_theory(nlev,z,zi,u10,v10,gravity) vs0%value = _ZERO_ ds%value = _ZERO_ -! wind direction wind_speed = sqrt(u10**2+v10**2) + + if (wind_speed .gt. _ZERO_) then + +! wind direction xcomp = u10 / wind_speed ycomp = v10 / wind_speed @@ -66,6 +69,8 @@ subroutine stokes_drift_theory(nlev,z,zi,u10,v10,gravity) vs0%value = us0_to_u10 * v10 ds%value = stokes_srf(0)*abs(zi(0))/max(SMALL, sqrt(us0%value**2.+vs0%value**2.)) + end if + return end subroutine stokes_drift_theory diff --git a/src/turbulence/CMakeLists.txt b/src/turbulence/CMakeLists.txt index a793328e..7a317702 100644 --- a/src/turbulence/CMakeLists.txt +++ b/src/turbulence/CMakeLists.txt @@ -5,6 +5,7 @@ add_library(turbulence cmue_b.F90 cmue_c.F90 cmue_d.F90 + cmue_d_h15.F90 cmue_ma.F90 cmue_sg.F90 compute_cpsi3.F90 diff --git a/src/turbulence/alpha_mnb.F90 b/src/turbulence/alpha_mnb.F90 index df5d583e..69aaa540 100644 --- a/src/turbulence/alpha_mnb.F90 +++ b/src/turbulence/alpha_mnb.F90 @@ -5,7 +5,7 @@ ! !ROUTINE: Update dimensionless alpha's\label{sec:alpha} ! ! !INTERFACE: - subroutine alpha_mnb(nlev,NN,SS) + subroutine alpha_mnb(nlev,NN,SS, SSCSTK, SSSTK) ! ! !DESCRIPTION: ! This subroutine updates the dimensionless numbers $\alpha_M$, $\alpha_N$, @@ -22,11 +22,18 @@ subroutine alpha_mnb(nlev,NN,SS) ! !USES: use turbulence, only: tke,eps,kb use turbulence, only: as,an,at + use turbulence, only: av, aw IMPLICIT NONE ! ! !INPUT PARAMETERS: integer, intent(in) :: nlev REALTYPE, intent(in) :: NN(0:nlev),SS(0:nlev) + +! Stokes-Eulerian cross-shear (1/s^2) + REALTYPE, intent(in), optional :: SSCSTK(0:nlev) + +! Stokes shear squared (1/s^2) + REALTYPE, intent(in), optional :: SSSTK (0:nlev) ! ! !REVISION HISTORY: ! Original author(s): Lars Umlauf @@ -35,21 +42,31 @@ subroutine alpha_mnb(nlev,NN,SS) !----------------------------------------------------------------------- ! !LOCAL VARIABLES: integer :: i - REALTYPE :: tau2 + REALTYPE :: tau2(0:nlev) !----------------------------------------------------------------------- !BOC - do i=0,nlev - tau2 = tke(i)*tke(i) / ( eps(i)*eps(i) ) - as(i) = tau2 * SS(i) - an(i) = tau2 * NN(i) - at(i) = tke(i)/eps(i) * kb(i)/eps(i) + do i=0,nlev + tau2(i) = tke(i)*tke(i) / ( eps(i)*eps(i) ) + as(i) = tau2(i) * SS(i) + an(i) = tau2(i) * NN(i) + at(i) = tke(i)/eps(i) * kb(i)/eps(i) + +! clip negative values + as(i) = max(as(i),1.e-10*_ONE_) + at(i) = max(at(i),1.e-10*_ONE_) + end do + + if (present(SSCSTK) .and. present(SSSTK)) then + do i=0,nlev + av(i) = tau2(i) * SSCSTK(i) + aw(i) = tau2(i) * SSSTK(i) -! clip negative values - as(i) = max(as(i),1.e-10*_ONE_) - at(i) = max(at(i),1.e-10*_ONE_) - end do +! clip negative values + aw(i) = max(aw(i),1.e-10*_ONE_) + end do + end if return end subroutine alpha_mnb diff --git a/src/turbulence/cmue_d_h15.F90 b/src/turbulence/cmue_d_h15.F90 new file mode 100644 index 00000000..4ada3064 --- /dev/null +++ b/src/turbulence/cmue_d_h15.F90 @@ -0,0 +1,312 @@ +#include"cppdefs.h" +!----------------------------------------------------------------------- +! +! !ROUTINE: The Langmuir turbulence quasi-equilibrium stability functions after Harcourt (2015)\label{sec:cmueDH15} +! +! !INTERFACE: + subroutine cmue_d_h15(nlev) +! +! !DESCRIPTION: +! +! Old Description from GTOM: This subroutine updates the explicit solution of +! \eq{bijVertical} and \eq{giVertical} under the same assumptions +! as those discussed in \sect{sec:cmueC}. Now, however, an additional +! equilibrium assumption is invoked. With the help of \eq{PeVertical}, +! one can write the equilibrium condition for the TKE as +! \begin{equation} +! \label{quasiEquilibrium} +! \dfrac{P+G}{\epsilon} = +! \hat{c}_\mu(\alpha_M,\alpha_N) \alpha_M +! - \hat{c}'_\mu(\alpha_M,\alpha_N) \alpha_N = 1 +! \comma +! \end{equation} +! where \eq{alphaIdentities} has been used. This is an implicit relation +! to determine $\alpha_M$ as a function of $\alpha_N$. +! With the definitions given in \sect{sec:cmueC}, it turns out that +! $\alpha_M(\alpha_N)$ is a quadratic polynomial that is easily solved. +! The resulting value for $\alpha_M$ is substituted into the stability +! functions described in \sect{sec:cmueC}. For negative $\alpha_N$ +! (convection) the shear number $\alpha_M$ computed in this way may +! become negative. The value of $\alpha_N$ is limited such that this +! does not happen, see \cite{UmlaufBurchard2005a}. +! This Langmuir turbulence version includes the CL Vortex force in +! the algebraic models as well as in the vortex production of TKE and L or epsilon +! +! !USES: + use turbulence, only: an,as,at +! nondimensional forcing functions for Eulerian shear dot Stokes shear, and Stokes shear squared: +! Also, surface proximity function SPF=(1-fzs), goes to zero at surface as tanh(0.25*z/l_S) where l_S +! the vortex-production-weighted dissipation length scale + use turbulence, only: av, aw, SPF + use turbulence, only: tke, L + use turbulence, only: cmue1,cmue2 + use turbulence, only: cmue3 + use turbulence, only: sq, sl, sq_var, sl_var + use turbulence, only: length_lim + + IMPLICIT NONE +! +! !INPUT PARAMETERS: + +! number of vertical layers + integer, intent(in) :: nlev + +! !DEFINED PARAMETERS: + REALTYPE, parameter :: small = 1.0D-8 + +! +! !REVISION HISTORY: +! Original author(s): Lars Umlauf +! Converted by Ramsey Harcourt, last updated 31 July 2018. This version uses the Harcourt(2015) +! stability functions from the quasi-equilibrium Second Moment closure (SMC) with Craik-Leibovich +! terms, but it has been further modified by replacing the crude limiters applied +! individually to Gh, Gv and Gs in Harcourt(2015) under unstable/positive production conditions +! with a combinations of limitations on (L/q)^2 applied consistently across Gm, Gs, Gh, Gv, +! as a function of Gh and Gv input to the ARSM. This ARSM also applies the Galperin limit to +! L going into the ARSM (algebraic) diagnosis of Sm, Sh, Ss, regardless of whether it/s +! being enforced within the dyanamic model. +! +! Recomend running with e3=5 & length.lim=false, but e3=1.2 & length.lim=true is similar +! When length.lim=false, length scale or at least L/q is still limited within the ARSM for +! Stability fcns cmue1,cmue2,cmue3. length.lim=false allows the elevated length scale within +! the mixed layer to impact the transition zone, while restraining the Stability functions +! to the stability-limited length scale regime. +! +!EOP +!----------------------------------------------------------------------- +! !LOCAL VARIABLES: +! + integer, parameter :: rk = kind(_ONE_) + integer :: i + REALTYPE :: Gv, Gs, Gh, Gm + REALTYPE :: Sm, Ss, Sh + + REALTYPE, parameter :: my_A1 = 0.92D0 + REALTYPE, parameter :: my_A2 = 0.74D0 + REALTYPE, parameter :: my_B1 = 16.6D0 + REALTYPE, parameter :: my_B2 = 10.1D0 + REALTYPE, parameter :: my_C1 = 0.08D0 + REALTYPE, parameter :: my_C2 = 0.7D0 + REALTYPE, parameter :: my_C3 = 0.2D0 + REALTYPE, parameter :: h15_Ghmin = -0.28D0 + REALTYPE, parameter :: h15_Ghoff = 0.003D0 + REALTYPE, parameter :: h15_Gvoff = 0.006D0 + REALTYPE, parameter :: h15_Sxmax = 2.12D0 + + REALTYPE :: h15_Shn0, h15_Shnh, h15_Shns, h15_Shnv + REALTYPE :: h15_Shdah, h15_Shdav, h15_Shdbh + REALTYPE :: h15_Shdv, h15_Shdvh, h15_Shdvv + REALTYPE :: h15_Ssn0, h15_Ssdh, h15_Ssdv + REALTYPE :: h15_Smn0, h15_SmnhSh, h15_SmnsSs + REALTYPE :: h15_Smdh, h15_Smdv + + REALTYPE :: tmp0,tmp1,tmp2,tmp3,tmp4 + REALTYPE :: Ghcrit, Gvcrit + +!----------------------------------------------------------------------- +! These constants above & below could all be set or computed elsewhere in advance, +! subject to adjustments in A's, B's & C's. Just sticking them all in here for now. + + h15_Shn0=my_A2*(1.D0-6.D0*my_A1/my_B1) + h15_Shnh=-9.D0*my_A1*my_A2*( my_A2*(1.D0-6.D0*my_A1/my_B1)) + h15_Shns=9.D0*my_A1*my_A2*(1.D0-6.D0*my_A1/my_B1)* & + (2.D0*my_A1+my_A2) + h15_Shnv=9.D0*my_A1*my_A2* & + (my_A2*(1.D0-6.D0*my_A1/my_B1-3.D0*my_C1) & + -2.D0*my_A1*(1.D0-6.D0*my_A1/my_B1+3.D0*my_C1)) + h15_Shdah=-9.D0*my_A1*my_A2 + h15_Shdav=-36.D0*my_A1*my_A1 + h15_Shdbh=-3.D0*my_A2*(6.D0*my_A1+my_B2*(1.D0-my_C3)) + h15_Shdv=-9.D0*my_A2*my_A2*(1.D0-my_C2) + h15_Shdvh=-162.D0*my_A1*my_A1*my_A2*(2.D0*my_A1+(2.D0-my_C2)*my_A2) + h15_Shdvv=324.D0*my_A1*my_A1*my_A2*my_A2*(1.D0-my_C2) + h15_Ssn0=my_A1*(1.D0-6.D0*my_A1/my_B1) + h15_Ssdh=-9.D0*my_A1*my_A2 + h15_Ssdv=-9.D0*my_A1*my_A1 + h15_Smn0 = my_A1*(1.D0-6.D0*my_A1/my_B1-3.D0*my_C1) + h15_SmnhSh = 9.D0*my_A1*(2.D0*my_A1+my_A2*(1.D0-my_C2)) + h15_SmnsSs = 27.D0*my_A1*my_A1 + h15_Smdh = -9.D0*my_A1*my_A2 + h15_Smdv = -36.D0*my_A1*my_A1 + + do i=1,nlev-1 +! convert nondimensional forcing functions to q2-q2l formulation, +! at least until this is rederived in k-epsilon formulation + tmp0 = 4.D0/my_B1**2. + Gh = -tmp0*an(i) + Gm = tmp0*as(i) + Gv = tmp0*av(i) + Gs = tmp0*aw(i) + +! Limitation on the stable side +! if length scale is not strictly limited by stratification, limit it here within the ARSM: + + if (.not.(length_lim)) then + tmp1=1.D0 + + tmp2=h15_Ghmin/min(h15_Ghmin,Gh) + + tmp1=min(tmp1,tmp2) + + if (tmp1.lt.1.D0) then + Gh=Gh*tmp1 + Gv=Gv*tmp1 + Gs=Gs*tmp1 + endif + endif + +! Limitation on the unstable side +! Before SPF: Limit to less than l/q corresponding to equilibrium +! tke eqn soln for (l/q) in limit of gs=gm=0 +! + +! Find ratio to rescale to critical Sh & l/q with offset of -Gvoff and -Ghoff from critical to equilibrium curve for Gm=Gs=0 +! Substituting Gh=R*Gh+Ghoff, Gv=R*Gv+Gvoff to find R such that R*Gh, R*Gv are offset by [-Ghoff, -Gvoff] from critical [Gh, Gv] + + tmp0=2.D0 + + if(Gv.gt._ZERO_) then + +! Coefficient of linear R with contributions from offsets in terms ~Gx*Gy + tmp1= (h15_Shdah+h15_Shdbh)*Gh+(h15_Shdav+h15_Shdv)*Gv + tmp1=tmp1+(h15_Shdah*h15_Ghoff+h15_Shdav*h15_Gvoff)*(h15_Shdbh*Gh) + & + (h15_Shdvh*h15_Ghoff+h15_Shdvv*h15_Gvoff)*Gv + tmp1=tmp1+(h15_Shdah*Gh+h15_Shdav*Gv)*(h15_Shdbh*h15_Ghoff) + & + (h15_Shdvh*Gh+h15_Shdvv*Gv)*h15_Gvoff +! Coefficient of R^2 + tmp2=(h15_Shdah*Gh+h15_Shdav*Gv)*(h15_Shdbh*Gh) + & + (h15_Shdvh*Gh+h15_Shdvv*Gv)*Gv +! Constant term, something a bit smaller than 1 because of offset + tmp4=1.D0+(h15_Shdah+h15_Shdbh)*h15_Ghoff+(h15_Shdav+h15_Shdv)*h15_Gvoff & + +(h15_Shdah*h15_Ghoff+h15_Shdav*h15_Gvoff)*(h15_Shdbh*h15_Ghoff) + & + (h15_Shdvh*h15_Ghoff+h15_Shdvv*h15_Gvoff)*h15_Gvoff + +! Solve for ratio rescaling to critical a*r^2+b*r+c=0; c=1, b=tmp1, a=tmp2 +! Check determinant + + tmp3=tmp1*tmp1-4.D0*tmp2*tmp4 + +! We need the smallest positive root R + if ((tmp3.ge._ZERO_).and.(tmp2.lt._ZERO_)) then + tmp3=(-tmp1+sqrt(tmp3))/(2.D0*tmp2) + elseif ((tmp3.ge._ZERO_).and.(tmp3.gt._ZERO_)) then + tmp3=(-tmp1-sqrt(tmp3))/(2.D0*tmp2) + else +! there is no root i.e. direction of [Gh, Gv] is in stable sector. Set tmp3>1 for no rescaling + tmp3=2.D0 + endif + + if ((tmp3.gt.0.D0).and.(tmp3.lt.1.D0)) then + tmp0=min(tmp0,tmp3) + endif + + endif + +! Apply SPF *after* limiting l/q to equilibrium in limit of gm=gs=0 +! because tke equilibrium does not involve the near-surface pressure-strain +! terms effected by SPF. Apply SPF *before* the limitation on monotonicity of +! Sh/Gh because that might actually be affected by pressure-strain terms. + Gv = Gv*SPF(i) + Gs = Gs*SPF(i)**2. + +! Second limitation after Umlauf & Burchard (2005), Eq 47, approximating -d( Sh/Gh )/d(Gh) > 0 + if(Gh.gt._ZERO_) then + +! This curve is approximated on the Gs=0 plane using the critical Den{Sh}=0 curve but with +! Sh(Gh->2*Gh) when Gh>0, no offsets needed. Bit of a coincidence that hasn't yet been worked out. +! tmp1 is coefficient (b) of linear R with contributions from offsets in terms ~Gx*Gy + tmp1= 2.D0*(h15_Shdah+h15_Shdbh)*Gh+(h15_Shdav+h15_Shdv)*Gv +! Coefficient (a) of R^2 + tmp2=(2.D0*h15_Shdah*Gh+h15_Shdav*Gv)*(2.D0*h15_Shdbh*Gh) + & + (2.D0*h15_Shdvh*Gh+h15_Shdvv*Gv)*Gv +! Constant term is c=1 + tmp4=1.D0 + +! Solve for ratio rescaling to critical a*r^2+b*r+c=0; c=1, b=tmp1, a=tmp2 +! Check determinant + + tmp3=tmp1*tmp1-4.D0*tmp2*tmp4 + + if ((tmp3.ge._ZERO_).and.(tmp2.lt._ZERO_)) then + tmp3=(-tmp1+sqrt(tmp3))/(2.D0*tmp2) + elseif ((tmp3.ge._ZERO_).and.(tmp3.gt._ZERO_)) then + tmp3=(-tmp1-sqrt(tmp3))/(2.D0*tmp2) + else +! there is no root i.e. direction of [Gh, Gv] is in stable sector say tmp3>1 for no rescaling + tmp3=2.D0 + endif + + if ((tmp3.gt.0.D0).and.(tmp3.lt.1.D0)) then + tmp0=min(tmp0,tmp3) + endif + + endif + +! only rescale if R*Gh is less than Gh and same sign, i.e. tmp0 between 0 & 1 + if ((tmp0.gt.0.D0).and.(tmp0.lt.1.D0)) then + Gh = tmp0*Gh + Gm = tmp0*Gm + Gv = tmp0*Gv + Gs = tmp0*Gs + endif + + tmp1=(h15_Shn0+h15_Shnh*Gh+h15_Shns*Gs+h15_Shnv*Gv) + if (tmp1.lt._ZERO_) then + Sh=small + else + tmp2=(1.D0+h15_Shdah*Gh+h15_Shdav*Gv)* & + (1.D0+h15_Shdbh*Gh)+ & + (h15_Shdv+h15_Shdvh*Gh+h15_Shdvv*Gv)*Gv + if (tmp2.le._ZERO_) then + Sh=h15_Sxmax + else + Sh=min(max(small,tmp1/tmp2),h15_Sxmax) + endif + endif + + tmp2=(1.D0+h15_Ssdh*Gh+h15_Ssdv*Gv) + if (tmp2.lt._ZERO_) then + Ss=h15_Sxmax + else + Ss=min(max(small,h15_Ssn0/tmp2),h15_Sxmax) + endif + + tmp1 = (h15_Smn0+h15_SmnhSh*Gh*Sh+h15_SmnsSs*Gs*Ss) +! avoid singularity in Sm from denominator that would cancel -- still need to derive the expression but for now just sidestep + if ((tmp1.lt.small).and.(tmp1.ge._ZERO_)) then + Gh=Gh+small + Gv=Gv+small + tmp1 = (h15_Smn0+h15_SmnhSh*Gh*Sh+h15_SmnsSs*Gs*Ss) + elseif((tmp1.gt.-small).and.(tmp1.lt._ZERO_)) then + Gh=Gh-small + Gv=Gv-small + tmp1 = (h15_Smn0+h15_SmnhSh*Gh*Sh+h15_SmnsSs*Gs*Ss) + endif + + if (tmp1.lt._ZERO_) then + Sm=small + else + tmp2=(1.D0+h15_Smdh*Gh+h15_Smdv*Gv) + if (tmp2.le._ZERO_) then + Sm=h15_Sxmax + else + Sm=min(max(small,tmp1/tmp2),h15_Sxmax) + endif + endif + + Ss=Ss*SPF(i) + + cmue1(i) = sqrt(2.D0)*Sm + cmue2(i) = sqrt(2.D0)*Sh + cmue3(i) = sqrt(2.D0)*Ss + + sq_var(i) = sqrt(sq**2 + (0.41_rk*Sh)**2) + sl_var(i) = sqrt(sl**2 + (0.41_rk*Sh)**2) + + end do + + return + end subroutine cmue_d_h15 + + +!----------------------------------------------------------------------- diff --git a/src/turbulence/lengthscaleeq.F90 b/src/turbulence/lengthscaleeq.F90 index cafe153b..0ac79236 100644 --- a/src/turbulence/lengthscaleeq.F90 +++ b/src/turbulence/lengthscaleeq.F90 @@ -72,7 +72,7 @@ subroutine lengthscaleeq(nlev,dt,depth,u_taus,u_taub,z0s,z0b,h,NN,SS) use turbulence, only: kappa,e1,e2,e3,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 + use turbulence, only: sl_var use util, only: Dirichlet,Neumann IMPLICIT NONE @@ -156,7 +156,7 @@ subroutine lengthscaleeq(nlev,dt,depth,u_taus,u_taub,z0s,z0b,h,NN,SS) do i=1,nlev-1 ! compute diffusivity - avh(i) = sl*sqrt(2.*tkeo(i))*L(i) + 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) diff --git a/src/turbulence/production.F90 b/src/turbulence/production.F90 index 44ac37cc..0a0cfc17 100644 --- a/src/turbulence/production.F90 +++ b/src/turbulence/production.F90 @@ -5,7 +5,7 @@ ! !ROUTINE: Update turbulence production\label{sec:production} ! ! !INTERFACE: - subroutine production(nlev,NN,SS,xP, SSCSTK) + subroutine production(nlev,NN,SS,xP, SSCSTK, SSSTK) ! ! !DESCRIPTION: ! This subroutine calculates the production terms of turbulent kinetic @@ -52,7 +52,7 @@ subroutine production(nlev,NN,SS,xP, SSCSTK) ! ! !USES: use turbulence, only: P,B,Pb, PSTK - use turbulence, only: num,nuh + use turbulence, only: num,nuh, nucl use turbulence, only: alpha,iw_model IMPLICIT NONE ! @@ -73,6 +73,9 @@ subroutine production(nlev,NN,SS,xP, SSCSTK) ! Stokes-Eulerian cross-shear (1/s^2) REALTYPE, intent(in), optional :: SSCSTK(0:nlev) + +! Stokes shear squared (1/s^2) + REALTYPE, intent(in), optional :: SSSTK (0:nlev) ! ! !REVISION HISTORY: ! Original author(s): Karsten Bolding, Hans Burchard @@ -103,11 +106,22 @@ subroutine production(nlev,NN,SS,xP, SSCSTK) enddo endif +! P is -du/dz for q2l production with e1, +! PSTK is -dus/dz for q2l prodcution with e6, +! see (6) in Harcourt2015. +! - = num*(du/dz) + nucl*(dus/dz), see (7) in Harcourt2015. if ( PRESENT(SSCSTK) ) then do i=0,nlev + P(i) = P(i) + nucl(i)*SSCSTK(i) PSTK(i) = num (i)*SSCSTK(i) end do end if + if ( PRESENT(SSSTK) ) then + if ( .not. PRESENT(SSCSTK) ) PSTK = _ZERO_ + do i=0,nlev + PSTK(i) = PSTK(i) + nucl(i)*SSSTK(i) + end do + end if return end subroutine production diff --git a/src/turbulence/q2over2eq.F90 b/src/turbulence/q2over2eq.F90 index 0880ba2a..40fa2b3a 100644 --- a/src/turbulence/q2over2eq.F90 +++ b/src/turbulence/q2over2eq.F90 @@ -47,7 +47,7 @@ subroutine q2over2eq(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS) use turbulence, only: P,B, 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 + use turbulence, only: sq_var use util, only: Dirichlet,Neumann IMPLICIT NONE @@ -100,7 +100,7 @@ subroutine q2over2eq(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS) do i=1,nlev-1 ! compute diffusivity - avh(i) = sq*sqrt( 2.*tke(i) )*L(i) + avh(i) = sq_var(i) * sqrt( 2.*tke(i) )*L(i) ! compute production terms in q^2/2-equation prod = P(i) + PSTK(i) diff --git a/src/turbulence/turbulence.F90 b/src/turbulence/turbulence.F90 index 41346c10..c29ce1f0 100644 --- a/src/turbulence/turbulence.F90 +++ b/src/turbulence/turbulence.F90 @@ -71,6 +71,11 @@ module turbulence REALTYPE, public, dimension(:), allocatable, target :: nuh REALTYPE, public, dimension(:), allocatable :: nus +! turbulent eddy coefficient for momentum flux down Stokes gradient +! in second moment closures with Craik-Leibovich vortex force in the +! algebraic Reynolds stress and flux models. + REALTYPE, public, dimension(:), allocatable :: nucl + ! non-local fluxes of momentum REALTYPE, public, dimension(:), allocatable :: gamu,gamv @@ -79,7 +84,10 @@ module turbulence REALTYPE, public, dimension(:), allocatable :: gamb,gamh,gams ! non-dimensional stability functions - REALTYPE, public, dimension(:), allocatable :: cmue1,cmue2 + REALTYPE, public, dimension(:), allocatable :: cmue1,cmue2, cmue3 + +! spatially varying sq and sl for q2over2 and lengthscale + REALTYPE, public, dimension(:), allocatable :: sq_var, sl_var ! non-dimensional counter-gradient term REALTYPE, public, dimension(:), allocatable :: gam @@ -87,6 +95,13 @@ module turbulence ! alpha_M, alpha_N, and alpha_B REALTYPE, public, dimension(:), allocatable :: as,an,at +! alpha_V, alpha_W +! dimensionless Stokes-Eulerian cross-shear and Stokes shear^2 + REALTYPE, public, dimension(:), allocatable :: av,aw + +! surface proximity function + REALTYPE, public, dimension(:), allocatable :: SPF + ! time scale ratio r REALTYPE, public, dimension(:), allocatable :: r @@ -258,6 +273,7 @@ module turbulence integer, parameter :: quasi_Eq=1 integer, parameter :: weak_Eq_Kb_Eq=2 integer, parameter :: weak_Eq_Kb=3 + integer, parameter :: quasi_Eq_H15=4 ! method to solve equation for k_b integer, parameter :: kb_algebraic=1 @@ -740,7 +756,7 @@ subroutine init_turbulence_yaml(branch) twig => branch%get_child('scnd', 'second-order model') call twig%get(scnd_method, 'method', 'method', & - options=(/option(quasi_eq, 'quasi-equilibrium', 'quasi_eq'), option(weak_eq_kb_eq, 'weak equilibrium with algebraic buoyancy variance', 'weak_eq_kb_eq')/), default=weak_eq_kb_eq) + options=(/option(quasi_eq, 'quasi-equilibrium', 'quasi_eq'), option(weak_eq_kb_eq, 'weak equilibrium with algebraic buoyancy variance', 'weak_eq_kb_eq'),option(quasi_eq_h15, 'quasi-equilibrium with Langmuir (Harcourt, 2015)', 'quasi_eq_h15')/), default=weak_eq_kb_eq) call twig%get(kb_method, 'kb_method', 'equation for buoyancy variance', & options=(/option(kb_algebraic, 'algebraic', 'algebraic'), option(kb_dynamic, 'prognostic', 'prognostic')/), default=1, display=display_advanced) call twig%get(epsb_method, 'epsb_method', 'equation for variance destruction', & @@ -926,6 +942,10 @@ subroutine post_init_turbulence(nlev) if (rc /= 0) stop 'init_turbulence: Error allocating (nus)' nus = 1.0D-6 + allocate(nucl(0:nlev),stat=rc) + if (rc /= 0) stop 'init_turbulence: Error allocating (nucl)' + nucl = _ZERO_ + allocate(gamu(0:nlev),stat=rc) if (rc /= 0) stop 'init_turbulence: Error allocating (gamu)' gamu = _ZERO_ @@ -954,6 +974,18 @@ subroutine post_init_turbulence(nlev) if (rc /= 0) stop 'init_turbulence: Error allocating (cmue2)' cmue2 = _ZERO_ + allocate(cmue3(0:nlev),stat=rc) + if (rc /= 0) stop 'init_turbulence: Error allocating (cmue3)' + cmue3 = _ZERO_ + + allocate(sq_var(0:nlev),stat=rc) + if (rc /= 0) stop 'init_turbulence: Error allocating (sq_var)' + sq_var = sq + + allocate(sl_var(0:nlev),stat=rc) + if (rc /= 0) stop 'init_turbulence: Error allocating (sl_var)' + sl_var = sl + allocate(gam(0:nlev),stat=rc) if (rc /= 0) stop 'init_turbulence: Error allocating (gam)' gam = _ZERO_ @@ -970,6 +1002,18 @@ subroutine post_init_turbulence(nlev) if (rc /= 0) stop 'init_turbulence: Error allocating (at)' at = _ZERO_ + allocate(av(0:nlev),stat=rc) + if (rc /= 0) stop 'init_turbulence: Error allocating (av)' + av = _ZERO_ + + allocate(aw(0:nlev),stat=rc) + if (rc /= 0) stop 'init_turbulence: Error allocating (aw)' + aw = _ZERO_ + + allocate(SPF(0:nlev),stat=rc) + if (rc /= 0) stop 'init_turbulence: Error allocating (SPF)' + SPF = _ONE_ + allocate(r(0:nlev),stat=rc) if (rc /= 0) stop 'init_turbulence: Error allocating (r)' r = _ZERO_ @@ -2402,7 +2446,7 @@ end subroutine report_model ! ! !INTERFACE: subroutine do_turbulence(nlev,dt,depth,u_taus,u_taub,z0s,z0b,h, & - NN,SS,xP, SSCSTK) + NN,SS,xP, SSCSTK, SSSTK) ! ! !DESCRIPTION: This routine is the central point of the ! turbulence scheme. It determines the order, in which @@ -2458,13 +2502,22 @@ subroutine do_turbulence(nlev,dt,depth,u_taus,u_taub,z0s,z0b,h, & IMPLICIT NONE interface - subroutine production(nlev,NN,SS,xP, SSCSTK) + subroutine production(nlev,NN,SS,xP, SSCSTK, SSSTK) integer, intent(in) :: nlev REALTYPE, intent(in) :: NN(0:nlev) REALTYPE, intent(in) :: SS(0:nlev) REALTYPE, intent(in), optional :: xP(0:nlev) REALTYPE, intent(in), optional :: SSCSTK(0:nlev) + REALTYPE, intent(in), optional :: SSSTK (0:nlev) end subroutine production + + subroutine alpha_mnb(nlev,NN,SS, SSCSTK, SSSTK) + integer, intent(in) :: nlev + REALTYPE, intent(in) :: NN(0:nlev) + REALTYPE, intent(in) :: SS(0:nlev) + REALTYPE, intent(in), optional :: SSCSTK(0:nlev) + REALTYPE, intent(in), optional :: SSSTK (0:nlev) + end subroutine alpha_mnb end interface ! @@ -2503,6 +2556,9 @@ end subroutine production ! Stokes-Eulerian cross-shear (1/s^2) REALTYPE, intent(in), optional :: SSCSTK(0:nlev) + +! Stokes shear squared (1/s^2) + REALTYPE, intent(in), optional :: SSSTK (0:nlev) ! ! !REVISION HISTORY: ! Original author(s): Karsten Bolding, Hans Burchard, @@ -2531,10 +2587,10 @@ end subroutine production if ( PRESENT(xP) ) then ! with seagrass turbulence - call production(nlev,NN,SS,xP, SSCSTK=SSCSTK) + call production(nlev,NN,SS,xP, SSCSTK=SSCSTK, SSSTK=SSSTK) else ! without - call production(nlev,NN,SS, SSCSTK=SSCSTK) + call production(nlev,NN,SS, SSCSTK=SSCSTK, SSSTK=SSSTK) end if call alpha_mnb(nlev,NN,SS) @@ -2552,10 +2608,10 @@ end subroutine production if ( PRESENT(xP) ) then ! with seagrass turbulence - call production(nlev,NN,SS,xP, SSCSTK=SSCSTK) + call production(nlev,NN,SS,xP, SSCSTK=SSCSTK, SSSTK=SSSTK) else ! without - call production(nlev,NN,SS, SSCSTK=SSCSTK) + call production(nlev,NN,SS, SSCSTK=SSCSTK, SSSTK=SSSTK) end if select case(scnd_method) @@ -2589,6 +2645,21 @@ end subroutine production STDERR 'Program execution stopped ...' stop 'turbulence.F90' + case (quasi_Eq_H15) + ! quasi-equilibrium model + ! SMC model of Langmuir turbulence (Harcourt, 2015) +! KK-TODO: test whether SSCSTK and SSSTK are present? + + call alpha_mnb(nlev,NN,SS, SSCSTK=SSCSTK, SSSTK=SSSTK) + call surface_proximity_function(nlev,depth,h) + call cmue_d_h15(nlev) + call do_tke(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS) + call do_kb(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS) + call do_lengthscale(nlev,dt,depth,u_taus,u_taub,z0s,z0b,h,NN,SS) + call do_epsb(nlev,dt,u_taus,u_taub,z0s,z0b,h,NN,SS) + call alpha_mnb(nlev,NN,SS, SSCSTK=SSCSTK, SSSTK=SSSTK) + call kolpran(nlev) + case default STDERR 'Not a valid method for second-order model' @@ -2890,6 +2961,8 @@ subroutine kolpran(nlev) nuh(i) = cmue2(i)*x ! salinity nus(i) = cmue2(i)*x +! momentum down Stokes gradient + nucl(i) = cmue3(i)*x end do return @@ -3083,6 +3156,78 @@ subroutine compute_cm0(turb_method,stab_method,scnd_method) end subroutine compute_cm0 !EOC +!----------------------------------------------------------------------- +!BOP +! +! !IROUTINE: Compute the surface proximity function +! +! !INTERFACE: + subroutine surface_proximity_function(nlev,depth,h) +! +! !DESCRIPTION: +! Compute the surface proximity function (SPF) used the second moment +! closure model of Langmuir turbulence \citep{Harcourt2015}. SPF is defined +! as $\mathrm{SPF} = 1-f_z^S$, where +! \begin{equation} +! f_z^S=1+\tanh\left(C_l^S z/l^S\right), +! \end{equation} +! with $C_l^S=0.25$. The length scale $l^S$ is an average of the dissipation +! length scale weighted by positive CL vortex production, +! \begin{equation} +! l^S=\langle l \rangle_{P^S>0}\equiv\frac{\int l P^S_+ dz}{\int P^S_+ dz} +! \end{equation} +! where $P^S_+=\max(0,P^S)$. +! +! !USES: + IMPLICIT NONE +! +! !INPUT PARAMETERS: + integer, intent(in) :: nlev + REALTYPE, intent(in) :: depth + REALTYPE, intent(in) :: h(0:nlev) +! +! !REVISION HISTORY: +! Original author(s): Qing Li +! +!EOP +! +! !LOCAL VARIABLES: + integer :: i + REALTYPE :: lS, znrm, dz + REALTYPE :: fzS(0:nlev), zz(0:nlev) + REALTYPE, parameter :: ClS = 0.25 +! +!----------------------------------------------------------------------- +!BOC + + ! find the length scale + lS = _ZERO_ + znrm = _ZERO_ + do i=1,nlev-1 + dz = 0.5*(h(i)+h(i+1)) + znrm = znrm + max(_ZERO_,PSTK(i))*dz + lS = lS + L(i)*max(_ZERO_,PSTK(i))*dz + end do + znrm = znrm + 0.5*max(_ZERO_,PSTK(nlev))*h(nlev) + lS = lS + 0.5*L(nlev)*max(_ZERO_,PSTK(nlev))*h(nlev) + lS = lS/max(SMALL,znrm) + + ! get depth at cell interface + zz(0) = -depth + do i=1,nlev + zz(i) = zz(i-1) + h(i) + end do + + fzS = _ZERO_ + do i=0,nlev + ! compute fzS + fzS(i) = _ONE_ + tanh(ClS*zz(i)/lS) + ! surface proximity function + SPF(i) = _ONE_ - fzS(i) + end do + + end subroutine surface_proximity_function +!EOC !----------------------------------------------------------------------- !BOP @@ -3874,6 +4019,7 @@ subroutine clean_turbulence() if (allocated(num)) deallocate(num) if (allocated(nuh)) deallocate(nuh) if (allocated(nus)) deallocate(nus) + if (allocated(nucl)) deallocate(nucl) if (allocated(gamu)) deallocate(gamu) if (allocated(gamv)) deallocate(gamv) if (allocated(gamb)) deallocate(gamb) @@ -3881,10 +4027,14 @@ subroutine clean_turbulence() if (allocated(gams)) deallocate(gams) if (allocated(cmue1)) deallocate(cmue1) if (allocated(cmue2)) deallocate(cmue2) + if (allocated(cmue3)) deallocate(cmue3) if (allocated(gam)) deallocate(gam) if (allocated(an)) deallocate(an) if (allocated(as)) deallocate(as) if (allocated(at)) deallocate(at) + if (allocated(av)) deallocate(av) + if (allocated(aw)) deallocate(aw) + if (allocated(SPF)) deallocate(SPF) if (allocated(r)) deallocate(r) if (allocated(Rig)) deallocate(Rig) if (allocated(xRf)) deallocate(xRf) @@ -3931,12 +4081,14 @@ subroutine print_state_turbulence() LEVEL2 'tkeo',tkeo LEVEL2 'kb,epsb',kb,epsb LEVEL2 'P,B,Pb,PSTK',P,B,Pb, PSTK - LEVEL2 'num,nuh,nus',num,nuh,nus + LEVEL2 'num,nuh,nus,nucl',num,nuh,nus, nucl LEVEL2 'gamu,gamv',gamu,gamv LEVEL2 'gamb,gamh,gams',gamb,gamh,gams - LEVEL2 'cmue1,cmue2',cmue1,cmue2 + LEVEL2 'cmue1,cmue2,cmue3',cmue1,cmue2, cmue3 LEVEL2 'gam',gam LEVEL2 'as,an,at',as,an,at + LEVEL2 'av aw', av, aw + LEVEL2 'SPF', SPF LEVEL2 'r',r LEVEL2 'Rig',Rig LEVEL2 'xRf',xRf