! Author: J. Pérez-Losada
! e-mail: joaquim.perez@udg.edu
! This program simulates the Neolithic dispersal along the SEA route (blue lines in Figs. 1b and 3)
!    
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Output files:
!  Archaelogical_Time_Sea_Model.dat (times of first arrival of farmers at the given archaelogical sites)
!  Fraction_Sea_Model.dat           (fraction of farmers with Haplo-K at given time and distance)
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
program main
!USE MSFLIB                                          ! Library 

implicit none
integer:: i,j,x,y,z,random_t                        ! Service variables
integer:: row_sea,col_sea                           ! row #, column # of the starting point of the dispersal
real*4 aleatori, decimal, jum_suma                  ! Service random variables
real*8 :: r_p                                       ! Control variable dispersal
character(len=35) :: sitename                       ! Name of archaelogical site                     
integer :: ARQ_FRA,ARQ_ITA,ARQ_POR                  ! Node # of Arqueological sites 
integer :: ARQ_NES,ARQ_NAV,ARQ_CYP     
integer :: CRT_05,CRT_06,CRT_07                     ! Control stwiches to find time of arrival to a site 1D
integer :: CRT_08,CRT_09,CRT_13  

real*8,dimension (-1:301,-1:301):: rr_hg            ! Net Growth Population factor Hunter-Gatherers per generation
real*8,dimension (-1:301,-1:301):: a_p              ! Matrix Farmers without Haplo-K Population time t
real*8,dimension (-1:301,-1:301):: a_n              ! Matrix Farmers with Haplo-K Population time t
real*8,dimension (-1:301,-1:301):: a_hg             ! Matrix Hunter-Gatherers Population
real*8,dimension (-1:301,-1:301):: b_p              ! Matrix Farmers Population time t+1
real*8,dimension (-1:301,-1:301):: b_n              ! Matrix Farmers Population time t+1
real*8,dimension (-1:301,-1:301):: eta              ! Matrix mixing factor between Farmers and Hunter-Gatherers 
real*8 :: couples_hn                                ! Number of hybrid couples farmers without Haplo-K
real*8 :: couples_hp                                ! Number of hybrid couples farmers with Haplo-K
real*8 :: farmers                                   ! Total number of farmers
real*8 :: fraction                                  ! Fraction of farmers with HK at starting point
real*8 :: pe_p                                      ! Persitency Farmers without Haplo-K 
real*8 :: pe_n                                      ! Persitency Farmers with Haplo-K
real*8 :: de_f                                      ! Saturation Density of Farmers (number/km2)
real*8 :: de_hg                                     ! Saturation Density of Hunter-Gatherers (number/km2)
real*8 :: eta_value                                 ! Interbreeding intensity between Farmers and Hunter-Gatherers
real*8 :: rr_hg_value                               ! Net Growth Population factor Hunter-Gatherers per generation 

integer:: t                                         ! Time (gen) of each simulation
integer:: gen_time                                  ! Generation (years) 
integer:: n_runs                                    ! Number of simulations

real*8,parameter:: jum_p     = 0.155D0              ! Fraction jumping each direction without Haplo-K 
real*8,parameter:: jum_n     = 0.155D0              ! Fraction jumping each direction with Haplo-K
real*8 :: jum_p_c                                   ! Fraction jumping each direction for corners without Haplo-K 
real*8 :: jum_n_c                                   ! Fraction jumping each direction for corners with Haplo-K
real*8 :: jum_p_s                                   ! Fraction jumping each direction for sides without Haplo-K 
real*8 :: jum_n_s                                   ! Fraction jumping each direction for sides with Haplo-K

integer :: dis_p                                    ! Grid size (km) of Farmer's dispersal 

