!     .../Projects/2025_Perez_real_numbers/coast_Fort/14_de_12_origin_point_O_Figs_S3_S4a/2026_Figs_S3_S4a.f90

! Author: J. Fort

! THIS IS THE SIMPLEST CODE USED IN THIS PAPER (because it does not distinguish between farmers with and without haplogroup K). 
! THIS CODE HAS BEEN USED TO OBTAIN FIGS. S3 AND S4a.

! This code is a modification of the code by J. Fort used in Fort, AAS (2022). The major change is coupled
! logistic reproduction here (versus uncoupled logistics in AAS 2022).

! To obtain Figs. S3 and S4a we do not need to distinguish between farmers with and without haplogroup K.
! Thus we consider only the total number of farmers, i.e. the sum of Eqs. (S29) and (S30) for reproduction 
! and the sum of Eqs. (S33) and (S34) for interbreeding. 

integer::i,j,t,x,y,xx,yy

! For a distance of 22 km between nodes (as used below), to reach 3,000 km in both directions we need only 273 x 273 nodes (because 6,000/22=272.7). 
! However, we use 2,999 x 2,999 nodes, just in case we want to run simulations with smaller distances between nodes.
real,dimension (2999,2999)::N,N2,M,M2

real,dimension (2999)::rfront
real dis,pe,rN,Nmax,rP,Mmax,eta

! t = number of iterations = final number of generations:
t=120
! For t=120 generations, this code runs in 2 minutes on a personal computer with 16 GB RAM and 2.60 GHz.

! eta = interbreeding intensity:
eta=0.8
												    
! dis = DISTANCE BETWEEN NODES (km) in this code = length of inland jumps (km) = d in the main paper and Supp. Info.:
dis=22
! Thus in this code, an inland jump has the length of 1 node or 22 km.

! i = length of coastal jumps divided by length of inland jumps, i.e. DELTA/d in the paper & Supp. Info. (i must be an integer number):
i=1
! In fact we only need the case i=1 to obtain Figs. S3 and S4a, because they correspond to the case DELTA=d (eastern Mediterranean), so a coastal jump also has the length of 1 node or 22 km.

! rN = Neolithic initial growth rate, in gen^-1.
! Neolithic net fecundity = RoN = 2.45, so rN=LOG(2.45)=0.896 gen^-1 /32 yr = 0.028 yr^-1 (Isern, Fort & Perez-Losada, J. Stat. Mechs.: Theor. & Exp. 2008). 
! In Fortran, LOG is the natural logarithm:
rN=LOG(2.45)	  

! rP = Mesolithic initial growth rate, directly in gen^-1 (Fort, Pujol & Cavalli-Sforza, Cambridge Archaeol. J. 2004): 
rP=0.59

! Nmax = Neolithic carrying capacity (Currat & Excoffier, Proc. R. Soc. B 2005):
Nmax=1.28

! Mmax = Mesolithic carrying capacity (Currat & Excoffier, Proc. R. Soc. B 2005):
Mmax=0.064

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

! Initially there are farmers only at 1 cell at the center of the coast (simplest case; others would yield the same speed):
N=0
N(1,1500)=Nmax
N2=0

! Initially there are HGs everywhere except at the cell initially occupied by farmers:
M=Mmax
M(1,1500)=0
M2=0

! The coast is the left side of the grid, i.e. the line x=1 or (1,y).
! We save the initial Neolithic profile along the coast:
	OPEN (5,FILE='NEOLPROF0.DAT',STATUS='UNKNOWN')
  	do y=1,2999
	WRITE(5,*) y,N(1,y)	
	end do
	CLOSE(5)

! We save the initial Mesolithic profile along the coast:
	OPEN (5,FILE='MESOLPROF0.DAT',STATUS='UNKNOWN')
  	do y=1,2999
	WRITE(5,*) y,M(1,y)	
	end do
	CLOSE(5)

