From ec7d12260ef572bfb0d14ee7fe48bb551997da81 Mon Sep 17 00:00:00 2001 From: Grzegorz Kowal Date: Wed, 16 Nov 2022 13:41:43 -0300 Subject: [PATCH] FORCING: Adjust forcing modes for domain smaller than unity. Signed-off-by: Grzegorz Kowal --- sources/forcing.F90 | 175 ++++++++++++++++++++++++++++---------------- 1 file changed, 111 insertions(+), 64 deletions(-) diff --git a/sources/forcing.F90 b/sources/forcing.F90 index c0129ff..7404c23 100644 --- a/sources/forcing.F90 +++ b/sources/forcing.F90 @@ -169,10 +169,11 @@ module forcing ! subroutine initialize_forcing(verbose, status) - use constants , only : pi2 - use helpers , only : print_message - use parameters, only : get_parameter - use random , only : randuni, randnorz + use constants , only : pi2 + use coordinates, only : xlen, ylen, zlen + use helpers , only : print_message + 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,53 +399,63 @@ module forcing ! driving modes ! #if NDIMS == 3 - do i = 1, kmax - do j = -kmax, kmax - do k = -kmax, kmax - k2 = i**2 + j**2 + k**2 + 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 + nmodes = nmodes + 1 + kcount(k2) = kcount(k2) + 2 + end if + end do + end do + end if + end do + 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 nmodes = nmodes + 1 kcount(k2) = kcount(k2) + 2 end if end do - end do + end if end do - do j = 1, kmax - do k = -kmax, kmax - k2 = j**2 + k**2 + 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 do - end do - do k = 1, kmax - 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 do #else /* NDIMS == 3 */ - do i = 1, kmax - do j = -kmax, kmax - k2 = i**2 + j**2 + 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 + nmodes = nmodes + 1 + kcount(k2) = kcount(k2) + 2 + end if + end do + end if + end do + 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 do - end do - do j = 1, kmax - 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 do #endif /* NDIMS == 3 */ @@ -435,61 +472,71 @@ module forcing ! l = 0 #if NDIMS == 3 - do i = 1, kmax - do j = -kmax, kmax - do k = -kmax, kmax - k2 = i**2 + j**2 + k**2 + 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 + l = l + 1 + kmodes(l,1) = i + kmodes(l,2) = j + kmodes(l,3) = k + end if + end do + end do + end if + end do + 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 l = l + 1 - kmodes(l,1) = i + kmodes(l,1) = 0 kmodes(l,2) = j kmodes(l,3) = k end if end do - end do + end if end do - do j = 1, kmax - do k = -kmax, kmax - k2 = j**2 + k**2 + do k = 0, kzm, kzs + if (k > 0) then + k2 = k**2 kv2 = 1.0d+00 * k2 if (kl2 <= kv2 .and. kv2 <= ku2) then l = l + 1 kmodes(l,1) = 0 - kmodes(l,2) = j + kmodes(l,2) = 0 kmodes(l,3) = k end if - end do - end do - do k = 1, kmax - k2 = k**2 - kv2 = 1.0d+00 * k2 - if (kl2 <= kv2 .and. kv2 <= ku2) then - l = l + 1 - kmodes(l,1) = 0 - kmodes(l,2) = 0 - kmodes(l,3) = k end if end do #else /* NDIMS == 3 */ - do i = 1, kmax - do j = -kmax, kmax - k2 = i**2 + j**2 + 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 + l = l + 1 + kmodes(l,1) = i + kmodes(l,2) = j + end if + end do + end if + end do + do j = 0, kym, kys + if (j > 0) then + k2 = j**2 kv2 = 1.0d+00 * k2 if (kl2 <= kv2 .and. kv2 <= ku2) then l = l + 1 - kmodes(l,1) = i + kmodes(l,1) = 0 kmodes(l,2) = j end if - end do - end do - do j = 1, kmax - k2 = j**2 - kv2 = 1.0d+00 * k2 - if (kl2 <= kv2 .and. kv2 <= ku2) then - l = l + 1 - kmodes(l,1) = 0 - kmodes(l,2) = j end if end do #endif /* NDIMS == 3 */