! .../Projects/2025_Perez_real_numbers/coast_Fort/14__de_13__origin_point_O__22x22km/Figs_S3_S4a.f90
! Author: J. Fort
! This code is a modification of the code by J. Fort used in Fort, AAS (2022). The major change is coupled
! logistic reproduction here (vs. indep. logistics in AAS 2022) 
!
!---------------------------------------------------------------------------------------
! Per exectuar amb "Developer studio 1994-95" (Fortran 90): Build Compile..., Build Build..., Build Execute...
!
! Per a que indiqui la línia on hi ha error: F4
!
! Per sortir i poder executar-ne un altre: File Close Worksplace (si no l'he executat: File Close) 
!			   
! Si dóna error "compiler is out of heap space in pass 2"-> possibles solucions:
!	1/ passar arxiu *.f90 de subcarpeta a carpeta -> ruta més curta
!	2/ tornar a fer copy paste del programa vell des d'una altra carpeta i canvis

! Si diu: "... error LNK2005: _main already defined in ...[nom d'un altre programa].obj
!          ... fatal error LNK 1169: one or more multiply defined symbols found",
!	és perquè hi ha 2 programes oberts-> fer File Close Worksplace i sortir del fortran.

! Si diu: "LINK : fatal error LNK1168: cannot open xxx.exe for writing",
!	crec que és pq hi ha una versió prèvia del mateix programa oberta ->
!	fer File Close Worksplace i sortir del fortran.
!----------------------------------------------------------------------------------------
integer::i,j,y,t,z,x,n,Cyprus,CAnatolia,DAnatolia,EAnatolia,Crete,GGreece,HGreece,SItaly,CItaly,Arene,Languedoc,Catalonia,Valencia,Andalusia,SPortugal,CPortugal
real,dimension (2999,2999)::k,h,a,b
real,dimension (2999)::rfront
real dis,pe,rr,am,rp,bm,f,ga

! Iterations = final number of generations:
t=120

! Horizontal cultural transmission (or vertical, for ga=1):
f=0
ga=1
												    
! Distance between nodes = length of jumps:
dis=22

! Length of coastal jumps (it must be an integer multiple of dis) divided by dis:
n=1

! Ro=2.45, so rr=LOG(2.45)0.896 gen^-1 /32 yr = 0.028 yr^-1 (Isern et al. 2008). 
! In Fortran, LOG is the natural logarithm.
rr=LOG(2.45)	  

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

! Neolithic carrying capacity (Courrant+Excoffier 2005):
am=1.28

! Mesolithic carrying capacity (Courrant+Excoffier 2005):
bm=0.064

! Persistency (Fort, PRE 2007):
pe=0.38

! Initially farmers at 1 coastal cell (simplest case; others would yield the same speed):
a=0
a(1,1500)=1.28
b=0

! Initially HGs everywhere except at the cell initially occupied by farmers:
k=0.064
k(1,1500)=0
h=0

! We save the initial Neolithic profile along the coast:
	OPEN (5,FILE='NEOLPROF0.DAT',STATUS='UNKNOWN')
  	do zz=1,2999
	WRITE(5,*) zz,a(1,zz)	
	end do
	CLOSE(5)

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

!================================================================================
! 1) INTEGREM EN EL TEMPS:
do j=1,t
!================================================================================
! 2) REPRODUCTION: COUPLED LOGISTIC TERMS FROM LAPOLICE, WILLIAMS & HUBER, Nature Comm. 2025:
	do z=1,2999
	do x=1,2999
	!menyspreem els punts amb densitat quasi nula (si no triga molt!):
		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,2999
	do x=1,2999

	    !menyspreem els punts amb densitat quasi nula (si no, dóna 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:

!   Faig primer l'interior (do i=2,2999 enlloc de i=1,2999) [la costa és la recta (i,y) amb i=1]:
	do y=1,2999
	do i=2,2999
!-----------------------
!   Si els valors inicials son aprox 0 passem a un nou punt (si no triga molt!!):
		if (a(i,y)<1E-10) then
		goto 100
		endif
!-----------------------
		do z=-1,1
		do x=-1,1	

!   a(i,y) és la densitat al punt d'origen dels salts:

		if (x==0.and.z==0) then
		b(i+x,y+z)=b(i+x,y+z)+(pe)*a(i,y)
		goto 80
		endif
		
		if(x==0.or.z==0) then
		b(i+x,y+z)=b(i+x,y+z)+((1-pe)/4)*a(i,y)
		goto 80
		endif
		
		80 continue
		end do
		end do		

100	end do
	end do

!***********************
!   COSTA: la costa és la recta (x,y) amb x=1:
	do y=2,2999
!-----------------------
!   Si els valors inicials son aprox 0 passem a un nou punt (si no triga molt!!):
		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+x,y+z)=b(1+x,y+z)+(pe)*a(1,y)
		goto 120
		endif