!================================================================================
! 1) ITERATIONS IN TIME (each iteration is one generation):
do j=1,t
!================================================================================
! 2) REPRODUCTION: COUPLED LOGISTIC TERMS FROM LAPOLICE, WILLIAMS & HUBER, Nature Comm. 2025:
	do y=1,2999
	do x=1,2999
	!We disregard almost empty nodes (otherwise the code runs very slowly):
		if (N(x,y)<1D-10) then
			else
			N2(x,y) = N(x,y) + rN * N(x,y) * (1-  ((N(x,y)+M(x,y))/Nmax) )
		endif
		! 2 lines above, I have used N2(x,y) rather than N(x,y) because otherwise N(x,y) would be different in the right-hand side there and 3 lines below.

		if (M(x,y)<1D-10) then															   
			else
		    M(x,y) = M(x,y) + rP * M(x,y) * (1-  ((N(x,y)+M(x,y))/Mmax) ) 
		end if	

		N(x,y) = N2(x,y)
	
	end do
	end do
	N2=0
!=================================================================================
! 3) INTERBREEDING:
	M2=M
	do y=1,2999
	do x=1,2999

	    ! We disregard almost empty nodes (otherwise, floating-point error):
		if (M(x,y)<1E-10) then
		else

		! 6 lines above I have written M2=M, i.e. M2(x,y)=M(x,y). The reason is that we use M2(x,y) rather than M(x,y) in the next 2 eqs., because
		! otherwise the value of M(x,y) in their right-hand sides would be different since the first equation changes the value of M(x,y):
			M(x,y)=M2(x,y)-eta*M2(x,y)*N(x,y)/(N(x,y)+M2(x,y))
			N(x,y)=N(x,y)+eta*M2(x,y)*N(x,y)/(N(x,y)+M2(x,y))
		end if
 
		if (M(x,y)<0) then 
			M(x,y)=0 
		end if


		if (N(x,y)>1.28) then 
			N(x,y)=1.28 
		end if	
	
	end do
	end do

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

!	First the inland nodes: do x=2,2999 instead of x=1,2999 because the coast is the side x=1, i.e. (1,y):
	do y=1,2999
	do x=2,2999
!-----------------------
	! We disregard almost empty nodes because there are many of them and otherwise, the code runs too slowly:
		if (N(x,y)<1E-10) then
		goto 100
		endif
!-----------------------
		! This is used above line 80 to compute jumps in 4 directions (up, down, left and right):
		do yy=-1,1
		do xx=-1,1	

!   N(x,y) is the farmer population density at the node from where farmers jump:
!   N2(x,y) is the farmer population density at the node into which farmers jump:
 
		! Jump of zero length:
		if (xx==0.and.yy==0) then
		N2(x,y)=N2(x,y)+(pe)*N(x,y)
		goto 80
		endif
		
		! In Fortran "or" is "or/and" but the "and" case has been already considered in the previous lines and jumped to line 80.
		! Jumps in 4 directions (up, down, left and right):
		if(xx==0.or.yy==0) then
		N2(x+xx,y+yy)=N2(x+xx,y+yy)+((1-pe)/4)*N(x,y)
		goto 80
		endif
		
		80 continue
		end do
		end do		

100	end do
	end do

!***********************
!	COAST: the coast is the vertical side, i.e. x=1 or (1,y):
	do y=2,2999
!-----------------------
		!We disregard almost empty nodes because there are many of them and otherwise, the code runs too slowly:
		if (N(1,y)<1E-10) then
		goto 150
		endif
