Revision d0a9ad7ca4459a1cc221b7bf1d1d311733400f0a authored by jaros007 on 15 March 2021, 11:09:16 UTC, committed by GitHub on 15 March 2021, 11:09:16 UTC
1 parent 2cf6940
dilution_var_gamma_alpha.f90
PROGRAM many_resources ! population-based consumer-resource model, dimension-dependent
IMPLICIT NONE ! exponents in alpha constraints, explicit resource dynamics
! with repletion of resources and dilution of populations
DOUBLE PRECISION, ALLOCATABLE :: X(:,:),C(:),X0(:),BR(:),RC(:),GA(:)
DOUBLE PRECISION, ALLOCATABLE :: MC(:),S(:),PA(:),PAN(:),DX(:)
LOGICAL, ALLOCATABLE :: LV(:)
INTEGER N,K,K1,K2,K3,I,IMAX,I0,IL,IT,IW,NH,NC,J1
DOUBLE PRECISION T,DT,TMAX,TMEAS,RAN3,EPS,MUT,BRT,W,TVF,RR,R,SPEED
DOUBLE PRECISION DTM,TV,aa,comp, DIST,XP,GW,GS,SA,POPMIN,SQ3,DCMEAS
DOUBLE PRECISION dmin,plotcount,A1,A2,A3,AE,SQ2,SQ6,SQ15,PIN,DISPCH
DOUBLE PRECISION DTSEED, DTMERGE,TSEED,TMERGE,PT,BRW,DCC,DCMIN,MU,MK
DOUBLE PRECISION PTOT,KK,RM,CTOT,CMIN,CIN
N=3 ! SPACE DIMENSIONALITY
IMAX=500 ! MAX NUMBER OF POPULATIONS
EPS=0.005D0 ! SPREAD AROUND THE INITIAL POSITION
MUT=0.005D0 ! MUTATION AMPLITUDE
MU=0.1D0 ! FRACTION OF ANCESTOR THAT GOES INTO THE OFFSPRING
PIN=1.D-3 ! INITIAL AND POST-DILUTION POPULATION
KK=1.D0 ! PARAMETER IN MONOD FUNCTION
CMIN=1.D-8 ! MINIMAL RESOURCE BEFORE REPLETION
CIN=0.1d0 ! CONCENTRATION OF EACH RESOURCE AFTER REPLETION
DCMIN=MUT ! MERGING DISTANCE
DCMEAS=MUT ! MERGING DISTANCE
I0=1 ! INITIAL NUMBER OF POPULATIONS
POPMIN=1.D-6 ! D-5
DT=1.D-2 ! TIME STEP
TMAX=4000000.D0 ! END TIME
DTSEED=1.D0 ! TIME TO SEED A NEW CLUSTER
DTMERGE=100.D0*DTSEED !TIME TO MERGE CLUSTERS IN SEEDING PER CLUSTER
DTM=2000.D0 ! MEASURE TIME
ALLOCATE(X(N,IMAX),C(N),BR(IMAX),LV(IMAX),S(N))
ALLOCATE(X0(N),RC(N),MC(IMAX),PA(IMAX),DX(N),GA(N))
! GA=(/1.1D0,1.1D0,1.1D0/) ! EXPONENTS IN THE CONTRAINT
GA=(/.9D0,.9D0,.9D0/) ! EXPONENTS IN THE CONSTRAINT
open(unit=12,file='rand.dat') ! READING THE RANDOM SEED
read(12,*)iw
! iw=-87615
close(12)
SQ2=DSQRT(2.D0)
SQ6=DSQRT(6.D0)
SQ15=DSQRT(1.5D0)
SQ3=DSQRT(3.D0)
X0=0.D0
LV=.FALSE.
S=1.D0 ! SETTING THE RESOURCES
MK=0.1D0!S(1)/4 ! THE DECAY COEFFICIENT FOR RESOURCES
X0=(/0.9D0,0.05D0,0.05D0/)
I=I0
DO K=1,3
X(K,I)=X0(K)
END DO
LV(I)=.TRUE.
PA(I)=PIN
IT=I0
DO I=1,I0 ! PUTTING EVERYBODY ON CURVED SIMPLEX
A1=0.D0
DO K=1,N
A1=A1+X(K,I)**GA(K)
END DO
DO K=1,N
X(K,I)=X(K,I)/A1**(1.D0/GA(K))
END DO
END DO
!!$ open(unit=49,file='initial_config.dat')
!!$ DO K=1,IT
!!$ IF(LV(K)) WRITE(49,*) (X(2,K)-X(1,K))/SQ2,X(3,K)*SQ15
!!$ END DO
!!$ close(49)
!!$it=4000 ! reading the previous output, it should be adjusted
!!$ open(unit=35, file='points.dat') !reading a snapshot
!!$! read(35,*)t
!!$ DO K=1,IT
!!$ LV(K)=.TRUE.
!!$ read (35,*) (X(K1,K),K1=1,N)
!!$ END DO
!!$ close (35)
IL=IT ! TOTAL NUMBER OF INDIVIDUALS=NUMBER OF LIVE INDIVIDUALS
T=0.D0
TMEAS=DTM
TSEED=0.D0
TMERGE=DTMERGE-DTSEED
plotcount=0
C=1.D0 ! INITIAL RESOURCES
OPEN(UNIT=10,FILE='fvar_pop.dat')
OPEN(UNIT=40,FILE='total_var.dat')
DO WHILE (T.LE.TMAX.AND.IT.LT.IMAX-1) !!!! BEGINNING OF UPDATE !!!!
BR=0.D0 ! BIRTH RATES
BRW=0.D0
DO K=1,N ! RESOURCE CONSUMPTION
RC(K)=C(K)/(C(K)+KK)
END DO
DO I=1,IT ! BIRTH RATE
IF(LV(I)) THEN
DO K=1,N
BR(I)=BR(I)+X(K,I)*RC(K)
END DO
BRW=BRW+BR(I)*PA(I)
END IF
END DO
IF (T.GE.TSEED) THEN ! SPLITTING EVENT
TSEED=T+DTSEED
RR=ran3(iw)
GS=0.D0
K=0
DO WHILE (GS.LE.RR) ! CHOOSING THE FATHER K
K=K+1 ! PROPORTIONAL TO ITS BIRTH RATE
IF (LV(K)) GS=GS+BR(K)*PA(K)/BRW
END DO
I=1
DO WHILE (LV(I)) ! FINDING A SPARE SPOT I
I=I+1
END DO
IF(I.GT.IT) IT=I ! ADDING A NEWBORN TO THE RIGHT END OF THE LIST
IL=IL+1 ! INCREMENTING THE NUMBER OF LIVE CLUSTERS
LV(I)=.TRUE. ! ASSINGING A SEAT
PA(I)=PA(K)*mu ! SPLITTING THE ANCESTOR
PA(K)=PA(K)*(1.d0-mu)
BR(I)=BR(K) ! TEMPORARILY THE SAME BIRTH RATE
A2=0.D0
DO K1=1,N ! COMPUTING DX
DX(K1) = MUT*(2*ran3(iw)-1.D0)
X(K1,I)=X(K1,K)+DX(K1)
X(K1,I)=DMAX1(1.D-8,X(K1,I))
X(K1,I)=DMIN1(1.D0,X(K1,I))
A2=A2+X(K1,I)**GA(K1)
END DO
DO K2=1,N ! PUTTING X ON CONSTRAINT SURFACE
X(K2,I)=DMIN1(1.D0,X(K2,I)/A2**(1.D0/GA(K2)))
END DO
END IF
PTOT=0.D0
DO I=1,IT ! UPDATE OF POPULATION
IF(LV(I)) THEN
PA(I)=PA(I)+BR(I)*PA(I)*DT ! NEW POPULATION
PTOT=PTOT+PA(I)
IF (PA(I).LE.POPMIN) THEN
PA(I)=0.D0
! WRITE(*,*)'EXTINCTION'
LV(I)=.FALSE.
IL=IL-1
END IF
END IF
END DO
CTOT=0.D0
DO K=1,N ! UPDATE OF RESOURCES
RM=0.D0
DO I=1,IT
IF(LV(I)) RM=RM+PA(I)*X(K,I)
END DO
C(K)=C(K)-RM*RC(K)*DT
CTOT=CTOT+C(K)
END DO
IF (CTOT.LE.CMIN) THEN ! DILUTION AND REPLETING RESOURCES
C=CIN
PA=PIN*PA/PTOT
END IF
T=T+DT
IF(T.GE.TMERGE.OR.T.GT.TMAX) THEN ! MERGING CLOSE CLUSTERS
TMERGE=T+DTMERGE
DO K1=1,IT
IF(LV(K1)) THEN ! removing coinciding clusters
DO K2=1,K1-1
IF(LV(K2).AND.LV(K1)) THEN
DCC=0.D0
DO J1=1,N
DCC=DCC+(X(J1,K1)-X(J1,K2))**2
END DO
DCC=DSQRT(DCC)
IF (DCC.LE.DCMIN) THEN !LIQUIDATING CLUSTER KC
DO J1=1,N ! PRESERVING CENTRE OF MASS POSITION
X(J1,K2)=(X(J1,K1)*PA(K1)+X(J1,K2)*PA(K2))/(PA(K1)+PA(K2))
END DO
LV(K1)=.FALSE.
PA(K2)=PA(K2)+PA(K1)
PA(K1)=0.D0
IL=IL-1
END IF
END IF
END DO
END IF
END DO
END IF ! END OF MERGING
IF(T.GE.TMEAS) THEN ! MEASUREMENT
TMEAS=T+DTM
DO K1=1,IT !merging once again, now with larger cutoff
IF(LV(K1)) THEN ! removing coinciding clusters
DO K2=1,K1-1
IF(LV(K2).AND.LV(K1)) THEN
DCC=0.D0
DO J1=1,N
DCC=DCC+(X(J1,K1)-X(J1,K2))**2
END DO
DCC=DSQRT(DCC)
IF (DCC.LE.DCMEAS) THEN !LIQUIDATING CLUSTER KC
DO J1=1,N ! PRESERVING CENTRE OF MASS POSITION
X(J1,K2)=(X(J1,K1)*PA(K1)+X(J1,K2)*PA(K2))/(PA(K1)+PA(K2))
END DO
LV(K1)=.FALSE.
PA(K2)=PA(K2)+PA(K1)
PA(K1)=0.D0
IL=IL-1
END IF
END IF
END DO
END IF
END DO
WRITE(10,*) (plotcount*100,K1=1,3) ! WRITING TO THE MOVIE FILE
plotcount=plotcount+1
DO K=1,IT
IF(LV(K)) THEN
WRITE(10,*) (X(2,K)**GA(2)-X(1,K)**GA(1))/DSQRT(2.D0),X(3,K)**GA(3)*DSQRT(1.5D0),PA(K)
END IF
END DO
WRITE(40,*) T, IL
WRITE(*,*) T, IL
!!$ open(unit=75, file='final_points_ran.dat') !writing last snapshot,
!!$ DO K=1,IT
!!$ IF(LV(K)) WRITE(75,*) (X(2,K)-X(1,K))/DSQRT(2.D0),X(3,K)*DSQRT(1.5D0),PA(K)
!!$ END DO
!!$ close (75)
END IF
END DO
CLOSE(10)
CLOSE(40)
open(unit=12, file='rand.dat')
write(12,*)nint(-ran3(iw)*1.d+6)
close(12)
write(*,*) 'survivors, it, il', it, il
END PROGRAM many_resources
!==========================================================================
! ran3 from Numerical Recipes, 2nd edition
!--------------------------------------------------------------------------
function ran3(idum)
!==========================================================================
! implicit real*4(m)
integer idum
integer mbig,seed,mz
! real mbig,seed,mz
double precision ran3,fac
! parameter (mbig=4000000.,mseed=1618033.,mz=0.,fac=2.5e-7)
! parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1.e-9)
parameter (mbig=1000000000,mseed=161803398,mz=0,fac=1./mbig)
integer i,iff,ii,inext,inextp,k
integer mj,mk,ma(55)
! real mj,mk,ma(55)
save iff,inext,inextp,ma
data iff /0/
1 if(idum.lt.0.or.iff.eq.0)then
iff=1
mj=mseed-iabs(idum)
mj=mod(mj,mbig)
ma(55)=mj
mk=1
do 11 i=1,54
ii=mod(21*i,55)
ma(ii)=mk
mk=mj-mk
if(mk.lt.mz)mk=mk+mbig
mj=ma(ii)
11 continue
do 13 k=1,4
do 12 i=1,55
ma(i)=ma(i)-ma(1+mod(i+30,55))
if(ma(i).lt.mz)ma(i)=ma(i)+mbig
12 continue
13 continue
inext=0
inextp=31
idum=1
endif
inext=inext+1
if(inext.eq.56)inext=1
inextp=inextp+1
if(inextp.eq.56)inextp=1
mj=ma(inext)-ma(inextp)
if(mj.lt.mz)mj=mj+mbig
ma(inext)=mj
ran3=mj*fac
if(ran3.le.0.or.ran3.ge.1) goto 1
return
end function ran3
Computing file changes ...