!   Dels punts de la costa salta població a 3 punts (no a 4 com abans de línia 100):
		if(x==0.and.z==-1) then
		b(1+x,y+n*z)=b(1+x,y+n*z)+((1-pe)/3)*a(1,y)
		goto 120
		endif

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

		! this is an inland jump, so its lenght is z, not n*z:
		if(x==1.and.z==0) then
		b(1+x,y+z)=b(1+x,y+z)+((1-pe)/3)*a(1,y)
		goto 120
		endif
! no incloc el punt (x=-1.and.z==0) pq (1+x,y+z)=(0,y) és mar. De fet, incloure'l potser
! donaria igual. El que sí importa és que als altres 3 punts he posat 1/3, no 1/4.
		
		120 continue
		end do
		end do
	150	end do	
		
!**********************

! guardem la població calculada a la matriu a, que serà fixa al proper bucle:
a=b

! reiniciem la matriu b, que anirà variant durant el proper bucle:
b=0

!===============================================================================
! 6) MESOLITHIC DISPERSAL:
!Aquestes 2 línies, les explico més avall:
h=k
k=0

!La costa és la recta (i,y) amb i=1: 
	do y=1,2999
	do i=2,2999
! Va lentíssim (en 1 hora no fa la 1ª iteració=generació)->als punts saturats no calculem:
		if (0.062-h(i,y)<1E-10) then
		goto 200
		endif
!-----------------------
		do z=-1,1
		do x=-1,1
		
		if (x==0.and.z==0) then
		k(i+x,y+z)=k(i+x,y+z)+(pe)*h(i,y)
		goto 180
		endif

! Si calculéssim com a la dispersió neolítica, calculariem als punts saturats, però
! aniria massa lent (veure el comentari anterior): en 1 hora no faria la 1a iteració.
! Per tant, aquí calculem directament la densitat al punt de DESTÍ dels salts=k(i,y),
! i per això he posat més amunt h=k i k=0.
!Explicació afegida (MOLT MÉS CLARA) 12/2/2021 (també al programa del PNAS 2012, d'on he tret la dispersió paleolítica):
!Població neolítica: hi ha molts nodes buits. Calcular a ells seria massa lent. Com no en SURT població, els ignorem: calculem la població que SURT dels nodes no buits. NOTA: No podem calcular la població que ARRIBA a cada node ignorant els buits (pq als buits sí n'hi arriba, i s'ha de calcular).
!Població paleolítica: hi ha molts nodes plens. Calcular a ells seria massa lent. No podem calcular la població que SURT de cada node (com amb la població neolítica) ignorant els plens, pq surt població dels plens (i en va cap a nodes no plens, que s'ha de calcular). Podem calcular la població que ARRIBA a cada node ignorant els plens? Sí, pq als plens es perd població per dispersió (si hi ha nodes no plens propers), però amb la reproducció es torna a saturar abans de la propera dispersió . En efecte, si el front és una línia de nodes saturats, com a màxim cada node perd 1/4=0.25 de Pmax, però amb la reproducció és multiplica per rp=0.59 i, per tant, torna a quedar saturat (sense tenir en compte l'efecte de la transmissió cultural, que ve després de la reproducció i per tant no modifica que s'ha saturat).	
		if(x==0.or.z==0) then
		k(i,y)=k(i,y)+((1-pe)/4)*h(i+x,y+z)
		goto 180
		endif
		
		180 continue
		end do
		end do		

200	end do
	end do
	
!***********************
!   COSTA: la costa és la recta (x,y) amb x=1:
	do y=1,2999
!-----------------------
! Com abans de línia 180, tampoc a la costa calculem als punts saturats:
		if (0.062-h(i,y)<1E-10) then
		goto 250
		endif
!-----------------------
		do z=-1,1
		do x=-1,1	
		if (x==0.and.z==0) then
		k(1+x,y+z)=k(1+x,y+z)+(pe)*h(1,y)
		goto 220
		endif

!   Dels punts de la costa salta població a 3 punts (no a 4 com abans de línia 200):
		if(x==0.and.z==-1) then
		k(1+x,y+n*z)=k(1+x,y+n*z)+((1-pe)/3)*h(1,y)
		goto 220
		endif

		if(x==0.and.z==1) then
		k(1+x,y+n*z)=k(1+x,y+n*z)+((1-pe)/3)*h(1,y)
		goto 220
		endif

		! this is an inland jump, so its lenght is z, not n*z:
		if(x==1.and.z==0) then
		k(1+x,y+z)=k(1+x,y+z)+((1-pe)/3)*h(1,y)
		goto 220
		endif
! no incloc el punt (x=-1.and.z==0) pq (1+x,y+z)=(0,y) és mar. De fet, incloure'l potser
! donaria igual. El que sí importa és que als altres 3 punts he posat 1/3, no 1/4.

		220 continue
		end do
		end do		
	250	end do

!*********************

k=h

!===============================================================================
! 7) GUARDEM ALGUNS RESULTATS ABANS DE PASSAR AL SEGÜENT INSTANT DE TEMPS:
!--------------------------------------------------------------------------------------------- 
! 7a) Arxiu que guarda la POSICIÓ del front EN DIRECCIÓ VERTICAL en funció del temps en generacions: 

