      PROGRAM Sim_pop
         IMPLICIT NONE
         integer, parameter :: nvm = 70  ! maximum number of connections for each (a,b)
         INTEGER :: day,a,b, k, aa, bb, n, nvis, nvmu
         REAL*8 :: ar, br
         INTEGER, DIMENSION(:,:), ALLOCATABLE :: nvir
         INTEGER, DIMENSION(:,:,:), allocatable :: r1,r2
         LOGICAL, DIMENSION(:,:), ALLOCATABLE :: pop,  pop_new
         INTEGER(kind=2), DIMENSION(:,:), ALLOCATABLE :: R0, immuni, degenza
         LOGICAL, DIMENSION(:,:), ALLOCATABLE :: guarito
         INTEGER :: days, t0, startc
         INTEGER, DIMENSION(1:2) :: r
         REAL*8, DIMENSION(1:2) :: p
         REAL*8:: cont_rate, r0_print
         REAL*8 :: c, raggio, angolo, std
         INTEGER :: uscite, usc, quarantena, immuno_day
         INTEGER :: contagiati, infetti
         INTEGER, DIMENSION(1:2) :: z
         REAL*8, DIMENSION(1:2) :: rr
         INTEGER, ALLOCATABLE, DIMENSION(:) :: x,y
         INTEGER :: xmin, xmax, ymin, ymax, day_free, nv
         LOGICAL :: separa,found,added
         INTEGER(kind=1), ALLOCATABLE, DIMENSION(:,:) :: reg0
         INTEGER(kind=1), DIMENSION(8) :: error
         REAL :: tic, toc
         !
         CALL CPU_TIME(tic)
         n = 150
         z = [4,4]
         day_free = 100
         days = 60
         immuno_day = 15
         cont_rate = 0.01d0
         usc = 15
         std = DBLE(n)*dsqrt(2.d0)/2.d0
         t0 = 4
         quarantena = 15
         CALL random_seed(PUT = [23,42])
!
         ALLOCATE(pop(1:n, 1:n), STAT = error(1))
         ALLOCATE(R0(1:n,1:n), STAT = error(2))
         ALLOCATE(pop_new(1:n, 1:n), STAT = error(3))
         ALLOCATE(degenza(1:n, 1:n), STAT = error(4))
         ALLOCATE(guarito(1:n, 1:n), STAT = error(5))
         ALLOCATE(immuni(1:n, 1:n), STAT = error(6)) 
         ALLOCATE(reg0(1:t0,1:2), STAT = error(7))     
         allocate(nvir(1:n,1:n), r1(1:nvm,1:n,1:n), r2(1:nvm,1:n,1:n),stat=error(8))
         IF (sum(error).eq.0) THEN
             WRITE(*,*) 'No Error in allocation'
         ELSE
             WRITE(*,*) 'Allocation Errors for these items:'
             DO k=1,8
                 if(error(k) /= 0)WRITE(*,'(4x,2I4)',advance='no') k,error(k)
             END DO
             print *
             STOP
         END IF
         startc = 0
         pop(:,:) = .FALSE.
         R0(:,:) = 0
         pop_new(:,:) = .FALSE.
         degenza(:,:) = -1
         immuni(:,:) = 0
         guarito(:,:) = .FALSE.
         separa = .FALSE.
         nvir = 0
!
         OPEN(1000, FILE = 'res.txt')        
!
         DO k=1,2
             IF (Z(k).gt.1) THEN
                separa =.TRUE.
                IF (k.eq.1) THEN
                    ALLOCATE(x(1:z(k)))
                    CALL linspace(DBLE(1),DBLE(n),z(k),x)
                ELSE
                    ALLOCATE(y(1:z(k)))
                    CALL linspace(DBLE(1),DBLE(n),z(k),y)
                END IF
             END IF
         END DO