!-----------------------
		! This is used below to compute jumps:
		do yy=-1,1
		do xx=-1,1
		
		! Jump of zero length:	
		if (xx==0.and.yy==0) then
		N2(1,y)=N2(1,y)+(pe)*N(1,y)
		goto 120
		endif

		! From coastal nodes, farmers jump to 3 nodes (not to 4 as above line 100). The 3 jumps are down, up & right: 
		! There is no jump to the left, i.e. (xx=-1.and.yy==0) because then (1+xx,y+yy)=(0,y) is on the sea.
		! The length of coastal jumps is i nodes (as already stated at the beginning of this code).

		! Jump downwards:
		if(xx==0.and.yy==-1) then
		N2(1,y-i)=N2(1,y-i)+((1-pe)/3)*N(1,y)
		goto 120
		endif

		! Jump upwards:
		if(xx==0.and.yy==1) then
		N2(1,y+i)=N2(1,y+i)+((1-pe)/3)*N(1,y)
		goto 120
		endif

		!The jump to the right is an inland jump, so its lenght is 1 node (not i nodes as along the coast):
		if(xx==1.and.yy==0) then
		N2(1+xx,y)=N2(1+xx,y)+((1-pe)/3)*N(1,y)
		goto 120
		endif
		
		120 continue
		end do
		end do
	150	end do	
		
!**********************

! We keep the farmer population density in matrix N, which will be fixed during the next iteration:
N=N2

! We re-start matrix N2, which will keep varying during the next iteration:
N2=0


!===============================================================================
! 5) WE KEEP SOME RESULTS BEFORE JUMPING TO THE NEXT ITERATION OR GENERATION:

! 5a) FRONT POSITION in the coast (vertical direction) as a function of time in generations (this file yields Fig. 3c):
OPEN (2,ACCESS='APPEND',FILE='FRONT.DAT',STATUS='UNKNOWN')

! I identify the last coastal node with farmer density >0.1, and find out the position of the front (farmer density = 0.1) by fitting a straight line:
! The coast is the side x=1, i.e. (1,y):

	do y=2999,50,-1						    
		! In the next line I use y<2250·22=49,500 km, which is more than enough because it is much larger that the necessary distances (about 2,500 km, see Fig. S3c):
		if (y>2250) then
		goto 400
		endif	
	if (N(1,y)>0.1) then
		rfront(j)=dis*(y+(N(1,y)-0.1)/(N(1,y)-N(1,y+1))-1500)
		exit
	end if
	400 continue
	end do
IF (.not.(j<1)) WRITE (2,*)  j,rfront(j)
close(2)

!---------------------------------------------------------------------------------------------	
! 5b) FRONT PROFILES ALONG THE COAST: 
! Neolithic coastal profile at time=60 generations (this file and the next one yield Fig. S3a):
	if (j==60) then
	OPEN (4,FILE='NEOLPROF60.DAT',STATUS='UNKNOWN')
   	do y=1,2999
	WRITE(4,*) y,N(1,y)
	end do
	CLOSE(4)
	end if

! Neolithic coastal profile at time=90 generations (this file and the previous one yield Fig. S3a):
	if (j==90) then
	OPEN (5,FILE='NEOLPROF90.DAT',STATUS='UNKNOWN')
   	do y=1,2999
	WRITE(5,*) y,N(1,y)
	end do
	CLOSE(5)
	end if

! Mesolithic coastal profile at time=60 generations  (this file and the next one yield Fig. S3b):
	if (j==60) then
	OPEN (6,FILE='MESOLPROF60.DAT',STATUS='UNKNOWN')
  	do y=1,2999
	WRITE(6,*) y,M(1,y)
	end do
	CLOSE(6)
	end if

! Mesolithic coastal profile at time=90 generations  (this file and the previous one yield Fig. S3b):
	if (j==90) then
	OPEN (7,FILE='MESOLPROF90.DAT',STATUS='UNKNOWN')
  	do y=1,2999
	WRITE(7,*) y,M(1,y)
	end do
	CLOSE(7)
	end if

!===============================================================================
! WE SHOW THE ITERATION NUMBER (I.E. THE TIME MEASURED IN GENERATIONS) ON THE SCREEN AND JUMP TO THE NEXT ONE:
print *,j
end do

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

print '(" Files NEOLPROF60.DAT and NEOLPROF90.DAT yield Fig. S3a, files MESOLPROF60.DAT and MESOLPROF90.DAT yield Fig. S3b, and file FRONT.dat yields Fig. S3c. ")'

end program