FORCING: Adjust forcing modes for domain smaller than unity.
Signed-off-by: Grzegorz Kowal <grzegorz@amuncode.org>
This commit is contained in:
parent
c06bce0ee0
commit
ec7d12260e
@ -170,8 +170,9 @@ module forcing
|
||||
subroutine initialize_forcing(verbose, status)
|
||||
|
||||
use constants , only : pi2
|
||||
use coordinates, only : xlen, ylen, zlen
|
||||
use helpers , only : print_message
|
||||
use parameters, only : get_parameter
|
||||
use parameters , only : get_parameter
|
||||
use random , only : randuni, randnorz
|
||||
|
||||
implicit none
|
||||
@ -186,6 +187,11 @@ module forcing
|
||||
integer :: i, j, l, k2
|
||||
#if NDIMS == 3
|
||||
integer :: k = 0
|
||||
#endif /* NDIMS == 3 */
|
||||
integer :: kxs, kxm
|
||||
integer :: kys, kym
|
||||
#if NDIMS == 3
|
||||
integer :: kzs, kzm
|
||||
#endif /* NDIMS == 3 */
|
||||
real(kind=8) :: kl2, ku2, kv2, kv
|
||||
real(kind=8) :: fa, fi, uu
|
||||
@ -313,6 +319,14 @@ module forcing
|
||||
end if
|
||||
kf = 2.0d+00
|
||||
end if
|
||||
if (kf * min(xlen, ylen, zlen) < 1.0d+00) then
|
||||
if (verbose) then
|
||||
write(*,*)
|
||||
write(*,wfmt) "'driving_wavenumber' is too small for the domain!"
|
||||
kf = 1.0d+00 / min(xlen, ylen, zlen)
|
||||
write(*,wfmt) "Resetting it to the minimum allowed value!"
|
||||
end if
|
||||
end if
|
||||
if (kc <= 0.0d+00) then
|
||||
if (verbose) then
|
||||
write(*,*)
|
||||
@ -358,6 +372,19 @@ module forcing
|
||||
!
|
||||
kmax = ceiling(ku)
|
||||
|
||||
! get the wavenumber step and the maximum value due to the domain size
|
||||
!
|
||||
kxs = max(1, int(1.0d+00 / xlen))
|
||||
kys = max(1, int(1.0d+00 / ylen))
|
||||
#if NDIMS == 3
|
||||
kzs = max(1, int(1.0d+00 / zlen))
|
||||
#endif /* NDIMS == 3 */
|
||||
kxm = kmax - mod(kmax, kxs)
|
||||
kym = kmax - mod(kmax, kys)
|
||||
#if NDIMS == 3
|
||||
kzm = kmax - mod(kmax, kzs)
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
! allocate arrays to count the modes with the same wave number
|
||||
!
|
||||
allocate(kcount(NDIMS * kmax**2), stat = status)
|
||||
@ -372,9 +399,10 @@ module forcing
|
||||
! driving modes
|
||||
!
|
||||
#if NDIMS == 3
|
||||
do i = 1, kmax
|
||||
do j = -kmax, kmax
|
||||
do k = -kmax, kmax
|
||||
do i = 0, kxm, kxs
|
||||
if (i > 0) then
|
||||
do j = - kym, kym, kys
|
||||
do k = - kzm, kzm, kzs
|
||||
k2 = i**2 + j**2 + k**2
|
||||
kv2 = 1.0d+00 * k2
|
||||
if (kl2 <= kv2 .and. kv2 <= ku2) then
|
||||
@ -383,9 +411,11 @@ module forcing
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
do j = 1, kmax
|
||||
do k = -kmax, kmax
|
||||
do j = 0, kym, kys
|
||||
if (j > 0) then
|
||||
do k = - kzm, kzm, kzs
|
||||
k2 = j**2 + k**2
|
||||
kv2 = 1.0d+00 * k2
|
||||
if (kl2 <= kv2 .and. kv2 <= ku2) then
|
||||
@ -393,18 +423,22 @@ module forcing
|
||||
kcount(k2) = kcount(k2) + 2
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
do k = 1, kmax
|
||||
do k = 0, kzm, kzs
|
||||
if (k > 0) then
|
||||
k2 = k**2
|
||||
kv2 = 1.0d+00 * k2
|
||||
if (kl2 <= kv2 .and. kv2 <= ku2) then
|
||||
nmodes = nmodes + 1
|
||||
kcount(k2) = kcount(k2) + 2
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
#else /* NDIMS == 3 */
|
||||
do i = 1, kmax
|
||||
do j = -kmax, kmax
|
||||
do i = 0, kxm, kxs
|
||||
if (i > 0) then
|
||||
do j = - kym, kym, kys
|
||||
k2 = i**2 + j**2
|
||||
kv2 = 1.0d+00 * k2
|
||||
if (kl2 <= kv2 .and. kv2 <= ku2) then
|
||||
@ -412,14 +446,17 @@ module forcing
|
||||
kcount(k2) = kcount(k2) + 2
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
do j = 1, kmax
|
||||
do j = 0, kym, kys
|
||||
if (j > 0) then
|
||||
k2 = j**2
|
||||
kv2 = 1.0d+00 * k2
|
||||
if (kl2 <= kv2 .and. kv2 <= ku2) then
|
||||
nmodes = nmodes + 1
|
||||
kcount(k2) = kcount(k2) + 2
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
@ -435,9 +472,10 @@ module forcing
|
||||
!
|
||||
l = 0
|
||||
#if NDIMS == 3
|
||||
do i = 1, kmax
|
||||
do j = -kmax, kmax
|
||||
do k = -kmax, kmax
|
||||
do i = 0, kxm, kxs
|
||||
if (i > 0) then
|
||||
do j = - kym, kym, kys
|
||||
do k = - kzm, kzm, kzs
|
||||
k2 = i**2 + j**2 + k**2
|
||||
kv2 = 1.0d+00 * k2
|
||||
if (kl2 <= kv2 .and. kv2 <= ku2) then
|
||||
@ -448,9 +486,11 @@ module forcing
|
||||
end if
|
||||
end do
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
do j = 1, kmax
|
||||
do k = -kmax, kmax
|
||||
do j = 0, kym, kys
|
||||
if (j > 0) then
|
||||
do k = - kzm, kzm, kzs
|
||||
k2 = j**2 + k**2
|
||||
kv2 = 1.0d+00 * k2
|
||||
if (kl2 <= kv2 .and. kv2 <= ku2) then
|
||||
@ -460,8 +500,10 @@ module forcing
|
||||
kmodes(l,3) = k
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
do k = 1, kmax
|
||||
do k = 0, kzm, kzs
|
||||
if (k > 0) then
|
||||
k2 = k**2
|
||||
kv2 = 1.0d+00 * k2
|
||||
if (kl2 <= kv2 .and. kv2 <= ku2) then
|
||||
@ -470,10 +512,12 @@ module forcing
|
||||
kmodes(l,2) = 0
|
||||
kmodes(l,3) = k
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
#else /* NDIMS == 3 */
|
||||
do i = 1, kmax
|
||||
do j = -kmax, kmax
|
||||
do i = 0, kxm, kxs
|
||||
if (i > 0) then
|
||||
do j = - kym, kym, kys
|
||||
k2 = i**2 + j**2
|
||||
kv2 = 1.0d+00 * k2
|
||||
if (kl2 <= kv2 .and. kv2 <= ku2) then
|
||||
@ -482,8 +526,10 @@ module forcing
|
||||
kmodes(l,2) = j
|
||||
end if
|
||||
end do
|
||||
end if
|
||||
end do
|
||||
do j = 1, kmax
|
||||
do j = 0, kym, kys
|
||||
if (j > 0) then
|
||||
k2 = j**2
|
||||
kv2 = 1.0d+00 * k2
|
||||
if (kl2 <= kv2 .and. kv2 <= ku2) then
|
||||
@ -491,6 +537,7 @@ module forcing
|
||||
kmodes(l,1) = 0
|
||||
kmodes(l,2) = j
|
||||
end if
|
||||
end if
|
||||
end do
|
||||
#endif /* NDIMS == 3 */
|
||||
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user