real*8  :: a_fmax                                   ! Maximum population of Farmers (individuals/cell)
real*8  :: a_hgmax                                  ! Maximum population of Hunter-Gatherers (individuals/cell)
integer :: randjump                                 ! Threshold population number of farmers
integer :: randjump_s                               ! Threshold population number of farmers (side). Only for Sea Model
real*8  :: rr_p                                     ! Net Growth Population factor Farmers without Haplo-K per generation  
real*8  :: rr_n                                     ! Net Growth Population factor Farmers with Haplo-K per generation
real*8 arrival                                      ! Arrival condition, expressed as number of farmers when it is at least 0.10*a_fmax

r_p = 0.00D0                                        ! Service variable dispersal
 
t = 185                ! Time (number of generations)
n_runs = 1             ! Number of runs (they yield different results because we use random numbers).
gen_time = 32          ! Number of years per generation
eta_value = 0.06D0     ! Interbreeding intensity between Farmers and Hunter-Gatherers  
rr_hg_value = 1.0D0    ! Net Growth Population factor Hunter-Gatherers
rr_p = 2.45D0          ! Net Growth Population factor Farmers without Haplo-K 
rr_n = 2.45D0          ! Net Growth Population factor Farmers with Haplo-K
pe_p = 0.38D0          ! Persitency Farmers without Haplo-K 
pe_n = 0.38D0          ! Persitency Farmers with Haplo-K
de_f = 1.28D0          ! Saturation Density of Farmers (number/km2)
de_hg = 0.064D0        ! Saturation Density of Hunter-Gatherers (number/km2)
fraction = 0.523D0     ! Fraction of farmers with HK at starting point (same value as that obtained by trial and error in the INLAND program for this value of eta).
randjump   = 1897      ! Threshold population number of farmers (Inside). It is a function of Grid size (km) of Farmer's dispersal (1/4). See Secs. S3-B and S5-A.
randjump_s = 1423      ! Threshold population number of farmers (Side). It is a function of Grid size (km) of Farmer's dispersal (1/3). See Secs. S3-B and S5-A.
dis_p = 70             ! Grid size (km) of Farmer's dispersal 											
 
ARQ_FRA = 60           ! Node # of Arqueological Site southern France (from Data S4)
ARQ_NES = 62           ! Node # of Arqueological Site NE Spain
ARQ_NAV = 69           ! Node # of Arqueological Site Navarre
ARQ_POR = 82           ! Node # of Arqueological Site Portugal
ARQ_ITA = 31           ! Node # of Arqueological Site Italy
ARQ_CYP = 10           ! Node # of Arqueological Site Cyprus
  
row_sea = 150          ! Row Node # of starting point Sea Model (X=150 nodes). X increases to the LEFT.
col_sea = 0            ! Column Node # of starting point Sea Model (Y=0 nodes). Y increaes upwards.
 
 open (201,file='Fraction_Sea_Model.dat',status='UNKNOWN', access='APPEND')
   write(201,*) " time(gen)   dist(km)  %HK   region               eta "
 close(201)

 open (556,file='Archaelogical_Time_Sea_Model.dat',status='UNKNOWN', access='APPEND')
   write(556,*) "dist(km) time(gen)    region "
 close(556)
 
do random_t = 1,n_runs  
    
!===================================================================================
! Initial Conditions

 CRT_05 = 0
 CRT_06 = 0
 CRT_07 = 0
 CRT_08 = 0
 CRT_09 = 0
 CRT_13 = 0

 a_fmax    = NINT(de_f *dis_p*dis_p)
 a_hgmax   = NINT(de_hg*dis_p*dis_p) 
 arrival   = NINT(0.1d0*a_fmax)
 
 jum_p_c   = (1-pe_p)/2.0D0       ! Fraction jumping each direction for corners without Haplo-K 
 jum_n_c   = (1-pe_n)/2.0D0       ! Fraction jumping each direction for corners with Haplo-K
 jum_p_s   = (1-pe_p)/3.0D0       ! Fraction jumping each direction for sides without Haplo-K 
 jum_n_s   = (1-pe_n)/3.0D0       ! Fraction jumping each direction for sides with Haplo-K

 
 a_hg  = a_hgmax 
 eta   = eta_value
 rr_hg = rr_hg_value
 
 a_n   = 0.00D0
 a_p   = 0.00D0
 b_n   = 0.00D0
 b_p   = 0.00D0
 
 a_n (row_sea,col_sea) = dnint(fraction *(a_fmax))
 a_p (row_sea,col_sea) = a_fmax - a_n(row_sea,col_sea)   
 a_hg(row_sea,col_sea) = 0.00D0
   
