!----------------------------------
MODULE DECLARE
IMPLICIT NONE
INTEGER, PARAMETER :: sp = selected_real_kind(6, 37)
INTEGER, PARAMETER :: dp = selected_real_kind(15, 307)
CHARACTER (LEN=90) :: system_name, print_format, solver, errormsg
INTEGER :: number_of_particles, number_of_walls, ios, step, total_steps, step_to_print, temporary_step
LOGICAL :: gravity_is_on = .FALSE.
REAL(KIND=DP) :: delta_time_to_print,delta_time,final_time,time,elasticity_modulus,poisson_coefficient,density
REAL(KIND=DP), DIMENSION(:), SAVE, ALLOCATABLE :: radius, mass
REAL(KIND=DP), DIMENSION(:,:), SAVE, ALLOCATABLE :: force_pre_iter, position_pre_iter, velocity_pre_iter, wall_point, wall_normal, &
    & position, preceding_position, velocity, preceding_velocity, force, preceding_force
REAL(KIND=DP), PARAMETER :: pi = 3.14159265358979_dp
REAL(KIND=DP), DIMENSION(3), PARAMETER :: gravity = (/ 0.0_dp, -9.7803267714_dp, 0.0_dp /)
END MODULE DECLARE
!----------------------------------
MODULE OUTPUT
IMPLICIT NONE
CONTAINS
SUBROUTINE dispatch_for_print
USE DECLARE
IMPLICIT NONE
SELECT CASE (print_format)
CASE ("csv_format")
    CALL print_csv
CASE ("vtk_ascii_format")
    CALL print_vtk_ascii
CASE ("vtk_binary_format")
    CALL print_vtk_binary
CASE DEFAULT
    print *, "Command not recognized; please re-enter"
END SELECT
END SUBROUTINE dispatch_for_print
!-----------------------------
function itoa(i) result(res)
  character(:),ALLOCATABLE :: res
  INTEGER,intent(in) :: i
  character(range(i)+2) :: tmp
  write(tmp,'(i0)') i
  res = trim(tmp)
END function
!-----------------------------
SUBROUTINE print_vtk_binary
USE DECLARE
IMPLICIT NONE
INTEGER :: print_vtk=0, print_number=0, i
OPEN(newunit=print_vtk, access='stream', form='unformatted', FILE='results'//itoa(print_number)//'.vtk', &
    & STATUS='replace', convert='big_endian', IOSTAT= ios,iomsg=errormsg)
IF (ios/=0) THEN
    print *, "error opening print_vtk_binary", print_vtk, errormsg
    STOP
ENDIF
WRITE(print_vtk)'# vtk DataFile Version 3.0',NEW_LINE('a')
WRITE(print_vtk)TRIM(system_name),NEW_LINE('a')
WRITE(print_vtk)'BINARY',NEW_LINE('a')
WRITE(print_vtk)'DATASET UNSTRUCTURED_GRID',NEW_LINE('a')
WRITE(print_vtk)'POINTS '//itoa(number_of_particles)//' double',NEW_LINE('a')
WRITE(print_vtk)(real(position(i,1),kind=8),real(position(i,2),kind=8),real(position(i,3),kind=8),i=1,number_of_particles)
WRITE(print_vtk)'POINT_DATA '//itoa(number_of_particles),NEW_LINE('a')
WRITE(print_vtk)'SCALARS Particle_radius double',NEW_LINE('a')
WRITE(print_vtk)'LOOKUP_TABLE default',NEW_LINE('a')
WRITE(print_vtk)real(radius(number_of_particles),kind=8)
WRITE(print_vtk)'VECTORS Particle_velocity double',NEW_LINE('a')
WRITE(print_vtk)(real(velocity(i,1),kind=8),real(velocity(i,2),kind=8),real(velocity(i,3),kind=8),i=1,number_of_particles)
CLOSE(unit=print_vtk)
print_number=print_number+1
END SUBROUTINE print_vtk_binary
!------------------------------------
SUBROUTINE print_vtk_ascii
USE DECLARE
IMPLICIT NONE
INTEGER :: print_vtk=0,print_number=0,i=0
OPEN(newunit=print_vtk, action='write', FILE='results'//itoa(print_number)//'.vtk', STATUS='replace',form="formatted", & 
    & IOSTAT= ios,iomsg=errormsg)
