module params
      implicit none

      !
      ! Define the shorthand for the common double.
      !
      integer,parameter :: rkd=selected_real_kind(10,40)

      !
      ! Define pi.
      !
      integer,parameter :: pi=4.D0*ATAN(1.D0)

      !
      ! Define the simulation box boundaries.
      ! This is an array of the form (width,height). The box is centered on (0,0).
      ! If a particle leaves the box, it is lost.
      !
      real(kind=rkd) :: sbox(2)
      !
      ! Define the basic parameters for the particle/individual interactions.
      !
      ! For the gravity we have the gravitational constant and the masses.
      !
      real(kind=rkd) :: gconst
      real(kind=rkd) :: gpars(6)

      !
      ! For the Langevian dynamics we have the gamma and the random force.
      !
      real(kind=rkd) :: lpars(6,2)

      ! 
      ! For the spring potential, or Hooke's Law, used in bonds we have
      ! the spring constant, equilibrium distance and cut off radius.
      ! They depend on the pairing of the particles.
      !
      real(kind=rkd) :: spars(6,6,3)

      !
      ! The Lennard-Jones potential has three parameters: epsilon, phi and
      ! the cut off radius. They are pairing dependent.
      !
      real(kind=rkd) :: ljpars(6,6,3)

      !
      ! The parameters for the particles. Define npp and npra for the total number of the
      ! population and the resource particles in the simulation,
      ! nnx, nny, nfx, nfy for the population particles,
      ! and nr, na for the resource and the alcohol particles.
      ! The actual particle data is stored in the arrays popg, resg and alcg,
      ! containing the population, resources and alcohol particles, respectively.
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      ! OLD INFORMATION (to be removed):                                           !
      ! The resource and the alcohol particles are in their own arrays as they are !
      ! more stable than the population particles and require less manipulation.   !
      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
      ! The vital rates, or actually probabilities, are stored in the array vitg
      ! which contains the particle type (1:4), age (0:101), and the probailities with
      ! modifiers (death rate, resource, alcohol, birth rate, resource, alcohol).
      !
      integer :: npp, npra, nnx, nny, nfx, nfy, nr, na
      real(kind=rkd),allocatable :: popg(:,:),resg(:,:),alcg(:,:),vitg(:,:,:)
      real(kind=rkd) :: agebprob=1.0D+00, idcount=0.0D+00

      !
      ! The size of the simulation area (x,y).
      !
      real(kind=rkd) :: simarea(2)

      !
      ! The arrays to set and store the random seed(s).
      !
      integer :: seed_size=1
      integer,allocatable :: seed_array(:)

      !
      ! Time integration parameters
      !
      real(kind=rkd) :: tmin,tmax,dt,tminus,quench_t,quench_s=-1.0D+00,vitcc
      real(kind=rkd),allocatable :: tavrg(:)
      integer :: outcc,tavrgcc
      contains
              !
              ! The function to form the particle identifier for the output
              ! This is here as the output depends on the definition of the parameters.
              !
              function form_ident(p)
                      real(kind=rkd) :: p(21)
                      character(len=4) :: form_ident

                      if (int(p(1))==3 .and. int(p(21))==0) then
                              form_ident="NxS0"
                      else if (int(p(1))==3 .and. int(p(21))/=0) then
                              form_ident="NxB0"
                      else if (int(p(1))==4 .and. int(p(21))==0) then
                              form_ident="NyS0"
                      else if (int(p(1))==4 .and. int(p(21))/=0) then
                              form_ident="NyB0"
                      else if (int(p(1))==5 .and. int(p(21))==0) then
                              form_ident="FxS0"
                      else if (int(p(1))==5 .and. int(p(21))/=0) then
                              form_ident="FxB0"
                      else if (int(p(1))==6 .and. int(p(21))==0) then
                              form_ident="FyS0"
                      else if (int(p(1))==6 .and. int(p(21))/=0) then
                              form_ident="FyB0"
                      else if (int(p(1))==1) then
                              form_ident="Res0"
                      else if (int(p(1))==2) then
                              form_ident="Alc0"
                      else
                              form_ident="0000"
                      end if
              end function form_ident
              !
              ! Initialize all the parameters with ones.
              !
              subroutine primitive_hard_init()
                      real(kind=rkd) :: tmp(2)
                      integer :: i,j,k


                      !
                      ! Set the basic individual particle parameters.
                      !
                      gconst=1.0D+00
                      gpars=1.0D+00
                      spars=1.0D+00
                      ljpars=1.5D+00

                      nnx=1!00
                      nny=1!00
                      nfx=0!1!00
                      nfy=0!1!00
                      nr=0!1!00
                      na=0!1!00
                      npp=nnx+nny+nfx+nfy
                      npra=nr+na

                      !
                      ! Set the size (x,y) of the simulation area proportional to the size
                      ! of the population. The sides of the area are at least 100x100 units. 
                      !
                      simarea(1:2)=max(100,npp)

                      !
                      ! For each individual/particle we define
                      ! in Stage 1:
                      ! type, pos_x, pos_y, vel_x, vel_y,
                      ! force_x, force_y, force_x_0, force_Y_0,
                      ! bonded_pair, pressure tensor contributions xx, yy, xy,
                      ! and kinetic energy.
                      ! Stage 2:
                      ! age
                      !
                      allocate(popg(npp,1:21))
                      allocate(resg(nr,1:11))
                      allocate(alcg(na,1:11))

                      !
                      ! Each individual particle is assigned a type.
                      ! We assign random postions for the individual particles.
                      ! The initial velocities and forces are all set to zero.
                      !
                      !
                      ! Population nx
                      !
                      j=1
                      i=1
                      if (nnx > 0) then
                              do i=j,nnx
                              popg(i,1)=3.0D+00
                              call random_number(tmp)
                              popg(i,2:3)=tmp*simarea
                              popg(i,4:10)=0.0D+00
                              end do
                      end if
                      j=i
                      
                      !
                      ! Population ny
                      !
                      if (nny > 0) then
                              do i=j,j+nny-1
                              popg(i,1)=4.0D+00
                              call random_number(tmp)
                              popg(i,2:3)=tmp*simarea
                              popg(i,4:10)=0.0D+00
                              end do
                      end if
                      j=i

                      !
                      ! Population fx
                      !
                      if (nfx > 0) then
                              do i=j,j+nfx-1
                              popg(i,1)=5.0D+00
                              call random_number(tmp)
                              popg(i,2:3)=tmp*simarea
                              popg(i,4:10)=0.0D+00
                              end do
                      end if
                      j=i

                      !
                      ! Population fy
                      !
                      if (nfy > 0) then
                              do i=j,j+nfy-1
                              popg(i,1)=6.0D+00
                              call random_number(tmp)
                              popg(i,2:3)=tmp*simarea
                              popg(i,4:10)=0.0D+00
                              end do
                      end if
                      popg(1,10)=2.0D+00
                      popg(2,10)=1.0D+00

                      !
                      ! Each resource and alcohol particle is assigned a type.
                      ! We assign random postions the particles.
                      ! The initial velocities and forces are all set to zero.
                      !
                      !
                      ! The resource particles
                      !
                      if (nr > 0) then
                              do i=1,nr
                              resg(i,1)=1.0D+00
                              call random_number(tmp)
                              resg(i,2:3)=tmp*simarea
                              resg(i,4:10)=0.0D+00
                              end do
                      end if

                      !
                      ! The alcohol particles
                      !
                      if (na > 0) then
                              do i=1,na
                              alcg(i,1)=2.0D+00
                              call random_number(tmp)
                              alcg(i,2:3)=tmp*simarea
                              alcg(i,4:10)=0.0D+00
                              end do
                      end if

                      !
                      ! Initialize the vital rates.
                      ! These include particle (1) type, (2) age specific annual probability of death,
                      ! age specific death rate modifiers by (3) general resources and (4) alcohol reources,
                      ! age specific (5) annual probability of reproduction, and the rate modifiers by
                      ! (6) general resources and (7) alcohol resources.
                      !
                      allocate(vitg(1:4,0:101,1:6))

                      do i=3,6
                      do j=1,100
                      do k=1,7
                      vitg(i,j,k)=0.01 ! NOTE: This is a dummy array, and the true values must be read from file.
                                       ! The statistics are from OSF. Death rates are based on the mean death rates
                                       ! in the years 2022 ans 2012.
                      end do
                      end do
                      end do
              end subroutine primitive_hard_init

              !
              ! Initialize the parameters with hardcoded values that
              ! depend on the particle pairing.
              !
              subroutine selective_hard_init()
                      integer :: i,j,k

                      !
                      ! Allocate the dummy population.
                      !
                      allocate(popg(npp,1:21))
                      !allocate(resg(nr,1:10))
                      !allocate(alcg(na,1:10))

                      !
                      ! Allocate the dummy vital rates.
                      !
                      allocate(vitg(1:4,0:101,1:6))

                      do i=3,6
                      do j=1,100
                      do k=1,7
                      vitg(i,j,k)=0.01 ! NOTE: This is a dummy array, and the true values must be read from file.
                                       ! The statistics are from OSF. Death rates are based on the mean death rates
                                       ! in the years 2022 ans 2012.
                      end do
                      end do
                      end do

                      !
                      ! Gravitational constant
                      !
                      gconst=10.0D+00
                      !
                      ! The mass, the whole array
                      !
                      gpars(1)=500.0D+00 ! rp, the resource
                      gpars(2)=200.0D+00 ! ap, the alcohol
                      gpars(3)=4.0D+00 ! nx
                      gpars(4)=4.0D+00 ! ny
                      gpars(5)=3.0D+00 ! fx
                      gpars(6)=3.0D+00 ! fy

                      !
                      ! The spring constant, lower triangle
                      !
                      spars(1,1,1)=1.0D+00 ! rp - rp
                      spars(1,2,1)=1.0D+00 ! rp - ap
                      spars(1,3,1)=1.0D+00 ! rp - nx
                      spars(1,4,1)=1.0D+00 ! rp - ny
                      spars(1,5,1)=1.0D+00 ! rp - fx
                      spars(1,6,1)=1.0D+00 ! rp - fy
                      spars(2,2,1)=1.0D+00 ! ap - ap
                      spars(2,3,1)=1.0D+00 ! ap - nx
                      spars(2,4,1)=1.0D+00 ! ap - ny
                      spars(2,5,1)=1.0D+00 ! ap - fx
                      spars(2,6,1)=1.0D+00 ! ap - fy
                      spars(3,3,1)=1.0D+00 ! nx - nx
                      spars(3,4,1)=1.0D+00 ! nx - ny
                      spars(3,5,1)=1.0D+00 ! nx - fx
                      spars(3,6,1)=1.0D+00 ! nx - fy
                      spars(4,4,1)=1.0D+00 ! ny - ny
                      spars(4,5,1)=1.0D+00 ! ny - fx
                      spars(4,6,1)=1.0D+00 ! ny - fy
                      spars(5,5,1)=1.0D+00 ! fx - fx
                      spars(5,6,1)=1.0D+00 ! fx - fy
                      spars(6,6,1)=1.0D+00 ! fy - fy
                      !
                      ! The equilibrium constant, lower triangle
                      !
                      spars(1,1,2)=2.0D+00 ! rp - rp
                      spars(1,2,2)=2.0D+00 ! rp - ap
                      spars(1,3,2)=2.0D+00 ! rp - nx
                      spars(1,4,2)=2.0D+00 ! rp - ny
                      spars(1,5,2)=2.0D+00 ! rp - fx
                      spars(1,6,2)=2.0D+00 ! rp - fy
                      spars(2,2,2)=2.0D+00 ! ap - ap
                      spars(2,3,2)=2.0D+00 ! ap - nx
                      spars(2,4,2)=2.0D+00 ! ap - ny
                      spars(2,5,2)=2.0D+00 ! ap - fx
                      spars(2,6,2)=2.0D+00 ! ap - fy
                      spars(3,3,2)=2.0D+00 ! nx - nx
                      spars(3,4,2)=2.0D+00 ! nx - ny
                      spars(3,5,2)=2.0D+00 ! nx - fx
                      spars(3,6,2)=2.0D+00 ! nx - fy
                      spars(4,4,2)=2.0D+00 ! ny - ny
                      spars(4,5,2)=2.0D+00 ! ny - fx
                      spars(4,6,2)=2.0D+00 ! ny - fy
                      spars(5,5,2)=2.0D+00 ! fx - fx
                      spars(5,6,2)=2.0D+00 ! fx - fy
                      spars(6,6,2)=2.0D+00 ! fy - fy
                      !
                      ! The cut off radius, lower triangle
                      !
                      spars(1,1,3)=5.0D+00 ! rp - rp
                      spars(1,2,3)=5.0D+00 ! rp - ap
                      spars(1,3,3)=5.0D+00 ! rp - nx
                      spars(1,4,3)=5.0D+00 ! rp - ny
                      spars(1,5,3)=5.0D+00 ! rp - fx
                      spars(1,6,3)=5.0D+00 ! rp - fy
                      spars(2,2,3)=5.0D+00 ! ap - ap
                      spars(2,3,3)=5.0D+00 ! ap - nx
                      spars(2,4,3)=5.0D+00 ! ap - ny
                      spars(2,5,3)=5.0D+00 ! ap - fx
                      spars(2,6,3)=5.0D+00 ! ap - fy
                      spars(3,3,3)=5.0D+00 ! nx - nx
                      spars(3,4,3)=5.0D+00 ! nx - ny
                      spars(3,5,3)=5.0D+00 ! nx - fx
                      spars(3,6,3)=5.0D+00 ! nx - fy
                      spars(4,4,3)=5.0D+00 ! ny - ny
                      spars(4,5,3)=5.0D+00 ! ny - fx
                      spars(4,6,3)=5.0D+00 ! ny - fy
                      spars(5,5,3)=5.0D+00 ! fx - fx
                      spars(5,6,3)=5.0D+00 ! fx - fy
                      spars(6,6,3)=5.0D+00 ! fy - fy
                      !
                      ! Fill out the upper triangle by mirroring.
                      !
                      do i=2,6
                      do j=1,i-1
                      do k=1,3
                      spars(i,j,k)=spars(j,i,k)
                      end do
                      end do
                      end do

                      !
                      ! Epsilon for the Lennard-Jones potential, lower triangle
                      !
                      ljpars(1,1,1)=4.0D+00 ! rp - rp
                      ljpars(1,2,1)=4.0D+00 ! rp - ap
                      ljpars(1,3,1)=3.0D+00 ! rp - nx
                      ljpars(1,4,1)=3.0D+00 ! rp - ny
                      ljpars(1,5,1)=1.0D+00 ! rp - fx
                      ljpars(1,6,1)=1.0D+00 ! rp - fy
                      ljpars(2,2,1)=4.0D+00 ! ap - ap
                      ljpars(2,3,1)=1.0D+00 ! ap - nx
                      ljpars(2,4,1)=1.0D+00 ! ap - ny
                      ljpars(2,5,1)=6.0D+00 ! ap - fx
                      ljpars(2,6,1)=6.0D+00 ! ap - fy
                      ljpars(3,3,1)=1.0D+00 ! nx - nx
                      ljpars(3,4,1)=4.0D+00 ! nx - ny
                      ljpars(3,5,1)=1.0D+00 ! nx - fx
                      ljpars(3,6,1)=5.0D+00 ! nx - fy
                      ljpars(4,4,1)=1.0D+00 ! ny - ny
                      ljpars(4,5,1)=5.0D+00 ! ny - fx
                      ljpars(4,6,1)=1.0D+00 ! ny - fy
                      ljpars(5,5,1)=1.0D+00 ! fx - fx
                      ljpars(5,6,1)=6.0D+00 ! fx - fy
                      ljpars(6,6,1)=1.0D+00 ! fy - fy
                      !
                      ! Phi for the Lennard-Jones potential, lower triangle
                      !
                      ljpars(1,1,2)=40.0D-01 ! rp - rp
                      ljpars(1,2,2)=40.0D-01 ! rp - ap
                      ljpars(1,3,2)=16.0D-01 ! rp - nx
                      ljpars(1,4,2)=16.0D-01 ! rp - ny
                      ljpars(1,5,2)=20.0D-01 ! rp - fx
                      ljpars(1,6,2)=20.0D-01 ! rp - fy
                      ljpars(2,2,2)=32.0D-01 ! ap - ap
                      ljpars(2,3,2)=20.0D-01 ! ap - nx
                      ljpars(2,4,2)=20.0D-01 ! ap - ny
                      ljpars(2,5,2)=8.0D-01 ! ap - fx
                      ljpars(2,6,2)=8.0D-01 ! ap - fy
                      ljpars(3,3,2)=16.0D-01 ! nx - nx
                      ljpars(3,4,2)=8.0D-01 ! nx - ny
                      ljpars(3,5,2)=32.0D-01 ! nx - fx
                      ljpars(3,6,2)=8.0D-01 ! nx - fy
                      ljpars(4,4,2)=16.0D-01 ! ny - ny
                      ljpars(4,5,2)=6.0D-01 ! ny - fx
                      ljpars(4,6,2)=16.0D-01 ! ny - fy
                      ljpars(5,5,2)=6.0D-01 ! fx - fx
                      ljpars(5,6,2)=4.0D-01 ! fx - fy
                      ljpars(6,6,2)=6.0D-01 ! fy - fy
                      !
                      ! The cut off radius, lower triangle
                      !
                      ljpars(1,1,3)=20.0D+00 ! rp - rp
                      ljpars(1,2,3)=20.0D+00 ! rp - ap
                      ljpars(1,3,3)=20.0D+00 ! rp - nx
                      ljpars(1,4,3)=20.0D+00 ! rp - ny
                      ljpars(1,5,3)=20.0D+00 ! rp - fx
                      ljpars(1,6,3)=20.0D+00 ! rp - fy
                      ljpars(2,2,3)=20.0D+00 ! ap - ap
                      ljpars(2,3,3)=20.0D+00 ! ap - nx
                      ljpars(2,4,3)=20.0D+00 ! ap - ny
                      ljpars(2,5,3)=20.0D+00 ! ap - fx
                      ljpars(2,6,3)=20.0D+00 ! ap - fy
                      ljpars(3,3,3)=20.0D+00 ! nx - nx
                      ljpars(3,4,3)=20.0D+00 ! nx - ny
                      ljpars(3,5,3)=20.0D+00 ! nx - fx
                      ljpars(3,6,3)=20.0D+00 ! nx - fy
                      ljpars(4,4,3)=20.0D+00 ! ny - ny
                      ljpars(4,5,3)=20.0D+00 ! ny - fx
                      ljpars(4,6,3)=20.0D+00 ! ny - fy
                      ljpars(5,5,3)=20.0D+00 ! fx - fx
                      ljpars(5,6,3)=20.0D+00 ! fx - fy
                      ljpars(6,6,3)=20.0D+00 ! fy - fy
                      !
                      ! Fill out the upper triangle by mirroring.
                      !
                      do i=2,6
                      do j=1,i-1
                      do k=1,3
                      ljpars(i,j,k)=ljpars(j,i,k)
                      end do
                      end do
                      end do
              end subroutine selective_hard_init

              !
              ! Initialize the parameters with hardcoded values.
              !
              subroutine hard_init(n)
                      integer :: n

                      !
                      ! Initialize the pseudo RNG.
                      !
                      allocate(seed_array(seed_size))
                      call random_seed()

                      !
                      ! Choose which initialization scheme we use.
                      !
                      if (n==1) then
                              call primitive_hard_init()
                      else
                              call selective_hard_init()
                      end if
              end subroutine hard_init

              subroutine cheat_pop(n,x,y)
                      !
                      ! Cheat in initial positioning for benchmarking
                      ! purposes. Set the first n particles in known
                      ! positions. If x or y is non zero, then adjust
                      ! the x and/or y coordinates of the particles.
                      ! 
                      integer :: n,x,y,i

                      if (n > 0 .AND. (x > 0 .OR. y > 0)) then
                        popg(1:n,2:3) = 50.0D+00
                        if (x > 0) then
                          do i=1,n
                          popg(i,2) = popg(i,2)+real(i*2.0D+00,kind=rkd)
                          end do
                        end if
                        if (y > 0) then
                          do i=1,n
                          popg(i,y) = popg(i,3)+real(i*2.0D+00,kind=rkd)
                          end do
                        end if
                      end if
              end subroutine cheat_pop

              subroutine read_time_integration_params()
                      !
                      ! Here we read from the file the parameters
                      ! for the time integration. These are
                      ! start time (FOR FUTURE USE, continue simulation)
                      ! time delta
                      ! the simulation end time
                      ! the quenching time
                      ! the time step multiplier during the quenching 
                      ! frequency of data output (every nth iteration)
                      ! the time average time window length
                      ! the frequency for probing deaths and births
                      !
                      character(len=50) :: paramin_name

                      open(20,file='paramint.txt', action='READ')
                      read(20,*) paramin_name, sbox(1:2)
                      read(20,*) paramin_name, tmin 
                      read(20,*) paramin_name, dt 
                      read(20,*) paramin_name, tmax
                      read(20,*) paramin_name, quench_t
                      read(20,*) paramin_name, quench_s
                      read(20,*) paramin_name, outcc
                      read(20,*) paramin_name, tavrgcc
                      read(20,*) paramin_name, vitcc
                      close(20)

                      allocate(tavrg(tavrgcc))
              end subroutine read_time_integration_params

              subroutine read_force_params()
                      !
                      ! Here we read from the file the parameters
                      ! for the force calculations.
                      !
                      character(len=50) :: paramin_name

                      open(20,file='paramfor.txt', action='READ')
                      read(20,*) paramin_name, gconst
                      read(20,*) paramin_name, gpars(1)
                      read(20,*) paramin_name, gpars(2)
                      read(20,*) paramin_name, gpars(3)
                      read(20,*) paramin_name, gpars(4)
                      read(20,*) paramin_name, gpars(5)
                      read(20,*) paramin_name, gpars(6)
                      read(20,*) paramin_name, lpars(1,1:2)
                      read(20,*) paramin_name, lpars(2,1:2)
                      read(20,*) paramin_name, lpars(3,1:2)
                      read(20,*) paramin_name, lpars(4,1:2)
                      read(20,*) paramin_name, lpars(5,1:2)
                      read(20,*) paramin_name, lpars(6,1:2)
                      read(20,*) paramin_name, spars(1,1,1)
                      read(20,*) paramin_name, spars(1,2,1)
                      read(20,*) paramin_name, spars(1,3,1)
                      read(20,*) paramin_name, spars(1,4,1)
                      read(20,*) paramin_name, spars(1,5,1)
                      read(20,*) paramin_name, spars(1,6,1)
                      read(20,*) paramin_name, spars(2,2,1)
                      read(20,*) paramin_name, spars(2,3,1)
                      read(20,*) paramin_name, spars(2,4,1)
                      read(20,*) paramin_name, spars(2,5,1)
                      read(20,*) paramin_name, spars(2,6,1)
                      read(20,*) paramin_name, spars(3,3,1)
                      read(20,*) paramin_name, spars(3,4,1)
                      read(20,*) paramin_name, spars(3,5,1)
                      read(20,*) paramin_name, spars(3,6,1)
                      read(20,*) paramin_name, spars(4,4,1)
                      read(20,*) paramin_name, spars(4,5,1)
                      read(20,*) paramin_name, spars(4,6,1)
                      read(20,*) paramin_name, spars(5,5,1)
                      read(20,*) paramin_name, spars(5,6,1)
                      read(20,*) paramin_name, spars(6,6,1)
                      read(20,*) paramin_name, spars(1,1,2)
                      read(20,*) paramin_name, spars(1,2,2)
                      read(20,*) paramin_name, spars(1,3,2)
                      read(20,*) paramin_name, spars(1,4,2)
                      read(20,*) paramin_name, spars(1,5,2)
                      read(20,*) paramin_name, spars(1,6,2)
                      read(20,*) paramin_name, spars(2,2,2)
                      read(20,*) paramin_name, spars(2,3,2)
                      read(20,*) paramin_name, spars(2,4,2)
                      read(20,*) paramin_name, spars(2,5,2)
                      read(20,*) paramin_name, spars(2,6,2)
                      read(20,*) paramin_name, spars(3,3,2)
                      read(20,*) paramin_name, spars(3,4,2)
                      read(20,*) paramin_name, spars(3,5,2)
                      read(20,*) paramin_name, spars(3,6,2)
                      read(20,*) paramin_name, spars(4,4,2)
                      read(20,*) paramin_name, spars(4,5,2)
                      read(20,*) paramin_name, spars(4,6,2)
                      read(20,*) paramin_name, spars(5,5,2)
                      read(20,*) paramin_name, spars(5,6,2)
                      read(20,*) paramin_name, spars(6,6,2)
                      read(20,*) paramin_name, spars(1,1,3)
                      read(20,*) paramin_name, spars(1,2,3)
                      read(20,*) paramin_name, spars(1,3,3)
                      read(20,*) paramin_name, spars(1,4,3)
                      read(20,*) paramin_name, spars(1,5,3)
                      read(20,*) paramin_name, spars(1,6,3)
                      read(20,*) paramin_name, spars(2,2,3)
                      read(20,*) paramin_name, spars(2,3,3)
                      read(20,*) paramin_name, spars(2,4,3)
                      read(20,*) paramin_name, spars(2,5,3)
                      read(20,*) paramin_name, spars(2,6,3)
                      read(20,*) paramin_name, spars(3,3,3)
                      read(20,*) paramin_name, spars(3,4,3)
                      read(20,*) paramin_name, spars(3,5,3)
                      read(20,*) paramin_name, spars(3,6,3)
                      read(20,*) paramin_name, spars(4,4,3)
                      read(20,*) paramin_name, spars(4,5,3)
                      read(20,*) paramin_name, spars(4,6,3)
                      read(20,*) paramin_name, spars(5,5,3)
                      read(20,*) paramin_name, spars(5,6,3)
                      read(20,*) paramin_name, spars(6,6,3)
                      read(20,*) paramin_name, ljpars(1,1,1)
                      read(20,*) paramin_name, ljpars(1,2,1)
                      read(20,*) paramin_name, ljpars(1,3,1)
                      read(20,*) paramin_name, ljpars(1,4,1)
                      read(20,*) paramin_name, ljpars(1,5,1)
                      read(20,*) paramin_name, ljpars(1,6,1)
                      read(20,*) paramin_name, ljpars(2,2,1)
                      read(20,*) paramin_name, ljpars(2,3,1)
                      read(20,*) paramin_name, ljpars(2,4,1)
                      read(20,*) paramin_name, ljpars(2,5,1)
                      read(20,*) paramin_name, ljpars(2,6,1)
                      read(20,*) paramin_name, ljpars(3,3,1)
                      read(20,*) paramin_name, ljpars(3,4,1)
                      read(20,*) paramin_name, ljpars(3,5,1)
                      read(20,*) paramin_name, ljpars(3,6,1)
                      read(20,*) paramin_name, ljpars(4,4,1)
                      read(20,*) paramin_name, ljpars(4,5,1)
                      read(20,*) paramin_name, ljpars(4,6,1)
                      read(20,*) paramin_name, ljpars(5,5,1)
                      read(20,*) paramin_name, ljpars(5,6,1)
                      read(20,*) paramin_name, ljpars(6,6,1)
                      read(20,*) paramin_name, ljpars(1,1,2)
                      read(20,*) paramin_name, ljpars(1,2,2)
                      read(20,*) paramin_name, ljpars(1,3,2)
                      read(20,*) paramin_name, ljpars(1,4,2)
                      read(20,*) paramin_name, ljpars(1,5,2)
                      read(20,*) paramin_name, ljpars(1,6,2)
                      read(20,*) paramin_name, ljpars(2,2,2)
                      read(20,*) paramin_name, ljpars(2,3,2)
                      read(20,*) paramin_name, ljpars(2,4,2)
                      read(20,*) paramin_name, ljpars(2,5,2)
                      read(20,*) paramin_name, ljpars(2,6,2)
                      read(20,*) paramin_name, ljpars(3,3,2)
                      read(20,*) paramin_name, ljpars(3,4,2)
                      read(20,*) paramin_name, ljpars(3,5,2)
                      read(20,*) paramin_name, ljpars(3,6,2)
                      read(20,*) paramin_name, ljpars(4,4,2)
                      read(20,*) paramin_name, ljpars(4,5,2)
                      read(20,*) paramin_name, ljpars(4,6,2)
                      read(20,*) paramin_name, ljpars(5,5,2)
                      read(20,*) paramin_name, ljpars(5,6,2)
                      read(20,*) paramin_name, ljpars(6,6,2)
                      read(20,*) paramin_name, ljpars(1,1,3)
                      read(20,*) paramin_name, ljpars(1,2,3)
                      read(20,*) paramin_name, ljpars(1,3,3)
                      read(20,*) paramin_name, ljpars(1,4,3)
                      read(20,*) paramin_name, ljpars(1,5,3)
                      read(20,*) paramin_name, ljpars(1,6,3)
                      read(20,*) paramin_name, ljpars(2,2,3)
                      read(20,*) paramin_name, ljpars(2,3,3)
                      read(20,*) paramin_name, ljpars(2,4,3)
                      read(20,*) paramin_name, ljpars(2,5,3)
                      read(20,*) paramin_name, ljpars(2,6,3)
                      read(20,*) paramin_name, ljpars(3,3,3)
                      read(20,*) paramin_name, ljpars(3,4,3)
                      read(20,*) paramin_name, ljpars(3,5,3)
                      read(20,*) paramin_name, ljpars(3,6,3)
                      read(20,*) paramin_name, ljpars(4,4,3)
                      read(20,*) paramin_name, ljpars(4,5,3)
                      read(20,*) paramin_name, ljpars(4,6,3)
                      read(20,*) paramin_name, ljpars(5,5,3)
                      read(20,*) paramin_name, ljpars(5,6,3)
                      read(20,*) paramin_name, ljpars(6,6,3)
                      close(20)
              end subroutine read_force_params

              subroutine read_popg()
                      integer :: n,i
                      real(kind=rkd) :: rnd1
                      !
                      ! Open the file and read in the number of
                      ! individual particles. Delete the old popg array
                      ! and allocate space for the new one.
                      !
                      open(20,file='popg.txt', action='READ')
                      read(20,*) npp
                      deallocate(popg)
                      allocate(popg(1:npp,1:21))
                      !
                      ! Read the particles.
                      !
                      ! The data consist:
                      ! 1 Type
                      ! 2,3 Position X,Y
                      ! 4,5 Velocity X,Y
                      ! 6,7 Force X,Y
                      ! 8,9 Force history X,Y
                      ! 10 ID
                      ! 11,12,13 Virial components
                      ! 14 Kinetic energy
                      ! 15 Age
                      ! 16 cumulative minimum distance to Res
                      ! 17 cumulative minimum distance to Alc
                      ! 18 cumulative minimum distance to Ny
                      ! 19 cumulative minimum distance to Fy
                      ! 20 age at giving birth
                      ! 21 Bond
                      do i=1,npp
                      read(20,*) popg(i,1:21)
                      end do
                      !
                      ! Close the file.
                      !
                      close(20)

                      ! Retrieve the current running particle ID.
                      idcount=maxval(popg(:,10))

                      do i=1,npp
                      ! If IDs are not found, assign ID to individuals.
                      if (popg(i,1)>2.0D+00 .and. popg(i,1)<7.0D+00 .and. popg(i,10)<1.0D+00) then
                              idcount=idcount+1.0D+00
                              popg(i,10)=idcount
                      end if
                      ! Additionally, add a random amount to the age to avoid catastrophic
                      ! simultaneous bond breaking when ages reach 18.
                      if (popg(i,1)>2.0 .and. popg(i,1)<7.0) then
                              call random_number(rnd1)
                              popg(i,15)=popg(i,15)+rnd1
                      end if
                      ! If the average minimum distances are 0, they are not true values.
                      ! Mark them as initializing values.
                      if (popg(i,16)<=0.0D+00) then
                              popg(i,16)=-1000.0D+00
                      end if
                      if (popg(i,17)<=0.0D+00) then
                              popg(i,17)=-1000.0D+00
                      end if
                      if (popg(i,18)<=0.0D+00) then
                              popg(i,18)=-1000.0D+00
                      end if
                      if (popg(i,19)<=0.0D+00) then
                              popg(i,19)=-1000.0D+00
                      end if
                      end do
              end subroutine read_popg

              subroutine read_vitg()
                      integer :: i,j
                      real :: r1,r2
                      !
                      ! Open the file and read the vital rate parameters.
                      !
                      open(20,file='vitg.txt', action='READ')
                      deallocate(vitg)
                      allocate(vitg(3:6,0:101,1:6))
                      !
                      ! Read the probabilities and the modifiers.
                      !
                      do i=3,6
                      do j=0,101
                      read(20,*) r1, r2, vitg(i,j,1:6)
                      end do
                      end do
                      !
                      ! Close the file.
                      !
                      close(20)
                      ! Scale the probabilities of reproduction.
                      vitg(:,:,4)=vitg(:,:,4)/maxval(vitg(:,:,4))
              end subroutine read_vitg
end module params