!=================================================================================
! Time Integration:
 !call SEED(RND$TIMESEED)
 call RANDOM_SEED()
 
 Print '(" Processing.... ")'

 do j=1,t     
   
!=================================================================================
! Farmers without haplogroup K, DISPERSAL:
do y=0,300
do i=0,300
     
 if (y.eq.col_sea) then      
    if(a_p(i,col_sea) + a_n(i,col_sea) .le. randjump_s) then  
     b_p(i,col_sea) = b_p(i,col_sea) + a_p(i,col_sea)
     goto 100  
   endif
 else        
    if (a_p(i,y) + a_n(i,y) .le. randjump) then   
     b_p(i,y) = b_p(i,y) + a_p(i,y)
     goto 100  
   endif
 endif
 
 ! coordinates (i,y): i increases to the LEFT and y increases upwards. 
 ! (i,y) here is (X,Y) in Data S2, where X increases to the LEFT and Y increases upwards.

 select case (y)
     case(0)
       select case(i)
         case(0)   ! Corner SE: (i,y)=(0,0). 
            b_p(i,y) = (b_p(i,y)+(pe_p)   *(a_p(i,y))) 
               
            b_p(i+1,y) = (b_p(1,y)  +(jum_p_c)*(a_p(i,y)))  
            b_p(i,y+1) = (b_p(i,y+1)+(jum_p_c)*(a_p(i,y)))  
                
         case(300)  ! Corner SW: (i,y)=(300,0).
            b_p(i,y) = (b_p(i,y)+(pe_p)   *(a_p(i,y))) 
                
            b_p(i,y+1) = (b_p(i,y+1)+(jum_p_c)*(a_p(i,y)))  
            b_p(i-1,y) = (b_p(i-1,y)+(jum_p_c)*(a_p(i,y)))  
                
         case default  ! Side S: y=0.
            b_p(i,y)   =  (b_p(i,y)  +(pe_p)   *(a_p(i,y))) 
              
     !       b_p(i,y)   =  (b_p(i,y)  +(jum_p_s)*(a_p(i+1,y)))  
     !       b_p(i,y)   =  (b_p(i,y)  +(jum_p_s)*(a_p(i-1,y)))
     !       b_p(i,y)   =  (b_p(i,y)  +(jum_p)  *(a_p(i,y+1)))
                
            b_p(i+1,y) =  (b_p(i+1,y)+(jum_p_s)*(a_p(i,y)))  
            b_p(i-1,y) =  (b_p(i-1,y)+(jum_p_s)*(a_p(i,y))) 
            b_p(i,y+1) =  (b_p(i,y+1)+(jum_p_s)*(a_p(i,y)))  
         end select

     case(300)
       select case(i)
         case(0) ! Corner NE: (i,y)=(0,300).
           b_p(i,y) =  (b_p(i,y)+(pe_p)   *(a_p(i,y)))
                
           b_p(i+1,y) =  (b_p(i+1,y)+(jum_p_c)*(a_p(i,y)))  
           b_p(i,y-1) =  (b_p(i,y-1)+(jum_p_c)*(a_p(i,y)))   
                
         case(300) ! Corner NW: (i,y)=(300,300).            
           b_p(i,y) =  (b_p(i,y)+(pe_p)   *(a_p(i,y)))
     
           b_p(i-1,y) =  (b_p(i-1,y)+(jum_p_c)*(a_p(i,y)))  
           b_p(i,y-1) =  (b_p(i,y-1)+(jum_p_c)*(a_p(i,y)))   
     
         case default  ! Side N: y=300.  
           b_p(i,y)   =  (b_p(i,y)  +(pe_p)   *(a_p(i,y)))
                
     !      b_p(i,y)   =  (b_p(i,y)  +(jum_p_s)*(a_p(i+1,y)))
     !      b_p(i,y)   =  (b_p(i,y)  +(jum_p_s)*(a_p(i-1,y)))
     !      b_p(i,y)   =  (b_p(i,y)  +(jum_p)  *(a_p(i,y-1)))
                
           b_p(i+1,y) =  (b_p(i+1,y)+(jum_p_s)*(a_p(i,y)))  
           b_p(i-1,y) =  (b_p(i-1,y)+(jum_p_s)*(a_p(i,y))) 
           b_p(i,y-1) =  (b_p(i,y-1)+(jum_p_s)*(a_p(i,y)))  
       end select
         
     case default
       select case(i)
         case(0)    ! Side E: i=0.
           b_p(i,y)   =   (b_p(i,y)  +(pe_p)   *(a_p(i,y)))
                
    !       b_p(i,y)   =   (b_p(i,y)  +(jum_p_s)*(a_p(i,y-1))) 
    !       b_p(i,y)   =   (b_p(i,y)  +(jum_p_s)*(a_p(i,y+1))) 
    !       b_p(i,y)   =   (b_p(i,y)  +(jum_p)  *(a_p(i+1,y))) 
                
           b_p(i,y-1) =  (b_p(i,y-1)+(jum_p_s)*(a_p(i,y))) 
           b_p(i,y+1) =  (b_p(i,y+1)+(jum_p_s)*(a_p(i,y))) 
           b_p(i+1,y) =  (b_p(i+1,y)+(jum_p_s)*(a_p(i,y))) 
                
         case(300)  ! Side W: i=300.
           b_p(i,y)   =  (b_p(i,y)  +(pe_p)   *(a_p(i,y)))
               
    !       b_p(i,y)   =  (b_p(i,y)  +(jum_p_s)*(a_p(i,y-1)))
    !       b_p(i,y)   =  (b_p(i,y)  +(jum_p_s)*(a_p(i,y+1)))
    !       b_p(i,y)   =  (b_p(i,y)  +(jum_p)  *(a_p(i-1,y)))
                 
           b_p(i,y-1) =  (b_p(i,y-1)+(jum_p_s)*(a_p(i,y))) 
           b_p(i,y+1) =  (b_p(i,y+1)+(jum_p_s)*(a_p(i,y))) 
           b_p(i-1,y) =  (b_p(i-1,y)+(jum_p_s)*(a_p(i,y)))
                						   
         case default    
           if (a_p(i,y) + a_n(i,y) .le. randjump) then   
             b_p(i,y) = b_p(i,y) + a_p(i,y)
             goto 804  
           endif
    
           decimal  = jum_p*a_p(i,y) - dint(jum_p*a_p(i,y))
           jum_suma = 0.0d0
    
           do z=-1,1
	       do x=-1,1
    	     r_p= dis_p*sqrt(real(x)**2+real(z)**2)
             
             if (r_p==0) then
               goto 802
             endif

  		     if (x==0.or.z==0) then
               call RANDOM_NUMBER(aleatori)                                
                 if(aleatori.gt.decimal) then
                    b_p(i+x,y+z) = (b_p(i+x,y+z)+dint(jum_p*a_p(i,y)))
                    jum_suma = jum_suma + dint(jum_p*a_p(i,y))
                 else
                    b_p(i+x,y+z) = (b_p(i+x,y+z)+(1+dint(jum_p*a_p(i,y))))
                    jum_suma = jum_suma + dint(jum_p*a_p(i,y)) + 1
                 endif
                 goto 802
             end if