!
         DO WHILE (startc.lt.t0)
            CALL random_number(p)
            ar = (1.d0 + p(1) * (DBLE(n)))
            br = (1.d0 + p(2) * (DBLE(n)))
            a = floor(ar)
            b = floor(br)
            IF (.not. pop(a,b)) THEN
               found = .false.
			   IF (separa) THEN
				   DO aa=1,z(1)
					   if (a.eq.x(aa)) THEN
						   found = .true.
						   EXIT
					   END IF
				   END DO
				   IF (.not. found) THEN
					DO bb=1,z(2)
						if (b.eq.y(bb)) THEN
						  found = .true.
						  EXIT
						END IF
					END DO
				   END IF
			   END IF
               IF (.not.found) THEN
                pop(a,b) = .TRUE.
                degenza(a,b) = -1
                startc = startc + 1
                reg0(startc,:) = [a,b]
               END IF
            END IF
         END DO
         pop_new = pop
         print *,'At the beginning, the indices (a, b) where POP(a,b) = .T. are:'
         call prt_pop(pop,n)
         print *
         CALL scrivi_res(0,pop,n)
         DO day = 1,days
            infetti = 0
            DO a=1,n
                DO b=1,n
                    IF (pop(a,b)) THEN
                        infetti = infetti + 1
                    END IF
                END DO
            END DO
            contagiati = SUM(SUM(R0, DIM =1), DIM = 1)
            IF (infetti.ne.0) THEN
                r0_print = DBLE(contagiati)/DBLE(infetti)
            ELSE
                r0_print = 0.d0
            END IF
            WRITE(*,9000)  day, infetti, contagiati, r0_print
            R0(:,:) = 0
            nvir = 0
            r1 = 0
            r2 = 0
            DO a=1,n
               DO b=1,n
                  uscite = 0
                  IF ((separa).and.(day.le.day_free)) THEN
                     CALL confina(a,z(1),x,xmin,xmax)   
                     CALL confina(b,z(2),y,ymin,ymax)  
                  ELSE
                    xmin = 0
                    xmax = n+1
                    ymin = 0
                    ymax = n+1
                  END IF
                  DO WHILE (uscite.le.usc)
                     CALL random_norm(raggio)
                     raggio = dabs(raggio) * std * 0.38532d0
                     CALL random_number(c)
                     angolo = c * 2*3.14159265358979d0
                     rr(1) = (DBLE(a) + raggio*cos(angolo))
                     rr(2) = (DBLE(b) + raggio*sin(angolo))
                     r(1)  = idnint(rr(1))
                     r(2)  = idnint(rr(2))
                     IF (a.eq.r(1).and.b.eq.r(2)) THEN
                         CYCLE
                     END IF
                     ! Nessuno può uscire dalla sua "cella" dentro pop
                     IF ((r(1).gt.xmin.and.r(1).lt.xmax).and.(r(2).gt.ymin.and.r(2).lt.ymax)) THEN
						   CALL search(nvm, n, nvir, r1, r2, a, b, r, found)
                           IF (.not. found) THEN
                               uscite = uscite + 1
                               IF (pop(a,b).eqv.pop(r(1),r(2))) THEN ! entrambi infettati o no
                                CYCLE
                               END IF
                               CALL random_number(c)
                               IF (c.le.cont_rate) THEN
                                  ! Calcolo R0 e propago virur								   								 
                                     IF ((.not. pop(a,b)).and.(immuni(a,b).eq.0)) THEN
									        pop_new(a,b) = .TRUE.
                                            R0(r(1),r(2)) = R0(r(1),r(2)) + 1
                                     ELSEIF ((.not. pop(r(1), r(2))).and.(immuni(r(1), r(2)).eq.0)) THEN
											pop_new(r(1),r(2)) = .TRUE.
                                            R0(a,b) = R0(a,b) + 1
                                     END IF
                               END IF
                           END IF
                       END IF
                  END DO
               END DO
            END DO
            pop = pop_new
            CALL scrivi_res(day,pop,n)
            ! verifico le guarigioni
            DO aa=1,n
               DO bb=1,n
                  IF (pop(aa,bb)) THEN
                     degenza(aa,bb) = degenza(aa,bb) + 1
                  END IF
                  IF (guarito(aa,bb)) THEN
                    immuni(aa,bb) = immuni(aa,bb) -1
                  END IF
                  IF (immuni(aa,bb).le.0) THEN
                    guarito(aa,bb) = .FALSE.
                    immuni(aa,bb) = 0
                  END IF
                  IF ((pop(aa,bb)).and.(degenza(aa,bb).ge.quarantena)) THEN
                     degenza(aa,bb) = -1
                     pop(aa,bb) = .FALSE.
                     guarito(aa,bb) = .TRUE.
                     immuni(aa,bb) = immuno_day
                  END IF
               END DO
            END DO
            pop_new = pop
            print '(1x,A,2x,i3,2x,A)','At end of day ',day, &
               ' the indices (a, b) where POP(a,b) = .T. are:'
            call PRT_POP(pop,n)
            print *
        END DO
         print '(/,1x,A,I4,A,I4)', 'End of simulation. Infected = ', &
            infetti, ', Contagious = ',contagiati
         print '(1x,A,2x,I3,A)',' After ',days,' days'
         print '(/,A,2x,i4,/)','NVM_max used: ',nvmu
         CALL CPU_TIME(toc)
         print '("Time = ",f12.3," minutes.")',(toc-tic)/60.d0
         CLOSE(10)