IF (ios/=0) THEN
    print *, "error opening print_vtk_ascii", print_vtk, errormsg
    STOP
ENDIF
WRITE(print_vtk,"(A)")'# vtk DataFile Version 3.0'
WRITE(print_vtk,"(A)")TRIM(system_name)
WRITE(print_vtk,"(A)")'ASCII'
WRITE(print_vtk,"(A)")'DATASET UNSTRUCTURED_GRID'
WRITE(print_vtk,"(A,I0,A)")'POINTS '//itoa(number_of_particles)//' double'
WRITE(print_vtk,"(E17.10,E17.10,E17.10)")(preceding_position(i,1),preceding_position(i,2),preceding_position(i,3),i=1,number_of_particles)
WRITE(print_vtk,"(A,I0)")'POINT_DATA '//itoa(number_of_particles)
WRITE(print_vtk,"(A)")'SCALARS Particle_radius double'
WRITE(print_vtk,"(A)")'LOOKUP_TABLE default'
WRITE(print_vtk,"(E17.10)")(radius(i),i=1,number_of_particles)
WRITE(print_vtk,"(A)")'VECTORS Particle_velocity double'
WRITE(print_vtk,"(E17.10,E17.10,E17.10)")(preceding_velocity(i,1),preceding_velocity(i,2),preceding_velocity(i,3),i=1,number_of_particles)
CLOSE(unit=print_vtk)
print_number=print_number+1
END SUBROUTINE print_vtk_ascii
!-------------------------
SUBROUTINE print_csv
USE DECLARE
IMPLICIT NONE
INTEGER :: csv=0
CHARACTER (LEN=8) :: date
CHARACTER (LEN=10) :: timet
CHARACTER (LEN=5) :: zone
IF (time==0) THEN
    CALL DATE_AND_TIME (date,timet,zone)
    OPEN(newunit=csv, file='file'//date//' '//timet//zone//'.csv', action='write', status='replace', decimal='comma', & 
        & IOSTAT= ios,iomsg=errormsg)
    IF (ios/=0) THEN
        print *, "error opening .csv", csv
        STOP
    ENDIF
    WRITE(csv,*) preceding_velocity
ELSE IF (time==final_time) THEN
    WRITE(csv,*) preceding_velocity
    CLOSE(unit=csv)
    RETURN
ELSE
    WRITE(csv,*) preceding_velocity