802            continue              
           end do  ! z
           end do  ! x
    
515      continue
    
           if (b_p(i,y)+a_p(i,y).lt.jum_suma) then       
             call RANDOM_NUMBER(aleatori)        
           if(aleatori.le.0.25) then
             b_p(i,y-1) = b_p(i,y-1)-1
             jum_suma   = jum_suma  -1
           else if (aleatori.gt.0.25.and.aleatori.le.0.50) then
             b_p(i+1,y) = b_p(i+1,y)-1
             jum_suma   = jum_suma  -1
           else if (aleatori.gt.0.50 .and.aleatori.le.0.75) then
             b_p(i,y+1) = b_p(i,y+1)-1
             jum_suma   = jum_suma  -1
           else
             b_p(i-1,y) = b_p(i-1,y)-1    
             jum_suma   = jum_suma  -1
           endif
           goto 515
           else      
             b_p(i,y) = b_p(i,y)+a_p(i,y)-jum_suma
           endif 
       
804      continue   
            
       end select
 end select ! case y
 
100	end do
end do	

!=================================================================================
! Farmers with haplogoup K, DISPERSAL:
do y=0,300
do i=0,300
     
  if (y.eq.col_sea) then      
    if(a_n(i,col_sea) + a_p(i,col_sea) .le. randjump_s) then  
     b_n(i,col_sea) = b_n(i,col_sea) + a_n(i,col_sea)
     goto 1000  
   endif
  else   
    if (a_n(i,y) + a_p(i,y) .le. randjump) then  
     b_n(i,y) = b_n(i,y) + a_n(i,y)
     goto 1000  
    endif
  endif 

