program Infest
  implicit none
  !--- Model parameters:
  integer, parameter :: itmax = 1000 !Days(Weeks) to simulate.
  integer, parameter :: Ip = 500 !Number of population.
  integer, parameter :: Ic =  10 !Number of affected people at t=0.
  integer, parameter :: Im =  60 !Medical capacity.
  integer, parameter :: Nds =  3 !Days(Weeks) to start the disease.
  integer, parameter :: Ndr =  6 !Days(Weeks) to get recovered.
  double precision, parameter :: pmax = 0.80d0 !Probability of new infection in meeting of "affect-normal" pair.
  double precision, parameter :: fd = 0.6d0 !1.d0 !"Social-distance" factor: if fd = 0.d0, people are never affected, except the initial (t=0) affection.
  double precision, parameter :: pdie = 0.25d0 !Death probability without medical cares.
  double precision, parameter :: pcol = 0.12d0 !Death probability, although under medical cares.
  !--- Local variables, changable in the process:
  integer :: it, k, l, nn1, nn2
  integer :: n0, n1, nd, nr !Number of "normal", "affected", "dead", and "recovered".
  integer, dimension(:), allocatable :: nn !Labels of people.
  integer, dimension(:), allocatable :: itnn !Time counter from when affected.
  double precision :: p = 0.d0
  write(6,*) "Number of population =", Ip, " thousands"
  write(6,*) "Number of affected at t=0 =", Ic, " thousands"
  write(6,*) "       Time :     Normal,   Affected,     Dead,   Recovered."
  10001 allocate(nn(Ip), itnn(Ip))
  nn = 0; itnn = 0
  do l = 1, Ic
     nn(l) = 1
  end do
  !--- Start the calculation.
  do 16911 it = 0, itmax
     if((n0+nd+nr).eq.Ip) exit
     if(it.eq.0) go to 71530
     !--- New affection.
     do l = 1, Ip
        nn1 = nn(l)
        do k = l+1, Ip
           nn2 = nn(k)
           !--- For "affect-affect" pair, nothing happens.
           if(nn1.eq.1 .and. nn2.eq.nn1) exit
           !--- For "normal-normal" pair, nothing happens.
           if(nn1.eq.0 .and. nn2.eq.nn1) exit
           !--- If normal meets affected,
           if((min(nn1,nn2).eq.0) .and. (max(nn1,nn2).eq.1)) then
              call random_at_date_and_time(1, p)
              if (p.lt.pmax*fd) then
                 nn(l) = 1; nn(k) = 1 !Newly affected person.
              end if
              exit
           end if
        end do
     end do
     !--- Record the time of each person from when affected.
     do l = 1, Ip
        if(nn(l).eq.1) itnn(l) = itnn(l)+1
     end do
     !--- Die or survive, or recover?
     k = 0
     do l = 1, Ip
        if(nn(l).eq.1) k = k+1 !--- Count the number of affected people to compare with the medical capacity.
        if     (itnn(l).ge.Ndr) then
           nn(l) = -1 !------------------This one gets recovered!
        elseif (itnn(l).ge.Nds) then
           call random_at_date_and_time(1, p)
           if(k .le. Im) then
              if(p.lt.pcol) nn(l) = -9 !----This one dies in hospital...
           else
              if(p.lt.pdie) nn(l) = -9 !----This one dies without cares...
           end if
        end if
     end do
     !--- Result at this time.
     71530 n0=0; n1=0; nd=0; nr=0
     do l = 1, Ip
        if(nn(l).eq. 0) n0 = n0+1
        if(nn(l).eq. 1) n1 = n1+1
        if(nn(l).eq.-9) nd = nd+1
        if(nn(l).eq.-1) nr = nr+1
     end do
     write(6,*) it,":", n0, n1, nd, nr
  16911 end do
  !--- End of the calculation.
  20001 deallocate(nn, itnn)
  write(6,*) "        Day :     Normal,   Affected,     Dead,   Recovered."
  write(6,*) ""
  write(6,*) " Number of population =", Ip
  write(6,*) " Number of deaths =", nd
  write(6,*) ""
  write(6,*) " Social-distance factor =", fd, ":"
  write(6,*) "   when this factor is ZERO, no new infections happen,"
  write(6,*) "   except the initial (t=0) cases."
  99999 continue
  write(6,*) ""
  write(6,*) "This program terminated safely. (^_^)b"
  write(6,*) ""
end program Infest

!--- Subroutine to get the random double-precision x(n) =< 1.0d0.
   subroutine random_at_date_and_time(n,x)
   implicit none
   integer, intent(IN) :: n
   double precision, intent(OUT) :: x(n)
   integer :: nsize, md(8)
   integer, allocatable :: iseed(:)
   double precision :: rnd(n)
   call random_seed(size=nsize) ; allocate(iseed(nsize))
   call random_seed(get=iseed) !; write(6,*) iseed
   call date_and_time(values = md) !only md(8) is used here.
   iseed = iseed*md(8)
   call random_seed(put=iseed)
   call random_number(rnd)
   deallocate(iseed)
   x = rnd
   end subroutine random_at_date_and_time