gkowal/amun-code
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'master' into reconnection

Browse Source
This commit is contained in:
Grzegorz Kowal 2022-11-16 15:52:56 -03:00
commit 909d0c4d5f
1 changed files with 111 additions and 64 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

175
sources/forcing.F90
View File

@ -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 */