! coordinates (i,y): i increases to the LEFT and y increases upwards. 
! (i,y) here is (X,Y) in Data S2, where X increases to the LEFT and Y increases upwards.
 
  select case (y)
    case(0)
      select case(i)
        case(0)     ! Corner SE: (i,y)=(0,0).
          b_n(i,y) = (b_n(i,y)+(pe_n)   *(a_n(i,y))) 
                
          b_n(i+1,y) = (b_n(i+1,y)+(jum_n_c)*(a_n(i,y)))  
          b_n(i,y+1) = (b_n(i,y+1)+(jum_n_c)*(a_n(i,y)))  
                
        case(300)     ! Corner SW: (i,y)=(300,0).
          b_n(i,y) = (b_n(i,y)+(pe_n)   *(a_n(i,y))) 
              
          b_n(i,y+1) = (b_n(i,y+1)+(jum_n_c)*(a_n(i,y)))  
          b_n(i-1,y) = (b_n(i-1,y)+(jum_n_c)*(a_n(i,y)))  
                
        case default  ! Side S: y=0.
          b_n(i,y)   =  (b_n(i,y)  +(pe_n)   *(a_n(i,y))) 
                
     !     b_n(i,y)   =  (b_n(i,y)  +(jum_n_s)*(a_n(i+1,y)))  
     !     b_n(i,y)   =  (b_n(i,y)  +(jum_n_s)*(a_n(i-1,y)))
     !     b_n(i,y)   =  (b_n(i,y)  +(jum_n)*(a_n(i,y+1)))
                
          b_n(i+1,y) =  (b_n(i+1,y)+(jum_n_s)*(a_n(i,y)))  
          b_n(i-1,y) =  (b_n(i-1,y)+(jum_n_s)*(a_n(i,y))) 
          b_n(i,y+1) =  (b_n(i,y+1)+(jum_n_s)*(a_n(i,y)))  
      end select
        
    case(300)
      select case(i)
        case(0)       ! Corner NE: (i,y)=(0,300).
          b_n(i,y) =  (b_n(i,y)+(pe_n)   *(a_n(i,y)))
                
          b_n(i+1,y) =  (b_n(1,y)+(jum_n_c)*(a_n(i,y)))  
          b_n(i,y-1) =  (b_n(i,y-1)+(jum_n_c)*(a_n(i,y)))   
                
        case(300)      ! Corner NW: (i,y)=(300,300).            
          b_n(i,y) =  (b_n(i,y)+(pe_n)   *(a_n(i,y)))
     
          b_n(i-1,y) =  (b_n(i-1,y)+(jum_n_c)*(a_n(i,y)))  
          b_n(i,y-1) =  (b_n(i,y-1)+(jum_n_c)*(a_n(i,y)))   
     
        case default    ! Side N: y=0.  
          b_n(i,y)   =  (b_n(i,y)  +(pe_n)   *(a_n(i,y)))
                
    !      b_n(i,y)   =  (b_n(i,y)  +(jum_n_s)*(a_n(i+1,y)))
    !      b_n(i,y)   =  (b_n(i,y)  +(jum_n_s)*(a_n(i-1,y)))
    !      b_n(i,y)   =  (b_n(i,y)  +(jum_n)*(a_n(i,y-1)))
               
          b_n(i+1,y) =  (b_n(i+1,y)+(jum_n_s)*(a_n(i,y)))  
          b_n(i-1,y) =  (b_n(i-1,y)+(jum_n_s)*(a_n(i,y))) 
          b_n(i,y-1) =  (b_n(i,y-1)+(jum_n_s)*(a_n(i,y)))  
      end select
        
    case default
      select case(i)
        case(0)       ! Side E: i=0.
          b_n(i,y)   =   (b_n(i,y)  +(pe_n)   *(a_n(i,y)))
                
    !      b_n(i,y)   =   (b_n(i,y)  +(jum_n_s)*(a_n(i,y-1))) 
    !      b_n(i,y)   =   (b_n(i,y)  +(jum_n_s)*(a_n(i,y+1))) 
    !      b_n(i,y)   =   (b_n(i,y)  +(jum_n)  *(a_n(i+1,y))) 
                
          b_n(i,y-1) =  (b_n(i,y-1)+(jum_n_s)*(a_n(i,y))) 
          b_n(i,y+1) =  (b_n(i,y+1)+(jum_n_s)*(a_n(i,y))) 
          b_n(i+1,y) =  (b_n(i+1,y)+(jum_n_s)*(a_n(i,y))) 
                
        case(300)      ! Side W: i=300.
          b_n(i,y)   =  (b_n(i,y)  +(pe_n)   *(a_n(i,y)))
                
    !      b_n(i,y)   =  (b_n(i,y)  +(jum_n_s)*(a_n(i,y-1)))
    !      b_n(i,y)   =  (b_n(i,y)  +(jum_n_s)*(a_n(i,y+1)))
    !      b_n(i,y)   =  (b_n(i,y)  +(jum_n)  *(a_n(i-1,y)))
                 
          b_n(i,y-1) =  (b_n(i,y-1)+(jum_n_s)*(a_n(i,y))) 
          b_n(i,y+1) =  (b_n(i,y+1)+(jum_n_s)*(a_n(i,y))) 
          b_n(i-1,y) =  (b_n(i-1,y)+(jum_n_s)*(a_n(i,y))) 
                
        case default    
          if (a_n(i,y) + a_p(i,y) .le. randjump) then  
            b_n(i,y) = b_n(i,y) + a_n(i,y)
            goto 904  
          endif
    
          decimal  = jum_n*a_n(i,y) - dint(jum_n*a_n(i,y))
          jum_suma = 0.0d0
    
          do z=-1,1
		  do x=-1,1
            r_p=dis_p*sqrt(real(x)**2+real(z)**2)
	       
		    if (r_p==0) then
			  goto 820
		    endif

  			if (x==0.or.z==0) then
              call RANDOM_NUMBER(aleatori)                               
              if(aleatori.gt.decimal) then
                 b_n(i+x,y+z) = (b_n(i+x,y+z)+ dint(jum_n*a_n(i,y)))
                 jum_suma     = jum_suma     + dint(jum_n*a_n(i,y))
              else
                 b_n(i+x,y+z) = (b_n(i+x,y+z)+ (1+dint(jum_n*a_n(i,y))))
                 jum_suma     = jum_suma + 1 + dint(jum_n*a_n(i,y))
              endif    
			  goto 820
			end if