!
9000 FORMAT(' Beginning day = ',I4,' Infected = ',I6, &
    ' Contagious = ',I6,' R0 = ', F5.3,/)
      STOP
      END PROGRAM
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Normally distributed random number generation
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! 
      SUBROUTINE random_norm(c)
         IMPLICIT NONE
         ReAL*8, INTENT(OUT) :: c
         REAL*8, DIMENSION(1:2) :: u
         REAL*8 :: pi = 3.14159265358979d0
!
         CALL random_number(u)
         C = dsqrt(-2.d0*dlog(u(1)))*dcos(2.d0*pi*u(2))
!
         RETURN
    END SUBROUTINE
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Linpsace
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!    
      SUBROUTINE linspace(a,b,n,yi)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: n
      REAL*8, INTENT(IN) :: a,b
      REAL*8, DIMENSION(1:n) :: y
      INTEGER, INTENT(OUT), DIMENSION(1:n) :: yi
      INTEGER :: j
      REAL*8 :: p
!
      p = (b-a)/DBLE(n-1)
      DO j=1,n
        y(j) = a + p*DBLE(j-1)
        yi(j) = idnint(y(j))
      END DO
!
      RETURN
      END SUBROUTINE  
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
! Confina
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!       
      SUBROUTINE confina(a,n,x,xmin,xmax)
      IMPLICIT NONE
      INTEGER, INTENT(IN) :: a, n
      INTEGER, DIMENSION(1:n), INTENT(IN) :: x
      INTEGER, INTENT(OUT) :: xmin, xmax
      INTEGER :: jl,jr,j
!CNXS Binary search instead of linear search, for sorted X(:)      
      jl=0; jr=n
      do while(jr-jl > 1)
         j = (jl+jr)/2
         if (x(j) >= a) then
            jr = j
         else
            jl = j
         endif
      end do
      if(jl==0)jl=jl+1
      xmin = x(jl); xmax=x(jl+1)
      return
!CNXS original linear search follows, not used now      
      DO j=1,(n-1)
        IF(a.ge.x(j).and.a.le.x(j+1)) THEN
            xmin = x(j)
            xmax = x(j+1)
            EXIT
        END IF
      END DO
      if(j /= jl)then
         write(*,'(i4,4x,2I4,4x,2i4,4x,2i4)')a,jl,j,x(jl),x(jl+1),xmin,xmax
      endif
!
    RETURN
    END SUBROUTINE

	subroutine search(nvm, n, nvir, r1, r2, a, b, r, found)
	implicit none
	integer, intent(in)  :: nvm,n,a,b
	integer, dimension(1:2), intent(in) :: r
	logical, intent(out) :: found
	integer, dimension(1:nvm,1:n,1:n), intent(inout) :: r1,r2
	integer, dimension(1:n,1:n), intent(inout) :: nvir
	integer :: i,m
	!
	found=.false.
	m = nvir(a,b)
    ! controllo se non sono troppo vicino alla soglia
    IF ((m+1).eq.nvm) THEN
        WRITE(*,*) 'vnm too small: enlarge it'
        READ(*,*)
        STOP
    END IF
    IF (m.ne.0) THEN       
	    do i=1,m
	       if (r1(a,b,i).eq.r(1).and.r2(a,b,i).eq.r(2)) then
		      found = .true.
		      return
	       end if
        end do
    END IF
	m = m+1
	r1(a,b,m) = r(1)
	r2(a,b,m) = r(2)
	!...è reciproco
	m = nvir(r(1),r(2))
	m = m+1
	r1(r(1),r(2),m) = a
	r2(r(1),r(2),m) = b
	return
	end subroutine

subroutine prt_pop(pop,n)
integer, intent(in) :: n
logical, intent(in) :: pop(n,n)
integer k, a, b
   k = 0
   do b = 1,n
      do a = 1,n
         if(pop(a,b))then
            k = k+1
            if(mod(k,5) == 0)then
               write(*,40)a,b
            else
               write(*,40, advance='no')a,b
            endif
         endif
   40 format(4x,'(',i3,',',i3,')')                  
      end do
   end do
   return
    end subroutine prt_pop
    
    subroutine scrivi_res(nfile,pop,n)
    logical, intent(in), dimension(1:n, 1:n) :: pop
    integer, intent(in) :: n, nfile
    integer:: j
    integer(kind=2), dimension(1:n) :: vec
    character(len=4) :: str
    !
    write(str,'(I4.4)') nfile
    OPEN(nfile, FILE ='res_'//str//'.txt')
    Do j=1,n
        where (pop(j,:))
            vec = 1
        elsewhere
            vec = 0
        end where
        WRITE(1000,1000) vec
        WRITE(nfile, 1000) vec
    END do
    CLOSE(nfile)
1000 FORMAT(1000(1X,I2))
    return
    end subroutine
            

