! Author: J. Pérez-Losada
! e-mail: joaquim.perez@udg.edu
! This program simulates the Neolithic dispersal along the INLAND route (red lines in Figs. 1b, 3 and 4)
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Output files:
!  Archaeological_Time_Land_Model.dat (times of first arrival of farmers at the given archaeological sites)
!  Fraction_Land_Model.dat           (fraction of farmers with Haplo-K at given time and distance)
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
   
program main

implicit none
integer:: i,j,x,y,z,random_t         ! Service variables
integer:: row,col                    ! 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 region                     

integer :: DCRT_01,DCRT_02,DCRT_03   ! Control stwiches to find time of arrival to a site 1D
integer :: DCRT_04,DCRT_05,DCRT_06,DCRT_07  

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         ! Density of Farmers (number/km2)
real*8 :: de_hg        ! Density of Hunter-Gatherers (number/km2)
real*8 :: eta_value    ! Mixing coefficient btw 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
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 (please vary it until you obtain,e.g., 47.4% at genetic region 1 Upper Mesopotamia)
randjump = 970       ! Threshold population number of farmers (inland). It is a function of the Grid size below of Farmer's dispersal (1/4). See Secs. S3-B and S5-A.
dis_p = 50           ! Grid size (km) of Farmer's dispersal 
 
 row = 150      ! Row Node # of starting point Land Model
 col = 150      ! Column Node # of starting point Land Model

 open (201,file='Fraction_Land_Model.dat',status='UNKNOWN', access='APPEND')
   write(201,*) " time(gen)   dist(km)  %HK   region               eta "
 close(201)

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

 DCRT_01 = 0
 DCRT_02 = 0
 DCRT_03 = 0
 DCRT_04 = 0
 DCRT_05 = 0
 DCRT_06 = 0
 DCRT_07 = 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,col) = dnint(fraction *(a_fmax))
 a_p (row,col) = a_fmax - a_n(row,col)   
 a_hg(row,col) = 0.00D0
   
!=================================================================================
! Time Integration:

 call RANDOM_SEED()

 Print '(" Processing.... ")'

 do j=1,t     
     
!=================================================================================
! Farmers without haplogroup K, dispersal:
do y=0,300
do i=0,300
           
    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
 
 ! 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 haplogroup K, dispersal:
do y=0,300
do i=0,300
     
  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

