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