Newer
Older
MODULE TURB
use NUMERICS, only : PROGONKA, IND_STAB_FACT_DB

Victor Stepanenko
committed
use INOUT, only : CHECK_UNIT
use LAKE_DATATYPES, only : ireals, iintegers
use DRIVING_PARAMS, only : missing_value

Victor Stepanenko
committed
SUBROUTINE K_EPSILON(ix, iy, nx, ny, year, month, day, &
& hour, kor, a_veg, c_veg, h_veg, dt, &
& b0, tau_air, tau_wav, tau_i, tau_gr, roughness, fetch)
! KT_eq calculates eddy diffusivity (ED) in water coloumn following parameterizations:
! 1) empirical profile of ED;
! 2) Semi-empirical formulations for ED;
! 2) k-epsilon parameterization for ED, including Kolmogorov relation.
use COMPARAMS, only: &
& num_ompthr
use DRIVING_PARAMS , only : &
& nstep_keps, &
& turb_out, &
& Turbpar, &
& path, &
& stabfunc, &
& kepsbc, &
& kwe, &
& M, &
& omp, &
& eos, lindens, &
& zserout
use ARRAYS_TURB , only : &
& S_integr_positive, &
& S_integr_negative, &
& Gen_integr, &
& Seps_integr_positive, &
& Seps_integr_negative, &
& Geneps_integr, &
& epseps_integr, &
& eps_integr, &
& E_integr, &
& TKE_balance, &
& eps_balance, &
& Seps_integr_positive, &
& Seps_integr_negative, &
& Geneps_integr, &
& epseps_integr, &
& veg_sw, &
& E_it1, E_it2, &
& E1, E2, &
& eps_it1, eps_it2, &
& eps1, eps2, &
& k2_mid, k2, k2t, &
& E12, eps12, &
& row, row2, &
& k5_mid, k5, &
& L, &
& TF, &
& k3_mid, &
& WU_, WV_, &
& GAMT, GAMU, GAMV, &
& S, &
& F, &
& KT, &
& TKE_turb_trans, &
& Re, &
& C1aup , &
& k_turb_T_flux, &
& H_entrainment, &
& Buoyancy0, signwaveheight, &
& Eseiches, &
& TKE_budget_terms
use ARRAYS_GRID, only : &
& ddz, ddz2, ddz05, ddz052, ddz054, &
& dzeta_int, dzeta_05int, z_half
use ARRAYS_BATHYM, only : &
& dhw, dhw0, &
use ARRAYS_WATERSTATE, only : &
& Tw1, Tw2, &
& Sal1, Sal2, lamw, lamsal
use ARRAYS, only : &
& um, u2, u1, &
& vm, v2, v1, &
& nstep, &
& time, &
& uv
use PHYS_CONSTANTS, only: &
& row0, roa0, &
& lamw0, &
& cw, &
& g, &
& cw_m_row0, &
& roughness0, niu_wat, &

Victor Stepanenko
committed
& alsal, &
& Racr
use WATER_DENSITY, only: &
& water_dens_t_unesco, &
& water_dens_ts, &
& DDENS_DTEMP0, DDENS_DSAL0, &
& SET_DDENS_DTEMP, SET_DDENS_DSAL
use INOUT_PARAMETERS, only : &
& lake_subr_unit_min, &
& lake_subr_unit_max
use NUMERICS, only : STEP
use SEICHES_PARAM, only: SEICHE_ENERGY, C_Deff, cnorm
implicit none
! Input variables
! Reals
real(kind=ireals), intent(in) :: dt
real(kind=ireals), intent(in) :: kor
real(kind=ireals), intent(in) :: h_veg, a_veg, c_veg
real(kind=ireals), intent(in) :: hour
real(kind=ireals), intent(in) :: b0
real(kind=ireals), intent(in) :: tau_air, tau_wav, tau_i, tau_gr
real(kind=ireals), intent(in) :: roughness
real(kind=ireals), intent(in) :: fetch
integer(kind=iintegers), intent(in) :: ix, iy
integer(kind=iintegers), intent(in) :: nx, ny
integer(kind=iintegers), intent(in) :: year, month, day
real(kind=ireals) :: C_eps1(M), C_eps2(M), C_eps3(M)
real(kind=ireals) :: CE(M), CEt(M)
real(kind=ireals) :: lam_E(M+1), lam_eps(M+1)
real(kind=ireals) :: a(vector_length)
real(kind=ireals) :: b(vector_length)
real(kind=ireals) :: c(vector_length)
real(kind=ireals) :: d(vector_length)
real(kind=ireals) :: AG(5)
real(kind=ireals) :: dt05

Victor Stepanenko
committed
real(kind=ireals) :: wr
real(kind=ireals) :: dhw_keps, dhw0_keps
real(kind=ireals) :: pi
real(kind=ireals) :: ufr
real(kind=ireals) :: dist_surf
real(kind=ireals) :: ACC, ACC2, ACCk
real(kind=ireals) :: dE_it, deps_it
real(kind=ireals) :: xx(1:2), yy(1:2), zz
real(kind=ireals) :: lm
real(kind=ireals) :: ext_lamw
real(kind=ireals) :: month_sec, day_sec, hour_sec
real(kind=ireals) :: al_it
real(kind=ireals) :: maxTKEinput, maxTKEinput_i, maxTKEinput_gr

Victor Stepanenko
committed
real(kind=ireals) :: FS, FTKES, TKES, DISS_TKES
real(kind=ireals) :: FB, FTKEB, TKEB, DISS_TKEB
real(kind=ireals) :: al, alft
real(kind=ireals) :: GAMUN
real(kind=ireals) :: urels, vrels
real(kind=ireals) :: ksurf1(1:2), ksurf2(1:2)
real(kind=ireals) :: Esurf1(1:2), Esurf2(1:2)
real(kind=ireals) :: ceps3
real(kind=ireals) :: vdamp
real(kind=ireals) :: dE_it_max, deps_it_max

