program test_calculate

use m_parameters

use m_calculate

implicit none

real*8, allocatable, dimension(:,:,:) :: aaa,bbb,ccc,ddd

integer :: ierr

real*8 :: xx,yy,zz,fxyz

integer :: i,j,k

real*4, dimension(2) :: startt
real*4, dimension(2) :: endt
real*4 :: result

nxp = 512
nyp = 256
nzp = 256

l_0 = 2.0
hyp=l_0*pi/REAL(nyp)
hxp=hyp
hzp=hyp

allocate(aaa(nxp, nyp, nzp), stat=ierr); aaa=0.0
allocate(bbb(nxp, nyp, nzp), stat=ierr); bbb=0.0
allocate(ccc(nxp, nyp, nzp), stat=ierr); ccc=0.0
allocate(ddd(nxp, nyp, nzp), stat=ierr); ddd=0.0


call m_calculate_init


write(*,*) "dimension - ", nxp, nyp, nzp
write(*,*) "spacing   - ", hxp, hyp, hzp

do k = 1,nzp
do j = 1,nyp
do i = 1,nxp

xx = i * hxp
yy = j * hyp
zz = k * hzp

aaa(i,j,k) = sin(1.1 * xx) * sin(3.0 * yy) * sin(2.0 * zz)

end do
end do
end do

write(*,*) "test ddx"

call ETIME(startt, result)


call ddx(bbb, aaa)
call ETIME(endt, result)
print *, "usert", endt(1) - startt(1)
print *, "systt", endt(2) - startt(2)

do k = 1,nzp
do j = 1,nyp
do i = nxp/4,3*nxp/4

xx = i * hxp
yy = j * hyp
zz = k * hzp
fxyz = (1.1 * cos(1.1 * xx) * sin(3.0 * yy) * sin(2.0 * zz))
ccc(i,j,k) = (bbb(i,j,k) - fxyz) / (fxyz)

if ((abs(ccc(i,j,k)) > 0.01 ) .and. (abs(bbb(i,j,k)) > 1.0e-14 )) write(*,*) bbb(i,j,k), fxyz

end do
end do
end do

write(*,*) "relerr_min_max", minval(ccc), maxval(ccc)



write(*,*) "test ddy"

call ETIME(startt, result)


call ddy(bbb, aaa)
call ETIME(endt, result)
print *, "usert", endt(1) - startt(1)
print *, "systt", endt(2) - startt(2)

do k = 1,nzp
do j = 1,nyp
do i = nxp/4,3*nxp/4

xx = i * hxp
yy = j * hyp
zz = k * hzp
fxyz = (3.0 * sin(1.1 * xx) * cos(3.0 * yy) * sin(2.0 * zz))
ccc(i,j,k) = (bbb(i,j,k) - fxyz) / (fxyz)

if ((abs(ccc(i,j,k)) > 0.01 ) .and. (abs(bbb(i,j,k)) > 1.0e-14 )) write(*,*) bbb(i,j,k), fxyz
end do
end do
end do

write(*,*) "relerr_min_max", minval(ccc), maxval(ccc)





write(*,*) "test ddz"

call ETIME(startt, result)


call ddz(bbb, aaa)
call ETIME(endt, result)
print *, "usert", endt(1) - startt(1)
print *, "systt", endt(2) - startt(2)

do k = 1,nzp
do j = 1,nyp
do i = nxp/4,3*nxp/4

xx = i * hxp
yy = j * hyp
zz = k * hzp
fxyz = (2.0 * sin(1.1 * xx) * sin(3.0 * yy) * cos(2.0 * zz))
ccc(i,j,k) = (bbb(i,j,k) - fxyz) / (fxyz)

if ((abs(ccc(i,j,k)) > 0.01 ) .and. (abs(bbb(i,j,k)) > 1.0e-14 )) write(*,*) bbb(i,j,k), fxyz
end do
end do
end do

write(*,*) "relerr_min_max", minval(ccc), maxval(ccc)





write(*,*) "test d2dx"

call ETIME(startt, result)


call d2dx(bbb, aaa)
call ETIME(endt, result)
print *, "usert", endt(1) - startt(1)
print *, "systt", endt(2) - startt(2)

do k = 1,nzp
do j = 1,nyp
do i = nxp/4,3*nxp/4

xx = i * hxp
yy = j * hyp
zz = k * hzp
fxyz = -(1.21 * sin(1.1 * xx) * sin(3.0 * yy) * sin(2.0 * zz))
ccc(i,j,k) = (bbb(i,j,k) - fxyz) / (fxyz)

if ((abs(ccc(i,j,k)) > 0.01 ) .and. (abs(fxyz) > 1.0e-14 )) write(*,*) bbb(i,j,k), fxyz

end do
end do
end do

write(*,*) "relerr_min_max", minval(ccc), maxval(ccc)



write(*,*) "test d2dy"

call ETIME(startt, result)


call d2dy(bbb, aaa)
call ETIME(endt, result)
print *, "usert", endt(1) - startt(1)
print *, "systt", endt(2) - startt(2)

do k = 1,nzp
do j = 1,nyp
do i = nxp/4,3*nxp/4

xx = i * hxp
yy = j * hyp
zz = k * hzp
fxyz = -(9.0 * sin(1.1 * xx) * sin(3.0 * yy) * sin(2.0 * zz))
ccc(i,j,k) = (bbb(i,j,k) - fxyz) / (fxyz)

if ((abs(ccc(i,j,k)) > 0.01 ) .and. (abs(fxyz) > 1.0e-14 )) write(*,*) bbb(i,j,k), fxyz
end do
end do
end do

write(*,*) "relerr_min_max", minval(ccc), maxval(ccc)


write(*,*) "test d2dz"

call ETIME(startt, result)


call d2dz(bbb, aaa)
call ETIME(endt, result)
print *, "usert", endt(1) - startt(1)
print *, "systt", endt(2) - startt(2)


do k = 1,nzp
do j = 1,nyp
do i = nxp/4,3*nxp/4

xx = i * hxp
yy = j * hyp
zz = k * hzp
fxyz = -(4.0 * sin(1.1 * xx) * sin(3.0 * yy) * sin(2.0 * zz))
ccc(i,j,k) = (bbb(i,j,k) - fxyz) / (fxyz)

if ((abs(ccc(i,j,k)) > 0.01 ) .and. (abs(fxyz) > 1.0e-14 )) write(*,*) bbb(i,j,k), fxyz
end do
end do
end do

write(*,*) "relerr_min_max", minval(ccc), maxval(ccc)





call m_calculate_finalize


deallocate(aaa)
deallocate(bbb)
deallocate(ccc)
deallocate(ddd)

end program