! ...\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