820         continue
            
		  end do  ! x
          end do  ! z
   
615     continue    
        
          if (b_n(i,y)+a_n(i,y).lt.jum_suma) then       
            call RANDOM_NUMBER(aleatori)        
          if(aleatori.le.0.25) then
            b_n(i,y-1) = b_n(i,y-1)-1
            jum_suma = jum_suma  -1
          else if (aleatori.gt.0.25.and.aleatori.le.0.50) then
            b_n(i+1,y) = b_n(i+1,y)-1
            jum_suma = jum_suma  -1
          else if (aleatori.gt.0.50 .and.aleatori.le.0.75) then
            b_n(i,y+1) = b_n(i,y+1)-1
            jum_suma = jum_suma  -1
          else
            b_n(i-1,y) = b_n(i-1,y)-1    
            jum_suma = jum_suma  -1
          endif
          goto 615
          else      
            b_n(i,y) = b_n(i,y)+a_n(i,y)-jum_suma
          endif 
  
904     continue         
            
      end select
  end select  ! y

1000 end do
 end do	

!=================================================================================

! Update time matrices
  a_p=b_p
  a_n=b_n

! Reset time matrices
  b_p=0.0D0
  b_n=0.0D0
  
!=================================================================================
! Arrival times to the oldest archaeological site in each region, for Fig. 1 (using the number of nodes from Data S4):

   open (555,file='Archaelogical_Time_Sea_Model.dat',status='UNKNOWN', access='APPEND')
   
   do i = 300,row_sea,-1
       
      if (i.eq.(row_sea+ARQ_CYP)) then 
       if(a_n(i,col_sea)+a_p(i,col_sea).gt.arrival.and.CRT_13.eq.0) then
            sitename = " 13 Cyprus "
            write(555,505) dis_p*(i-row_sea),j,sitename
            CRT_13 = 1
       endif     
      endif
  
      if (i.eq.(row_sea+ARQ_ITA)) then 
        if(a_n(i,col_sea)+a_p(i,col_sea).gt.arrival.and.CRT_09.eq.0) then
            sitename = " 9 Italy "
            write(555,505) dis_p*(i-row_sea),j,sitename
            CRT_09 = 1
        endif     
      endif
 
      if (i.eq.(row_sea+ARQ_FRA)) then 
        if(a_n(i,col_sea)+a_p(i,col_sea).gt.arrival.and.CRT_05.eq.0) then
            sitename = " 5 France "
            write(555,505) dis_p*(i-row_sea),j,sitename
            CRT_05 = 1
        endif
      endif 
  
      if (i.eq.(row_sea+ARQ_NES)) then 
        if(a_n(i,col_sea)+a_p(i,col_sea).gt.arrival.and.CRT_06.eq.0) then
            sitename = " 6 NE Spain "
            write(555,505) dis_p*(i-row_sea),j,sitename
            CRT_06 = 1
        endif
      endif 
  
      if (i.eq.(row_sea+ARQ_NAV)) then 
        if(a_n(i,col_sea)+a_p(i,col_sea).gt.arrival.and.CRT_07.eq.0) then
           sitename = " 7 Navarre "
           write(555,505) dis_p*(i-row_sea),j,sitename
           CRT_07 = 1
        endif
      endif
  
      if (i.eq.(row_sea+ARQ_POR)) then 
        if(a_n(i,col_sea)+a_p(i,col_sea).gt.arrival.and.CRT_08.eq.0) then
            sitename = " 8 Portugal "
            write(555,505) dis_p*(i-row_sea),j,sitename
            CRT_08 = 1
        endif
      endif
      
   enddo
   close(555)

