module m_terms

use m_openmpi
use m_parameters
use m_arrays
use m_calculate

implicit none

logical, parameter :: pass2_required=.false.
character (len = *), parameter :: output_header="x avg_c avg_divn avg_fsd"
real*8, allocatable, dimension(:) :: avg_c
real*8, allocatable, dimension(:) :: avg_fsd
real*8, allocatable, dimension(:) :: avg_divn

! - file_handles and mpi_infos
! integer(kind=MPI_INTEGER_KIND) :: {field_name}_fh
! integer(kind=MPI_INTEGER_KIND) :: {field_name}_info

! - buffer
! real*8, allocatable, dimension(:,:,:) :: {field_name}_export_array

real*8, allocatable, dimension(:,:,:) :: xyzbuffer0
real*8, allocatable, dimension(:,:,:) :: xyzbuffer1
real*8, allocatable, dimension(:,:,:) :: xyzbuffer2
real*8, allocatable, dimension(:,:,:) :: xyzbuffer3
real*8, allocatable, dimension(:,:,:) :: xyzbuffer4

contains

subroutine m_terms_init

integer :: ierr

allocate(avg_c(nxp), stat=ierr) ; avg_c = 0.
allocate(avg_fsd(nxp), stat=ierr) ; avg_fsd = 0.
allocate(avg_divn(nxp), stat=ierr) ; avg_divn = 0.

! init
! call MPI_INFO_CREATE({field_name}_info, mpi_err)
! call MPI_FILE_OPEN(MPI_COMM_TASK,{field_name}_fname,MPI_MODE_WRONLY+MPI_MODE_CREATE,{field_name}_info,{field_name}_fh,mpi_err)
! allocate({field_name}_export_array({1}), stat=ierr) ; {field_name}_export_array = 0.

allocate(xyzbuffer0(nxp,nyp,nzp), stat=ierr) ; xyzbuffer0 = 0.
allocate(xyzbuffer1(nxp,nyp,nzp), stat=ierr) ; xyzbuffer1 = 0.
allocate(xyzbuffer2(nxp,nyp,nzp), stat=ierr) ; xyzbuffer2 = 0.
allocate(xyzbuffer3(nxp,nyp,nzp), stat=ierr) ; xyzbuffer3 = 0.
allocate(xyzbuffer4(nxp,nyp,nzp), stat=ierr) ; xyzbuffer4 = 0.

end subroutine m_terms_init


subroutine m_terms_finalize

deallocate(avg_c)
deallocate(avg_fsd)
deallocate(avg_divn)

! finalize
! call MPI_FILE_CLOSE({field_name}_fh, mpi_err)
! call MPI_INFO_FREE({field_name}_info, mpi_err)
! deallocate({field_name}_export_array)

deallocate(xyzbuffer0)
deallocate(xyzbuffer1)
deallocate(xyzbuffer2)
deallocate(xyzbuffer3)
deallocate(xyzbuffer4)

end subroutine m_terms_finalize


subroutine m_terms_calculate_pass1

integer :: i, j, k


! ( 1.0 - y(i,j,k) )
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
xyzbuffer3(i,j,k) = ( 1.0 - y(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) + xyzbuffer3(i,j,k)
end do
end do
end do

call ddx ( xyzbuffer2, xyzbuffer3 )
call ddy ( xyzbuffer1, xyzbuffer3 )
call ddz ( xyzbuffer0, xyzbuffer3 )

! ( 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))) ) ) )
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
xyzbuffer3(i,j,k) = ( sqrt ( ( ( ((xyzbuffer2(i,j,k))*(xyzbuffer2(i,j,k))) + ((xyzbuffer1(i,j,k))*(xyzbuffer1(i,j,k))) ) + ((xyzbuffer0(i,j,k))*(xyzbuffer0(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) + xyzbuffer3(i,j,k)
end do
end do
end do


! ( ( - ddx_c(i,j,k) ) / fsd(i,j,k) )
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
xyzbuffer4(i,j,k) = ( ( - xyzbuffer2(i,j,k) ) / xyzbuffer3(i,j,k) )
end do
end do
end do

call ddx ( xyzbuffer2, xyzbuffer4 )

! ( ( - ddy_c(i,j,k) ) / fsd(i,j,k) )
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
xyzbuffer4(i,j,k) = ( ( - xyzbuffer1(i,j,k) ) / xyzbuffer3(i,j,k) )
end do
end do
end do


! ( ( - ddz_c(i,j,k) ) / fsd(i,j,k) )
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
xyzbuffer1(i,j,k) = ( ( - xyzbuffer0(i,j,k) ) / xyzbuffer3(i,j,k) )
end do
end do
end do

call ddy ( xyzbuffer3, xyzbuffer4 )
call ddz ( xyzbuffer0, xyzbuffer1 )

! ( ( ddx_nx(i,j,k) + ddy_ny(i,j,k) ) + ddz_nz(i,j,k) )
do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
xyzbuffer4(i,j,k) = ( ( xyzbuffer2(i,j,k) + xyzbuffer3(i,j,k) ) + xyzbuffer0(i,j,k) )
end do
end do
end do

! copy to array for export
! do k = 1, nzp
! do j = 1, nyp
! do i = {}, {}
! {}(i,j,k) = {}(i,j,k)
! end do
! end do
! end do

! after calculate
! count = {x1-x0+1} * ny * nz
! offset = time_count * count
! call MPI_FILE_WRITE_AT({field_name}_fh, offset, {field_name}_export_buffer, count, MPI_REAL8, mpi_status, mpi_err)


do k = 1, nzp
do j = 1, nyp
do i = 1, nxp
avg_divn(i) = avg_divn(i) + xyzbuffer4(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_divn, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)

avg_divn = avg_divn  / denum


call MPI_ALLREDUCE(MPI_IN_PLACE, avg_fsd, nxp, MPI_REAL8, MPI_SUM, MPI_COMM_TASK, mpi_err)

avg_fsd = avg_fsd  / denum


end subroutine m_terms_average_pass1


subroutine m_terms_calculate_pass2

integer :: i, j, k



end subroutine m_terms_calculate_pass2


subroutine m_terms_average_pass2 (nfiles)

integer :: nfiles
real*8 :: denum

denum=real(nfiles*nyp*nzp)



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,'(4e20.10)') real(i)*hxp, avg_c(i), avg_divn(i), avg_fsd(i)
end do

close (200)


end subroutine m_terms_write_result


end module m_terms