module io

    use elements
    use parameters
    use atoms
    use box

    implicit none

    integer :: outfilenum = 0
    character(len=100) :: outfiles(10)

   public 
   contains

    subroutine get_out_file(filename)

        implicit none

        character(len=100), intent(in) :: filename
        character(len=100) :: temp_outfile
        character(len=1) :: overwrite
        logical :: file_exists

        !If no filename is provided then this function is called with none and prompts user input
        if (filename=='none') then 
            print *, "Please specify a filename or extension to output to:"
            read(*,*) temp_outfile
        else
            temp_outfile = filename
        end if

        !Infinite loop which only exists if user provides valid filetype
        overwrite = 'r'
        do while(.true.)

            !Check to see if file exists, if it does then ask user if they would like to overwrite the file
            inquire(file=trim(temp_outfile), exist=file_exists)
            if (file_exists) then 
                if (overwrite == 'r') print *, "File ", trim(temp_outfile), " already exists. Would you like to overwrite? (Y/N)"
                read(*,*) overwrite
                if((scan(overwrite, "n") > 0).or.(scan(overwrite, "N") > 0)) then 
                    print *, "Please specify a new filename with extension:"
                    read(*,*) temp_outfile
                else if((scan(overwrite, "y") > 0).or.(scan(overwrite, "Y") > 0)) then 
                    continue
                else
                    print *, "Please pick either y or n"
                    read(*,*) overwrite
                end if

            end if

            if (scan(temp_outfile,'.',.true.) == 0) then 
                print *, "No extension included on filename, please type a full filename that includes an extension."
                read(*,*) temp_outfile
                cycle
            end if
            select case(temp_outfile(scan(temp_outfile,'.',.true.)+1:))
            case('xyz', 'lmp', 'vtk', 'mb')
                outfilenum=outfilenum+1
                outfiles(outfilenum) =  temp_outfile
                exit
            case default    
                print *, "File type: ", trim(temp_outfile(scan(temp_outfile,'.',.true.):)), "not currently accepted. ", &
                         "please input a filename with extension from following list: xyz, lmp, vtk."
                read(*,*) temp_outfile

            end select
        end do

    end subroutine get_out_file 


    subroutine write_out
        !This subroutine loops over alll of the outfile types defined and calls the correct writing subroutine

        integer :: i 

        !Find max esize which will be needed later
        call set_max_esize

        do i = 1, outfilenum
            !Pull out the extension of the file and call the correct write subroutine
            select case(trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))))
            case('xyz')
                call write_xyz(outfiles(i))
            case('lmp')
                call write_lmp(outfiles(i))
            case('vtk')
                call write_vtk(outfiles(i))
            case('mb')
                call write_mb(outfiles(i))
            case default
                print *, "The extension ", trim(adjustl(outfiles(i)(scan(outfiles(i),'.',.true.)+1:))), &
                         " is not accepted for writing. Please select from: xyz and try again"
                stop 

            end select
        end do 
    end subroutine write_out



    subroutine write_xyz(file)
        !This is the simplest visualization subroutine, it writes out all nodal positions and atom positions to an xyz file
        character(len=100), intent(in) :: file 

        integer :: node_num, i, inod, ibasis

        open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind')

        !Calculate total node number
        node_num=0
        do i = 1, ele_num
            node_num = node_num + basisnum(lat_ele(i))*ng_node(lat_ele(i))
        end do

        !Write total number of atoms + elements
        write(11, '(i16)') node_num+atom_num
         
        !Write comment line
        write(11, '(a)') "#Node + atom file created using cacmb"

        !Write nodal positions
        do i = 1, ele_num
            do inod = 1, ng_node(lat_ele(i))
                do ibasis = 1, basisnum(lat_ele(i))
                        write(11, '(a, 3f23.15)') basis_type(ibasis,lat_ele(i)), r_node(:,ibasis,inod,i)
                end do 
            end do 
        end do 

        !Write atom positions
        do i = 1, atom_num
            write(11, '(a, 3f23.15)') type_atom(i), r_atom(:,i) 
        end do

        !Finish writing 
        close(11)
    end subroutine write_xyz

    subroutine write_lmp(file)
        !This subroutine writes out a .lmp style dump file
        character(len=100), intent(in) :: file
        integer :: write_num, i, iatom, type_interp(max_basisnum*max_esize**3)
        real(kind=dp) :: r_interp(3, max_basisnum*max_esize**3), mass

        open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind')

        !Comment line
        write(11, '(a)') '# lmp file made with cacmb'
        write(11, '(a)')
        !Calculate total atom number
        write_num = atom_num
        do i = 1,ele_num
            if(type_ele(i) == 'fcc') write_num = write_num + size_ele(i)**3
        end do
        !Write total number of atoms + elements
        write(11, '(i16, a)') write_num, ' atoms'
        !Write number of atom types 
        write(11, '(i16, a)') atom_types, ' atom types'

        write(11,'(a)') ' '
        !Write box bd
        write(11, '(2f23.15, a)') box_bd(1:2), ' xlo xhi'
        write(11, '(2f23.15, a)') box_bd(3:4), ' ylo yhi'
        write(11, '(2f23.15, a)') box_bd(5:6), ' zlo zhi'

        !Masses
        write(11, '(a)') 'Masses'

        write(11, '(a)') ' '
        do i =1, atom_types
            call atommass(type_to_name(i),mass)
            write(11, '(i16, f23.15)') i, mass
        end do
        write(11, '(a)') ' '

        !Write atom positions
        write(11, '(a)') 'Atoms'
        write(11, '(a)') ' '
        do i = 1, atom_num
            write(11, '(2i16, 3f23.15)') i, type_atom(i), r_atom(:,i) 
        end do

        !Write refined element atomic positions
        do i = 1, ele_num
            call interpolate_atoms(type_ele(i), size_ele(i), lat_ele(i), r_node(:,:,:,i), type_interp, r_interp)
            select case(trim(adjustl(type_ele(i))))
            case('fcc')
                do iatom = 1, basisnum(lat_ele(i))*size_ele(i)**3
                    write(11, '(2i16, 3f23.15)') atom_num+iatom, type_interp(iatom), r_interp(:,iatom)
                end do
            end select
        end do
    end subroutine write_lmp

    subroutine write_vtk(file)
        !This subroutine writes out a vtk style dump file
        integer :: i, j, inod, ibasis
        character(len=100), intent(in) :: file

