553 lines
17 KiB
Fortran
553 lines
17 KiB
Fortran
subroutine rhs_velocity
|
|
|
|
use m_openmpi
|
|
use m_io
|
|
use m_parameters
|
|
use m_fields
|
|
use m_work
|
|
use x_fftw
|
|
use m_les
|
|
|
|
implicit none
|
|
|
|
integer :: i, j, k, n, nx3
|
|
real*8 :: t1(0:6), rtmp, wnum2, rnx3
|
|
|
|
|
|
|
|
! The IFFT of velocities has been done earlier in rhs_scalars
|
|
! the velocities were kept in wrk1...wrk3, intact.
|
|
|
|
!!$ ! putting the velocity field in the wrk array
|
|
!!$ wrk(:,:,:,1:3) = fields(:,:,:,1:3)
|
|
!!$ ! performing IFFT to convert them to the X-space
|
|
!!$ call xFFT3d(-1,1)
|
|
!!$ call xFFT3d(-1,2)
|
|
!!$ call xFFT3d(-1,3)
|
|
|
|
|
|
!-------------------------------------------------------------------------
|
|
! getting the Courant number (on the master process only)
|
|
wrk(:,:,:,4) = abs(wrk(:,:,:,1)) + abs(wrk(:,:,:,2)) + abs(wrk(:,:,:,3))
|
|
rtmp = maxval(wrk(1:nx,:,:,4))
|
|
call MPI_REDUCE(rtmp,courant,1,MPI_REAL8,MPI_MAX,0,MPI_COMM_TASK,mpi_err)
|
|
if (variable_dt) then
|
|
count = 1
|
|
call MPI_BCAST(courant,count,MPI_REAL8,0,MPI_COMM_TASK,mpi_err)
|
|
end if
|
|
courant = courant * dt / dx
|
|
|
|
!-------------------------------------------------------------------------
|
|
|
|
!--------------------------------------------------------------------------------
|
|
! Calculating the right-hand side for the velocities
|
|
!
|
|
! There are two options available: the standard 2/3 rule (dealias=0) and
|
|
! combination of phase shift and truncation (dealias=1). The latter retains
|
|
! more modes but requires more calculations thus slowing down the simulation.
|
|
! These are treated separately in two different "if" blocks. This is done in
|
|
! order not to complicate the logic. Also this way both blocks can be
|
|
! optimized separately.
|
|
!--------------------------------------------------------------------------------
|
|
|
|
two_thirds_rule: if (dealias.eq.0) then
|
|
|
|
! getting all 6 products of velocities
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx
|
|
t1(1) = wrk(i,j,k,1) * wrk(i,j,k,1)
|
|
t1(2) = wrk(i,j,k,1) * wrk(i,j,k,2)
|
|
t1(3) = wrk(i,j,k,1) * wrk(i,j,k,3)
|
|
t1(4) = wrk(i,j,k,2) * wrk(i,j,k,2)
|
|
t1(5) = wrk(i,j,k,2) * wrk(i,j,k,3)
|
|
t1(6) = wrk(i,j,k,3) * wrk(i,j,k,3)
|
|
|
|
do n = 1,6
|
|
wrk(i,j,k,n) = t1(n)
|
|
end do
|
|
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
|
|
! converting the products to the Fourier space
|
|
do n = 1,6
|
|
call xFFT3d(1,n)
|
|
end do
|
|
|
|
! Building the RHS.
|
|
! First, put into wrk arrays the convectove terms (that will be multiplyed by "i"
|
|
! later) and the factor that corresponds to the diffusion
|
|
|
|
! Do not forget that in Fourier space the indicies are (ix, iz, iy)
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx+2
|
|
|
|
t1(1) = - ( akx(i) * wrk(i,j,k,1) + aky(k) * wrk(i,j,k,2) + akz(j) * wrk(i,j,k,3) )
|
|
t1(2) = - ( akx(i) * wrk(i,j,k,2) + aky(k) * wrk(i,j,k,4) + akz(j) * wrk(i,j,k,5) )
|
|
t1(3) = - ( akx(i) * wrk(i,j,k,3) + aky(k) * wrk(i,j,k,5) + akz(j) * wrk(i,j,k,6) )
|
|
|
|
|
|
t1(4) = - nu * ( akx(i)**2 + aky(k)**2 + akz(j)**2 )
|
|
|
|
do n = 1,4
|
|
wrk(i,j,k,n) = t1(n)
|
|
end do
|
|
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
! now take the actual fields from fields(:,:,:,:) and calculate the RHSs
|
|
|
|
! at this moment the contains of wrk(:,:,:,1:3) are the convective terms in the RHS
|
|
! which are not yet multiplied by "i"
|
|
! wrk(:,:,:,4) contains the Laplace operator in Fourier space. To get the diffusion term
|
|
! we need to take wrk(:,:,:,4) and multiply it by the velocity
|
|
|
|
t1(6) = real(kmax,8)
|
|
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx+1,2
|
|
|
|
! If the dealiasing option is 2/3-rule (dealias=0) then we retain the modes
|
|
! inside the cube described by $| k_i | \leq k_{max}$, $i=1,2,3$.
|
|
! The rest of the modes is purged
|
|
|
|
if (ialias(i,j,k) .gt. 0) then
|
|
! setting the Fourier components to zero
|
|
wrk(i ,j,k,1:3) = zip
|
|
wrk(i+1,j,k,1:3) = zip
|
|
|
|
else
|
|
|
|
! RHS for u, v and w
|
|
do n = 1,3
|
|
! taking the convective term, multiply it by "i"
|
|
! (see how it's done in x_fftw.f90)
|
|
! and adding the diffusion term
|
|
rtmp = - wrk(i+1,j,k,n) + wrk(i ,j,k,4) * fields(i ,j,k,n)
|
|
wrk(i+1,j,k,n) = wrk(i ,j,k,n) + wrk(i+1,j,k,4) * fields(i+1,j,k,n)
|
|
wrk(i ,j,k,n) = rtmp
|
|
end do
|
|
|
|
end if
|
|
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
end if two_thirds_rule
|
|
|
|
!--------------------------------------------------------------------------------
|
|
! The second option (dealias=1). All pairwise products of velocities are
|
|
! dealiased using one phase shift of (dx/2,dy/2,dz/2).
|
|
!--------------------------------------------------------------------------------
|
|
phase_shifting: if (dealias.eq.1) then
|
|
|
|
! work parameters
|
|
wrk(:,:,:,0) = zip
|
|
|
|
! getting all 6 products of velocities
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx
|
|
|
|
t1(1) = wrk(i,j,k,1) * wrk(i,j,k,1)
|
|
t1(2) = wrk(i,j,k,1) * wrk(i,j,k,2)
|
|
t1(3) = wrk(i,j,k,1) * wrk(i,j,k,3)
|
|
t1(4) = wrk(i,j,k,2) * wrk(i,j,k,2)
|
|
t1(5) = wrk(i,j,k,2) * wrk(i,j,k,3)
|
|
t1(6) = wrk(i,j,k,3) * wrk(i,j,k,3)
|
|
do n = 1,6
|
|
wrk(i,j,k,n) = t1(n)
|
|
end do
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
! converting the products to the Fourier space
|
|
do n = 1,6
|
|
call xFFT3d(1,n)
|
|
end do
|
|
|
|
! Building the RHS.
|
|
! First, put into wrk arrays the convectove terms (that will be multiplyed by "i"
|
|
! later) and the factor that corresponds to the diffusion
|
|
|
|
! Do not forget that in Fourier space the indicies are (ix, iz, iy)
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx+2
|
|
|
|
t1(1) = - ( akx(i) * wrk(i,j,k,1) + aky(k) * wrk(i,j,k,2) + akz(j) * wrk(i,j,k,3) )
|
|
t1(2) = - ( akx(i) * wrk(i,j,k,2) + aky(k) * wrk(i,j,k,4) + akz(j) * wrk(i,j,k,5) )
|
|
t1(3) = - ( akx(i) * wrk(i,j,k,3) + aky(k) * wrk(i,j,k,5) + akz(j) * wrk(i,j,k,6) )
|
|
! putting a factor from the diffusion term into t1(4) (and later in wrk4)
|
|
t1(4) = - nu * ( akx(i)**2 + aky(k)**2 + akz(j)**2 )
|
|
do n = 1,4
|
|
wrk(i,j,k,n) = t1(n)
|
|
end do
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
! now use the actual fields from fields(:,:,:,:) to calculate the RHSs
|
|
|
|
! at this moment the contains of wrk(:,:,:,1:3) are the convective terms in the RHS
|
|
! which are not yet multiplied by "i"
|
|
|
|
! wrk(:,:,:,4) contains the Laplace operator in Fourier space. To get the diffusion term
|
|
! we need to take wrk(:,:,:,4) and multiply it by the velocity
|
|
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx+1,2
|
|
|
|
! If the dealiasing option is (dealias=1) then we retain the modes
|
|
! for which no more than one component of the k-vector is larger than nx/3.
|
|
! The rest of the modes is purged.
|
|
|
|
if (ialias(i,j,k) .gt. 1) then
|
|
! setting the Fourier components to zero
|
|
wrk(i ,j,k,1:3) = zip
|
|
wrk(i+1,j,k,1:3) = zip
|
|
else
|
|
! RHS for u, v and w
|
|
do n = 1,3
|
|
! taking the HALF of the convective term, multiply it by "i"
|
|
! and adding the diffusion term
|
|
rtmp = - 0.5d0 * wrk(i+1,j,k,n) + wrk(i ,j,k,4) * fields(i ,j,k,n)
|
|
wrk(i+1,j,k,n) = 0.5d0 * wrk(i ,j,k,n) + wrk(i+1,j,k,4) * fields(i+1,j,k,n)
|
|
wrk(i ,j,k,n) = rtmp
|
|
end do
|
|
end if
|
|
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
|
|
!--------------------------------------------------------------------------------
|
|
! Second part of the phase shifting technique
|
|
!--------------------------------------------------------------------------------
|
|
|
|
! since wrk1...3 are taken by parts of RHS constructed earlier, we can use
|
|
! only wrk0 and wrk4...6.
|
|
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx+1,2
|
|
|
|
! computing sines and cosines for the phase shift of dx/2,dy/2,dz/2
|
|
! and putting them into wrk0
|
|
wrk(i ,j,k,0) = cos(half*(akx(i )+aky(k)+akz(j))*dx)
|
|
wrk(i+1,j,k,0) = sin(half*(akx(i+1)+aky(k)+akz(j))*dx)
|
|
|
|
! wrk4 will have phase-shifted u
|
|
wrk(i ,j,k,4) = fields(i ,j,k,1) * wrk(i,j,k,0) - fields(i+1,j,k,1) * wrk(i+1,j,k,0)
|
|
wrk(i+1,j,k,4) = fields(i+1,j,k,1) * wrk(i,j,k,0) + fields(i ,j,k,1) * wrk(i+1,j,k,0)
|
|
|
|
! wrk5 will have phase-shifted v
|
|
wrk(i ,j,k,5) = fields(i ,j,k,2) * wrk(i,j,k,0) - fields(i+1,j,k,2) * wrk(i+1,j,k,0)
|
|
wrk(i+1,j,k,5) = fields(i+1,j,k,2) * wrk(i,j,k,0) + fields(i ,j,k,2) * wrk(i+1,j,k,0)
|
|
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
! transforming u+ and v+ into X-space
|
|
call xFFT3d(-1,4)
|
|
call xFFT3d(-1,5)
|
|
|
|
! now wrk4 and wrk5 contain u+ and v+
|
|
|
|
! getting (u+)*(u+) in real space, converting it to Fourier space,
|
|
! phase shifting back and adding -0.5*(the results) to the RHS for u
|
|
wrk(:,:,:,6) = wrk(:,:,:,4)**2
|
|
call xFFT3d(1,6)
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx+1,2
|
|
rtmp = wrk(i ,j,k,6) * wrk(i,j,k,0) + wrk(i+1,j,k,6) * wrk(i+1,j,k,0)
|
|
wrk(i+1,j,k,6) = wrk(i+1,j,k,6) * wrk(i,j,k,0) - wrk(i ,j,k,6) * wrk(i+1,j,k,0)
|
|
wrk(i ,j,k,6) = rtmp
|
|
end do
|
|
end do
|
|
end do
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx+1,2
|
|
wrk(i ,j,k,1) = wrk(i ,j,k,1) + 0.5d0 * akx(i+1) * wrk(i+1,j,k,6)
|
|
wrk(i+1,j,k,1) = wrk(i+1,j,k,1) - 0.5d0 * akx(i ) * wrk(i ,j,k,6)
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
! getting (u+)*(v+) in real space, converting it to Fourier space,
|
|
! phase shifting back and adding -0.5*(the results) to the RHSs for u and v
|
|
wrk(:,:,:,6) = wrk(:,:,:,4)*wrk(:,:,:,5)
|
|
call xFFT3d(1,6)
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx+1,2
|
|
rtmp = wrk(i ,j,k,6) * wrk(i,j,k,0) + wrk(i+1,j,k,6) * wrk(i+1,j,k,0)
|
|
wrk(i+1,j,k,6) = wrk(i+1,j,k,6) * wrk(i,j,k,0) - wrk(i ,j,k,6) * wrk(i+1,j,k,0)
|
|
wrk(i ,j,k,6) = rtmp
|
|
end do
|
|
end do
|
|
end do
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx+1,2
|
|
wrk(i ,j,k,1) = wrk(i ,j,k,1) + 0.5d0 * aky(k) * wrk(i+1,j,k,6)
|
|
wrk(i+1,j,k,1) = wrk(i+1,j,k,1) - 0.5d0 * aky(k) * wrk(i ,j,k,6)
|
|
|
|
wrk(i ,j,k,2) = wrk(i ,j,k,2) + 0.5d0 * akx(i+1) * wrk(i+1,j,k,6)
|
|
wrk(i+1,j,k,2) = wrk(i+1,j,k,2) - 0.5d0 * akx(i ) * wrk(i ,j,k,6)
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
! getting (v+)*(v+) in real space, converting it to Fourier space,
|
|
! phase shifting back and adding -0.5*(the results) to the RHS for v
|
|
wrk(:,:,:,6) = wrk(:,:,:,5)**2
|
|
call xFFT3d(1,6)
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx+1,2
|
|
rtmp = wrk(i ,j,k,6) * wrk(i,j,k,0) + wrk(i+1,j,k,6) * wrk(i+1,j,k,0)
|
|
wrk(i+1,j,k,6) = wrk(i+1,j,k,6) * wrk(i,j,k,0) - wrk(i ,j,k,6) * wrk(i+1,j,k,0)
|
|
wrk(i ,j,k,6) = rtmp
|
|
end do
|
|
end do
|
|
end do
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx+1,2
|
|
wrk(i ,j,k,2) = wrk(i ,j,k,2) + 0.5d0 * aky(k) * wrk(i+1,j,k,6)
|
|
wrk(i+1,j,k,2) = wrk(i+1,j,k,2) - 0.5d0 * aky(k) * wrk(i ,j,k,6)
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
! now get the (w+) in wrk6
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx+1,2
|
|
! wrk6 will have phase-shifted w
|
|
wrk(i ,j,k,6) = fields(i ,j,k,3) * wrk(i,j,k,0) - fields(i+1,j,k,3) * wrk(i+1,j,k,0)
|
|
wrk(i+1,j,k,6) = fields(i+1,j,k,3) * wrk(i,j,k,0) + fields(i ,j,k,3) * wrk(i+1,j,k,0)
|
|
end do
|
|
end do
|
|
end do
|
|
! transforming w+ into X-space
|
|
call xFFT3d(-1,6)
|
|
|
|
! at this point wrk4..6 contain (u+), (v+) and (w+) in real space.
|
|
! the combinations that we have not dealt with are: uw, vw and ww.
|
|
! we'll deal with all three of them at once.
|
|
|
|
! first get all three of these in wrk4...6 and
|
|
wrk(:,:,:,4) = wrk(:,:,:,4) * wrk(:,:,:,6)
|
|
wrk(:,:,:,5) = wrk(:,:,:,5) * wrk(:,:,:,6)
|
|
wrk(:,:,:,6) = wrk(:,:,:,6)**2
|
|
|
|
! transform them into Fourier space
|
|
call xFFT3d(1,4)
|
|
call xFFT3d(1,5)
|
|
call xFFT3d(1,6)
|
|
|
|
! phase shift back to origianl grid and add to corresponding RHSs
|
|
do n = 4,6
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx+1,2
|
|
rtmp = wrk(i ,j,k,n) * wrk(i,j,k,0) + wrk(i+1,j,k,n) * wrk(i+1,j,k,0)
|
|
wrk(i+1,j,k,n) = wrk(i+1,j,k,n) * wrk(i,j,k,0) - wrk(i ,j,k,n) * wrk(i+1,j,k,0)
|
|
wrk(i ,j,k,n) = rtmp
|
|
end do
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
! adding to corresponding RHSs
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx+1,2
|
|
|
|
! If the dealiasing option is (dealias=1) then we retain the modes
|
|
! for which no more than one component of the k-vector is larger than nx/3.
|
|
! The rest of the modes is purged.
|
|
|
|
if (ialias(i,j,k) .lt. 2) then
|
|
|
|
wrk(i ,j,k,1) = wrk(i ,j,k,1) + 0.5d0 * akz(j) * wrk(i+1,j,k,4)
|
|
wrk(i+1,j,k,1) = wrk(i+1,j,k,1) - 0.5d0 * akz(j) * wrk(i ,j,k,4)
|
|
|
|
wrk(i ,j,k,2) = wrk(i ,j,k,2) + 0.5d0 * akz(j) * wrk(i+1,j,k,5)
|
|
wrk(i+1,j,k,2) = wrk(i+1,j,k,2) - 0.5d0 * akz(j) * wrk(i ,j,k,5)
|
|
|
|
wrk(i ,j,k,3) = wrk(i ,j,k,3) + 0.5d0 * &
|
|
(akx(i+1)*wrk(i+1,j,k,4) + aky(k)*wrk(i+1,j,k,5) + akz(j)*wrk(i+1,j,k,6))
|
|
wrk(i+1,j,k,3) = wrk(i+1,j,k,3) - 0.5d0 * &
|
|
(akx(i )*wrk(i ,j,k,4) + aky(k)*wrk(i ,j,k,5) + akz(j)*wrk(i ,j,k,6))
|
|
|
|
else
|
|
wrk(i:i+1,j,k,1) = zip
|
|
wrk(i:i+1,j,k,2) = zip
|
|
wrk(i:i+1,j,k,3) = zip
|
|
end if
|
|
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
end if phase_shifting
|
|
|
|
! if performing large eddy simulations, call LES subroutine to augment
|
|
! the right hand side for velocioties
|
|
les_active: if (les) then
|
|
call les_rhs_velocity
|
|
end if les_active
|
|
|
|
return
|
|
end subroutine rhs_velocity
|
|
|
|
|
|
!================================================================================
|
|
!================================================================================
|
|
!================================================================================
|
|
!================================================================================
|
|
!================================================================================
|
|
!================================================================================
|
|
!================================================================================
|
|
|
|
subroutine test_rhs_velocity
|
|
|
|
use m_openmpi
|
|
use m_io
|
|
use m_parameters
|
|
use m_fields
|
|
use m_work
|
|
use x_fftw
|
|
|
|
implicit none
|
|
|
|
integer :: i,j,k, n
|
|
real*8 :: a,b,c, x,y,z
|
|
|
|
|
|
! defining very particular velocities so the RHS can be computed analytically
|
|
|
|
if (task.eq.'hydro') then
|
|
|
|
|
|
write(out,*) 'inside.'
|
|
call flush(out)
|
|
|
|
|
|
a = 1.d0
|
|
b = 1.d0
|
|
c = 1.d0
|
|
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx
|
|
|
|
x = dx*real(i-1)
|
|
y = dx*real(j-1)
|
|
z = dx*real(myid*nz + k-1)
|
|
|
|
wrk(i,j,k,1) = sin(a * x)
|
|
wrk(i,j,k,2) = sin(b * y)
|
|
wrk(i,j,k,3) = sin(c * z)
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
write(out,*) 'did work'
|
|
call flush(out)
|
|
|
|
|
|
|
|
do n = 1,3
|
|
call xFFT3d(1,n)
|
|
fields(:,:,:,n) = wrk(:,:,:,n)
|
|
end do
|
|
|
|
|
|
write(out,*) 'did FFTs'
|
|
call flush(out)
|
|
|
|
nu = 0.d0
|
|
|
|
call rhs_velocity
|
|
|
|
write(out,*) 'got rhs'
|
|
call flush(out)
|
|
|
|
do n = 1,3
|
|
call xFFT3d(-1,n)
|
|
end do
|
|
|
|
write(out,*) 'did FFTs'
|
|
call flush(out)
|
|
|
|
|
|
|
|
do k = 1,nz
|
|
do j = 1,ny
|
|
do i = 1,nx
|
|
|
|
x = dx*real(i-1)
|
|
y = dx*real(j-1)
|
|
z = dx*real(myid*nz + k-1)
|
|
|
|
! checking u
|
|
wrk(i,j,k,4) = -sin(a*x) * ( two*a*cos(a*x) + b*cos(b*y) + c*cos(c*z) + nu*a**2)
|
|
|
|
! checking v
|
|
wrk(i,j,k,5) = -sin(b*y) * ( two*b*cos(b*y) + a*cos(a*x) + c*cos(c*z) + nu*b**2)
|
|
|
|
! checking w
|
|
wrk(i,j,k,6) = -sin(c*z) * ( two*c*cos(c*z) + b*cos(b*y) + a*cos(a*x) + nu*c**2)
|
|
|
|
|
|
end do
|
|
end do
|
|
end do
|
|
|
|
!!$ do k = 1,nz
|
|
!!$ write(out,"(3e15.6)") wrk(1,1,k,3),wrk(1,1,k,5),wrk(1,1,k,4)
|
|
!!$ end do
|
|
|
|
wrk(:,:,:,0) = &
|
|
abs(wrk(:,:,:,1) - wrk(:,:,:,4)) + &
|
|
abs(wrk(:,:,:,2) - wrk(:,:,:,5)) + &
|
|
abs(wrk(:,:,:,3) - wrk(:,:,:,6))
|
|
|
|
print *,'Maximum error is ',maxval(wrk(1:nx,:,:,0))
|
|
|
|
|
|
tmp4(:,:,:) = wrk(1:nx,:,:,1) - wrk(1:nx,:,:,4)
|
|
fname = 'e1.arr'
|
|
call write_tmp4
|
|
|
|
tmp4(:,:,:) = wrk(1:nx,:,:,2) - wrk(1:nx,:,:,5)
|
|
fname = 'e2.arr'
|
|
call write_tmp4
|
|
|
|
tmp4(:,:,:) = wrk(1:nx,:,:,3) - wrk(1:nx,:,:,6)
|
|
fname = 'e3.arr'
|
|
call write_tmp4
|
|
|
|
|
|
|
|
end if
|
|
return
|
|
end subroutine test_rhs_velocity
|