incomp-flame-post/code/m_terms.f90
2019-05-03 02:03:08 +09:00

959 lines
20 KiB
Fortran

module m_terms
use m_openmpi
use m_parameters
use m_arrays
use m_calculate
implicit none
character (len = *), parameter :: output_header="x avg_c avg_u avg_v avg_w avg_fsd avg_lapc avg_ddx_c avg_ddy_c avg_ddz_c avg_d2dx_c avg_d2dy_c avg_d2dz_c fsd_avg_u fsd_avg_v fsd_avg_w fsd_avg_nx fsd_avg_ny fsd_avg_nz fsd_avg_t1 fsd_avg_t2 fsd_avg_t3 fsd_avg_t4 fsd_avg_t5 fsd_avg_t6 fsd_avg_t7 fsd_avg_t8 fsd_avg_absk fsd_avg_divn fsd_avg_sd"
real*8, allocatable, dimension(:) :: fsd_avg_divn
real*8, allocatable, dimension(:,:,:) :: t3_fsd
real*8, allocatable, dimension(:,:,:) :: ddx_nx
real*8, allocatable, dimension(:) :: fsd_avg_ny
real*8, allocatable, dimension(:) :: fsd_avg_nx
real*8, allocatable, dimension(:,:,:) :: t2_fsd
real*8, allocatable, dimension(:) :: fsd_avg_nz
real*8, allocatable, dimension(:,:,:) :: ddz_c
real*8, allocatable, dimension(:,:,:) :: vn
real*8, allocatable, dimension(:) :: avg_d2dz_c
real*8, allocatable, dimension(:) :: avg_d2dx_c
real*8, allocatable, dimension(:,:,:) :: ddy_c
real*8, allocatable, dimension(:,:,:) :: ddx_c
real*8, allocatable, dimension(:,:,:) :: t8_fsd
real*8, allocatable, dimension(:) :: fsd_avg_t3
real*8, allocatable, dimension(:) :: fsd_avg_t2
real*8, allocatable, dimension(:) :: fsd_avg_t1
real*8, allocatable, dimension(:) :: fsd_avg_t7
real*8, allocatable, dimension(:) :: fsd_avg_v
real*8, allocatable, dimension(:) :: fsd_avg_t5
real*8, allocatable, dimension(:,:,:) :: ddy_ny
real*8, allocatable, dimension(:,:,:) :: wrate
real*8, allocatable, dimension(:) :: fsd_avg_t8
real*8, allocatable, dimension(:) :: avg_fsd
real*8, allocatable, dimension(:) :: fsd_avg_absk
real*8, allocatable, dimension(:,:,:) :: nx
real*8, allocatable, dimension(:,:,:) :: ny
real*8, allocatable, dimension(:,:,:) :: nz
real*8, allocatable, dimension(:,:,:) :: ddz_nz
real*8, allocatable, dimension(:) :: fsd_avg_sd
real*8, allocatable, dimension(:,:,:) :: lapc
real*8, allocatable, dimension(:) :: avg_ddy_c
real*8, allocatable, dimension(:,:,:) :: divn
real*8, allocatable, dimension(:,:,:) :: t7_fsd
real*8, allocatable, dimension(:) :: avg_ddx_c
real*8, allocatable, dimension(:) :: avg_lapc
real*8, allocatable, dimension(:,:,:) :: t6_fsd
real*8, allocatable, dimension(:,:,:) :: absk
real*8, allocatable, dimension(:) :: avg_ddz_c
real*8, allocatable, dimension(:) :: fsd_avg_t6
real*8, allocatable, dimension(:,:,:) :: c
real*8, allocatable, dimension(:) :: avg_c
real*8, allocatable, dimension(:,:,:) :: fsd
real*8, allocatable, dimension(:,:,:) :: t5_fsd
real*8, allocatable, dimension(:,:,:) :: d2dy_c
real*8, allocatable, dimension(:) :: avg_d2dy_c
real*8, allocatable, dimension(:) :: fsd_avg_t4
real*8, allocatable, dimension(:,:,:) :: d2dx_c
real*8, allocatable, dimension(:,:,:) :: t1
real*8, allocatable, dimension(:) :: avg_u
real*8, allocatable, dimension(:) :: avg_v
real*8, allocatable, dimension(:) :: avg_w
real*8, allocatable, dimension(:,:,:) :: d2dz_c
real*8, allocatable, dimension(:) :: fsd_avg_w
real*8, allocatable, dimension(:) :: fsd_avg_u
real*8, allocatable, dimension(:,:,:) :: t4_fsd
real*8, allocatable, dimension(:,:,:) :: sd
contains
subroutine m_terms_init
integer :: ierr
allocate(fsd_avg_divn(nxp), stat=ierr) ; fsd_avg_divn = 0.
allocate(t3_fsd(nxp,nyp,nzp), stat=ierr) ; t3_fsd = 0.
allocate(ddx_nx(nxp,nyp,nzp), stat=ierr) ; ddx_nx = 0.
allocate(fsd_avg_ny(nxp), stat=ierr) ; fsd_avg_ny = 0.
allocate(fsd_avg_nx(nxp), stat=ierr) ; fsd_avg_nx = 0.
allocate(t2_fsd(nxp,nyp,nzp), stat=ierr) ; t2_fsd = 0.
allocate(fsd_avg_nz(nxp), stat=ierr) ; fsd_avg_nz = 0.
allocate(ddz_c(nxp,nyp,nzp), stat=ierr) ; ddz_c = 0.
allocate(vn(nxp,nyp,nzp), stat=ierr) ; vn = 0.
allocate(avg_d2dz_c(nxp), stat=ierr) ; avg_d2dz_c = 0.
allocate(avg_d2dx_c(nxp), stat=ierr) ; avg_d2dx_c = 0.
allocate(ddy_c(nxp,nyp,nzp), stat=ierr) ; ddy_c = 0.
allocate(ddx_c(nxp,nyp,nzp), stat=ierr) ; ddx_c = 0.
allocate(t8_fsd(nxp,nyp,nzp), stat=ierr) ; t8_fsd = 0.
allocate(fsd_avg_t3(nxp), stat=ierr) ; fsd_avg_t3 = 0.
allocate(fsd_avg_t2(nxp), stat=ierr) ; fsd_avg_t2 = 0.
allocate(fsd_avg_t1(nxp), stat=ierr) ; fsd_avg_t1 = 0.
allocate(fsd_avg_t7(nxp), stat=ierr) ; fsd_avg_t7 = 0.
allocate(fsd_avg_v(nxp), stat=ierr) ; fsd_avg_v = 0.
allocate(fsd_avg_t5(nxp), stat=ierr) ; fsd_avg_t5 = 0.
allocate(ddy_ny(nxp,nyp,nzp), stat=ierr) ; ddy_ny = 0.
allocate(wrate(nxp,nyp,nzp), stat=ierr) ; wrate = 0.
allocate(fsd_avg_t8(nxp), stat=ierr) ; fsd_avg_t8 = 0.
allocate(avg_fsd(nxp), stat=ierr) ; avg_fsd = 0.
allocate(fsd_avg_absk(nxp), stat=ierr) ; fsd_avg_absk = 0.
allocate(nx(nxp,nyp,nzp), stat=ierr) ; nx = 0.
allocate(ny(nxp,nyp,nzp), stat=ierr) ; ny = 0.
allocate(nz(nxp,nyp,nzp), stat=ierr) ; nz = 0.
allocate(ddz_nz(nxp,nyp,nzp), stat=ierr) ; ddz_nz = 0.
allocate(fsd_avg_sd(nxp), stat=ierr) ; fsd_avg_sd = 0.
allocate(lapc(nxp,nyp,nzp), stat=ierr) ; lapc = 0.
allocate(avg_ddy_c(nxp), stat=ierr) ; avg_ddy_c = 0.
allocate(divn(nxp,nyp,nzp), stat=ierr) ; divn = 0.
allocate(t7_fsd(nxp,nyp,nzp), stat=ierr) ; t7_fsd = 0.
allocate(avg_ddx_c(nxp), stat=ierr) ; avg_ddx_c = 0.
allocate(avg_lapc(nxp), stat=ierr) ; avg_lapc = 0.
allocate(t6_fsd(nxp,nyp,nzp), stat=ierr) ; t6_fsd = 0.
allocate(absk(nxp,nyp,nzp), stat=ierr) ; absk = 0.
allocate(avg_ddz_c(nxp), stat=ierr) ; avg_ddz_c = 0.
allocate(fsd_avg_t6(nxp), stat=ierr) ; fsd_avg_t6 = 0.
allocate(c(nxp,nyp,nzp), stat=ierr) ; c = 0.
allocate(avg_c(nxp), stat=ierr) ; avg_c = 0.
allocate(fsd(nxp,nyp,nzp), stat=ierr) ; fsd = 0.
allocate(t5_fsd(nxp,nyp,nzp), stat=ierr) ; t5_fsd = 0.
allocate(d2dy_c(nxp,nyp,nzp), stat=ierr) ; d2dy_c = 0.
allocate(avg_d2dy_c(nxp), stat=ierr) ; avg_d2dy_c = 0.
allocate(fsd_avg_t4(nxp), stat=ierr) ; fsd_avg_t4 = 0.
allocate(d2dx_c(nxp,nyp,nzp), stat=ierr) ; d2dx_c = 0.
allocate(t1(nxp,nyp,nzp), stat=ierr) ; t1 = 0.
allocate(avg_u(nxp), stat=ierr) ; avg_u = 0.
allocate(avg_v(nxp), stat=ierr) ; avg_v = 0.
allocate(avg_w(nxp), stat=ierr) ; avg_w = 0.
allocate(d2dz_c(nxp,nyp,nzp), stat=ierr) ; d2dz_c = 0.
allocate(fsd_avg_w(nxp), stat=ierr) ; fsd_avg_w = 0.
allocate(fsd_avg_u(nxp), stat=ierr) ; fsd_avg_u = 0.
allocate(t4_fsd(nxp,nyp,nzp), stat=ierr) ; t4_fsd = 0.
allocate(sd(nxp,nyp,nzp), stat=ierr) ; sd = 0.
end subroutine m_terms_init
subroutine m_terms_finalize
deallocate(fsd_avg_divn)
deallocate(t3_fsd)
deallocate(ddx_nx)
deallocate(fsd_avg_ny)
deallocate(fsd_avg_nx)
deallocate(t2_fsd)
deallocate(fsd_avg_nz)
deallocate(ddz_c)
deallocate(vn)
deallocate(avg_d2dz_c)
deallocate(avg_d2dx_c)
deallocate(ddy_c)
deallocate(ddx_c)
deallocate(t8_fsd)
deallocate(fsd_avg_t3)
deallocate(fsd_avg_t2)
deallocate(fsd_avg_t1)
deallocate(fsd_avg_t7)
deallocate(fsd_avg_v)
deallocate(fsd_avg_t5)
deallocate(ddy_ny)
deallocate(wrate)
deallocate(fsd_avg_t8)
deallocate(avg_fsd)
deallocate(fsd_avg_absk)
deallocate(nx)
deallocate(ny)
deallocate(nz)
deallocate(ddz_nz)
deallocate(fsd_avg_sd)
deallocate(lapc)
deallocate(avg_ddy_c)
deallocate(divn)
deallocate(t7_fsd)
deallocate(avg_ddx_c)
deallocate(avg_lapc)
deallocate(t6_fsd)
deallocate(absk)
deallocate(avg_ddz_c)
deallocate(fsd_avg_t6)
deallocate(c)
deallocate(avg_c)
deallocate(fsd)
deallocate(t5_fsd)
deallocate(d2dy_c)
deallocate(avg_d2dy_c)
deallocate(fsd_avg_t4)
deallocate(d2dx_c)
deallocate(t1)
deallocate(avg_u)
deallocate(avg_v)
deallocate(avg_w)
deallocate(d2dz_c)
deallocate(fsd_avg_w)
deallocate(fsd_avg_u)
deallocate(t4_fsd)
deallocate(sd)
end subroutine m_terms_finalize
subroutine m_terms_calculate_pass1
integer :: i, j, k
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
c(i,j,k) = ( 1.0 - y(i,j,k) )
end do
end do
end do
call ddz ( ddz_c, c )
call ddy ( ddy_c, c )
call ddx ( ddx_c, c )
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
wrate(i,j,k) = ( rxn_rate ( c(i,j,k) ) )
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
avg_ddy_c(i) = avg_ddy_c(i) + ddy_c(i,j,k)
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
avg_ddx_c(i) = avg_ddx_c(i) + ddx_c(i,j,k)
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
avg_ddz_c(i) = avg_ddz_c(i) + ddz_c(i,j,k)
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
avg_c(i) = avg_c(i) + c(i,j,k)
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
fsd(i,j,k) = ( sqrt ( ( ( ((ddx_c(i,j,k))*(ddx_c(i,j,k))) + ((ddy_c(i,j,k))*(ddy_c(i,j,k))) ) + ((ddz_c(i,j,k))*(ddz_c(i,j,k))) ) ) )
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
fsd_avg_v(i) = fsd_avg_v(i) + v(i,j,k) * fsd(i,j,k)
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
avg_fsd(i) = avg_fsd(i) + fsd(i,j,k)
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
nx(i,j,k) = ( ( - ddx_c(i,j,k) ) / fsd(i,j,k) )
end do
end do
end do
call ddx ( ddx_nx, nx )
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
fsd_avg_nx(i) = fsd_avg_nx(i) + nx(i,j,k) * fsd(i,j,k)
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
ny(i,j,k) = ( ( - ddy_c(i,j,k) ) / fsd(i,j,k) )
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
fsd_avg_ny(i) = fsd_avg_ny(i) + ny(i,j,k) * fsd(i,j,k)
end do
end do
end do
call ddy ( ddy_ny, ny )
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
nz(i,j,k) = ( ( - ddz_c(i,j,k) ) / fsd(i,j,k) )
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
fsd_avg_nz(i) = fsd_avg_nz(i) + nz(i,j,k) * fsd(i,j,k)
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
vn(i,j,k) = ( ( ( u(i,j,k) * nx(i,j,k) ) + ( v(i,j,k) * ny(i,j,k) ) ) + ( w(i,j,k) * nz(i,j,k) ) )
end do
end do
end do
call ddz ( ddz_nz, nz )
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
divn(i,j,k) = ( ( ddx_nx(i,j,k) + ddy_ny(i,j,k) ) + ddz_nz(i,j,k) )
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
fsd_avg_divn(i) = fsd_avg_divn(i) + divn(i,j,k) * fsd(i,j,k)
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
absk(i,j,k) = ( dabs ( divn(i,j,k) ) )
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
fsd_avg_absk(i) = fsd_avg_absk(i) + absk(i,j,k) * fsd(i,j,k)
end do
end do
end do
call d2dy ( d2dy_c, c )
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
avg_d2dy_c(i) = avg_d2dy_c(i) + d2dy_c(i,j,k)
end do
end do
end do
call d2dx ( d2dx_c, c )
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
avg_d2dx_c(i) = avg_d2dx_c(i) + d2dx_c(i,j,k)
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
avg_u(i) = avg_u(i) + u(i,j,k)
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
avg_v(i) = avg_v(i) + v(i,j,k)
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
avg_w(i) = avg_w(i) + w(i,j,k)
end do
end do
end do
call d2dz ( d2dz_c, c )
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
avg_d2dz_c(i) = avg_d2dz_c(i) + d2dz_c(i,j,k)
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
lapc(i,j,k) = ( ( d2dx_c(i,j,k) + d2dy_c(i,j,k) ) + d2dz_c(i,j,k) )
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
avg_lapc(i) = avg_lapc(i) + lapc(i,j,k)
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
fsd_avg_w(i) = fsd_avg_w(i) + w(i,j,k) * fsd(i,j,k)
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
fsd_avg_u(i) = fsd_avg_u(i) + u(i,j,k) * fsd(i,j,k)
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
sd(i,j,k) = ( ( ( ( ( d2dx_c(i,j,k) + d2dy_c(i,j,k) ) + d2dz_c(i,j,k) ) * rod ) + wrate(i,j,k) ) / fsd(i,j,k) )
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
fsd_avg_sd(i) = fsd_avg_sd(i) + sd(i,j,k) * fsd(i,j,k)
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
t1(i,j,k) = ( ( vn(i,j,k) + sd(i,j,k) ) * nx(i,j,k) )
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
fsd_avg_t1(i) = fsd_avg_t1(i) + t1(i,j,k) * fsd(i,j,k)
end do
end do
end do
end subroutine m_terms_calculate_pass1
subroutine m_terms_average_pass1 (nfiles)
integer :: nfiles
real*8 :: denum
denum=real(nfiles*nyp*nzp)
call MPI_ALLREDUCE(MPI_IN_PLACE, avg_c, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
avg_c = avg_c / denum
call MPI_ALLREDUCE(MPI_IN_PLACE, avg_u, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
avg_u = avg_u / denum
call MPI_ALLREDUCE(MPI_IN_PLACE, avg_v, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
avg_v = avg_v / denum
call MPI_ALLREDUCE(MPI_IN_PLACE, avg_w, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
avg_w = avg_w / denum
call MPI_ALLREDUCE(MPI_IN_PLACE, avg_fsd, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
avg_fsd = avg_fsd / denum
call MPI_ALLREDUCE(MPI_IN_PLACE, avg_lapc, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
avg_lapc = avg_lapc / denum
call MPI_ALLREDUCE(MPI_IN_PLACE, avg_ddx_c, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
avg_ddx_c = avg_ddx_c / denum
call MPI_ALLREDUCE(MPI_IN_PLACE, avg_ddy_c, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
avg_ddy_c = avg_ddy_c / denum
call MPI_ALLREDUCE(MPI_IN_PLACE, avg_ddz_c, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
avg_ddz_c = avg_ddz_c / denum
call MPI_ALLREDUCE(MPI_IN_PLACE, avg_d2dx_c, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
avg_d2dx_c = avg_d2dx_c / denum
call MPI_ALLREDUCE(MPI_IN_PLACE, avg_d2dy_c, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
avg_d2dy_c = avg_d2dy_c / denum
call MPI_ALLREDUCE(MPI_IN_PLACE, avg_d2dz_c, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
avg_d2dz_c = avg_d2dz_c / denum
call MPI_ALLREDUCE(MPI_IN_PLACE, fsd_avg_u, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
fsd_avg_u = fsd_avg_u / denum / avg_fsd
call MPI_ALLREDUCE(MPI_IN_PLACE, fsd_avg_v, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
fsd_avg_v = fsd_avg_v / denum / avg_fsd
call MPI_ALLREDUCE(MPI_IN_PLACE, fsd_avg_w, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
fsd_avg_w = fsd_avg_w / denum / avg_fsd
call MPI_ALLREDUCE(MPI_IN_PLACE, fsd_avg_nx, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
fsd_avg_nx = fsd_avg_nx / denum / avg_fsd
call MPI_ALLREDUCE(MPI_IN_PLACE, fsd_avg_ny, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
fsd_avg_ny = fsd_avg_ny / denum / avg_fsd
call MPI_ALLREDUCE(MPI_IN_PLACE, fsd_avg_nz, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
fsd_avg_nz = fsd_avg_nz / denum / avg_fsd
call MPI_ALLREDUCE(MPI_IN_PLACE, fsd_avg_t1, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
fsd_avg_t1 = fsd_avg_t1 / denum / avg_fsd
call MPI_ALLREDUCE(MPI_IN_PLACE, fsd_avg_absk, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
fsd_avg_absk = fsd_avg_absk / denum / avg_fsd
call MPI_ALLREDUCE(MPI_IN_PLACE, fsd_avg_divn, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
fsd_avg_divn = fsd_avg_divn / denum / avg_fsd
call MPI_ALLREDUCE(MPI_IN_PLACE, fsd_avg_sd, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
fsd_avg_sd = fsd_avg_sd / denum / avg_fsd
end subroutine m_terms_average_pass1
subroutine m_terms_calculate_pass2
integer :: i, j, k
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
c(i,j,k) = ( 1.0 - y(i,j,k) )
end do
end do
end do
call ddz ( ddz_c, c )
call ddy ( ddy_c, c )
call ddx ( ddx_c, c )
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
wrate(i,j,k) = ( rxn_rate ( c(i,j,k) ) )
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
fsd(i,j,k) = ( sqrt ( ( ( ((ddx_c(i,j,k))*(ddx_c(i,j,k))) + ((ddy_c(i,j,k))*(ddy_c(i,j,k))) ) + ((ddz_c(i,j,k))*(ddz_c(i,j,k))) ) ) )
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
nx(i,j,k) = ( ( - ddx_c(i,j,k) ) / fsd(i,j,k) )
end do
end do
end do
call ddx ( ddx_nx, nx )
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
ny(i,j,k) = ( ( - ddy_c(i,j,k) ) / fsd(i,j,k) )
end do
end do
end do
call ddy ( ddy_ny, ny )
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
nz(i,j,k) = ( ( - ddz_c(i,j,k) ) / fsd(i,j,k) )
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
vn(i,j,k) = ( ( ( u(i,j,k) * nx(i,j,k) ) + ( v(i,j,k) * ny(i,j,k) ) ) + ( w(i,j,k) * nz(i,j,k) ) )
end do
end do
end do
call ddz ( ddz_nz, nz )
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
divn(i,j,k) = ( ( ddx_nx(i,j,k) + ddy_ny(i,j,k) ) + ddz_nz(i,j,k) )
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
t7_fsd(i,j,k) = (((divn(i,j,k) - fsd_avg_divn(i)))*((divn(i,j,k) - fsd_avg_divn(i))))
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
fsd_avg_t7(i) = fsd_avg_t7(i) + t7_fsd(i,j,k) * fsd(i,j,k)
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
t6_fsd(i,j,k) = (((nx(i,j,k) - fsd_avg_nx(i)))*((nx(i,j,k) - fsd_avg_nx(i))))
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
fsd_avg_t6(i) = fsd_avg_t6(i) + t6_fsd(i,j,k) * fsd(i,j,k)
end do
end do
end do
call d2dy ( d2dy_c, c )
call d2dx ( d2dx_c, c )
call d2dz ( d2dz_c, c )
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
absk(i,j,k) = ( dabs ( divn(i,j,k) ) )
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
t8_fsd(i,j,k) = (((absk(i,j,k) - fsd_avg_absk(i)))*((absk(i,j,k) - fsd_avg_absk(i))))
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
fsd_avg_t8(i) = fsd_avg_t8(i) + t8_fsd(i,j,k) * fsd(i,j,k)
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
sd(i,j,k) = ( ( ( ( ( d2dx_c(i,j,k) + d2dy_c(i,j,k) ) + d2dz_c(i,j,k) ) * rod ) + wrate(i,j,k) ) / fsd(i,j,k) )
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
t1(i,j,k) = ( ( vn(i,j,k) + sd(i,j,k) ) * nx(i,j,k) )
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
t2_fsd(i,j,k) = ( (t1(i,j,k) - fsd_avg_t1(i)) * (nx(i,j,k) - fsd_avg_nx(i)) )
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
fsd_avg_t2(i) = fsd_avg_t2(i) + t2_fsd(i,j,k) * fsd(i,j,k)
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
t5_fsd(i,j,k) = (((t1(i,j,k) - fsd_avg_t1(i)))*((t1(i,j,k) - fsd_avg_t1(i))))
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
fsd_avg_t5(i) = fsd_avg_t5(i) + t5_fsd(i,j,k) * fsd(i,j,k)
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
t3_fsd(i,j,k) = ( (t1(i,j,k) - fsd_avg_t1(i)) * (divn(i,j,k) - fsd_avg_divn(i)) )
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
fsd_avg_t3(i) = fsd_avg_t3(i) + t3_fsd(i,j,k) * fsd(i,j,k)
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
t4_fsd(i,j,k) = ( (t1(i,j,k) - fsd_avg_t1(i)) * (absk(i,j,k) - fsd_avg_absk(i)) )
end do
end do
end do
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
fsd_avg_t4(i) = fsd_avg_t4(i) + t4_fsd(i,j,k) * fsd(i,j,k)
end do
end do
end do
end subroutine m_terms_calculate_pass2
subroutine m_terms_average_pass2 (nfiles)
integer :: nfiles
real*8 :: denum
denum=real(nfiles*nyp*nzp)
call MPI_ALLREDUCE(MPI_IN_PLACE, fsd_avg_t2, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
fsd_avg_t2 = fsd_avg_t2 / denum / avg_fsd
call MPI_ALLREDUCE(MPI_IN_PLACE, fsd_avg_t3, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
fsd_avg_t3 = fsd_avg_t3 / denum / avg_fsd
call MPI_ALLREDUCE(MPI_IN_PLACE, fsd_avg_t4, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
fsd_avg_t4 = fsd_avg_t4 / denum / avg_fsd
call MPI_ALLREDUCE(MPI_IN_PLACE, fsd_avg_t5, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
fsd_avg_t5 = fsd_avg_t5 / denum / avg_fsd
call MPI_ALLREDUCE(MPI_IN_PLACE, fsd_avg_t6, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
fsd_avg_t6 = fsd_avg_t6 / denum / avg_fsd
call MPI_ALLREDUCE(MPI_IN_PLACE, fsd_avg_t7, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
fsd_avg_t7 = fsd_avg_t7 / denum / avg_fsd
call MPI_ALLREDUCE(MPI_IN_PLACE, fsd_avg_t8, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)
fsd_avg_t8 = fsd_avg_t8 / denum / avg_fsd
end subroutine m_terms_average_pass2
subroutine m_terms_write_result
integer :: i
open (200, file="qEdge_X.dat")
write (200,*) output_header
do i=1,nxp
write (200,'(30e20.10)') real(i)*hxp, avg_c(i), avg_u(i), avg_v(i), avg_w(i), avg_fsd(i), avg_lapc(i), avg_ddx_c(i), avg_ddy_c(i), avg_ddz_c(i), avg_d2dx_c(i), avg_d2dy_c(i), avg_d2dz_c(i), fsd_avg_u(i), fsd_avg_v(i), fsd_avg_w(i), fsd_avg_nx(i), fsd_avg_ny(i), fsd_avg_nz(i), fsd_avg_t1(i), fsd_avg_t2(i), fsd_avg_t3(i), fsd_avg_t4(i), fsd_avg_t5(i), fsd_avg_t6(i), fsd_avg_t7(i), fsd_avg_t8(i), fsd_avg_absk(i), fsd_avg_divn(i), fsd_avg_sd(i)
end do
close (200)
end subroutine m_terms_write_result
end module m_terms