! .../Projects/2025_Perez_real_numbers/coast_Fort/37a_de_28_Fig_4__eta_0/2026_xx_xx__Fig_3.f90
! Author: J. Fort
! This code is a modification of the code Fort, AAS (2022). The major change is coupled	logistic 
! reproduction here (vs. indep. logistics in AAS 2022). Genetics as in Sci.Rep.2017 & Nat.Com.2025.
!----------------------------------------------------------------------------------------
integer::i,j,y,t,z,x,n,m,Cyprus,CAnatolia,DAnatolia,EAnatolia,Crete,GGreece
real,dimension (6999,6999)::k,h,a,an,ap,bn,bp
real,dimension (200)::rfront
real,dimension (6999)::East,Eastm,previous,ee
real dis,pe,rr,am,rp,bm,f,ga,fraction

! Iterations = final number of generations: t=172 because Data S2, region 5, 171 gen.
t=172

! Horizontal (or vertical, for ga=1 and eta=f) cultural transmission parameters:
f=0
ga=1
												    
! distance between nodes (km):
dis=1

! Ro=2.45 as in Fort, PNAS 2017. In fortran LOG is the natural logarithm, so 
! a=rr=ln(Ro)=ln(2.45)=0.896 gen^-1 · 1 gen/32 yr =0.028 yr^-1:
rr=LOG(2.45)	  

! Mesolithic initial growth rate in gen^-1 (Fort, Pujol & Cavalli-Sforza 2004): 
rp=0.59

! Neolithic satutarion density (Currat+Excoffier 2005):
am=1.28
! Mesolithic saturation density (Currat+Excoffier 2005):
bm=0.064

! Persistency (Fort et al., Phys. Rev. E 2007):
pe=0.38

!================================================================================
! INITIAL CONDITIONS:
fraction= 0.474   ! Fraction of farmers with haplogroup K at point O.
! Initally only farmers at point O (y=100); an=K farmers, ap=non-K farmers, k=HGs:
an=0
ap=0
a=0
an(1,100)=1.28*fraction
ap(1,100)=1.28-an(1,100)
a(1,100)=an(1,100)+ap(1,100)
bn=0
bp=0

! Initially HGs everywhere except at point O:
k=0.064
k(1,100)=0
h=0

! We keep the initial Neolithic profile:
	OPEN (5,FILE='NEOPROF0.DAT',STATUS='UNKNOWN')
  	do zz=1,6999
	WRITE(5,*) zz,a(1,zz)	
	end do
	CLOSE(5)

! We keep the initial Mesolithic profile:
	OPEN (5,FILE='MESOPROF0.DAT',STATUS='UNKNOWN')
  	do zz=1,6999
	WRITE(5,*) zz,k(1,zz)	
	end do
	CLOSE(5)

! This deletes the file RFRONT.DAT:
OPEN (2,FILE='RFRONT.DAT',STATUS='UNKNOWN')
WRITE (2,*) ' '
CLOSE (2)

!================================================================================
! 1) INTEGRATION IN TIME:
do j=1,t
!================================================================================

! The dispersal in the WESTERN Mediterranean begins as close as possible to the intersection of the 2 regresssions (121.5 gen and 1411 km from point O).  
! We know from running this code that at t=121 gen rfront(121)=1394 km, i.e. the EASTERN front is at 1394 km from point O or y=1394+100=1494 km.
! We use the %K from the EASTERN expansion as initial %K for the WESTERN expansion, but there are no farmers at y=1494 km because 1394/18=77.44 is not
! an integer. The closest node to y=1494 km (1394 km from point O) to which farmers arrive is y=1504 km (1404 km from point O) because 1404/18=78 is
! an integer. So we have to evaluate the %K at y=1504 km and apply it at y=1494 km, immediately after all calculations for the EASTERN expansion and 
! t=121 gen, so t=122 gen:
		if (j==122) then
		!print *, an(1,1504)  ! result: 1.990595E-02, different from zero-->OK
		!print *, ap(1,1504)  ! result: 2.212932E-02, different from zero-->OK
		bcd=an(1,1504)
		efg=ap(1,1504)
		an=0
		ap=0
		a=0
		an(1,1494)=bcd*1.28/(bcd+efg)	! This line and next imply an+ap=1.28=am, and this is necessary to obtain rfront.dat and Fig1b.dat as for Fig. 1b.
		ap(1,1494)=efg*1.28/(bcd+efg)
		bn=0
		bp=0
		b=0
		k=0.064
		k(1,1494)=0
		h=0
		endif

! 2) REPRODUCTION: COUPLED LOGISTIC TERMS FROM LAPOLICE, WILLIAMS & HUBER, Nature Comm. 2025:
	do z=1,6999
	do x=1,6999
	! We neglect nodes with approx. zero population density (otherwise it is very slow):
		if (an(x,z)<1D-10) then
			else
			bn(x,z) = an(x,z) + rr * an(x,z) * (1-  ((an(x,z)+ap(x,z)+k(x,z))/am) )
		endif

 		if (ap(x,z)<1D-10) then
			else
			bp(x,z) = ap(x,z) + rr * ap(x,z) * (1-  ((an(x,z)+ap(x,z)+k(x,z))/am) )
		endif

		if (k(x,z)<1D-10) then															   
			else
		    k(x,z) = k(x,z) + rp * k(x,z) * (1-  ((an(x,z)+ap(x,z)+k(x,z))/bm) ) 
		end if	

		an(x,z) = bn(x,z)
		ap(x,z) = bp(x,z)

	end do
	end do
	bn=0
	bp=0
!=================================================================================
! 3) INTERBREEDING:
	h=k
	do z=1,6999
	do x=1,6999

	! We neglect nodes with approx. zero Mesol. popul. density (otherwise, floating-point error!):
		if (k(x,z)<1E-10) then
		else
			k(x,z)=h(x,z)-f*h(x,z)*(an(x,z)+ap(x,z))/(an(x,z)+ap(x,z)+ga*h(x,z))
			ap(x,z)=ap(x,z)+f*h(x,z)*(an(x,z)+ap(x,z))/(an(x,z)+ap(x,z)+ga*h(x,z))
		end if
 
		if (k(x,z)<0) then 
			k(x,z)=0 
		end if

		if (a(x,z)>1.28) then 
			a(x,z)=1.28 
		end if	
	
	end do
	end do

	h=0
!================================================================================
! 4) NEOLITHIC DISPERSAL: 

! The dispersal in the WESTERN Mediterranean begins at 122 generations, see also line 200 below:
		if (j>121) then
		goto 200
		endif

! ======DISPERSAL IN THE EASTERN MEDITERRANEAN:==================================    
! length of coastal jumps (km) in the EASTERN Mediterranean:
 		n=18
! length of inland jumps (km) in the EASTERN Mediterranean:
 		m=18

! INLAND jumps: i=2,6999 (not i=1,6999) because the coast is the straight line (i,y) with i=1:
	do y=1,6999
	do i=2,6999
!-----------------------
! If initial values are approx. 0, we jump to the next node (otherwise it is very slow!!):
		if (an(i,y)+ap(i,y)<1E-10) then
		goto 100
		endif
!-----------------------
		do z=-1,1
		do x=-1,1	
! a(i,y) is the density of farmers at the initial node of jumps:

		if (x==0.and.z==0) then
		bn(i,y)=bn(i,y)+(pe)*an(i,y)
		bp(i,y)=bp(i,y)+(pe)*ap(i,y)
		goto 80
		endif
		
		if(x==0.or.z==0) then
		bn(i+m*x,y+m*z)=bn(i+m*x,y+m*z)+((1-pe)/4)*an(i,y)
		bp(i+m*x,y+m*z)=bp(i+m*x,y+m*z)+((1-pe)/4)*ap(i,y)
		goto 80
		endif
		
		80 continue
		end do
		end do		

100	end do
	end do

!***********************
!   COASTAL jumps: the coast is the straight line (i,y) with i=1:
	do y=1,6999
!-----------------------
! If initial values are approx. 0, we jump to the next node (otherwise it is very slow!!):
		if (an(1,y)+ap(1,y)<1E-10) then
		goto 150
		endif
!-----------------------
		do z=-1,1
		do x=-1,1	
		if (x==0.and.z==0) then
		bn(1,y)=bn(1,y)+(pe)*an(1,y)
		bp(1,y)=bp(1,y)+(pe)*ap(1,y)
		goto 120
		endif

! From coastal nodes, farmers jump to 3 nodes (not to 4 as above line 100):  
		if(x==0.and.z==-1) then
		bn(1,y+n*z)=bn(1,y+n*z)+((1-pe)/3)*an(1,y)
		bp(1,y+n*z)=bp(1,y+n*z)+((1-pe)/3)*ap(1,y)
		goto 120
		endif

		if(x==0.and.z==1) then
		bn(1,y+n*z)=bn(1,y+n*z)+((1-pe)/3)*an(1,y)
		bp(1,y+n*z)=bp(1,y+n*z)+((1-pe)/3)*ap(1,y)
		goto 120
		endif

		! this is an inland jump, so its lenght is m, not n:
		if(x==1.and.z==0) then
		bn(1+m*x,y)=bn(1+m*x,y)+((1-pe)/3)*an(1,y)
		bp(1+m*x,y)=bp(1+m*x,y)+((1-pe)/3)*ap(1,y)
		goto 120
		endif
! Node (x=-1.and.z==0) is not included because (1+x,y+z)=(0,y) is sea (not coast). 
		
		120 continue
		end do
		end do
	150	end do	
	
	goto 370

! ======DISPERSAL IN THE WESTERN MEDITERRANEAN:==============================
! length of coastal jumps (km) in the WESTERN Mediterranean:
	200 n=430	  
! length of inland jumps (km) in the WESTERN Mediterranean:
 		m=50

! The dispersal in the WESTERN Mediterranean begins as close as possible to the intersection of 
! the 2 regresssions. We know from running this code that at t=121 gen the EASTERN front is at 
! rfront=1394 km from point O, i.e. y=1394+100=1494 km. So we assume that the WESTERN front starts 
! at this position and 122 gen:

! INLAND jumps: i=2,6999 (not i=1,6999) because the coast is the straight line (i,y) with i=1:
	do y=1494,6999
	do i=2,6999
!-----------------------
! If initial values are approx. 0, we jump to the next node (otherwise it is very slow!!):
		if (an(i,y)+ap(i,y)<1E-10) then
		goto 300
		endif
!-----------------------
		do z=-1,1
		do x=-1,1	
! a(i,y) is the density of farmers at the initial node of jumps:

		if (x==0.and.z==0) then
		bn(i,y)=bn(i,y)+(pe)*an(i,y)
		bp(i,y)=bp(i,y)+(pe)*ap(i,y)
		goto 280
		endif
		
		if(x==0.or.z==0) then
		bn(i+m*x,y+m*z)=bn(i+m*x,y+m*z)+((1-pe)/4)*an(i,y)
		bp(i+m*x,y+m*z)=bp(i+m*x,y+m*z)+((1-pe)/4)*ap(i,y)
		goto 280
		endif
		
		280 continue
		end do
		end do		

300	end do
	end do

!***********************
!   COASTAL jumps: the coast is the straight line (i,y) with i=1:
	do y=1494,6999
!-----------------------
! If initial values are approx. 0, we jump to the next node (otherwise it is very slow!!):
		if (an(1,y)+ap(1,y)<1E-10) then
		goto 350
		endif
!-----------------------
		do z=-1,1
		do x=-1,1	
		if (x==0.and.z==0) then
		bn(1,y)=bn(1,y)+(pe)*an(1,y)
		bp(1,y)=bp(1,y)+(pe)*ap(1,y)
		goto 320
		endif

! From coastal nodes, farmers jump to 3 nodes (not to 4 as above line 100):  
		if(x==0.and.z==-1) then
		bn(1,y+n*z)=bn(1,y+n*z)+((1-pe)/3)*an(1,y)
		bp(1,y+n*z)=bp(1,y+n*z)+((1-pe)/3)*ap(1,y)
		goto 320
		endif

		if(x==0.and.z==1) then
		bn(1,y+n*z)=bn(1,y+n*z)+((1-pe)/3)*an(1,y)
		bp(1,y+n*z)=bp(1,y+n*z)+((1-pe)/3)*ap(1,y)
		goto 320
		endif

		! this is an inland jump, so its lenght is m, not n:
		if(x==1.and.z==0) then
		bn(1+m*x,y)=bn(1+m*x,y)+((1-pe)/3)*an(1,y)
		bp(1+m*x,y)=bp(1+m*x,y)+((1-pe)/3)*ap(1,y)
		goto 320
		endif
! Node (x=-1.and.z==0) is not included because (1+x,y+z)=(0,y) is sea (not coast). 
		
		320 continue
		end do
		end do
	350	end do	
		


	370 continue
		
!**********************

! We keep the Neolithic population in matrices a, to be used in the next iteration:
an=bn
ap=bp

! We reinitiate matrices b, they will change during the next iteration:
bn=0
bp=0

a=an+ap

!===============================================================================
! 5) WE KEEP SOME RESULTS BEFORE JUMPING TO THE NEXT TIME STEP:

!------------------
! 5a) FRONT POSITION along the vertical direction as a function of time: 
! the slope is the spread rate in km/generation (import with e.g. origin + linear fit)
OPEN (2,ACCESS='APPEND',FILE='RFRONT.DAT',STATUS='UNKNOWN')
! I idenfity the first node with rho<0.1, and find the point with rho=0.1 (by fitting a straight line):
! COAST=EDGE OF THE GRID =(1,zz) [I use zz because z causes error "cannot link"!]:
	do zz=6999,100,-1						    
		! the upper limit of coast (x,y)=(1,6999) fails, so I look at y<6250:
		if (zz>6250) then
		goto 400
		endif

	! zz+n appears below because the node zz+1 is always empty of farmers: they arrive only to coastal nodes separated n km.	
	if (a(1,zz)>0.1) then
		rfront(j)=dis*(zz+(a(1,zz)-0.1)*n/(a(1,zz)-a(1,zz+n))-100)
		exit
	end if
	400 continue
	end do
IF (.not.(j<1)) WRITE (2,*)  j,rfront(j)
close(2)

!------------------
! 5b) FRONT PROFILES: 
! The dispersal in the WESTERN Mediterranean begins as close as possible to the intersection of 
! the 2 regressions. We know from running this code that at t=121 gen the EASTERN front is at y=1394+100=1494 km
	
! Neolithic profile at generation 140:
	if (j==140) then
	OPEN (4,FILE='NEOPROF.DAT',STATUS='UNKNOWN',POSITION='APPEND')

	do zz=1,1493
	WRITE(4,*) zz,East(zz)
	end do

	do zz=1494,6999
	WRITE(4,*) zz,a(1,zz)
	end do

	CLOSE(4)
	end if

! Mesolithic profile at generation 140:
	if (j==140) then
	OPEN (5,FILE='MESOPROF.DAT',STATUS='UNKNOWN',POSITION='APPEND')

	do zz=1,1493
	WRITE(5,*) zz,Eastm(zz)
	end do

  	do zz=1494,6999
	WRITE(5,*) zz,k(1,zz)
	end do

	CLOSE(5)
	end if

!------------------
! 5c) Dates for Fig. 1b (using distances along the coast in km from Data S1):
	OPEN (6,FILE='FIG1b.DAT',STATUS='UNKNOWN',POSITION='APPEND')

 ! Eastern Mediterranean:

	if (j==1) then
	WRITE(6,*) 0,0
	endif

	! Below, no need to subtract 100 km (initial position of farmers) because already subtracted in rfront above.

	if (Cyprus==1) then
	goto 500
	endif
	if (rfront(j)>453) then
	WRITE(6,*) 453,j-1+(453-rfront(j-1))/(rfront(j)-rfront(j-1))
	Cyprus=1 
	endif
 	500 continue

	if (CAnatolia==1) then
	goto 503
	endif	
	if ((rfront(j))>583) then
 	WRITE(6,*) 583,j-1+(583-rfront(j-1))/(rfront(j)-rfront(j-1))
	CAnatolia=1
	endif 
 	503 continue

	if (DAnatolia==1) then
	goto 510
	endif	
	if ((rfront(j))>846) then
 	WRITE(6,*) 846,j-1+(846-rfront(j-1))/(rfront(j)-rfront(j-1))
	DAnatolia=1
	endif
 	510 continue
	
	if (EAnatolia==1) then
	goto 513
	endif	
	if ((rfront(j))>1122) then
 	WRITE(6,*) 1122,j-1+(1122-rfront(j-1))/(rfront(j)-rfront(j-1))
	EAnatolia=1
	endif
 	513 continue

	if (Crete==1) then
	goto 516
	endif	
	if ((rfront(j))>1172) then
 	WRITE(6,*) 1172,j-1+(1172-rfront(j-1))/(rfront(j)-rfront(j-1))
	Crete=1
	endif
 	516 continue

	if (GGreece==1) then
	goto 520
	endif	
	if ((rfront(j))>1231) then
 	WRITE(6,*) 1231,j-1+(1231-rfront(j-1))/(rfront(j)-rfront(j-1))
	GGreece=1
	endif
 	520 continue

	if (j==121) then
 	WRITE(6,*) NINT(rfront(121)),121
	endif

	! Western Mediterranean:
	! Using the distance of the oldest site per region does not lead to a straight line. I use one node every n km:

	690 continue

	do xx=4,14

	if (ee(xx)==1) then 
	goto 700
	endif

 	if (a(1,1494+n*xx)>0.1) then
	WRITE(6,*) 1494+n*xx-100,j-(a(1,1494+n*xx)-0.1)/(a(1,1494+n*xx)-previous(1494+n*xx))
	ee(xx)=1
	end if

	700 continue
	end do

	CLOSE(6)
!===============================================================================

!	This is used to obtain the times in Fig. 1b (otherwise we get no straight line in the West Med.):
		do yy=1,6999
		previous(yy)=a(1,yy)
		end do

!------------------
!5d) Numbers of generations, km from point O and %K for Fig. 3 (regions in Data S2):

OPEN (7,FILE='FIG3.DAT',STATUS='UNKNOWN',POSITION='APPEND')

!---------------------------------------------------------------------	
! We include some values of %K in Fig.3 between point 0 and region 1 because it is a curve, not a straight line.
! The time (j) is not important because %K changes <1% in 20 gen (Fig. S6c in Fort & Pérez, Nature Comm. 2025).
! We need a node with farmers: its distance from point O is an integer multiple of n. y is this distance +100km.

if (j==50) then
y = 17*n+100
write(7,*) j, y-100, 100*an(1,y)/(an(1,y)+ap(1,y))
endif

if (j==70) then
y = 34*n+100
write(7,*) j, y-100, 100*an(1,y)/(an(1,y)+ap(1,y))
endif

if (j==90) then
y = 50*n+100
write(7,*) j, y-100, 100*an(1,y)/(an(1,y)+ap(1,y))
endif

if (j==120) then
y = 67*n+100
write(7,*) j, y-100, 100*an(1,y)/(an(1,y)+ap(1,y))
endif

if (j==121) then
y = 77*n+100
write(7,*) j, y-100, 100*an(1,y)/(an(1,y)+ap(1,y))
endif

!---------------------------------------------------------------------	
        
selectcase(j)

case(1)		! 1 is the number of generations j; below, j=0 gen and i=0 km for region 1 = initial region (it contains point O):
	write(7,*)  0, 0, 47.4

!---------------------------------------------------------------------	
case(126)	! Region 2 = Greece and North Macedonia: 126 generations. In fact Data S2 says 116 generations, so in
			! principle this should be case(116) = 116 generations. However, then "Error M6101. Flotaing-point error"  
			! because the front has not yet arrived at 1670+100 km (Fig 1b or rfront.dat), i.e. an=0, ap=0 so %K=0/0.
			! It is reasonable that it has not arrived because the model is homogeneous.
			! But we know that at 126 gen it has arrived from rfront.dat. Using a larger value of time 
			! is OK because %K changes <1% in 20 generations (Fig. S6c in Fort & Pérez-Losada, Nature Comm. 2025).

    y = 1670+100	! 100 km is the initial position of farmers, i.e. the position of point O in this code. 
					! 1670 km is the distance from point O to region 2 (Data S2).

			! The node y=1670+100 km is empty of farmers because they jump 430 km from y=1494 km, so an=0, ap=0 so %K=0/0. 
			! We look for the closest node where farmers have arrived from a coastal jump:

 	abc=ABS(y-1494)

	do xx=1,14

	if (ABS(y-(1494+n*xx))<ABS(y-(1494+n*(xx-1)))) then 
	abe=1494+n*xx
	abc=ABS(y-(1494+n*xx))
	endif
	continue
	end do

    write(7,*) j, y-100, 100*an(1,abe)/(an(1,abe)+ap(1,abe))
!---------------------------------------------------------------------
       
case(135) ! Region 3 = Croatia: 135 is the number of generations j (from Data S2).	 
    y = 2527+100 
 	abc=ABS(y-1494)
	do xx=1,14

	if (ABS(y-(1494+n*xx))<ABS(y-(1494+n*(xx-1)))) then 
	abe=1494+n*xx
	abc=ABS(y-(1494+n*xx))
	endif
	continue
	end do

    write(7,*) j, y-100, 100*an(1,abe)/(an(1,abe)+ap(1,abe))
!---------------------------------------------------------------------
                
case(144) ! Region 4 = Italy: 144 is the number of generations j (from Data S2). 
    y = 3081+100 
 	abc=ABS(y-1494)
	do xx=1,14

	if (ABS(y-(1494+n*xx))<ABS(y-(1494+n*(xx-1)))) then 
	abe=1494+n*xx
	abc=ABS(y-(1494+n*xx))
	endif
	continue
	end do

    write(7,*) j, y-100, 100*an(1,abe)/(an(1,abe)+ap(1,abe))
!---------------------------------------------------------------------

case(171) ! Region 5 = Southern France: 171 is the number of generations j (from Data S2). 
    y = 4066+100 
 	abc=ABS(y-1494)
	do xx=1,14

	if (ABS(y-(1494+n*xx))<ABS(y-(1494+n*(xx-1)))) then 
	abe=1494+n*xx
	abc=ABS(y-(1494+n*xx))
	endif
	continue
	end do

    write(7,*) j, y-100, 100*an(1,abe)/(an(1,abe)+ap(1,abe))
!---------------------------------------------------------------------
    
case(150) ! Region 6 = Spain: 150 is the number of generations j (from Data S2).
    y = 4574+100 
 	abc=ABS(y-1494)
	do xx=1,14

	if (ABS(y-(1494+n*xx))<ABS(y-(1494+n*(xx-1)))) then 
	abe=1494+n*xx
	abc=ABS(y-(1494+n*xx))
	endif
	continue
	end do

    write(7,*) j, y-100, 100*an(1,abe)/(an(1,abe)+ap(1,abe))
!---------------------------------------------------------------------
                       
case(147) ! Region 7 = Portugal: 147 is the number of generations j (from Data S2).
    y = 5987+100 
 	abc=ABS(y-1494)
	do xx=1,14

	if (ABS(y-(1494+n*xx))<ABS(y-(1494+n*(xx-1)))) then 
	abe=1494+n*xx
	abc=ABS(y-(1494+n*xx))
	endif
	continue
	end do

    write(7,*) j, y-100, 100*an(1,abe)/(an(1,abe)+ap(1,abe))

endselect
close (7)	

!===============================================================================
! We keep the final profile in the Eastern Mediterranean:
		if (j==121) then

		do yy=1,6999
		East(yy)=a(1,yy)
		end do

		do yy=1,6999
		Eastm(yy)=k(1,yy)
		end do

		endif

! WE SHOW THE GENERATION NUMBER ON THE SCREEN AND JUMP TO THE NEXT ONE:
print *,j
end do

print '(" To obtain the values for Fig. 1b, import file FIG1b.dat ")'
print '("  ")'

print '(" EASTERN MEDITERRANEAN: To obtain the spread rate, import the file RFRONT.DAT, make a linear regression up to e.g. 100 generations and ")'
print '(" the spread rate in km/gen is the slope. In Fig. 2a we have fitted the regression from 51 to 100 generations")'
print '(" Multiply this spread rate by 32 to convert it into km/yr. ")'
print '("  ")'

print '(" WESTERN MEDITERRANEAN: To obtain the spread rate, import the file RFRONT.DAT, make a linear regression from e.g. 122 generations and ")'
print '(" the spread rate in km/gen is the slope.")'
print '(" Multiply this spread rate by 32 to convert it into km/yr. ")'

print '(" To see the profiles, import NEOPROF.DAT, etc.  ")'
print '("  ")'

print '(" Please find the results for Fig. 3 in file FIG3.dat ")'

end program