Victor Stepanenko
committed
real(kind=ireals) :: E_min, E_min_lamw0, E_min_Racr, eps_min, eps_min0
real(kind=ireals) :: Buoyancy1, dBdz0, z0, bvf2
real(kind=ireals) :: Eseiches_diss
real(kind=ireals) :: tf_integr
real(kind=ireals) :: GRADU, GRADV, GRADT, GAMVN
real(kind=ireals) :: TFR
real(kind=ireals) :: COGRTN, COGRUN, COGRVN
real(kind=ireals) :: DTDZH
real(kind=ireals), allocatable :: work(:,:)
real(kind=ireals), allocatable :: rhotemp(:), rhosal(:)

Victor Stepanenko
committed
integer(kind=iintegers), parameter :: maxiter = 35
integer(kind=iintegers) :: iter, iter_t
integer(kind=iintegers) :: nl
integer(kind=iintegers) :: i, j, k
integer(kind=iintegers) :: keps_coef
! integer(kind=iintegers) :: stabfunc
integer(kind=iintegers) :: time_deriv
integer(kind=iintegers) :: i_entrain, i_entrain_old
integer(kind=iintegers) :: nunit_ = lake_subr_unit_min
integer(kind=iintegers) :: grav_wave_Gill = 0
integer(kind=iintegers) :: seiches_Goudsmit = 0
! Logicals
logical, allocatable :: init_keps(:,:)
logical :: indstab
logical :: ind_bound
logical :: iterat
logical :: firstcall
logical :: cycle_keps
logical :: next_tstep_keps
logical :: smooth
logical :: perdam
logical :: semi_implicit_scheme
logical :: contrgrad
! Characters
character :: tp_name*10
character :: numt*1
data firstcall /.true./
! Externals
real(kind=ireals), external :: DZETA
SAVE
firstcall_if : if (firstcall) then
if (turb_out%par == 1) then