ENDIF
END SUBROUTINE print_csv
END MODULE OUTPUT
!----------------------------------
SUBROUTINE read_vtk_ascii
USE DECLARE
USE OUTPUT
IMPLICIT NONE
CHARACTER (LEN=90) :: input, dummy
INTEGER :: read_vtk=0,i=1,dimmy=0
print *, "enter input filename without .vtk"
read *, input
OPEN(newunit=read_vtk,action='read',FILE=trim(input)//'.vtk',STATUS='old',encoding='utf-8',IOSTAT= ios,iomsg=errormsg)
IF (ios/=0) THEN
    print *, "error opening read_vtk", read_vtk, errormsg
    STOP
ENDIF
!# vtk DataFile Version 3.0
READ(read_vtk,*)
!system name max 256 char
READ(read_vtk,*) solver, print_format, system_name, gravity_is_on
!ASCII/BINARY
READ(read_vtk,*)
!DATASET type
READ(read_vtk,*)
number_of_particles = 0
!POINTS n dataType / for unstructured grid
READ(read_vtk,*) dummy, number_of_particles, dummy
!ALLOCATE
ALLOCATE(force(number_of_particles,3),preceding_force(number_of_particles,3),radius(number_of_particles), & 
    & preceding_position(number_of_particles,3),preceding_velocity(number_of_particles,3),position(number_of_particles,3), &
    & velocity(number_of_particles,3), STAT=ios, ERRMSG=errormsg)
IF(ios /= 0) THEN
     print *, "Error trying to ALLOCATE array", errormsg
     STOP 
END IF
radius=0.0_dp
elasticity_modulus=0.0_dp
poisson_coefficient=0.0_dp
density=0.0_dp
preceding_position=0.0_dp
preceding_velocity=0.0_dp
preceding_force=0.0_dp
force=0.0_dp
position=0.0_dp
velocity=0.0_dp
!position of points
READ(read_vtk,*) (preceding_position(i,1),preceding_position(i,2),preceding_position(i,3),i=1,number_of_particles)
!FIELD FieldData numberarrays
READ(read_vtk,*)
!Time 1 1 double
READ(read_vtk,*)
time=0.0_dp
!real(time,kind=8)
READ(read_vtk,*) time
!system_specifications 1 3 double
READ(read_vtk,*)
delta_time_to_print=0.0_dp
delta_time=0.0_dp
final_time=0.0_dp
!dtprint, dt, final, phi, contact_damping_rate, tolerance_velocity, tolerance_position
READ(read_vtk,*) delta_time,delta_time_to_print,final_time
!particle_specifications 1 3 double
READ(read_vtk,*)
!elasticity, poisson, density, particle number is implied
READ(read_vtk,*) elasticity_modulus,poisson_coefficient,density
number_of_walls = 0
!wall_normal number_of_walls 3 double
READ(read_vtk,*)dummy,number_of_walls,dimmy,dummy
ALLOCATE(wall_normal(1:number_of_walls,3),wall_point(1:number_of_walls,3), STAT=ios, ERRMSG=errormsg)
IF(ios /= 0) THEN
     print *, "Error trying to ALLOCATE array", errormsg
     STOP 
END IF
wall_normal=0.0_dp
wall_point=0.0_dp
!wall_normal&point
READ(read_vtk,*)(wall_normal(i,1),wall_normal(i,2),wall_normal(i,3),wall_point(i,1),wall_point(i,2),wall_point(i,3),i=1,number_of_walls)
!POINT_DATA number_of_particles
READ(read_vtk,*)
!SCALARS Particle_radius double
READ(read_vtk,*)
!LOOKUP_TABLE default
READ(read_vtk,*)
!real(radius(number_of_particles),kind=8)
READ(read_vtk,*) (radius(i),i=1,number_of_particles)
!VECTORS Particle_velocity double
READ(read_vtk,*)
!(real(velocity(i,:),kind=8),i=1,number_of_particles)
READ(read_vtk,*) (preceding_velocity(i,1),preceding_velocity(i,2),preceding_velocity(i,3),i=1,number_of_particles)
CLOSE(UNIT=read_vtk)
ALLOCATE(mass(number_of_particles), STAT=ios, ERRMSG=errormsg)
IF(ios /= 0) THEN
     print *, "Error trying to ALLOCATE array", errormsg
     STOP 
END IF
DO i=1, number_of_particles
    mass(i) = density*1.33333333333333_dp*pi*(radius(i)**3)
END DO
total_steps = final_time/delta_time
step_to_print = delta_time_to_print/delta_time
CALL dispatch_for_print
END SUBROUTINE read_vtk_ascii
!------------------------------
SUBROUTINE EULER_EXPLICIT
USE DECLARE
USE OUTPUT
IMPLICIT NONE
REAL(KIND=DP) :: overlap_wall=0.0_dp,overlap=0.0_dp,effective_elasticity_modulus=0.0_dp,effective_mass=0.0_dp, &
    & effective_radius=0.0_dp,norm_position=0.0_dp
REAL(KIND=DP), DIMENSION(3) :: fcont=0,normal_dir=0
INTEGER :: i=0, j=0, k=0
DO step=0, total_steps
    time=step*delta_time
    force=0.0_dp
    preceding_force=0.0_dp
    position=0.0_dp
    velocity=0.0_dp
    DO i=1, number_of_particles
        DO k=1, number_of_walls
            norm_position = ABS(DOT_PRODUCT((wall_point(k,:)-preceding_position(i,:)),wall_normal(k,:)))/NORM2(wall_normal)
            overlap_wall = radius(i) - norm_position
            IF (overlap_wall>0.0_dp) THEN
                normal_dir = wall_normal(k,:)/NORM2(wall_normal(k,:))
                effective_elasticity_modulus = elasticity_modulus/(1_dp-poisson_coefficient*poisson_coefficient)
                effective_radius = radius(i)
                effective_mass = mass(i)
                fcont = 1.33333333333333_dp*effective_elasticity_modulus*SQRT(effective_radius)*(overlap_wall**1.5_dp)*normal_dir
                force(i,:) = force(i,:) + fcont
            END IF
        END DO
        DO j=i+1, number_of_particles
            norm_position = ABS(NORM2((preceding_position(i,:)-preceding_position(j,:))))
            overlap = radius(i) + radius(j) - norm_position
            IF (overlap>0.0_dp) THEN
                normal_dir = (preceding_position(i,:)-preceding_position(j,:))/norm_position
                effective_elasticity_modulus = elasticity_modulus*elasticity_modulus/(elasticity_modulus+elasticity_modulus- &
                & poisson_coefficient*poisson_coefficient*elasticity_modulus-poisson_coefficient*poisson_coefficient*elasticity_modulus)
                effective_radius = radius(i)*radius(j)/(radius(i)+radius(j))
                effective_mass = mass(i)*mass(j)/(mass(i)+mass(j))
                fcont = 1.3333333333333333_dp*effective_elasticity_modulus*SQRT(effective_radius)*(overlap**1.5_dp)*normal_dir
                force(i,:) = force(i,:) + fcont
                force(j,:) = force(j,:) - fcont
            END IF
        END DO
        IF (gravity_is_on == .TRUE. .AND. radius(i).NE.preceding_position(i,2)) THEN
            force(i,:) = force(i,:) + mass(i)*gravity
        END IF
        velocity(i,:) = preceding_velocity(i,:) + 1.0_dp/mass(i)*force(i,:)*delta_time
        position(i,:) = preceding_position(i,:) + velocity(i,:)*delta_time
    END DO
    preceding_force = force
    preceding_position = position
    preceding_velocity = velocity
    IF(temporary_step==step_to_print) THEN
        temporary_step = 0
        CALL dispatch_for_print
    END IF
    temporary_step = temporary_step + 1
END DO
DEALLOCATE(force,preceding_force,radius,preceding_position,preceding_velocity,position,velocity, STAT=ios, ERRMSG=errormsg)
IF(ios /= 0) THEN
     WRITE(*,*)"Error trying to DEALLOCATE array", errormsg
     STOP 
END IF
END SUBROUTINE EULER_EXPLICIT
!--------------------------------------
SUBROUTINE EULER_IMPLICIT
USE DECLARE
USE OUTPUT
IMPLICIT NONE
REAL(KIND=DP) :: position_sn=0.0_dp,velocity_sn=0.0_dp,velocity_normalizer=0.0_dp,position_normalizer=0.0_dp, &
    & velocity_residual=0.0_dp,position_residual =0.0_dp,tolerance_error_velocity=1.0D-6,tolerance_error_position=1.0D-6, &
    & error_position=0.0_dp,error_velocity=0.0_dp,overlap_wall=0.0_dp,overlap=0.0_dp,temporary_time=0.0_dp, & 
    & effective_elasticity_modulus=0.0_dp,effective_mass=0.0_dp,effective_radius=0.0_dp,norm_position=0.0_dp
REAL(KIND=DP), DIMENSION(3) :: fcont=0,normal_dir=0
INTEGER :: i=0,j=0,k=0,l=0,max_iterations=6
ALLOCATE(position_pre_iter(number_of_particles,3),velocity_pre_iter(number_of_particles,3),force_pre_iter(number_of_particles,3), &
    & STAT=ios, ERRMSG=errormsg)
IF(ios /= 0) THEN
     print *, "Error trying to ALLOCATE array", errormsg
     STOP 
END IF
position_pre_iter = 0.0_dp
velocity_pre_iter = 0.0_dp
force_pre_iter = 0.0_dp
DO step=0, total_steps
    time=step*delta_time
    DO l=1, max_iterations 
        if (l==1) THEN
            position = preceding_position
            velocity = preceding_velocity
            position_pre_iter = preceding_position
            velocity_pre_iter = preceding_velocity
        END if
        force = 0.0_dp
        DO i=1, number_of_particles
            DO k=1, number_of_walls
                norm_position = ABS(DOT_PRODUCT((wall_point(k,:)-position(k,:)),wall_normal(k,:)))/NORM2(wall_normal)
                overlap_wall = radius(i) - norm_position
                IF (overlap_wall>0.0_dp) THEN
                    normal_dir = wall_normal(k,:)/NORM2(wall_normal(k,:))
                    effective_elasticity_modulus = elasticity_modulus/(1_dp-poisson_coefficient*poisson_coefficient)
                    effective_radius = radius(i)
                    effective_mass = mass(i)
                    fcont = 1.33333333333333_dp*effective_elasticity_modulus*SQRT(effective_radius)*(overlap_wall**1.5_dp)*normal_dir
                    force(i,:) = force(i,:) + fcont
                END IF
            END DO
            DO j=1, number_of_particles
                IF (i>=j) CYCLE
                norm_position = NORM2((position(i,:)-position(j,:)))
                overlap = radius(i) + radius(j) - norm_position
                IF (overlap>0.0_dp) THEN
                    normal_dir = (position(i,:)-position(j,:))/norm_position
                    effective_elasticity_modulus = elasticity_modulus*elasticity_modulus/(elasticity_modulus+elasticity_modulus-poisson_coefficient*poisson_coefficient*elasticity_modulus-poisson_coefficient*poisson_coefficient*elasticity_modulus)
                    effective_radius = radius(i)*radius(j)/(radius(i)+radius(j))
                    effective_mass = mass(i)*mass(j)/(mass(i)+mass(j))
                    fcont = -1.3333333333333333_dp*effective_elasticity_modulus*SQRT(effective_radius)*(overlap**1.5_dp)*normal_dir
                    force(i,:) = force(i,:) + fcont
                    force(j,:) = force(j,:) - fcont
                END IF
            END DO
            IF (gravity_is_on .EQV. .TRUE. .AND. radius(i).NE.preceding_position(i,2)) THEN
                force(i,:) = force(i,:) + mass(i)*gravity
            END IF
            velocity(i,:) = preceding_velocity(i,:) + 0.5_dp*delta_time*(preceding_force(i,:)+force(i,:))/mass(i)
            position(i,:) = preceding_position(i,:) + 0.5_dp*delta_time*(preceding_velocity(i,:)+velocity(i,:))
        END DO
        position_residual = sum(norm2((position-position_pre_iter),DIM=2))
        velocity_residual = sum(norm2((velocity-velocity_pre_iter),DIM=2))
        position_normalizer = sum(norm2((position-preceding_position),DIM=2))
        velocity_normalizer = sum(norm2((velocity-preceding_velocity),DIM=2))
        position_sn = sum(norm2(position,DIM=2))
        velocity_sn = sum(norm2(velocity,DIM=2))
        IF (position_normalizer < 1.0D-14*position_sn) THEN
            position_normalizer = position_sn
            tolerance_error_position = 0.01_dp*tolerance_error_position
            IF (position_sn < 1.0D-14) THEN
                position_normalizer = 1.0_dp
                tolerance_error_position = 1.0D-15
            END IF
        END IF
        IF (velocity_normalizer < 1.0D-14*velocity_sn) THEN
            velocity_normalizer = velocity_sn
            tolerance_error_velocity = 0.01_dp*tolerance_error_velocity
            IF (velocity_sn < 1.0D-14) THEN
                velocity_normalizer = 1.0_dp
                tolerance_error_velocity = 1.0D-15
            END IF
        END IF
        IF (position_normalizer<=TINY(1.0_dp)) position_normalizer = 1.0_dp
        IF (velocity_normalizer<=TINY(1.0_dp)) velocity_normalizer = 1.0_dp
        error_position = position_residual/position_normalizer
        error_velocity = velocity_residual/velocity_normalizer
        IF (error_velocity<tolerance_error_velocity .and. error_position<tolerance_error_position) EXIT
        position_pre_iter = position
        velocity_pre_iter = velocity
        force_pre_iter = force
    END DO
    preceding_force = force
    preceding_position = position
    preceding_velocity = velocity
    IF(temporary_step==step_to_print) THEN
        temporary_step = 0
        CALL dispatch_for_print
    END IF
    temporary_step = temporary_step + 1
END DO
DEALLOCATE(position_pre_iter,velocity_pre_iter,force_pre_iter,force,preceding_force,radius,preceding_position,preceding_velocity, &
    & position,velocity, STAT=ios, ERRMSG=errormsg)
IF(ios /= 0) THEN
     WRITE(*,*)"Error trying to DEALLOCATE array", errormsg
     STOP 
END IF
END SUBROUTINE EULER_IMPLICIT
!--------------------------
PROGRAM MAIN
USE DECLARE
IMPLICIT NONE
EXTERNAL :: EULER_EXPLICIT, EULER_IMPLICIT, READ_VTK_ASCII
CALL READ_VTK_ASCII
SELECT CASE (solver)
CASE ("explicit")
    CALL EULER_EXPLICIT
CASE ("implicit")
    CALL EULER_IMPLICIT
CASE DEFAULT
    PRINT *, "Command not recognized; please re-enter"
END SELECT
END PROGRAM MAIN
!--------------------------