1       format('# vtk DataFile Version 4.0.1', / &
              'CAC output -- cg', / &
              'ASCII')
11      format('# vtk DataFile Version 4.0.1', / &
               'CACmb output -- atoms', / &
               'ASCII')
2       format('DATASET UNSTRUCTURED_GRID')
3       format('POINTS', i16, ' float')
4       format(/'CELLS', 2i16)
5       format(/'CELL_TYPES', i16)
12      format(/'CELL_DATA', i16)
16      format(/'POINT_DATA', i16)
17      format('SCALARS weight float', / &
               'LOOKUP_TABLE default')
18      format('SCALARS atom_type float', / &
               'LOOKUP_TABLE default')

20      format('SCALARS lattice_type float', /&
               'LOOKUP_TABLE default')

        !First we write the vtk file containing the atoms
        open(unit=11, file='atoms_'//trim(adjustl(file)), action='write', status='replace',position='rewind')

        write(11, 11)
        write(11, 2)
        write(11, 3) atom_num
        do i = 1, atom_num
            write(11, '(3f23.15)') r_atom(:,i)
        end do 
        write(11,4) atom_num, atom_num*2
        do i = 1, atom_num
            write(11, '(2i16)') 1, i-1
        end do
        write(11, 5) atom_num
        do i = 1, atom_num
            write(11, '(i16)') 1 
        end do
        write(11, 16) atom_num
        write(11, 18)
        do i = 1, atom_num
            write(11, '(i16)') type_atom(i)
        end do
        close(11)

        open(unit=11, file='cg_'//trim(adjustl(file)), action='write', status='replace',position='rewind')
        write(11,1)
        write(11,2)
        write(11,3) node_num
        do i = 1, ele_num
            do inod=1, ng_node(lat_ele(i))  
                do ibasis = 1, basisnum(lat_ele(i))
                    write(11, '(3f23.1)') sum(r_node(:,:,inod,i),2)/basisnum(lat_ele(i))
                end do 
            end do 
        end do
        write(11, 4) ele_num, ele_num + node_num
        do i =1, ele_num
            write(11, '(9i16)') ng_node(lat_ele(i)), (j, j = (i-1)*ng_node(lat_ele(i)), i*ng_node(lat_ele(i))-1)
        end do
        write(11,5) ele_num
        do i = 1, ele_num
            if(trim(adjustl(type_ele(i))) == 'fcc') write(11, '(i16)') 12
        end do
        write(11,12) ele_num
        write(11,20)
        do i = 1, ele_num
            write(11, '(i16)') lat_ele(i)
        end do
        close(11)
    end subroutine

    subroutine write_mb(file)

        !This subroutine writes the cacmb formatted file which provides necessary information for building models
        character(len=100), intent(in) :: file

        integer :: i, j, inod, ibasis

        !Open the .mb file for writing
        open(unit=11, file=trim(adjustl(file)), action='write', status='replace',position='rewind')

        !First write the box boundary information
        !Write the global box boundaries
        write(11,*) box_bd(:)
        !Write the number of sub_boxes in the system
        write(11,*) sub_box_num
        !For every subbox write the orientation and sub box boundary
        do i = 1, sub_box_num
            write(11,*) sub_box_ori(:,:,i)
            write(11,*) sub_box_bd(:,i)
        end do

        !Write the number of atom types in the current model and all of their names
        write(11,*) atom_types, (type_to_name(i), i=1, atom_types)
        !Write the number of lattice_types, basisnum and number of nodes for each lattice type
        write(11,*) lattice_types, (basisnum(i), i = 1, lattice_types), (ng_node(i), i = 1, lattice_types)
        !Now for every lattice type write the basis atom types
        write(11,*) ((basis_type(i,j), i = 1, basisnum(j)), j = 1, lattice_types)

        !Now write the numbers of elements and atoms
        write(11,*) atom_num, ele_num

        !Write out atoms first
        do i = 1, atom_num
            write(11,*) type_atom(i), r_atom(:,i)
        end do

        !Write out the elements, this is written in two stages, one line for the element and then 1 line for 
        !every basis at every node
        do i = 1, ele_num
            write(11, *) i, lat_ele(i), size_ele(i), type_ele(i)
            do inod = 1, ng_node(lat_ele(i))
                do ibasis =1, basisnum(lat_ele(i))
                    write(11,*) inod, ibasis, r_node(:, ibasis, inod, i)
                end do
            end do
        end do

    end subroutine write_mb
end module io