!=======================================
! Binary output of scalar field data
! Latest update: 2014.1.10
! Reference : http://wwwoa.ees.hokudai.ac.jp/people/numa/enshu/fileio.html
!=======================================
program scalar_field_binary_out
        implicit none
        real(8),dimension(:,:),allocatable :: point_coord  ! three-dimensional coordinates (x,y,z)
        real(8),dimension(:),allocatable :: field_data ! scalar function f(x,y,z)
        real(8) :: Lx,Ly,Lz
        real(8) :: cx,cy,cz
        real(8) :: dx,dy,dz
        integer :: nx,ny,nz
        real(8) :: x,y,z,r
        integer :: dim ! the number of spatial dimension
        integer :: size ! the number of data points
        integer :: reclen ! record length
        integer :: i,j,k,d,count

        write(*,'(a)') "********************************************"
        write(*,'(a)') "     Binary output of scalar field data     "
        write(*,'(a)') "********************************************"

        write(*,'(a)') "********************************************"
        write(*,'(a)') "               Initialization               "
        write(*,'(a)') ""

        dim = 3

        Lx = 5.0d0
        Ly = Lx
        Lz = Lx

        cx = Lx/2.0d0
        cy = Ly/2.0d0
        cz = Lz/2.0d0

        nx = 51
        ny = nx 
        nz = nx

        dx = Lx/(nx-1)
        dy = Ly/(ny-1)
        dz = Lz/(nz-1)

        size = nx*ny*nz

        allocate(point_coord(dim,size))
        allocate(field_data(size))

        point_coord(:,:) = 0.0d0
        field_data(:) = 0.0d0

        write(*,*) "dim =",3
        write(*,*) "Lx =",Lx,"Ly =",Ly,"Lz =",Lz
        write(*,*) "cx =",cx,"cy =",cy,"cz =",cz
        write(*,*) "dx =",dx,"dy =",dy,"dz =",dz
        write(*,*) "nx =",nx,"ny =",ny,"z =",nz
        write(*,*) "The number of data points :",size
        write(*,'(a)') ""
        write(*,'(a)') "            End of initialization           "
        write(*,'(a)') "********************************************"

        write(*,'(a)') "********************************************"
        write(*,'(a)') "  Generation of grid points and field data  "
        write(*,'(a)') ""

        count = 1
        z = 0.0d0
        do k=1,nz
                y = 0.0d0 
                do j=1,ny
                        x = 0.0d0
                        do i=1,nx
                                point_coord(1,count) = x ! set x
                                point_coord(2,count) = y ! set y
                                point_coord(3,count) = z ! set z
                                
                                r = sqrt((x-cx)*(x-cx)+(y-cy)*(y-cy)+(z-cz)*(z-cz)) ! distance from the central point
                                field_data(count) = exp(-r*r) ! Set scalar function value. The centeral point is (cx,cy,cz).
                                
                                count = count + 1
                                x = x + dx
                        end do
                        y = y + dy
                end do
                z = z + dz
        end do 

        write(*,'(a)') "           End of data generation           "
        write(*,'(a)') "********************************************"

        write(*,'(a)') "********************************************"
        write(*,'(a)') "                Binary output               "
        write(*,'(a)') ""
        
        INQUIRE(IOLENGTH=reclen) (point_coord(d,i),d=1,3),field_data(1)
        write(*,*) "Record length:",reclen,"Byte"

        open(10,file='scalar_field_3d.dat',form='unformatted',access='direct',recl=reclen)
        do i=1,size
            write(10,rec=i) (point_coord(d,i),d=1,3),field_data(i)
        end do
        close(10)

        write(*,*) ""
        write(*,'(a)') "            End of Binary output            "
        write(*,'(a)') "********************************************"
        
        deallocate(point_coord)
        deallocate(field_data)

        write(*,'(a)') "********************************************"
        write(*,'(a)') "              End of program                "
        write(*,'(a)') "********************************************"

        stop
end program scalar_field_binary_out