Victor Stepanenko
committed
call CHECK_UNIT(lake_subr_unit_min,lake_subr_unit_max,nunit_)
open (nunit_,file=path(1:len_trim(path))//'results/debug/err_keps_iter.dat', &

Victor Stepanenko
committed
write(unit=nunit_,fmt=*) 'Errors of iterative process in k-epsilon solver'
endif
allocate (init_keps(1:nx, 1:ny) )
init_keps(:,:) = .false.
month_sec = 30.*24.*60.*60.
day_sec = 24*60.*60.
hour_sec = 60.*60.

Victor Stepanenko
committed
! Parameters of numerical scheme
semi_implicit_scheme = .true. ! No iterations in k-epsilon, semi-implicit scheme
al_it = 0.2 !0.8
ACC = 1.d-20 !1.d-20
ACC2 = 1.d-20 !1.d-20
ACCk = 1.d-20 !1.d-20
knum = 0.
z0 = 1.d-2

Victor Stepanenko
committed
eps_min0 = 1.d-12 ! The range for dissipation rate values
! 10**(-11) - 10**(-6) m**2/s**3 (Wuest and Lorke, 2003)
E_min_lamw0 = sqrt(eps_min0*1.E-5*lamw0*(cw_m_row0*CEt0)**(-1) ) ! This expression ensures that
! eddy diffusivity in the decaying
! turbulence (E -> Emin, eps -> eps_min)
! is much less than the molecular diffusivity
dE_it_max = 1.d-11 !1.d-9
deps_it_max = 1.d-15 !1.d-12
smooth = .false.
perdam = .false.
vdamp = 1.d-3 !1.d-3
contrgrad = .false.

Victor Stepanenko
committed
! AG(1) = 0._ireals
! AG(2) = 0._ireals
! AG(3) = 0._ireals
! AG(4)=1. !0.
! AG(5)=1. !0.

Victor Stepanenko
committed
! keps_coef = 1 - standard empirical coefficients in epsilon-equation
! keps_coef = 2 - the coefficient by (Aupoix et al., 1989) is used
! stabfunc = 1 - stability functions CE and CEt are constants
! stabfunc = 2 - stability functions CE and CEt are dependent
! on shear and buoyancy (Canuto et al., 2001)
! stabfunc = 3 - stability functions CE and CEt are dependent
! only on buoyancy (Galperin et al., 1988)
! FRICTION: VEGETATION
if (h_veg>0) then
print*, 'Vegetation friction parameterization &
& is not operational currently: STOP'
STOP
endif
veg_sw = 0.
do i = 1, M+1
if (h1-dzeta_int(i)*h1 < h_veg) veg_sw(i) = 1.
enddo
endif firstcall_if

Victor Stepanenko
committed
if (l1 > 0._ireals) then
!The minimal values for E and epsilon are determined from two conditions:
!1.corresponding turbulent diffusion coefficient is much less than molecular coefficient
!2.free convection occurs in quiscent state under Ra>Racr
eps_min = 1.E-5*(lamw0/cw_m_row0)**2*Racr*niu_wat/(maxval(ddz)*h1)**4
E_min_Racr = eps_min*(maxval(ddz)*h1)**2*(CEt0*Racr*niu_wat*lamw0/cw_m_row0)**(-0.5)
!E_min = max(E_min_lamw0,E_min_Racr)
E_min = E_min_Racr

Victor Stepanenko
committed
else
eps_min = eps_min0
E_min = E_min_lamw0

Victor Stepanenko
committed
endif
dt05 = 0.5d0*dt
allocate (work(1:M+1,1:2))
allocate (rhotemp(1:M), rhosal(1:M))

Victor Stepanenko
committed
! Implementation of the Crank-Nicolson numerical scheme
! for k-epsilon parameterization using simple iterations to
! solve the set of nonlinear finite-difference equations
! k_turb_T_flux used here from the previous time step
! Mixed layer defined from the gradient of heat turbulent flux
!i_entrain = M
!do i = M-1, 1, -1
! if (abs(k_turb_T_flux(i+1)-k_turb_T_flux(i))/abs(k_turb_T_flux(i+1) + 1.d-20) &
! & > 0.2 ) then ! 0.2 - quite arbitrary value
! i_entrain = i
! exit
! endif
!enddo
! Seiche energy evolution, TKE production due to seiche dissipation
if (seiches_Goudsmit == 1) then
call SEICHE_ENERGY &

Victor Stepanenko
committed
& (M,bathymwater,ls,tau_air-tau_wav,sqrt(uwind*uwind+vwind*vwind), &
& vol,dt,Eseiches,Eseiches_diss)
!A part of seiche energy dissipation feeding TKE (Goudsmit et al. 2002, eg. (24))
do i = 1, M
bvf2 = max(AL*(row(i+1) - row(i))/(h1*ddz(i)),0._ireals)
Gen_seiches(i) = (1. - 10.*sqrt(C_Deff))*bvf2/(row0*cnorm*botar)*&
& (area_int(i) - area_int(i+1))/area_half(i)*Eseiches_diss
enddo
else
Gen_seiches(1:M) = 0.
endif
! Mixed layer defined from TKE gradient
i_entrain = minloc(E1(1:M),1)
!do i = 1, M
! if (E1(i) < 1.E-6) then
! i_entrain = i
! exit
! endif
!enddo
!do i = i_entrain, 1, -1
! if (E1(i)/E1(i_entrain) > 5.E+0_ireals) then
! i_entrain = i
! exit
! endif
!enddo
! Bulk Richardson number for mixed layer
Ri_bulk = g/row0*(row(i_entrain)-row(1))*(dzeta_int(i_entrain)*h1)/ &
& ( (u1(i_entrain) - u1(1))**2 + (v1(i_entrain) - v1(1))**2 + small_value)
! Significant wave height
if (kepsbc%par == 2 .or. kepsbc%par == 3) then
signwaveheight = H13(tau_air,fetch)
endif
dhw_keps = dhw / real(nstep_keps%par)
dhw0_keps = dhw0 / real(nstep_keps%par)
! This time interpolation ensures correspondence between mean energy
! viscous dissipation in Crank-Nickolson scheme for momentum equation and TKE shear production
! when semi-implicit scheme is used (semi_implicit_scheme = .true.)
do i = 1, M+1
Um(i) = 0.5d0 * (u2(i) + u1(i))
Vm(i) = 0.5d0 * (v2(i) + v1(i))
enddo

Victor Stepanenko
committed
! The TKE injection from waves breaking is limited
! by the wind kinetic energy input to waves development
maxTKEinput = tau_wav/row0*sqrt(uwind*uwind + vwind*vwind)

Victor Stepanenko
committed
! maxTKEinput_i = tau_i/row0*sqrt(Um(1)*Um(1)+Vm(1)*Vm(1))
! maxTKEinput_gr = tau_gr/row0*sqrt(Um(M+1)*Um(M+1)+Vm(M+1)*Vm(M+1))

Victor Stepanenko
committed
! DTDZH=(Tw1(2)-Tw1(1))/(h1*ddz(1))
next_tstep_keps = .false.
keps_mode: do while (.not.next_tstep_keps)
if (.not.init_keps(ix,iy)) then

Victor Stepanenko
committed
! Initialization mode of k-epsilon parameterization: time derivatives are neglected

Victor Stepanenko
committed
! Evolution mode of k-epsilon parameterization: full equations are solved
time_deriv = 1
endif
do i = 1, M+1
E_it1(i) = E1(i)
eps_it1(i) = eps1(i)
k2_mid(i) = 0._ireals
k2(i) = 0._ireals
k2t(i) = 0._ireals
enddo
iter_t = 0
cycle_keps = .true.
if (semi_implicit_scheme) then
iterat = .false.
else
iterat = .true.
endif
! Setting density derivatives on temperature and salinity

Victor Stepanenko
committed
rhotemp(i) = SET_DDENS_DTEMP(eos%par,lindens%par,0.5*(Tw1(i)+Tw1(i+1)), &
0.5*(Sal1(i)+Sal1(i+1)))
rhosal(i) = SET_DDENS_DSAL (eos%par,lindens%par,0.5*(Tw1(i)+Tw1(i+1)), &
0.5*(Sal1(i)+Sal1(i+1)))

Victor Stepanenko
committed
! The cycle implements iterations to
! solve the set of nonlinear finite-difference equations
! of k-epsilon parameterization

Victor Stepanenko
committed
do i = 1, M+1
E12(i) = 0.5d0*(E1(i) + E_it1(i))
eps12(i) = 0.5d0*(eps1(i) + eps_it1(i))
enddo

Victor Stepanenko
committed
! E12 = E_it1
! eps12 = eps_it1
! Calculating the stability functions

Victor Stepanenko
committed
lam_eps(i) = lam_eps0
lam_E (i) = lam_E0
CE (i) = CE0
CEt (i) = CEt0

Victor Stepanenko
committed
work(i,1) = E12(i)*E12(i)/((eps12(i) + ACC2)*(eps12(i) + ACC2))* &
& (rhotemp(i) * & ! Assuming Pr = Sc

Victor Stepanenko
committed
& ( Tw1(i+1) - Tw1(i) ) + &
& alsal*rhosal(i) * &

Victor Stepanenko
committed
& ( Sal1(i+1) - Sal1(i)) ) / &
& (h1*ddz(i)) *AL ! ~ bouyancy production of TKE
work(i,2) = E12(i)*E12(i)/((eps12(i) + ACC2)*(eps12(i) + ACC2))* &
& ( (U1(i+1) - U1(i))*(U1(i+1) - U1(i)) + &
& (V1(i+1) - V1(i))*(V1(i+1) - V1(i)) ) / &
& (h1*h1*ddz(i)*ddz(i)) ! ~ shear production of TKE
if (stabfunc%par == 2) then
CE(i) = CE_CANUTO (work(i,2), work(i,1))
CEt(i) = CEt_CANUTO(work(i,2), work(i,1))
elseif (stabfunc%par == 3) then
CE(i) = sqrt(2.d0)*CL_K_KL_MODEL*SMOMENT_GALPERIN(work(i,1))
CEt(i) = sqrt(2.d0)*CL_K_KL_MODEL*SHEAT_GALPERIN(work(i,1))
endif
lam_eps(i) = CE(i)/sigmaeps
lam_E (i) = CE(i)/sigmaE
endif
enddo

Victor Stepanenko
committed
! Calculating eddy diffusivities for kinetic energy and its dissipation
! in the middle between half levels

Victor Stepanenko
committed
! E_mid=0.5d0*(E12(i-1)+E12(i))
! eps_mid=0.5d0*(eps12(i-1)+eps12(i))
work(i,1) = INTERPOLATE_TO_INTEGER_POINT(E12(i-1),E12(i),ddz(i-1),ddz(i)) ! interpolated TKE
work(i,2) = INTERPOLATE_TO_INTEGER_POINT(eps12(i-1),eps12(i),ddz(i-1),ddz(i)) ! interpolated dissipation
xx(1) = INTERPOLATE_TO_INTEGER_POINT(lam_eps(i-1),lam_eps(i),ddz(i-1),ddz(i))
xx(2) = INTERPOLATE_TO_INTEGER_POINT(lam_E(i-1),lam_E(i),ddz(i-1),ddz(i))
k2_mid(i) = work(i,1)*work(i,1)/(work(i,2) + ACCk)
k2_mid(i) = k2_mid(i)*al_it + work(i,1)*work(i,1)/(work(i,2) + ACCk)*(1-al_it)
endif
k5_mid(i) = niu_wat + max(k2_mid(i)*xx(1), min_visc/sigmaeps)
k3_mid(i) = niu_wat + max(k2_mid(i)*xx(2), min_visc/sigmaE)
k2(i) = max(E12(i)*E12(i)/(eps12(i) + ACCk), min_visc/CE(i))
k2t(i) = max(E12(i)*E12(i)/(eps12(i) + ACCk), min_diff/CEt(i))
k2(i) = k2(i) *al_it + &
& max(E12(i)*E12(i)/(eps12(i) + ACCk), min_visc/CE(i))*(1-al_it)
k2t(i) = k2t(i)*al_it + &
& max(E12(i)*E12(i)/(eps12(i) + ACCk), min_diff/CEt(i))*(1-al_it)
k5(i) = niu_wat + k2(i)*lam_eps(i)
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
enddo
do i = 2, M-1
WU_(i) = (E_it1(i+1)-E_it1(i-1))*(U2(i)-U2(i-1))/ &
& (h1*h1*ddz052(i)*ddz(i-1) )
WV_(i) = (E_it1(i+1)-E_it1(i-1))*(V2(i)-V2(i-1))/ &
& (h1*h1*ddz052(i)*ddz(i-1) )
enddo
if (contrgrad) then
do i = 2, M
GRADT = (Tw1(i+1)-Tw1(i-1))/(h1*ddz2(i-1))
GRADU = (U2(i)-U2(i-1))/(h1*ddz(i-1) )
GRADV = (V2(i)-V2(i-1))/(h1*ddz(i-1) )
GAMUN = - CONUV*(WU_(i+1)-WU_(i-1))/(h1*ddz2(i-1))
GAMVN = - CONUV*(WV_(i+1)-WV_(i-1))/(h1*ddz2(i-1))
COGRTN = DTDZH * 0.1
TFR=E_it1(i)/(l(i)*l(i)+ACC) + acc
COGRUN = GAMUN/TFR
COGRVN = GAMVN/TFR
if(gradT < 0.) then
GAMT(i) = - AG(3)*cogrtn
GAMU(i) = - AG(1)*cogrun
GAMV(i) = - AG(2)*cogrvn
else
GAMT(i) = - AG(3)*(AL*GRADT**2+CON0*TFR*COGRTN) / &
& (AL*GRADT+CON0*TFR)
GAMU(i) = - AG(1)*(AL*(CON2*GRADT+(CON2-1.)*GAMT(i))*GRADU + &
& CON1*TFR*COGRUN)/(AL*GRADT+CON1*TFR)
GAMV(i) = - AG(2)*(AL*(CON2*GRADT+(CON2-1.)*GAMT(i))*GRADV + &
& CON1*TFR*COGRVN)/(AL*GRADT+CON1*TFR)
endif
enddo
GAMT(1)=0.
GAMU(1)=0.
GAMV(1)=0.
else
do i = 1, M+1
GAMT(i) = 0.
GAMU(i) = 0.
GAMV(i) = 0.
enddo
endif
do i = 1, M
Gen(i) = ((Um(i+1) - Um(i))*(Um(i+1) - Um(i) + h1*ddz(i)*GAMU(i)) + &
& (Vm(i+1) - Vm(i))*(Vm(i+1) - Vm(i) + h1*ddz(i)*GAMV(i))) / &
& (h1**2*ddz(i)**2)*CE(i)

Victor Stepanenko
committed
! S(i) = - ( 0.5d0*(row(i+1) + row2(i+1) - row(i) - row2(i)) / &
! & (h1*ddz(i)) + GAMT(i))*ALFT*AL*CEt(i)
S(i) = - ( 0.5d0*( rhotemp(i) * & ! Assuming Pr = Sc
& (Tw1(i+1) + Tw2(i+1) - Tw1(i) - Tw2(i)) + &
& alsal*rhosal(i) * &
& (Sal1(i+1) + Sal2(i+1) - Sal1(i) - Sal2(i)) ) / &
& (h1*ddz(i)) + GAMT(i) )*ALFT*AL*CEt(i)

Victor Stepanenko
committed
!S(i) = max(S(i),0._ireals)
! Adding shear production of turbulence by gravitational waves in stable stratification
Gen(i) = Gen(i) - grav_wave_Gill*grav_wave_Gill_const*STEP(-S(i))*S(i)*CE(i)/CEt(i)
Buoyancy1 = S(1)*k2t(1)
dBdz0 = (Buoyancy1 - Buoyancy0)/(0.5*h1*ddz(1))
endif

Victor Stepanenko
committed
! Solution of the equation for TKE (turbulent kinetic energy)
if_tkebc : if (kepsbc%par <= 4) then
!Neumann boundary conditions for TKE
!Boundary condition at the top

Victor Stepanenko
committed
!Open water
if (kepsbc%par < 4) then
FTKES = TKE_FLUX_SHEAR(tau_air,kwe%par,maxTKEinput,kepsbc%par)
elseif (kepsbc%par == 4) then ! The new boundary condition,
! implemented currently only for open water case
FTKES = TKE_FLUX_BUOY(Buoyancy0,roughness0,H_entrainment)
endif
elseif (l1 > 0._ireals .and. l1 <= L0) then

Victor Stepanenko
committed
!Fractional ice
FTKES = b0*TKE_FLUX_SHEAR(tau_i,0._ireals,maxTKEinput,1) + &
& (1.d0-b0)*TKE_FLUX_SHEAR(tau_air,kwe%par,maxTKEinput,kepsbc%par)
elseif (l1 > L0) then
!Complete ice cover
FTKES = TKE_FLUX_SHEAR(tau_i,0._ireals,maxTKEinput,1)

Victor Stepanenko
committed
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
xx(1) = 0.5*(area_int(2)/area_half(1)*k3_mid(2)*dt / &
& (h1**2*ddz(1)*ddz05(1)) + &
& time_deriv*dzeta_05int(1)*dhw_keps/(h1*ddz05(1) ) - &
& time_deriv*dhw0_keps/(h1*ddz05(1) ) )
yy(1) = dt/(h1*ddz(1))
b(1) = - xx(1)
c(1) = - (xx(1) + time_deriv*1.d0) - (abs(S(1)) - S(1))/(TF(1) + ACC)*dt05 - TF(1)*dt
d(1) = - time_deriv*E1(1) - xx(1)*(E1(2) - E1(1)) - &
& area_int(1)/area_half(1)*yy(1)*FTKES - &
& k2t(1)*(S(1) + abs(S(1)))*dt05 - dt*(k2(1)*Gen(1) + Gen_seiches(1))!+eps_it1(i)*dt
! Boundary condition at the bottom
FTKEB = - TKE_FLUX_SHEAR(tau_gr,0._ireals,maxTKEinput,1)
xx(2) = 0.5*(area_int(M)/area_half(M)*k3_mid(M)*dt / &
& (h1**2*ddz(M)*ddz05(M-1) ) - &
& time_deriv*dzeta_05int(M)*dhw_keps/(h1*ddz05(M-1) ) + &
& time_deriv*dhw0_keps/(h1*ddz05(M-1) ) )
yy(2) = dt/(h1*ddz(M))
a(M) = - xx(2)
c(M) = - (xx(2) + time_deriv*1.d0) - (abs(S(M)) - S(M))/(TF(M) + ACC)*dt05 - &
& TF(M)*dt
d(M) = - time_deriv*E1(M) - xx(2)*(E1(M-1) - E1(M)) + &
& area_int(M+1)/area_half(M)*yy(2)*FTKEB - &
& k2t(M)*(S(M) + abs(S(M)))*dt05 - dt*(k2(M)*Gen(M) + Gen_seiches(M)) !+eps_it1(i)*dt
elseif (kepsbc%par == 5) then
!Dirichlet boundary condition for TKE
!Boundary condition at the top

Victor Stepanenko
committed
!Open water
TKES = TKE_WALL_SHEAR(tau_air)
elseif (l1 > 0._ireals .and. l1 <= L0) then

Victor Stepanenko
committed
!Fractional ice
TKES = b0*TKE_WALL_SHEAR(tau_i) + (1.-b0)*TKE_WALL_SHEAR(tau_air)
elseif (l1 > L0) then
!Complete ice cover
TKES = TKE_WALL_SHEAR(tau_i)
endif
b(1) = 0.
c(1) = 1.
d(1) = TKES
! Boundary condition at the bottom
TKEB = TKE_WALL_SHEAR(tau_gr)
a(M) = 0.
c(M) = 1.
d(M) = TKEB
endif if_tkebc
! The coefficients of linear algebraic equations in the internal points of the mesh
do i = 2, M-1
a(i) = time_deriv*dzeta_05int(i)*dhw_keps/(ddz054(i)*h1) - &
& time_deriv*dhw0_keps/(ddz054(i)*h1) - &
& area_int(i)/area_half(i)*k3_mid(i)*dt/(h1**2*ddz(i)*ddz2(i-1) )
b(i) = -time_deriv*dzeta_05int(i)*dhw_keps/(ddz054(i)*h1) + &
& time_deriv*dhw0_keps/(ddz054(i)*h1) - &
& area_int(i+1)/area_half(i)*k3_mid(i+1)*dt/(h1**2*ddz(i)*ddz2(i) )
c(i) = a(i) + b(i) - time_deriv*1.d0-(abs(S(i))-S(i))/(TF(i)+ACC)*dt05 - TF(i)*dt
d(i) = -time_deriv*E1(i)-k2t(i)*(S(i)+abs(S(i)))*dt05 - dt*(k2(i)*Gen(i) + Gen_seiches(i))
d(i) = d(i)-dt*(h1**2*ddz(i)*ddz2(i))**(-1)* &
& (area_int(i+1)/area_half(i)*k3_mid(i+1)*(E1(i+1)-E1(i)))
d(i) = d(i)+dt*(h1**2*ddz(i)*ddz2(i-1))**(-1)* &
& (area_int(i)/area_half(i)*k3_mid(i)*(E1(i)-E1(i-1)))
d(i) = d(i)-time_deriv*dzeta_05int(i)*dhw_keps* &
& (ddz054(i)*h1)**(-1)*(E1(i+1)-E1(i-1))
d(i) = d(i)+time_deriv*dhw0_keps* &
& (ddz054(i)*h1)**(-1)*(E1(i+1)-E1(i-1))
enddo
ind_bound = .false.
call IND_STAB_FACT_DB(a,b,c,1,M,indstab,ind_bound)
if (indstab .eqv. .false.) then
do i = 2, M-1
a(i) = - k3_mid(i)*dt/(h1**2*ddz(i)*ddz2(i-1) )
b(i) = - k3_mid(i+1)*dt/(h1**2*ddz(i)*ddz2(i) )
enddo
endif
call PROGONKA (a,b,c,d,E_it2,1,M)
E_it2(M+1) = 0.
call CHECK_MIN_VALUE(E_it2, M, E_min)

Victor Stepanenko
committed
! C1 is given following Satyanarayana et al., 1999; Aupoix et al., 1989
yy(1) = 0.5*(1. + sign(1._ireals,S(i)))
zz = 0.5*(1. - sign(1._ireals,S(i)))
ceps3 = ceps3_stable*zz + ceps3_unstable*yy(1)

Victor Stepanenko
committed
! ceps3 = 1.14d0
! ceps3 < -0.4 should be used for stable stratification (Burchard, 2002)
! ceps3 = 1. ! following (Kochergin and Sklyar, 1992)
! ceps3 = 1.14d0 ! Baum and Capony (1992)
! Other values for ceps3:
! 0.<ceps3<0.29 for stable and ceps3 = 1.44 for unstable conditions (Rodi, 1987);
! ceps3 >= 1 (Kochergin and Sklyar, 1992)
C_eps1(i) = ceps1
C_eps2(i) = ceps2
C_eps3(i) = ceps3
enddo
elseif (keps_coef == 2) then
do i = 1, M
Re(i) = ACC + (2*E_it1(i)/3.)**2/(eps_it1(i)*niu_wat + ACC)
C1aup(i) = Co/(1.d0 + 0.69d0*(2.d0 - Co)/sqrt(Re(i)))
C_eps1(i) = C1aup(i)
C_eps2(i) = C1aup(i)
C_eps3(i) = C1aup(i)
enddo
endif

Victor Stepanenko
committed
! The solution of the equation for dissipation rate
! dist_surf is the distance from the surface,
! where the b.c. for epsilon are formulated
! dist_surf = 0._ireals !0.5d0*ddz(1)*h1
if_dissbc: if (kepsbc%par <= 4) then
!Neumann boundary conditions for TKE dissipation
!Boundary condition at the top

Victor Stepanenko
committed
if (kepsbc%par < 4) then
! FS = CE0**(0.75d0)*k5(1)*E_it1(1)**(1.5d0)/kar* &
! & (roughness0 + dist_surf)**(-2)
! FS = DISSIPATION_FLUX_SHEAR(tau_air,roughness0)
! Extrapolation of the turbulent diffusivity to the boundary level, assuming
! logarithmic profile
ksurf1(1) = max(niu_wat, CE(1)*k2(1)-kar*sqrt(tau_air/row0)*0.5*ddz(1)*h1)
! ksurf1(1) = kar*sqrt(tau_air/row0)*z0 !roughness0
! ksurf1(1) = CE(1)*k2(1)
Esurf1(1) = E_it1(1) + 0.5*FTKES*sigmaE*h1*ddz(1)/ksurf1(1)
FS = DISSIPATION_FLUX_SHEAR(tau_air,ksurf1(1),Esurf1(1),z0, &
& kwe%par,fetch,maxTKEinput,kepsbc%par,signwaveheight)
elseif (kepsbc%par == 4) then ! The new boundary condition,
! implemented currently only for open water case
Esurf1(1) = E_it1(1)
FS = DISSIPATION_FLUX_BUOY(Buoyancy0,roughness0,H_entrainment,dBdz0,Esurf1(1))
endif
elseif (l1 > 0._ireals .and. l1 <= L0) then

Victor Stepanenko
committed
! FS = CE0**(0.75d0)*k5(1)*E_it1(1)**(1.5d0)/kar* &
! & (0.0._ireals + dist_surf)**(-2) !0.01 m - roughness of ice
! FS = b0*DISSIPATION_FLUX_SHEAR(tau_i,1.d-3) + & ! 1.d-3 roughness of ice
! & (1.d0-b0)*DISSIPATION_FLUX_SHEAR(tau_air,roughness0)
! Extrapolation of the turbulent diffusivity to the boundary level, assuming
! logarithmic profile
ksurf1(1) = max(niu_wat, CE(1)*k2(1) - kar*sqrt(tau_air/row0)*0.5*ddz(1)*h1)
ksurf2(1) = max(niu_wat, CE(1)*k2(1) - kar*sqrt(tau_i/row0)*0.5*ddz(1)*h1)

Victor Stepanenko
committed
! ksurf1(1) = kar*sqrt(tau_air/row0)*z0 !*roughness0
! ksurf2(1) = kar*sqrt(tau_i/row0)*z0 !1.d-3 is the roughness of ice bottom

Victor Stepanenko
committed
! ksurf2(1) = CE(1)*k2(1)
Esurf1(1) = E_it1(1) + 0.5*FTKES*sigmaE*h1*ddz(1)/ksurf1(1)

Victor Stepanenko
committed
Esurf2(1) = E_it1(1) + 0.5*FTKES*sigmaE*h1*ddz(1)/ksurf2(1)
FS = b0*DISSIPATION_FLUX_SHEAR(tau_i,ksurf2(1),Esurf2(1),z0, &
& 0._ireals,fetch,maxTKEinput,1,signwaveheight) + &
& (1.d0-b0)*DISSIPATION_FLUX_SHEAR(tau_air,ksurf1(1),Esurf1(1),z0, &
& kwe%par,fetch,maxTKEinput,kepsbc%par,signwaveheight)

Victor Stepanenko
committed
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
elseif (l1 > L0) then
! FS = DISSIPATION_FLUX_SHEAR(tau_i,1.d-3)
! Extrapolation of the turbulent diffusivity to the boundary level, assuming
! logarithmic profile
ksurf2(1) = max(niu_wat, CE(1)*k2(1) - kar*sqrt(tau_i/row0)*0.5*ddz(1)*h1)
! ksurf2(1) = kar*sqrt(tau_i/row0)*z0 !1.d-3 is the roughness of ice bottom
! ksurf2(1) = CE(1)*k2(1)
Esurf2(1) = E_it1(1) + 0.5*FTKES*sigmaE*h1*ddz(1)/ksurf2(1)
FS = DISSIPATION_FLUX_SHEAR(tau_i,ksurf2(1),Esurf2(1),z0, &
& 0._ireals,fetch,maxTKEinput,1,signwaveheight)
endif
xx(1) = 0.5*(area_int(2)/area_half(1)*k5_mid(2)*dt / &
& (h1**2*ddz(1)*ddz05(1)) + &
& time_deriv*dzeta_05int(1)*dhw_keps/(h1*ddz05(1) ) - &
& time_deriv*dhw0_keps/(h1*ddz05(1) ) )
yy(1) = dt/(h1*ddz(1))
b(1) = - xx(1)
c(1) = - (xx(1) + time_deriv*1.d0) - C_eps2(1)*TF(1)*dt + &
& C_eps3(1)*(- abs(S(1)) + S(1))/(TF(1) + ACC)*dt05
d(1) = - time_deriv*eps1(1) - xx(1)*(eps1(2) - eps1(1)) - &
& area_int(1)/area_half(1)*yy(1)*FS - C_eps3(1) * &
& (abs(S(1)) + S(1))*TF(1)*k2t(1)*dt05 - &
& C_eps1(1)*TF(1)*dt*(k2(1)*Gen(1) + Gen_seiches(1)) !+eps_it1(i)*dt
! Boundary condition at the bottom
! dist_surf = 0._ireals ! ddz(M)/2*h1

Victor Stepanenko
committed
! FB = - CE0**(0.75d0)*k5(M)* &
! & E_it1(M)**(1.5d0)/kar*(0.0._ireals + dist_surf)**(-2)
! FB = - DISSIPATION_FLUX_SHEAR(tau_gr,1.d-3) ! 1.d-3 m - rougnness of bottom

Victor Stepanenko
committed
! Extrapolation of the turbulent diffusivity to the boundary level, assuming
! logarithmic profile
ksurf2(2) = max(niu_wat, CE(M)*k2(M) - kar*sqrt(tau_gr/row0)*0.5*ddz(M)*h1)

Victor Stepanenko
committed
! ksurf2(2) = kar*sqrt(tau_gr/row0)*z0 !1.d-3 is the roughness of bottom
! ksurf2(2) = CE(M)*k2(M)
Esurf2(2) = E_it1(M) - 0.5*FTKEB*sigmaE*h1*ddz(M)/ksurf2(2)
FB = - DISSIPATION_FLUX_SHEAR(tau_gr,ksurf2(2),Esurf2(2),z0, &
& 0._ireals,fetch,maxTKEinput,1,signwaveheight)
xx(2) = 0.5*(area_int(M)/area_half(M)*k5_mid(M)*dt / &
& (h1**2*ddz(M)*ddz05(M-1) ) - &
& time_deriv*dzeta_05int(M)*dhw_keps/(h1*ddz05(M-1) ) + &
& time_deriv*dhw0_keps/(h1*ddz05(M-1) ))
yy(2) = dt/(h1*ddz(M))
a(M) = - xx(2)
c(M) = - (xx(2) + time_deriv*1.d0) - C_eps2(M)*TF(M)*dt + &
& C_eps3(M)*( - abs(S(M)) + S(M))/(TF(M) + ACC)*dt05
d(M) = - time_deriv*eps1(M) - xx(2)*(eps1(M-1) - eps1(M)) + &
& area_int(M+1)/area_half(M)*yy(2)*FB - C_eps3(M) * &
& (abs(S(M)) + S(M))*TF(M)*k2t(M)*dt05 - &
& C_eps1(M)*TF(M)*dt*(k2(M)*Gen(M) + Gen_seiches(M))
elseif (kepsbc%par == 5) then
!Dirichlet boundary conditions for TKE dissipation
!Boundary condition at the top
dist_surf = 0.5*ddz(1)*h1

Victor Stepanenko
committed
!Open water
DISS_TKES = DISSIPATION_WALL_SHEAR(tau_air,Buoyancy0,dist_surf)
elseif (l1 > 0._ireals .and. l1 <= L0) then

Victor Stepanenko
committed
!Fractional ice
DISS_TKES = b0 *DISSIPATION_WALL_SHEAR(tau_i ,Buoyancy0,dist_surf) + &
& (1.-b0)*DISSIPATION_WALL_SHEAR(tau_air,Buoyancy0,dist_surf)
elseif (l1 > L0) then
!Complete ice cover
DISS_TKES = DISSIPATION_WALL_SHEAR(tau_air,Buoyancy0,dist_surf)
endif
b(1) = 0.
c(1) = 1.
d(1) = DISS_TKES
!Boundary condition at the bottom
dist_surf = 0.5*ddz(M)*h1
DISS_TKEB = DISSIPATION_WALL_SHEAR(tau_gr,0._ireals,dist_surf)
a(M) = 0.
c(M) = 1.
d(M) = DISS_TKEB
endif if_dissbc
do i = 2, M-1
a(i) = time_deriv*dzeta_05int(i)*dhw_keps/(ddz054(i)*h1) - &
& time_deriv*dhw0_keps/(ddz054(i)*h1) - &
& area_int(i)/area_half(i)*k5_mid(i)*dt/(h1**2*ddz(i)*ddz2(i-1) )
b(i) = -time_deriv*dzeta_05int(i)*dhw_keps/(ddz054(i)*h1) + &
& time_deriv*dhw0_keps/(ddz054(i)*h1) - &
& area_int(i+1)/area_half(i)*k5_mid(i+1)*dt/(h1**2*ddz(i)*ddz2(i) )
c(i) = a(i) + b(i) - time_deriv*1.d0 - C_eps2(i)*TF(i)*dt + &
& C_eps3(i)*(-abs(S(i))+S(i))/(TF(i)+ACC)*dt05
d(i) = - time_deriv*eps1(i) - C_eps3(i)*(abs(S(i)) + S(i))*TF(i)*k2t(i)*dt05 - &
& C_eps1(i)*TF(i)*dt*(k2(i)*Gen(i) + Gen_seiches(i))
d(i) = d(i)-dt*(h1**2*ddz(i)*ddz2(i))**(-1)* &
& (area_int(i+1)/area_half(i)*k5_mid(i+1)*(eps1(i+1)-eps1(i)))
d(i) = d(i)+dt*(h1**2*ddz(i)*ddz2(i-1))**(-1)* &
& (area_int(i)/area_half(i)*k5_mid(i)*(eps1(i)-eps1(i-1)))
d(i) = d(i)-time_deriv*dzeta_05int(i)*dhw_keps* &
& (ddz054(i)*h1)**(-1)*(eps1(i+1)-eps1(i-1))
d(i) = d(i)+time_deriv*dhw0_keps* &
& (ddz054(i)*h1)**(-1)*(eps1(i+1)-eps1(i-1))
enddo
ind_bound = .false.
call IND_STAB_FACT_DB(a,b,c,1,M,indstab,ind_bound)
if (indstab .eqv. .false.) then
do i = 2, M-1
a(i) = - k5_mid(i)*dt/(h1**2*ddz(i)*ddz2(i-1))
b(i) = - k5_mid(i+1)*dt/(h1**2*ddz(i)*ddz2(i))
enddo
endif
call PROGONKA (a,b,c,d,eps_it2,1,M)
eps_it2(M+1) = 0.
call CHECK_MIN_VALUE(eps_it2, M, eps_min)
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
if (iterat) then
if ((dE_it > dE_it_max .or. deps_it > deps_it_max) .and. iter_t < maxiter) then
do i = 1, M+1
E_it1(i) = E_it1(i) * al_it + E_it2(i) * (1-al_it)
eps_it1(i) = eps_it1(i) * al_it + eps_it2(i) * (1-al_it)
enddo
iter = iter + 1
iter_t = iter_t + 1
if (iter == 8) iter = 0
else
if (dE_it < 1.E-6 .and. deps_it < 1.E-7) then
do i = 1, M+1
E2(i) = E_it2(i)
eps2(i) = eps_it2(i)
enddo
iter = 0
iter_t = 0
nl = 2
cycle_keps = .false.
else
do i = 1, M+1
eps_it1(i) = eps1(i)
E_it1(i) = E1(i)
enddo
iterat = .false.
iter_t = 0
nl = 1
endif
endif
else
do i = 1, M+1
E2(i) = E_it2(i)
eps2(i) = eps_it2(i)
enddo

Victor Stepanenko
committed
! E2 = E_it1 * al_it + E_it2 * (1-al_it)
! eps2 = eps_it1 * al_it + eps_it2 * (1-al_it)
iter = 0
iter_t = 0
cycle_keps = .false.
endif

Victor Stepanenko
committed
! if (dE_it<1.E-10.or.deps_it<1.E-13) nl=3
enddo iterat_keps
if (.not.init_keps(ix,iy)) then

Victor Stepanenko
committed
! K-epsilon solver is implemented in initialization mode, i.e.
! to obtain initial profiles of E1 and eps1.
! Time derivatives in k-epslilon equations are neglected in this case.
do i = 1, M+1
E1(i) = E2(i)
eps1(i) = eps2(i)
enddo
init_keps(ix,iy) = .true.
else

Victor Stepanenko
committed
! The case when k-epsilon solver is implemented in evolution mode, i.e.
! E2 and epsilon2 are obtained at the next time step.
next_tstep_keps = .true.
endif
enddo keps_mode

Victor Stepanenko
committed
! Smoothing of the k and epsilon profiles
if (smooth) then
call VSMOP_LAKE(E2,E2,lbound(E2,1),ubound(E2,1),1,M,vdamp,perdam)
call CHECK_MIN_VALUE(E2, M, E_min)
call VSMOP_LAKE(eps2,eps2,lbound(eps2,1),ubound(eps2,1),1,M,vdamp,perdam)
call CHECK_MIN_VALUE(eps2, M, eps_min)
endif
do i = 1, M+1
E12(i) = 0.5d0*(E1(i) + E2(i))
eps12(i) = 0.5d0*(eps1(i) + eps2(i))
enddo

Victor Stepanenko
committed
if (stabfunc%par /= 1) then
do i = 1, M
work(i,1) = E2(i)*E2(i)/((eps2(i) + ACC2)*(eps2(i) + ACC2))* &
& (rhotemp(i) * & ! Assuming Pr = Sc

Victor Stepanenko
committed
& ( Tw2(i+1) - Tw2(i) ) + &
& alsal*rhosal(i) * &

Victor Stepanenko
committed
& ( Sal2(i+1) - Sal2(i)) ) / &
& (h1*ddz(i)) *AL ! ~ bouyancy production of TKE
work(i,2) = E2(i)*E2(i)/((eps2(i) + ACC2)*(eps2(i) + ACC2))* &
& ( (U2(i+1) - U2(i))*(U2(i+1) - U2(i)) + &
& (V2(i+1) - V2(i))*(V2(i+1) - V2(i)) ) / &
& (h1*h1*ddz(i)*ddz(i)) ! ~ shear production of TKE
if (stabfunc%par == 2) then
CEt(i) = CEt_CANUTO(work(i,2), work(i,1))
elseif (stabfunc%par == 3) then
CEt(i) = sqrt(2.d0)*CL_K_KL_MODEL*SHEAT_GALPERIN(work(i,1))

Victor Stepanenko
committed
endif
enddo
endif

Victor Stepanenko
committed
KT(i) = max(CEt(i)*E2(i)**2/(eps2(i) + ACCk),min_diff)*cw_m_row0 !E2, eps2

Victor Stepanenko
committed
! Diagnostic calculation