505 format (i6,2x,i6,6x,A20) 
 
!=================================================================================        
! Reproduction and interbreeding:    
     
	do z=0,300
	do x=0,300
                    
        couples_hn = 0.0D0
        couples_hp = 0.0D0
        farmers    = 0.0D0
        
        couples_hn = dint(eta(x,z)*a_hg(x,z)*a_n(x,z)/(a_hg(x,z)+a_p(x,z)+a_n(x,z)))
        couples_hp = dint(eta(x,z)*a_hg(x,z)*a_p(x,z)/(a_hg(x,z)+a_p(x,z)+a_n(x,z)))
          
        a_hg(x,z) = dint(rr_hg(x,z)*(a_hg(x,z) - couples_hn - couples_hp))
        a_p(x,z)  = dint(rr_p* (a_p(x,z)  + couples_hn + couples_hp))        
        a_n(x,z)  = dint(rr_n* (a_n(x,z)))
       
       farmers = dint(a_p(x,z)+a_n(x,z))
       if(farmers > a_fmax) then
         a_p(x,z) = dint(a_p(x,z)*a_fmax/farmers)
         a_n(x,z) = dint(a_n(x,z)*a_fmax/farmers)
       end if
               
      if (a_hg(x,z)<1D-10) a_hg(x,z)      = 0.00D0
      if (a_hg(x,z).gt.a_hgmax) a_hg(x,z) = a_hgmax  

    end do
    end do

