! .../Projects/2025_Perez_real_numbers/coast_Fort/30a_de_28_Fig_S4b_400x50km_eta_0/2026_xx_xx__Fig_S4b.f90
! Author: J. Fort
! This code is a modification of the code by Quim Fort used 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
real,dimension (6999,6999)::k,h,a,b
real,dimension (200)::rfront
real,dimension (6999)::previous,ee
real dis,pe,rr,am,rp,bm,f,ga

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

! cultural horizontal transmission (or vertical, with eta=f and ga=1):
f=0
ga=1
												    
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, directly in gen^-1 (Fort, Pujol & Cavalli-Sforza 2004): 
rp=0.59

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

! persistency:
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 save 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 save 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 file RFRONT.DAT:
OPEN (2,FILE='RFRONT.DAT',STATUS='UNKNOWN')
WRITE (2,*) ' '
CLOSE (2)

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

! 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 (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 approx. 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: 

! length of coastal jumps (nodes, 1 node=50 km):
 		n=8
! length of inland jumps (nodes, 1 node=50 km):
 		m=1

! 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 (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
!-----------------------
! if initial values are approx. 0, we jump to the next node (otherwise it 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+m*x,y+z)=(1-m,y) is sea (not coast) because 1-m<=0 (coast is i=1). 
		
		120 continue
		end do
		end do
	150	end do	
	
!**********************

! We keep the Neolithic population in matrix a, to be used in the next iteration:
a=b

! We reinitiate matrix b, it will change during the next iteration:
b=0

!===============================================================================
! 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 regresssions. We know from running this code that at t=121 gen the EASTERN front is at y=1394 km
	
! Neolithic profile at generation 120:
	if (j==120) then
	OPEN (4,FILE='NEOPROF.DAT',STATUS='UNKNOWN',POSITION='APPEND')

	do zz=1,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,6999
	WRITE(5,*) zz,k(1,zz)
	end do

	CLOSE(5)
	end if

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

	! 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:

	do xx=4,14

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

 	if (a(1,100+n*xx)>0.1) then
	WRITE(6,*) 100+n*xx-100,j-(a(1,100+n*xx)-0.1)/(a(1,100+n*xx)-previous(100+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 SHOW THE GENERATION NUMBER ON THE SCREEN AND JUMP TO THE NEXT ONE:
print *,j
end do

print '(" To obtain the values for Fig. S4b, import file FIG1b.dat ")'
print '("  ")'
print '(" To obtain the spread rate (Fig. 2b), import the file RFRONT.DAT and ")'
print '(" the spread rate in node /gen is the slope.")'
print '(" Divide this spread rate by 32 and multiply it by 50 to convert it into km/yr. ")'
print '("  ")'
print '(" To see the profiles, import NEOPROF.DAT, etc.  ")'
print '("  ")'

end program