38 real(mytype),
dimension(:),
intent(out) :: aa,cc
39 complex(mytype),
dimension(:,:,:),
intent(out) :: bb
42 complex(mytype) :: kx,ky
45 bb(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
79 elseif (
cds .eq. 2)
then 81 kx = (
wx1**(3.d0*l) - &
82 54.d0*
wx1**(2.d0*l) + &
86 54.d0*
wx1**(-2.d0*l) + &
87 wx1**(-3.d0*l) ) /576.d0/
dx2 89 ky = (
wy1**(3.d0*m) - &
90 54.d0*
wy1**(2.d0*m) + &
94 54.d0*
wy1**(-2.d0*m) + &
95 wy1**(-3.d0*m) ) /576.d0/
dy2 97 elseif (
cds .eq. 3)
then 99 kx = -(2.d0*
pi/
lx*l)**2.d0
100 ky = -(2.d0*
pi/
ly*m)**2.d0
104 bb(i,j,k)= kx + ky - aa(k) - cc(k)
107 if (k .eq.
kstr3)
then 108 bb(i,j,k)=bb(i,j,k)+aa(k)
111 if (k .eq.
kend3)
then 112 bb(i,j,k)=bb(i,j,k)+cc(k)
134 real(mytype),
dimension(:,:,:),
intent(in) :: u,v,w
135 real(mytype),
dimension(:,:,:),
intent(out) :: qp
138 real(mytype) :: dudx,dvdy,dwdz
146 if (
cds .eq. 1)
then 148 dudx = ( u(i,j,k)-u(i-1,j,k) )/
dx 149 dvdy = ( v(i,j,k)-v(i,j-1,k) )/
dy 151 elseif (
cds .eq. 2)
then 153 dudx = ( -u(i+1,j,k) + &
156 u(i-2,j,k) ) / 24.d0 /
dx 158 dvdy = ( -v(i,j+1,k) + &
161 v(i,j-2,k) ) / 24.d0 /
dy 166 dwdz = ( w(i,j,k)-w(i,j,k-1) )/
dz(k)
168 qp(i,j,k)=(dudx+dvdy+dwdz)/
dt 187 real(mytype),
dimension(:),
intent(out) :: aa,cc
188 complex(mytype),
dimension(:,:,:),
intent(out) :: bb
191 complex(mytype) :: kx,ky
194 bb(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
222 if (
cds .eq. 1)
then 228 elseif (
cds .eq. 2)
then 230 kx =
dt * ( -
wx1**(2.d0*l) + &
236 ky =
dt * ( -
wy1**(2.d0*m) + &
242 elseif (
cds .eq. 3)
then 244 kx = 0.5d0*
dt*
nu*( -(2.d0*
pi/
lx*l)**2.d0 )
245 ky = 0.5d0*
dt*
nu*( -(2.d0*
pi/
ly*m)**2.d0 )
252 bb(i,j,k)= - kx - ky - aa(k) - cc(k)
257 if (k .eq.
kstr3)
then 258 bb(i,j,k)=bb(i,j,k)-aa(k)
264 if (k .eq.
kend3)
then 267 bb(i,j,k)=bb(i,j,k)-cc(k)
270 bb(i,j,k)=bb(i,j,k)+cc(k)
293 real(mytype),
dimension(:),
intent(out) :: aa,cc
294 complex(mytype),
dimension(:,:,:),
intent(out) :: bb
297 complex(mytype) :: kx,ky
300 bb(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
328 if (
cds .eq. 1)
then 334 elseif (
cds .eq. 2)
then 336 kx =
dt * ( -
wx1**(2.d0*l) + &
342 ky =
dt * ( -
wy1**(2.d0*m) + &
348 elseif (
cds .eq. 3)
then 350 kx = 0.5d0*
dt*
nu*( -(2.d0*
pi/
lx*l)**2.d0 )
351 ky = 0.5d0*
dt*
nu*( -(2.d0*
pi/
ly*m)**2.d0 )
358 bb(i,j,k)= - kx - ky - aa(k) - cc(k)
361 if (k .eq.
kstr3)
then 364 if (k .eq.
kend3-1)
then 385 real(mytype),
dimension(:),
intent(out) :: aa,cc
386 complex(mytype),
dimension(:,:,:) ,
intent(out) :: bb
389 complex(mytype) :: kx,ky
392 bb(:,:,:)=cmplx(0.d0,0.d0,kind=mytype)
420 if (
cds .eq. 1)
then 426 elseif (
cds .eq. 2)
then 428 kx =
dt * ( -
wx1**(2.d0*l) + &
434 ky =
dt * ( -
wy1**(2.d0*m) + &
440 elseif (
cds .eq. 3)
then 450 bb(i,j,k)= - kx - ky - aa(k) - cc(k)
455 if (k .eq.
kstr3)
then 457 bb(i,j,k)=bb(i,j,k)-aa(k)
464 if (k .eq.
kend3)
then 466 bb(i,j,k)=bb(i,j,k)-cc(k)
486 use,
intrinsic :: iso_c_binding
495 include
'include/fftw3.f03' 498 real(mytype),
dimension(:,:,:),
intent(inout) :: aa,cc
499 complex(mytype),
dimension(:,:,:),
intent(inout) :: bb
502 real(mytype),
dimension(:,:,:),
intent(in) :: r_src
513 integer,
intent(in) :: zstagger,delnull,tbnd,ubnd
516 real(mytype),
dimension(:,:,:),
intent(out) :: r_var
519 complex(mytype) :: tmp
521 real(mytype) :: aa_t_zstr,cc_t_zend
522 real(mytype) :: aa_u_zstr,cc_u_zend
523 integer :: i,j,k,kmax3
528 c_var(:,:,:)=cmplx(0.d0,0.d0 , kind=mytype )
529 c_src(:,:,:)=cmplx(0.d0,0.d0 , kind=mytype )
547 if (zstagger .eq. 1)
then 549 elseif (zstagger .eq. 0)
then 552 print *,
'Invalid zstagger flag!' 558 if (delnull .eq. 1)
then 560 if (
myid .eq. 0)
then 563 c_src(1,1,1)=cmplx(0.d0,0.d0 , kind=mytype )
565 elseif (delnull .eq. 0)
then 568 print *,
'delnull flag error!' 572 if (tbnd .eq. 1)
then 574 if (
myid .eq. 0)
then 582 elseif (tbnd .eq. 0)
then 585 print *,
'tbnd flag error!' 590 if (ubnd .eq. 1)
then 592 if (
myid .eq. 0)
then 602 elseif (ubnd .eq. 0)
then 605 print *,
'ubnd flag error!' 615 tmp=aa(i,j,k)/bb(i,j,k-1)
616 bb(i,j,k)=bb(i,j,k)-cc(i,j,k-1)*tmp
631 do k=kmax3-1,
kstr3,-1
636 cc(i,j,k)*
c_var(i,j,k-
ghst+1) ) / bb(i,j,k)
650 r_var(i+
ghst,j+
ghst,k)=
real(c_var(i,j,k-ghst), kind=mytype )
integer, save, protected myid
integer, save, protected ochannel
integer, save, protected ny
subroutine get_u_eqn_coeff(aa, bb, cc)
subroutine get_p_eqn_src(u, v, w, qp)
complex(mytype), dimension(:,:,:), allocatable, save c_src
real(mytype), parameter pi
real(mytype), save, protected dt
complex(mytype), save, protected wy1
integer, save, protected jend3
real(mytype), save, protected nu
real(mytype), save, protected ttop
integer, save, protected nz
subroutine get_t_eqn_coeff(aa, bb, cc)
real(mytype), save, protected ly
real(mytype), dimension(:), allocatable, save, protected dz
subroutine poisson_solver_fft(aa, bb, cc, r_src, r_var, zstagger, delnull, tbnd, ubnd)
real(mytype), save, protected dy2
real(mytype), save, protected alpha
subroutine get_p_eqn_coeff(aa, bb, cc)
integer, save, protected istr3
real(mytype), save, protected dx
integer, save, protected jstr3
real(mytype), dimension(:), allocatable, save, protected dz_b
integer, save, protected i_offset
real(mytype), save, protected dx2
integer, save, protected j_offset
real(mytype), save, protected u_mrf
integer, save, protected cds
real(mytype), save, protected lx
real(mytype), dimension(:), allocatable, save, protected dz_t
real(mytype), save, protected dy
integer, save, protected nx
subroutine get_w_eqn_coeff(aa, bb, cc)
complex(mytype), save, protected wx1
real(mytype), save, protected rk_c
integer, save, protected iend3
real(mytype), save, protected tbot
complex(mytype), dimension(:,:,:), allocatable, save c_var
integer, save, protected kstr3
integer, save, protected kend3