! show generation time on screen:
 print *,j

!=================================================================================
! Fraction of farmers with haplogroup K in each region, for Fig. 2b (usig the number of generations and nodes from Data S2):

    open (200,file='Fraction_Sea_Model.dat',status='UNKNOWN', access='APPEND')

    selectcase(j)
    case(80)
        sitename = " 11 Greece and North Macedonia "
        i = 27
        write(200,383)  j, dis_p*(i), 100*a_n(i+row_sea,col_sea)/(a_n(i+row_sea,col_sea)+a_p(i+row_sea,col_sea)), sitename ,eta(i+row_sea,col_sea)
       
    case(98) 
        sitename = " 12 Croatia "
        i = 39  
        write(200,383)  j, dis_p*(i), 100*a_n(i+row_sea,col_sea)/(a_n(i+row_sea,col_sea)+a_p(i+row_sea,col_sea)), sitename ,eta(i+row_sea,col_sea)
                
    case(108) 
        sitename = " 13 Italy "
        i = 47 
        write(200,383)  j, dis_p*(i), 100*a_n(i+row_sea,col_sea)/(a_n(i+row_sea,col_sea)+a_p(i+row_sea,col_sea)), sitename ,eta(i+row_sea,col_sea)

    case(134) 
        sitename = " 14 Southern France "
        i = 61 
        write(200,383)  j, dis_p*(i), 100*a_n(i+row_sea,col_sea)/(a_n(i+row_sea,col_sea)+a_p(i+row_sea,col_sea)), sitename ,eta(i+row_sea,col_sea)
    
    case(114) 
        sitename = " 15 Spain "
        i = 69 
        write(200,383)  j, dis_p*(i), 100*a_n(i+row_sea,col_sea)/(a_n(i+row_sea,col_sea)+a_p(i+row_sea,col_sea)), sitename ,eta(i+row_sea,col_sea)
                       
    case(130) ! This should be case(111) = 111 generations according to table S2, but then "Error M6101. Flotaing-point error" because the front has not arrived (Fig 1a). 
	          ! It is reasonable that it has not arrived because the model is homogeneous. Using a larger value of time is OK because %K changes <1% (Fig. S6c).
        sitename = " 16 Portugal "
        i = 89 
        write(200,383)  j, dis_p*(i), 100*a_n(i+row_sea,col_sea)/(a_n(i+row_sea,col_sea)+a_p(i+row_sea,col_sea)), sitename ,eta(i+row_sea,col_sea)

    endselect
    close (200)
	    
383 format(i5,7x,i7,2x,f7.2,x,A20,f6.4)
  
 end do ! End time loop
!===============================================================================

enddo ! random_t
end program
