Skip to content
Snippets Groups Projects
Commit 34d85a2f authored by Victor Stepanenko's avatar Victor Stepanenko
Browse files

1.Corrections of input data for Kitinen river experiment 2. Dirichlet boundary...

1.Corrections of input data for Kitinen river experiment 2. Dirichlet boundary conditions and code reformatting in k-epsilon closure
parent a9cef725
No related branches found
No related tags found
No related merge requests found
This diff is collapsed.
......@@ -1048,13 +1048,13 @@
2018,5.9030018432,4,200,-999,-999,-999,-999,-999,-999,-999,-999
2018,-0.3163131904,4,200,-999,-999,-999,-999,-999,-999,-999,-999
2018,-3.9275283712,4,200,-999,-999,-999,-999,-999,-999,-999,-999
2018,-0.450061900799999,4,200,-999,-999,-999,-999,-999,-999,-999,-999
2018,-0.4500619008,4,200,-999,-999,-999,-999,-999,-999,-999,-999
2018,-3.6600309504,4,200,-999,-999,-999,-999,-999,-999,-999,-999
2018,-0.3163131904,4,200,-999,-999,-999,-999,-999,-999,-999,-999
2018,-2.3894182016,4,200,-999,-999,-999,-999,-999,-999,-999,-999
2018,-1.5869259392,4,200,-999,-999,-999,-999,-999,-999,-999,-999
2018,-1.185679808,4,200,-999,-999,-999,-999,-999,-999,-999,-999
2018,0.686802137600001,4,200,-999,-999,-999,-999,-999,-999,-999,-999
2018,0.6868021376,4,200,-999,-999,-999,-999,-999,-999,-999,-999
2018,5.7023787776,4,200,-999,-999,-999,-999,-999,-999,-999,-999
2018,4.4317660288,4,200,-999,-999,-999,-999,-999,-999,-999,-999
2018,5.6355044224,4,200,-999,-999,-999,-999,-999,-999,-999,-999
......
This diff is collapsed.
......@@ -210,9 +210,9 @@ init_T 1
# the first N values are polynomial coefficients, where N stands after
# 'effl_outflow' keyword
#
area_lake 4.E+4
area_lake 4.E+7
#area_lake 9.E+10
cellipt 1.
cellipt 1000.
lakeform 2
trib_inflow -9999.
effl_outflow 1
......
......@@ -131,7 +131,7 @@ lindens 0
nmeltpoint 1
Turbpar 2
stabfunc 1
kepsbc 1
kepsbc 5
botfric 2
soiltype 5
soil_depth 10.
......@@ -149,7 +149,7 @@ salsoil 0
horvisc 0
backdiff 0
dyn_pgrad 5
nManning 5.5E-2
nManning 5.2E-2
zero_model 0
thermokarst_meth_prod 0.
soil_meth_prod 1.
......
......@@ -6,7 +6,8 @@ use INOUT, only : CHECK_UNIT
use LAKE_DATATYPES, only : ireals, iintegers
contains
SUBROUTINE K_EPSILON(ix, iy, nx, ny, year, month, day, hour, kor, a_veg, c_veg, h_veg, dt, &
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:
......@@ -15,8 +16,6 @@ SUBROUTINE K_EPSILON(ix, iy, nx, ny, year, month, day, hour, kor, a_veg, c_veg,
! 2) Semi-empirical formulations for ED;
! 2) k-epsilon parameterization for ED, including Kolmogorov relation.
use OMP_LIB
use COMPARAMS, only: &
& num_ompthr
......@@ -159,12 +158,6 @@ real(kind=ireals) :: b(vector_length)
real(kind=ireals) :: c(vector_length)
real(kind=ireals) :: d(vector_length)
!real(kind=ireals) :: am(vector_length,2,2)
!real(kind=ireals) :: bm(vector_length,2,2)
!real(kind=ireals) :: cm(vector_length,2,2)
!real(kind=ireals) :: ym(vector_length,2)
!real(kind=ireals) :: dm(vector_length,2)
real(kind=ireals) :: AG(5)
real(kind=ireals) :: dt05
......@@ -182,8 +175,8 @@ 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
real(kind=ireals) :: FS, FTKES
real(kind=ireals) :: FB, FTKEB
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
......@@ -207,7 +200,7 @@ real(kind=ireals), allocatable :: work(:,:)
real(kind=ireals), allocatable :: rhotemp(:), rhosal(:)
! Integers
integer(kind=iintegers), parameter :: maxiter = 35 !35
integer(kind=iintegers), parameter :: maxiter = 35
integer(kind=iintegers) :: iter, iter_t
integer(kind=iintegers) :: nl
integer(kind=iintegers) :: i, j, k
......@@ -245,41 +238,14 @@ real(kind=ireals), external :: DZETA
SAVE
!call TIMEC(0)
firstcall_if : if (firstcall) then
!$OMP SINGLE
if (turb_out%par == 1) then
! write (numt,'(i1)') Turbpar%par
! open(2114,file=path(1:len_trim(path))//'results/'// &
! & 'err_progon.dat', status = 'unknown')
! open(2115,file=path(1:len_trim(path))// &
! & 'results/'//'E-eps.dat', &
! & status = 'unknown')
! open(2116,file=path(1:len_trim(path))// &
! & 'results/'//'E-eps1.dat', &
! & status = 'unknown')
! open (2117,file=path(1:len_trim(path))//'results/' &
! & //'turb'//numt//'.dat', status = 'unknown')
! open (2118,file=path(1:len_trim(path))//'results/' &
! & //'temp'//numt//'.dat', status = 'unknown')
! open (2119,file=path(1:len_trim(path))//'results/' &
! & //'uv'//numt//'.dat', status = 'unknown')
! open (2120,file=path(1:len_trim(path))//'results/' &
! & //'KT'//numt//'.dat', status = 'unknown')
! open (2121,file=path(1:len_trim(path))//'results/' &
! & //'dE1'//numt//'.dat', status = 'unknown')
! open (2122,file=path(1:len_trim(path))//'results/' &
! & //'dE2'//numt//'.dat', status = 'unknown')
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', &
& status = 'unknown')
write(unit=nunit_,fmt=*) 'Errors of iterative process in k-epsilon solver'
! open (unit=1, file='test.dat', status='unknown')
endif
allocate (init_keps(1:nx, 1:ny) )
......@@ -304,7 +270,8 @@ firstcall_if : if (firstcall) then
z0 = 1.d-2
eps_min = 1.d-12 ! The range for dissipation rate values 10**(-11) - 10**(-6) m**2/s**3 (Wuest and Lorke, 2003)
eps_min = 1.d-12 ! The range for dissipation rate values
! 10**(-11) - 10**(-6) m**2/s**3 (Wuest and Lorke, 2003)
E_min = sqrt(eps_min*1.E-5*lamw0*(cw_m_row0*CEt0)**(-1) ) ! This expression ensures that
! eddy diffusivity in the decaying
! turbulence (E -> Emin, eps -> eps_min)
......@@ -335,7 +302,6 @@ firstcall_if : if (firstcall) then
! stabfunc = 3 - stability functions CE and CEt are dependent
! only on buoyancy (Galperin et al., 1988)
! FRICTION: VEGETATION
if (h_veg>0) then
......@@ -349,21 +315,15 @@ firstcall_if : if (firstcall) then
if (h1-dzeta_int(i)*h1 < h_veg) veg_sw(i) = 1.
enddo
!$OMP END SINGLE
endif firstcall_if
dt05 = 0.5d0*dt
!$OMP SINGLE
allocate (work(1:M+1,1:2))
allocate (rhotemp(1:M+1), rhosal(1:M+1))
!$OMP END SINGLE
!$OMP MASTER
! Implementation of the Krank-Nikolson numerical scheme
! Implementation of the Crank-Nicolson numerical scheme
! for k-epsilon parameterization using simple iterations to
! solve the set of nonlinear finite-difference equations
......@@ -382,7 +342,8 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
! Seiche energy evolution, TKE production due to seiche dissipation
if (seiches_Goudsmit == 1) then
call SEICHE_ENERGY &
& (M,bathymwater,ls,tau_air-tau_wav,sqrt(uwind*uwind+vwind*vwind),vol,dt,Eseiches,Eseiches_diss)
& (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)
......@@ -410,14 +371,11 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
!enddo
H_entrainment = dzeta_05int(i_entrain) * h1
!$OMP END MASTER
!$OMP FLUSH (i_entrain,H_entrainment,work)
! 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)
......@@ -426,7 +384,6 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
dhw_keps = dhw / real(nstep_keps%par)
dhw0_keps = dhw0 / real(nstep_keps%par)
!$OMP DO
! 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.)
......@@ -456,7 +413,6 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
time_deriv = 1
endif
!$OMP DO
do i = 1, M+1
E_it1(i) = E1(i)
eps_it1(i) = eps1(i)
......@@ -465,7 +421,6 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
k2t(i) = 0._ireals
enddo
iter_t = 0
cycle_keps = .true.
......@@ -482,27 +437,21 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
rhosal (i) = SET_DDENS_DSAL(eos%par,lindens%par,Tw1(i),Sal1(i))
enddo
!call TIMEC(1)
iterat_keps: do while (cycle_keps)
! The cycle implements iterations to
! solve the set of nonlinear finite-difference equations
! of k-epsilon parameterization
!$OMP DO
do i = 1, M+1
E12(i) = 0.5d0*(E1(i) + E_it1(i))
eps12(i) = 0.5d0*(eps1(i) + eps_it1(i))
enddo
! E12 = E_it1
! eps12 = eps_it1
! Calculating the stability functions
!$OMP DO
do i = 1, M
if (stabfunc%par == 1) then
lam_eps(i) = lam_eps0
......@@ -532,11 +481,8 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
endif
enddo
! call TIMEC(2)
! Calculating eddy diffusivities for kinetic energy and its dissipation
! in the middle between half levels
!$OMP DO
do i = 2, M
! E_mid=0.5d0*(E12(i-1)+E12(i))
! eps_mid=0.5d0*(eps12(i-1)+eps12(i))
......@@ -552,8 +498,6 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
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)
enddo
! call TIMEC(3)
!$OMP DO
do i = 1, M
if (iter_t == 0) then
k2(i) = max(E12(i)*E12(i)/(eps12(i) + ACCk), min_visc/CE(i))
......@@ -568,18 +512,14 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
l(i) = E12(i)**(1.5d0)/(eps12(i) + ACC)
TF(i) = sqrt(E12(i))/(l(i) + ACC)
enddo
! call TIMEC(4)
!$OMP DO
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
! call TIMEC(5)
if (contrgrad) then
!$OMP SINGLE ! the loop is not parallelized because the option contrgrad is normally not used
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) )
......@@ -606,9 +546,7 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
GAMT(1)=0.
GAMU(1)=0.
GAMV(1)=0.
!$OMP END SINGLE
else
!$OMP DO
do i = 1, M+1
GAMT(i) = 0.
GAMU(i) = 0.
......@@ -616,10 +554,6 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
enddo
endif
! call TIMEC(6)
! zz = OMP_GET_WTIME()
!$OMP DO
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))) / &
......@@ -631,14 +565,11 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
& alsal*rhosal(i) * &
& (Sal1(i+1) + Sal2(i+1) - Sal1(i) - Sal2(i)) ) / &
& (h1*ddz(i)) + GAMT(i) )*ALFT*AL*CEt(i)
!S(i) = max(S(i),0._ireals)
F(i) = Gen(i) + S(i)
! 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)
enddo
! print*, OMP_GET_WTIME() - zz, OMP_GET_THREAD_NUM()
! Sensitivity test
! S(1) = min(S(1), 0._ireals)
if (iter_t == 0) then
Buoyancy1 = S(1)*k2t(1)
......@@ -647,10 +578,11 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
! 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
!$OMP SECTIONS
!$OMP SECTION
if (l1 == 0) then
!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,
......@@ -658,9 +590,11 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
FTKES = TKE_FLUX_BUOY(Buoyancy0,roughness0,H_entrainment)
endif
elseif (l1 > 0 .and. l1 <= L0) then
!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)
endif
xx(1) = 0.5*(area_int(2)/area_half(1)*k3_mid(2)*dt / &
......@@ -673,10 +607,7 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
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
!$OMP SECTION
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) ) - &
......@@ -689,13 +620,31 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
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
if (l1 == 0) then
!Open water
TKES = TKE_WALL_SHEAR(tau_air)
elseif (l1 > 0 .and. l1 <= L0) then
!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
!$OMP END SECTIONS
! call TIMEC(4)
! The coefficients of linear algebraic equations in the internal points of the mesh
!$OMP DO
do i = 2, M-1
a(i) = time_deriv*dzeta_05int(i)*dhw_keps/(ddz054(i)*h1) - &
& time_deriv*dhw0_keps/(ddz054(i)*h1) - &
......@@ -718,33 +667,25 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
& (ddz054(i)*h1)**(-1)*(E1(i+1)-E1(i-1))
enddo
ind_bound = .false.
!$OMP MASTER
call IND_STAB_FACT_DB(a,b,c,1,M,indstab,ind_bound)
!$OMP END MASTER
!$OMP FLUSH(indstab)
if (indstab .eqv. .false.) then
!$OMP DO
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
!$OMP MASTER
call PROGONKA (a,b,c,d,E_it2,1,M)
E_it2(M+1) = 0.
call CHECK_MIN_VALUE(E_it2, M, E_min)
dE_it = maxval(abs(E_it1-E_it2))
!$OMP END MASTER
!$OMP FLUSH (E_it2,dE_it)
! C1 is given following Satyanarayana et al., 1999; Aupoix et al., 1989
if (keps_coef == 1) then
!$OMP DO PRIVATE (yy,zz,ceps3)
do i = 1, M
yy(1) = 0.5*(1. + sign(1._ireals,S(i)))
zz = 0.5*(1. - sign(1._ireals,S(i)))
......@@ -761,7 +702,6 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
C_eps3(i) = ceps3
enddo
elseif (keps_coef == 2) then
!$OMP DO
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)))
......@@ -771,8 +711,6 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
enddo
endif
! call TIMEC(5)
! The solution of the equation for dissipation rate
! dist_surf is the distance from the surface,
......@@ -780,9 +718,9 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
! 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
!$OMP SECTIONS
!$OMP SECTION
if (l1 == 0) then
if (kepsbc%par < 4) then
! FS = CE0**(0.75d0)*k5(1)*E_it1(1)**(1.5d0)/kar* &
......@@ -820,7 +758,8 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
! ksurf2(1) = CE(1)*k2(1)
Esurf1(1) = E_it1(1) + 0.5*FTKES*sigmaE*h1*ddz(1)/ksurf1(1)
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) + &
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)
elseif (l1 > L0) then
......@@ -833,7 +772,8 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
! 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)
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)) + &
......@@ -849,7 +789,6 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
& C_eps1(1)*TF(1)*dt*(k2(1)*Gen(1) + Gen_seiches(1)) !+eps_it1(i)*dt
! Boundary condition at the bottom
!$OMP SECTION
! dist_surf = 0._ireals ! ddz(M)/2*h1
! FB = - CE0**(0.75d0)*k5(M)* &
......@@ -863,7 +802,8 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
! 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)
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) ) + &
......@@ -876,11 +816,32 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
& 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))
!$OMP END SECTIONS
! call TIMEC(6)
elseif (kepsbc%par == 5) then
!Dirichlet boundary conditions for TKE dissipation
!Boundary condition at the top
dist_surf = 0.5*ddz(1)*h1
if (l1 == 0) then
!Open water
DISS_TKES = DISSIPATION_WALL_SHEAR(tau_air,Buoyancy0,dist_surf)
elseif (l1 > 0 .and. l1 <= L0) then
!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
!$OMP DO
do i = 2, M-1
a(i) = time_deriv*dzeta_05int(i)*dhw_keps/(ddz054(i)*h1) - &
& time_deriv*dhw0_keps/(ddz054(i)*h1) - &
......@@ -905,48 +866,32 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
& (ddz054(i)*h1)**(-1)*(eps1(i+1)-eps1(i-1))
enddo
ind_bound = .false.
!$OMP MASTER
call IND_STAB_FACT_DB(a,b,c,1,M,indstab,ind_bound)
!$OMP END MASTER
!$OMP FLUSH (indstab)
if (indstab .eqv. .false.) then
!$OMP DO
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
!$OMP MASTER
call PROGONKA (a,b,c,d,eps_it2,1,M)
eps_it2(M+1) = 0.
call CHECK_MIN_VALUE(eps_it2, M, eps_min)
deps_it = maxval(abs(eps_it1-eps_it2))
!$OMP END MASTER
!$OMP FLUSH(eps_it2,deps_it)
! call TIMEC(7)
! write(*,*) 'iter_t = ', iter_t
if (iterat) then
if ((dE_it > dE_it_max .or. deps_it > deps_it_max) .and. iter_t < maxiter) then
!$OMP DO
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
!$OMP MASTER
iter = iter + 1
iter_t = iter_t + 1
if (iter == 8) iter = 0
!$OMP END MASTER
!$OMP FLUSH (iter,iter_t)
else
if (dE_it < 1.E-6 .and. deps_it < 1.E-7) then
!$OMP DO
do i = 1, M+1
E2(i) = E_it2(i)
eps2(i) = eps_it2(i)
......@@ -956,7 +901,6 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
nl = 2
cycle_keps = .false.
else
!$OMP DO
do i = 1, M+1
eps_it1(i) = eps1(i)
E_it1(i) = E1(i)
......@@ -967,7 +911,6 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
endif
endif
else
!$OMP DO
do i = 1, M+1
E2(i) = E_it2(i)
eps2(i) = eps_it2(i)
......@@ -982,14 +925,11 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
enddo iterat_keps
! call TIMEC(8)
if (.not.init_keps(ix,iy)) then
! 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.
!$OMP DO
do i = 1, M+1
E1(i) = E2(i)
eps1(i) = eps2(i)
......@@ -1005,17 +945,13 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
! Smoothing of the k and epsilon profiles
!$OMP MASTER
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
!$OMP END MASTER
!$OMP FLUSH (E2,eps2)
!$OMP DO
do i = 1, M+1
E12(i) = 0.5d0*(E1(i) + E2(i))
eps12(i) = 0.5d0*(eps1(i) + eps2(i))
......@@ -1041,14 +977,12 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
enddo
endif
!$OMP DO
do i = 1, M
KT(i) = lam_T*max(CEt(i)*E2(i)**2/(eps2(i) + ACCk),min_diff)*cw_m_row0 !E2, eps2
enddo
! Diagnostic calculation
!Turbulent transport of TKE
!$OMP DO
do i = 2, M-1
TKE_turb_trans(i) = 1.d0/(area_half(i)*h1*h1*ddz(i))* &
& (area_int(i+1)*k3_mid(i+1)*0.5*(E2(i+1) + E1(i+1) - E2(i ) - E1(i) ) / ddz05(i) - &
......@@ -1063,13 +997,11 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
& area_int (M) *k3_mid(M)*0.5 *(E2(M) + E1(M) - E2(M-1) - E1(M-1)) / ddz05(M-1) )
! Diagnostic calculations
!$OMP DO
do i = 1, M
S (i) = S (i)*k2t(i)
Gen(i) = Gen(i)*k2(i)
enddo
S_integr_positive = 0._ireals
S_integr_negative = 0._ireals
Gen_integr = 0._ireals
......@@ -1084,8 +1016,6 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
E_integr = 0._ireals
TF_integr = 0._ireals
!$OMP SINGLE
do i = 1, i_entrain !M
xx(1) = ddz(i)*h1
yy(1) = 0.5*(1. + sign(1._ireals,S(i)))
......@@ -1171,13 +1101,9 @@ allocate (rhotemp(1:M+1), rhosal(1:M+1))
deallocate (work)
deallocate (rhotemp, rhosal)
!$OMP END SINGLE
!$OMP DO
do i = 1, M+1
E1(i) = E2(i)
eps1(i) = eps2(i)
! write(*,*) 'hihi', i, OMP_GET_THREAD_NUM()
enddo
! write(*,*) 'K_EPSILON finished', nstep
......@@ -1683,18 +1609,71 @@ select case (nbc)
! xx = min(kwe*(momentum_flux/row0)**1.5,maxTKEinput)
Esurf = momentum_flux/row0/(CE0**0.5)*const4**(2./3.) !Esurf_in TKE at the surface
xx = kwe*(momentum_flux/row0)**1.5
DISSIPATION_FLUX_SHEAR = const2*(momentum_flux/row0)**0.5*kar*z0s*const4**(1./3.)* &
DISSIPATION_FLUX_SHEAR = const2*(momentum_flux/row0)**0.5* &
& kar*z0s*const4**(1./3.)* &
& (const3*xx + kar*Esurf**1.5) / (z0s*z0s + small_number)
case(3)
! Generalized b.c. for sheared flow with wave breaking (Burchard, 2002)
const4 = (1.5*sigmaE)**0.5*CE0**0.25*kwe
const5 = sigmaeps**(-1)*(a + const4)**(1./3.)*((a + const4)/(z0s + small_number) + alpham*const4)
const5 = sigmaeps**(-1)*(a + const4)**(1./3.)*((a + const4)/ &
& (z0s + small_number) + alpham*const4)
DISSIPATION_FLUX_SHEAR = const5*(momentum_flux/row0)**2
end select
if (firstcall) firstcall = .false.
END FUNCTION DISSIPATION_FLUX_SHEAR
!> The function TKE_WALL_SHEAR calculates the turbulent kinetic energy
!! at the surface in the case of shear lorarithmic flow, i.e. no buoyancy flux but wind stress
!! given at the surface.
FUNCTION TKE_WALL_SHEAR(momentum_flux)
use PHYS_CONSTANTS, only : &
& row0
use TURB_CONST, only : CE0
implicit none
real(kind=ireals) :: TKE_WALL_SHEAR
! Input variables
real(kind=ireals), intent(in) :: momentum_flux !> Momentum flux at the wall (surface), N/m**2
!Local variables
TKE_WALL_SHEAR = (momentum_flux/row0)/sqrt(CE0)
END FUNCTION TKE_WALL_SHEAR
!> The function DISSIPATION_WALL_SHEAR calculates the turbulent kinetic energy
!! dissipation rate at the surface in the case of shear logarithmic flow,
!! i.e. no buoyancy flux but wind stress given at the surface.
FUNCTION DISSIPATION_WALL_SHEAR(momentum_flux,buoyancy_flux,dist)
use PHYS_CONSTANTS, only : &
& row0, kappa
implicit none
real(kind=ireals) :: DISSIPATION_WALL_SHEAR
! Input variables
real(kind=ireals), intent(in) :: momentum_flux !> Momentum flux at the wall (surface), N/m**2
real(kind=ireals), intent(in) :: buoyancy_flux !> Bouyancy flux at the wall (surface), m**2/s**3
real(kind=ireals), intent(in) :: dist !> Distance from the wall (surface), m
! Local variables
!logical, save :: firstcall = .true.
!if (firstcall) then
!endif
DISSIPATION_WALL_SHEAR = (momentum_flux/row0)**1.5/(kappa*dist) + buoyancy_flux
!if (firstcall) firstcall = .false.
END FUNCTION DISSIPATION_WALL_SHEAR
SUBROUTINE ED_TEMP_HOSTETLER(wind,frict_vel,zref,phi180,levels,row,ed_temp,N_levels)
......@@ -1821,9 +1800,11 @@ do i = 1, N_levels
else
Richardson = &
& (- 1. + sqrt(1. + &
& max(Brunt_Vaisala2_const*Brunt_Vaisala2*kappa**2*(levels(N_levels+1) - level)**2/ &
& max(Brunt_Vaisala2_const*Brunt_Vaisala2*kappa**2* &
& (levels(N_levels+1) - level)**2/ &
& (w_star**2 + small_value), -1._ireals) ))/Ri_const2
ed_temp(i) = kappa*w_star*(levels(N_levels+1) - level)/(PrandtL0*(1.+Ri_const1*Richardson**2) )
ed_temp(i) = kappa*w_star*(levels(N_levels+1) - level)/ &
& (PrandtL0*(1.+Ri_const1*Richardson**2) )
endif
enddo
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment