! .../Projects/2025_Perez_real_numbers/coast_Fort/47_de_28_FIG_1b_18km_418km/2026_Figs1b_2.f90

! Author: J. Fort
! This code is a modification of the code in Fort, AAS (2022). The major change is coupled
! logistic reproduction here (vs. indep. logistics in AAS 2022). 
!
!----------------------------------------------------------------------------------------
integer::i,j,y,t,z,x,n,m,Cyprus,CAnatolia,DAnatolia,EAnatolia,Crete,GGreece
real,dimension (6999,6999)::k,h,a,b
real,dimension (200)::rfront
real,dimension (6999)::East,Eastm,previous,ee
real dis,pe,rr,am,rp,bm,f,ga

! Iterations = final number of generations. Suggestion: t=145
t=145

! Cultural horitzontal transmission (vertical if eta=f and ga=1):
f=0
ga=1
												    
! Distance between nodes (km):
dis=1

! rr=LOG(2.45) = 0.896 gen^-1 = 0.028 yr^-1 (in fortran LOG is the Neperian logarithm):
rr=LOG(2.45)	  

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

! Neolithic carrying capacity (Currat+Excoffier 2005):
am=1.28
! Mesolithic carrying capacity (Currat+Excoffier 2005):
bm=0.064

! Persistency (Fort, Pérez-Losada & Isern 2007):
pe=0.38

! Initally only farmers at point O (y=100):
a=0
a(1,100)=1.28
b=0

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

! We keep the initial profile of the Neolithic population density:
	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 profile of the Mesolithic population density:
	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. We know from running this code that at t=121 gen the EASTERN front is at rfont=1394 km 
! from point O, i.e. y=1394+100=1494 km. So we assume the WESTERN front starts at this positon and 122 gen:
		if (j==122) then
		a=0
		a(1,1494)=1.28
		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 almost zero population density (otherwise the code is very slow):
		if (a(x,z)<1D-10) then
			else
			b(x,z) = a(x,z) + rr * a(x,z) * (1-  ((a(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-  ((a(x,z)+k(x,z))/bm) ) 
		end if	

		a(x,z) = b(x,z)
	
	end do
	end do
	b=0
!=================================================================================
! 3) INTERBREEDING:
	h=k
	do z=1,6999
	do x=1,6999

	! We neglect nodes with almost zero population density (otherwise, floating-point error):
		if (k(x,z)<1E-10) then
		else
			k(x,z)=h(x,z)-f*h(x,z)*a(x,z)/(a(x,z)+ga*h(x,z))
			a(x,z)=a(x,z)+f*h(x,z)*a(x,z)/(a(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
!-----------------------
! We neglect nodes with almost zero population density (otherwise the code is very slow):
		if (a(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
		b(i,y)=b(i,y)+(pe)*a(i,y)
		goto 80
		endif
		
		if(x==0.or.z==0) then
		b(i+m*x,y+m*z)=b(i+m*x,y+m*z)+((1-pe)/4)*a(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
!-----------------------
	! We neglect nodes with almost zero population density (otherwise the code is very slow):
		if (a(1,y)<1E-10) then
		goto 150
		endif
!-----------------------
		do z=-1,1
		do x=-1,1	
		if (x==0.and.z==0) then
		b(1,y)=b(1,y)+(pe)*a(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
		b(1,y+n*z)=b(1,y+n*z)+((1-pe)/3)*a(1,y)
		goto 120
		endif

		if(x==0.and.z==1) then
		b(1,y+n*z)=b(1,y+n*z)+((1-pe)/3)*a(1,y)
		goto 120
		endif

		! this is an inland jump, so its lenght is m, not n:
		if(x==1.and.z==0) then
		b(1+m*x,y)=b(1+m*x,y)+((1-pe)/3)*a(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=418	  
! 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 rfont=1394 km
! from point O, i.e. y=1394+100=1494 km. So we assume that the WESTERN front starts ar 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
!-----------------------
! We neglect nodes with almost zero population density (otherwise the code is very slow):
		if (a(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
		b(i,y)=b(i,y)+(pe)*a(i,y)
		goto 280
		endif
		
		if(x==0.or.z==0) then
		b(i+m*x,y+m*z)=b(i+m*x,y+m*z)+((1-pe)/4)*a(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
!-----------------------
! We neglect nodes with almost zero population density (otherwise the code is very slow):
		if (a(1,y)<1E-10) then
		goto 350
		endif
!-----------------------
		do z=-1,1
		do x=-1,1	
		if (x==0.and.z==0) then
		b(1,y)=b(1,y)+(pe)*a(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
		b(1,y+n*z)=b(1,y+n*z)+((1-pe)/3)*a(1,y)
		goto 320
		endif

		if(x==0.and.z==1) then
		b(1,y+n*z)=b(1,y+n*z)+((1-pe)/3)*a(1,y)
		goto 320
		endif

		! this is an inland jump, so its lenght is m, not n:
		if(x==1.and.z==0) then
		b(1+m*x,y)=b(1+m*x,y)+((1-pe)/3)*a(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 calculated population density in matrix a, which will be fixed during the next iteration:
a=b

! We re-initiate matrix b, which will vary during the next iteration:
b=0

!===============================================================================
! 5) WE KEEP SOME RESULTS BEFORE JUMPINT TO THE NEXT TIME:
!!------------------
! 5a) This file keeps the front position along the vertical direction (coast) as a function of time. 
! The slope is the spread rate in km/generation (you can import it with origin + linear fit):
OPEN (2,ACCESS='APPEND',FILE='RFRONT.DAT',STATUS='UNKNOWN')
! We identify the first node with rho<0.1, and find the point with rho=0.1 (by linear interpolation).
! COAST= SIDE OF GRID: I change a(1500,x) into a(1,zz) and x into zz in 5 places. 
! [I use zz, not z, to avoid the error "cannot link"!]:
	do zz=6999,100,-1						    
		! I use this because in the upper corner, i.e. (x,y)=(1,6999), the code fails (the Mesolithic
		! density grows and it should not), so I look below the upper corner, i.e. 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 regresssions. 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 120:
	if (j==120) 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 120:
	if (j==120) 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

!===============================================================================
! 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 121 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 '("  ")'
print '(" WESTERN MEDITERRANEAN: use THE FILE FIG1b.dat, it considers only nodes with farmers and is thus more accurate than RFRONT.dat ")'
print '(" To see the profiles, import NEOPROF.DAT, etc.  ")'
print '("  ")'

end program