program main

  real(4), parameter :: PI = 3.141592

  integer :: imax
  integer :: jmax
  integer :: kmax

  real(4) :: X0, DelX
  real(4) :: Y0, DelY
  real(4) :: Z0, DelZ
  real(4), allocatable :: x_X(:)
  real(4), allocatable :: y_Y(:)
  real(4), allocatable :: xyz_Z(:,:,:)
  real(4), allocatable :: xyz_Val(:,:,:)

  integer       :: FU
  character(64) :: FileName = 'out/test2.vtk'

  integer :: i
  integer :: j
  integer :: k

  imax = 21
  jmax = 21
  kmax =  1
  kmax = 21

  allocate( x_X  ( imax ) )
  allocate( y_Y  ( jmax ) )
  allocate( xyz_Z( imax, jmax, kmax ) )
  allocate( xyz_Val( imax, jmax, kmax ) )

  DelX = 2.0*PI/real(imax) ; X0 = -DelX*imax/2
  DelY = 2.0*PI/real(jmax) ; Y0 = -DelY*imax/2
  DelZ = 2.0*PI/real(kmax) ; Z0 = 0.0
  do i = 1, imax
    x_X(i) = X0 + DelX * (i-1)
  end do
  do j = 1, jmax
    y_Y(j) = Y0 + DelY * (j-1)
  end do
  do k = 1, kmax
!!$    xyz_Z(:,:,k) = DelZ + Z0 + DelZ * (k-1)
    do j = 1, jmax
      do i = 1, imax
!!$        xyz_Z(i,j,k) = 0.1 * ( x_X(i)**2 + y_Y(j)**2 ) + DelZ * (k-1)
        xyz_Z(i,j,k) = cos( 2.0 * x_X(i) ) * cos( y_Y(j)/2.0 ) + DelZ * (k-1)
      end do
    end do
  end do
  do k = 1, kmax
    do j = 1, jmax
      do i = 1, imax
        xyz_Val(i,j,k) = cos(x_X(i))*cos(y_Y(j))*xyz_Z(i,j,k)
      end do
    end do
  end do



  FU = 50
  open( FU, file = FileName, status='unknown' )

  write( FU, '(a)' ) '# vtk DataFile Version 1.0'
  write( FU, '(a)' ) trim( FileName )
  write( FU, '(a)' ) 'ASCII'

  write( FU, '(a)' ) 'DATASET STRUCTURED_GRID'
  write( FU, '(a,3(a,i3))' ) 'DIMENSIONS', ' ', imax, ' ', jmax, ' ', kmax
  write( FU, '(a,i6,a)' ) 'POINTS ', imax*jmax*kmax, ' float'
  do k = 1, kmax
    do j = 1, jmax
      do i = 1, imax
        write( FU, '(3(f5.2,a))' ) x_X(i), ' ', y_Y(j), ' ', xyz_Z(i,j,k), ' '
      end do
    end do
  end do

  write( FU, '(a,i8)' ) 'POINT_DATA ', imax*jmax*kmax

  write( FU, '(a)' ) 'SCALARS intensity float'
  write( FU, '(a)' ) 'LOOKUP_TABLE default'

  do k = 1, kmax
    do j = 1, jmax
      do i = 1, imax
        write( FU, '(f5.2)' ) xyz_Val(i,j,k)
      end do
    end do
  end do

  close( FU )

end program main
