! ...\Projects\2021-30\2021_Fort_coast__2_species\06__costa_salts_llargs\salts_llargs.f90
!
integer::i,j,y,t,z,x,n
real,dimension (2999,2999)::k,h,a,b
real,dimension (2999)::rfront
real dis,pe,rr,am,m,rp,bm,mp,f,ga

!Final number of generations:
t=100

!Cultural tranmission parameters:
f=0
ga=1

! distance between nodes = length of inland jumps:
dis=50

! length of coastal jumps (it must be an integer multiple of 50 km) divided by 50 km:
n=100/50

! rr is the Neolithic growh rate a=LOG(Ro) [in fortran LOG is the Napierian logarithm], value from Isern et al., J Stat Mechs: Theory & Exp. (2008):
rr=LOG(2.45)	  

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

! maximum Neolithic population density (Currant+Excoffier 2005 + Fort et al NJP 2008):
am=1.28
! maximum Mesolithic population density (Currant+Excoffier 2005 + Fort et al NJP 2008):
bm=0.064

! persistence:
pe=0.38

! Initial Neolithic population density:
a=0
do z=1,50
a(1,z)=1.28
end do

b=0

! Initial Mesolithic population density:
k=0.064
!

h=0
m=0

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

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

! I delete the files with the front position along the VERTICAL (MAX) and DIAGONAL (MAXDIAG) directions:
OPEN (2,FILE='MAX.DAT',STATUS='UNKNOWN')
WRITE (2,*) ' '
CLOSE (2)

OPEN (2,FILE='MAXDIAG.DAT',STATUS='UNKNOWN')
WRITE (2,*) ' '
CLOSE (2)
!================================================================================
! 1) 
do j=1,t
!================================================================================
! 2) Logistic Neolithic reproduction:
	do z=1,2999
	do x=1,2999
	!we save a lot of time by neglecting nodes with neglegible population:
		if (a(x,z)<1E-10) then
		else
			m=am+a(x,z)*(exp(rr)-1)
			a(x,z)=(a(x,z)*am*exp(rr))/m
			m=0
		end if	
	end do
	end do
!================================================================================
! 3) Logistic Mesolithic reproduction:
	do z=1,2999
	do x=1,2999
	!we save a lot of time by neglecting nodes with neglegible population:
		if (k(x,z)<1E-10) then
		else
			mp=bm+k(x,z)*(exp(rp)-1)
			k(x,z)=(k(x,z)*bm*exp(rp))/mp
			mp=0
		end if	
	end do
	end do
!=================================================================================
! 4) Cultural transmission:
	h=k
	do z=1,2999
	do x=1,2999

		!we save a lot of time by neglecting nodes with neglegible population:
		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
!================================================================================
! 5) Neolithic dispersal: 

!   Inland (i=2,2999, not i=1,2999 beacuse the coast is i=1):
	do y=1,2999
	do i=2,2999
!-----------------------
		!we save a lot of time by neglecting nodes with neglegible population:
		if (a(i,y)<1E-10) then
		goto 100
		endif
!-----------------------
		do z=-1,1
		do x=-1,1	

!   a(i,y) is the Neolithic density at the origin point of jumps:
		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

!*********************************************************************************
!   The coast is the straight line (edge of the grid) x=1:
	do y=2,2999
!-----------------------
		!we save a lot of time by neglecting nodes with neglegible population:
		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

! From coastal points (nodes) there are jumps to 3 points (not 4 as above line 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
! I do not include the node (x=-1.and.z==0) because (1+x,y+z)=(0,y)  is on the sea.
		
		120 continue
		end do
		end do
	150	end do	
		
!*********************************************************************************

! We keep the computed Neolithic population in matrix a, which will be fixed in the next iteration:
a=b

! We delete the matrix b, which will change in the next iteration:
b=0

!===============================================================================
! 6) Mesolithic dispersal:
! These 2 lines are explained below:
h=k
k=0

! The coast is (i,y) with i=1: 
	do y=1,2999
	do i=2,2999
		!we save a lot of time by neglecting almost saturated nodes:
		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

! Computing as for the Neolithic population would be too slow (see the comment above).
! So here we compute directly at the FINAL node of the jump =k(i,y),
! This is why above I have written h=k i k=0.
		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
	
!*********************************************************************************
!   COAST: (x,y) with x=1:
	do y=1,2999
!-----------------------
! As before line 180, we do not compute at saturated nodes on the coast:
		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

		! Forward and backward jumps:
		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
! I do not include the point (x=-1.and.z==0) because (1+x,y+z)=(0,y) is on the sea.

		220 continue
		end do
		end do		
	250	end do

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

k=h

!===============================================================================
! 7) We keep some results before jumping to the next generation:
!!------------------
! 7a) Front position in the VERTICAL direction as a function of time:  
! The slope is the speed in km/generation 
OPEN (2,ACCESS='APPEND',FILE='MAX.DAT',STATUS='UNKNOWN')
! I find the frist node with rho<0.1, and find the node rho=0.1 (by computing a straight line):
	do zz=2999,50,-1						    
		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)))
		exit
	end if
	400 continue
	end do
IF (.not.(j<1)) WRITE (2,*)  j,rfront(j)
close(2)
!
! 7b) Front profiles at a given generation:	

! e.g. generation 90, Neolithic front:
	if (j==90) then
	OPEN (4,FILE='MAXPERFIL.DAT',STATUS='UNKNOWN')
   	do zz=1,2999
	WRITE(4,*) zz,a(1,zz)
	end do
	CLOSE(4)
	end if

! e.g. generation 90, Mesolithic front:
	if (j==90) then
	OPEN (5,FILE='MAXPERPAL.DAT',STATUS='UNKNOWN')
  	do zz=1,2999
	WRITE(5,*) zz,k(1,zz)
	end do
	CLOSE(5)
	end if

!===============================================================================
! We show the time (generation) on the screen, and jump to the next one:
print *,j
end do

print '(" You can import MAX.dat and the front speed is ")'
print '(" the slope in km/generation (not m/yr) ")'

end program