OPEN (2,ACCESS='APPEND',FILE='FRONT.DAT',STATUS='UNKNOWN')

! identifico el primer punt amb rho<0.1, i trobo punt de rho=0.1 (ajustant recta):
! COSTA=VORA XARXA: CANVIO a(1500,x) per a(1,zz) i x per zz a 5 llocs 
! [poso zz, no z, pq. z dóna error "cannot link"!]:

	do zz=2999,50,-1						    
		! A la cantonada superior de la costa, i.e. (x,y)=(1,2999) falla (creix el núm. de
		! neolítics), per tant miro més avall de la cantonada superior, a y<2250:
		if (zz>2250) then
		goto 400
		endif	
	if (a(1,zz)>0.1) then
		rfront(j)=dis*(zz+(a(1,zz)-0.1)/(a(1,zz)-a(1,zz+1))-1500)
		exit
	end if
	400 continue
	end do
IF (.not.(j<1)) WRITE (2,*)  j,rfront(j)
close(2)

! 7b) Arxius que guarden el PERFIL del front a un cert instant (per a importar amb origin): 
!---------------------------------------------------------------------------------------------	
! perfil neolític a la generació 90:
	if (j==90) then
	OPEN (4,FILE='NEOLPROF.DAT',STATUS='UNKNOWN')
   	do zz=1,2999
	WRITE(4,*) zz,a(1,zz)
	end do
	CLOSE(4)
	end if

! perfil paleolític a la generació 90:
	if (j==90) then
	OPEN (5,FILE='MESOLPROF.DAT',STATUS='UNKNOWN')
  	do zz=1,2999
	WRITE(5,*) zz,k(1,zz)
	end do
	CLOSE(5)
	end if
!---------------------------------------------------------------------------------------------
! 7c) Dates for Fig. 1 (using distances along the coast in km from Data S1):
	OPEN (6,FILE='FIG1.DAT',STATUS='UNKNOWN',POSITION='APPEND')

 ! Eastern Mediterranean:

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

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

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

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

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

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

	if (j==106) then
 	WRITE(6,*) INT(rfront(106)),106
	endif

	if (j==107) then
 	WRITE(6,*) INT(rfront(107)),107
	endif

	! Western Mediterranean:

	if (HGreece==1) then
	goto 530
	endif	
	if ((rfront(j))>1948) then
 	WRITE(6,*) 1948,j
	HGreece=1
	endif
 	530 continue

	if (SItaly==1) then
	goto 553
	endif	
	if ((rfront(j))>2166) then
 	WRITE(6,*) 2166,j
	SItaly=1
	endif
 	553 continue

	if (CItaly==1) then
	goto 554
	endif	
	if ((rfront(j))>3044) then
 	WRITE(6,*) 3044,j
	CItaly=1
	endif
 	554 continue

	if (Arene==1) then
	goto 555
	endif	
	if ((rfront(j))>3565) then
 	WRITE(6,*) 3565,j
	Arene=1
	endif
 	555 continue

	if (Languedoc==1) then
	goto 558
	endif	
	if ((rfront(j))>3961) then
 	WRITE(6,*) 3961,j
	Languedoc=1
	endif
 	558 continue

	if (Catalonia==1) then
	goto 561
	endif	
	if ((rfront(j))>4328) then
 	WRITE(6,*) 4328,j
	Catalonia=1
	endif
 	561 continue

	if (Valencia==1) then
	goto 566
	endif	
	if ((rfront(j))>4670) then
 	WRITE(6,*) 4670,j
	Valencia=1
	endif
 	566 continue

	if (Andalusia==1) then
	goto 577
	endif	
	if ((rfront(j))>5074) then
 	WRITE(6,*) 5074,j
	Andalusia=1
	endif
 	577 continue
	
	if (SPortugal==1) then
	goto 580
	endif	
	if ((rfront(j))>5651) then
 	WRITE(6,*) 5651,j
	SPortugal=1
	endif
 	580 continue

	if (CPortugal==1) then
	goto 588
	endif	
	if ((rfront(j))>5945) then
 	WRITE(6,*) 5945,j
	CPortugal=1
	endif
 	588 continue

	CLOSE(6)


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

print '(" You may now import the file FRONT.dat and the spread rate in ")'
print '(" km/generation (not km/yr) is the slope FOR THE LAST 50 GENERATIONS ")'
print '("  ")'
print '(" You may know the simulation time by comparing those of NEOLPROF0.DAT i FRONT.DAT ")'
print '("  ")'


end program