! 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 time to the oldest archaeological site in each region, for Fig.1 (using the number of nodes from Data S4): 

    open (556,file='Archaeological_Time_Land_Model.dat',status='UNKNOWN', access='APPEND')
    
    if(a_n(row+0,row-0)+a_p(row+0,row-0).gt.arrival.and.DCRT_01.eq.0) then
            sitename =  " A  Syria  " 
            write(556,505)  nint(sqrt(real((dis_p*(0))**2)+real((dis_p*(0))**2))),j,sitename
            DCRT_01 = 1
    endif    
      
    if(a_n(row+10,row+4)+a_p(row+10,row+4).gt.arrival.and.DCRT_02.eq.0) then
            sitename =  " B  Anatolia " 
            write(556,505)  nint(sqrt(real((dis_p*(10))**2)+real((dis_p*(4))**2))),j,sitename
            DCRT_02 = 1            
    endif    
    
    if(a_n(row+32,row+21)+a_p(row+37,row+16).gt.arrival.and.DCRT_03.eq.0) then
            sitename =  " C  Serbia " 
            write(556,505)  nint(sqrt(real((dis_p*(37))**2)+real((dis_p*(16))**2))),j,sitename
            DCRT_03 = 1
    endif      
    
    if(a_n(row+40,row+34)+a_p(row+39,row+34).gt.arrival.and.DCRT_04.eq.0) then
            sitename =  " D  Germany LBK " 
            write(556,505)  nint(sqrt(real((dis_p*(39))**2)+real((dis_p*(34))**2))),j,sitename
            DCRT_04 = 1
    endif      
    
    if(a_n(row+39,row+44)+a_p(row+39,row+44).gt.arrival.and.DCRT_06.eq.0) then
            sitename =  " Sweden " 
            write(556,505)  nint(sqrt(real((dis_p*(39))**2)+real((dis_p*(44))**2))),j,sitename
            DCRT_06 = 1
    endif
      
    if(a_n(row+62,row+34)+a_p(row+62,row+34).gt.arrival.and.DCRT_07.eq.0) then
            sitename =  " E Engl Wales Scot " 
            write(556,505)  nint(sqrt(real((dis_p*(62))**2)+real((dis_p*(34))**2))),j,sitename
            DCRT_07 = 1
    endif  

  close(556)

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. 2a (using the number of generations and nodes from Data S2):
   
   open (201,file='Fraction_Land_Model.dat',status='UNKNOWN', access='APPEND')

   selectcase(j) 
    case(30) 
       sitename = " 1  Northern Mesopotamia  "
            i =  3 
            z = -2
        write(201,383)  j, nint(sqrt(real((dis_p*(i))**2)+real((dis_p*(z))**2))), 100*a_n(z+row,i+col)/(a_n(z+row,i+col)+a_p(z+row,i+col)), sitename,eta(z+row,i+col)      
    
    case(52) 
       sitename = " 2  Central Anatolia "
            i =  4 
            z =  9 
        write(201,383)  j, nint(sqrt(real((dis_p*(i))**2)+real((dis_p*(z))**2))), 100*a_n(z+row,i+col)/(a_n(z+row,i+col)+a_p(z+row,i+col)), sitename,eta(z+row,i+col)      
    
    case(77)
         sitename = " 3  Western Anatolia "    
            i =  10 
            z =  16 
        write(201,383)  j, nint(sqrt(real((dis_p*(i))**2)+real((dis_p*(z))**2))), 100*a_n(z+row,i+col)/(a_n(z+row,i+col)+a_p(z+row,i+col)), sitename,eta(z+row,i+col)      
                
    case(92)   
        sitename = " 4   Bulgaria (except Malak Preslavets) "
            i =  16 
            z =  22 
        write(201,383)  j, nint(sqrt(real((dis_p*(i))**2)+real((dis_p*(z))**2))), 100*a_n(z+row,i+col)/(a_n(z+row,i+col)+a_p(z+row,i+col)), sitename,eta(z+row,i+col)      
    
    case(101) 
        sitename = " 5 Romania and  Serbia (except Iron Gates) "
            i =  20 
            z =  26 
        write(201,383)  j, nint(sqrt(real((dis_p*(i))**2)+real((dis_p*(z))**2))), 100*a_n(z+row,i+col)/(a_n(z+row,i+col)+a_p(z+row,i+col)), sitename,eta(z+row,i+col)      
    
    case(107) 
        sitename = " 6 Hungary "
            i =  25 
            z =  32 
        write(201,383)  j, nint(sqrt(real((dis_p*(i))**2)+real((dis_p*(z))**2))), 100*a_n(z+row,i+col)/(a_n(z+row,i+col)+a_p(z+row,i+col)), sitename,eta(z+row,i+col)      

    case(109) 
        sitename = " 7 Austria & Czech Republic "
            i =  29 
            z =  36 
        write(201,383)  j, nint(sqrt(real((dis_p*(i))**2)+real((dis_p*(z))**2))), 100*a_n(z+row,i+col)/(a_n(z+row,i+col)+a_p(z+row,i+col)), sitename,eta(z+row,i+col)      

    case(113) 
        sitename = " 8 Germany LBK "
            i =  33 
            z =  45 
        write(201,383)  j, nint(sqrt(real((dis_p*(i))**2)+real((dis_p*(z))**2))), 100*a_n(z+row,i+col)/(a_n(z+row,i+col)+a_p(z+row,i+col)), sitename,eta(z+row,i+col)      

    case(170) 
        sitename = " 9 Scandinavia "
            i =  47 
            z =  38 
        write(201,383)  j, nint(sqrt(real((dis_p*(i))**2)+real((dis_p*(z))**2))), 100*a_n(z+row,i+col)/(a_n(z+row,i+col)+a_p(z+row,i+col)), sitename,eta(z+row,i+col)      

    case(128) 
        sitename = " 10 Northern France "
            i =  28 
            z =  55 
        write(201,383)  j, nint(sqrt(real((dis_p*(i))**2)+real((dis_p*(z))**2))), 100*a_n(z+row,i+col)/(a_n(z+row,i+col)+a_p(z+row,i+col)), sitename,eta(z+row,i+col)      
       
    endselect
    close(201)
            
383 format(i5,7x,i7,2x,f7.2,x,A20,f6.4)
    
 end do ! End time loop
!===============================================================================

enddo ! random_t
end program

