C* *******   MUTATION AND RECOMBINATION SIMULATION *****
C           FIVE LOCUS BABY MODEL - MUTATION, SELECTION
c			AND RECOMBINATION 
c			TWO POPULATIONS ONE WITH DIFFERENT RATES OF RECOMBINATION
C			Each gene can be of one of three states, alleles, 0, 1, 2
c			1 is selectively neutral, 0 is deleterious, 2 is favored
c			The relative fitness of a five locus genotype is proportional to the sum
c			of the fitnesses of the five 
c			The relative fitness of a five locus genotype is proportional to the sum
c			of the fitnesses of the five 
C			Fit (I,J,K,L,M) = C+ (1-C)(I+J+K+L+M)/10
C			Recombination - choose recipient, chose donor, chose gene
c			Only a single gene can be moved at a time
c			Rate constant of recombination, X 
c			Probability of recombinatin X*NT*DT per time interval
C			At each time step there is a probability, PB*DT of a bottlneck.   If the
c			random number is less than the bottle neck probability, the bottle neck occurs.
c			For that bottle neck we go through all of the genotypes and any with frequencies
c			less than the bottle neck fraction, b (0 < b <=1) are set equal to 0.  
c			If in the bottle neck the population goes from NT to Nb, b=(Nt-Nb)/Nt. For example if
c			NT=1E8, and Nb = 1E2 then genotypes with frequencies < 1e-2  would be set to 0.  
c			Once all the genotypes go through the bottle neck, the frequencies are
c			normalized by dividing them by the sum of the frequencies. 
C		MODEL 4 - Initiate by choosing NC clones at random, 0<NC<20.
c		In this incarnation of the model, the two competing popuationas may have different
c		genotype distributions for the start.
C	Parameters
C	RANSEED - Integer - random number seed 
C	TMAX - Maximum number of time units
c	There are two clones, one that does and one that does not recombine,
c	1 and 2 respectively. The last index N is the clone, N=1 or N-2.
C	C Fraction of fitness constant  (0 <=C <=1)
c	Fitness of a clone of genotype I,J,K,L,M
c	F(I,J,K,L,M,1) = W1*(C+ (1-C)(I+J+K+L+M)/10)
C	F(I,J,K,L,M,2) = W1*(C+ (1-C)(I+J+K+L+M)/10)
C 	Where W1 and W2 are the relative fitnesses of the clones
c	P(I,J,K,L,M,N)' = P(I,J,K,L,M)* WN*(C+ (1-C)(I+J+K+L+M)/10)/NORM
C	NORM = is the average fitness - weighted average over all genotypes
c	For both clones.  Within a clone
c	Highest fitness 3,3,3,3,3, lowest, 1,1,1,1,1

C	Parameters
C	RANSEED - Integer - random number seed 
C	IGEN - Maximum number of generations
c	IPRT - PRINT INTERVAL
C	IRAN - RANDOM NUMBER RUNNER -
c	DT - STEP SIZE
C	XMU - MUTATION RATE per gene
C	X1 FRACTION OF MUTANTS 1, X2 FRACTION OF MUTANTS 2
C	C - FITNESS PARAMETER, C=1.0 no selection
c	X- RECOMBINATION PARAMETER, X*NT PROB. REPLACEMENT BY RECOMBINATION
C	NT -  TOTAL POPULATION SIZE Clone 1 and clone 2
C	IRNS - NUMBER OF INDEPENDENT RUNS
c	PB, B - Probability of a bottle neck and fraction surviving
C   XPX(I,J,K,L,M,N) -Frequency of genotpe
C 	In this verson of the model, N is the clone, 1 (sexual) 2(asexual)
C	The role of this model is to calculate how long it takes for loss of a population
c	To facilitate these simulations, when DT*N(I)*N(I)*x > 1, 
c	DT*N(I)*N(I)*x 	recombinants are generated at random. No decision is made
c	for recombination to occur.
C	Fitness determined as an exponential, 
c		F(I,J,K,L,M)= c +(1-c)(I^E+J^E+K^E+L*E+M^E)/(5*2^E)
   		PROGRAM CXG
		EXTERNAL RANDOM 
		INTEGER JX,JY,IRUNS,JQ,RUNTYPE
		DIMENSION XPX(4,4,4,4,4,3),F(4,4,4,4,4,3),X(3),W(3)
		DIMENSION NT(3),WBAR(3),WB(3,50,50),NTP(3,50,50)
        REAL RANDOM, MU,NT,NTT,NMIN
		INTEGER RANSEED,XNT,L,M
        COMMON /RANDSUB/ RANSEED
        RANSEED = MAX (RANSEED, 1)
        RANSEED = MIN (RANSEED,2796202)
		OPEN (2, FILE = 'CXGIN.TXT')
		OPEN (1, FILE = 'CXGOUT.TXT')
		READ(2,*) RANSEED,IRAN,IGEN,IPRT,IRNS,NC,RUNTYPE,NMIN,E
		READ(2,*) XN1,XN2,XMU,X(1),X(2),XM1,XM2,DT,C,PB,B,W(1),W(2)
		WRITE (1,*) '##### MUTATION RECOMBINATION AND SELECTION #####'
		WRITE (1,*) ' Two competing populations, with different xs'
		WRITE (1,*) ' RANDOM SELECTION OF NC CLONES TO START'
		WRITE (1,*) '  --- 1 AND 2 HAVE THE SAME INITIAL VARIABILITY'
		WRITE (1,*) ' FITNESS A FUNCTION OF THE SUM OF E-TH POWER'
		WRITE(1,*)  ' IN THIS SIMULATION 1 AND 2 CAN HAVE DIFFERENT'
		WRITE (1,*) '    INITIAL  GENOTYPE DISTRIBUTIONS'
		WRITE (1,*) ' '
		WRITE (1,*) '-- PARAMETERS AND INITIAL VALUES OF VARIABLES--'
		WRITE(1,*) 'SEED  IRAN IGEN  IPRT  RUNS  NC RUNTYPE NMIN  E'
		WRITE(1,10) RANSEED,IRAN,IGEN,IPRT,IRNS,NC,RUNTYPE,NMIN,E
		
10		FORMAT (7(I7,2X),1E8.3,2X,F6.2)
	 	WRITE (1,*) 'NT1  NT2  XMU X1 X2  XM1  XM2  DT C PB B W1 W2'
	 	WRITE (1,20) XN1,XN2,XMU,X(1),X(2),XM1,XM2,DT,C,PB,B,W(1),W
	 	WRITE (1,*) ' '
	 	WRITE (1,*) '   ---- OUTPUT -----'
	 	WRITE (1,*) 'TERM-TIME        N1       N2     WBAR1    WBAR2 '	 	
20		FORMAT (5(E8.2,2X),4(F8.3),2X,2(E8.3,2X),2(F6.2,2X))
		PRINT *, 'GOT THE DATA AND NOW TO DO MY THING'
C 	Run random numbers
		DO 2 I=1,IRAN
2		RAN=RANDOM()
		XM12=XM1+XM2
		NTT=XN1+XN2					


C		START RUN LOOP
		XP=IPRT
		IPX=XP/DT
		DO 550 IRN=1,IRNS
		NT(1)=XN1
		NT(2)=XN2
		IF (NT(1) .LT. NT(2)) DEN=NT(1)
		IF (NT(2) .LE. NT(1)) DEN=NT(2)
C	Calculate fitness of all genotypes and set initial frequencies to 0
		DO 100 I=1,3
		DO 100 J=1,3
		DO 100 K=1,3
		DO 100 L=1,3
		DO 100 M=1,3
		DO 100 N=1,2
		DN=5*(2**E)
		F(I,J,K,L,M,N) =W(N)*(C+((1-C)*(I**E+J**E+K**E+L**E+M**E))/DN)
		XPX(I,J,K,L,M,N)=0.0
100     CONTINUE
c	Set initial frequencies 
		IF (RUNTYPE .GE. 1) GO TO 25
		XPX(2,2,2,2,2,1) = 1.0
		XPX(2,2,2,2,2,2) = 1.0
		WBAR(1) =1.0
		WBAR(2) = 1.0
		GO TO 30	
C		NC clones will be selected at random 
c	 	Separately for 1 and 2
25		Q1=0.333333333333
		Q2=0.666666666667
c 		Select initial distribution for 1 and 2
		XNSUM=0
		DO 930 N=1,2
		DO 900 IN=1,NC
		RAN=RANDOM()
		IF (RAN .LE. Q1) I=1
		IF (RAN .GT. Q1 .AND. RAN .LE. Q2) I=2
		IF (RAN .GT. Q2) I=3
		RAN=RANDOM()
		IF (RAN .LE. Q1) J=1
		IF (RAN .GT. Q1 .AND. RAN .LE. Q2) J=2
		IF (RAN .GT. Q2) J=3
		RAN=RANDOM()
		IF (RAN .LE. Q1) K=1
		IF (RAN .GT. Q1 .AND .RAN .LE. Q2) K=2
		IF (RAN .GT. Q2) K=3
		RAN=RANDOM()
		IF (RAN .LE. Q1) L=1
		IF (RAN .GT. Q1 .AND. RAN .LE. Q2) L=2
		IF (RAN .GT. Q2) L=3
		RAN=RANDOM()
		IF (RAN .LE. Q1) M=1
		IF (RAN .GT. Q1 .AND. RAN .LE. Q2) M=2
		IF (RAN .GT. Q2) M=3
		RAN=RANDOM()
		XPX(I,J,K,L,M,N)=RAN +XPX(I,J,K,L,M,N)
		XNSUM=XNSUM+XPX(I,J,K,L,M,N)
C		WRITE (1,925) IRN,I,J,K,L,M, N, XPX(I,J,K,L,M,1),XNSUM	
925		FORMAT( I4,2X,5I4,2X,I4,F8.5,2X,F8.5)
900		CONTINUE
C		NORMALIZE AND CALACULATE INITIAL MEAN FITNESS
		WBARX=0
		SUM=0
		DO 910 I=1,3
		DO 910 J=1,3
		DO 910 K=1,3
		DO 910 L=1,3
		DO 910 M=1,3
		TX=XPX(I,J,K,L,M,N)
		IF (TX .GT. 0 .AND. TX .LT. 1/DEN) XPX(I,J,K,L,M,N) = 1/DEN
		XPX(I,J,K,L,M,N) =XPX(I,J,K,L,M,N)/XNSUM
		XPX(I,J,K,L,M,N)= XPX(I,J,K,L,M,N)		
		WBAR(N)=WBAR(N)+XPX(I,J,K,L,M,N)*F(I,J,K,L,M,N)	

910		CONTINUE
		IXGEN=0
930		CONTINUE
c		WRITE (1,*) 'STARTING RUN ', IRN	
C
c 		WRITE (1,*) 'RUN ',IRN,IXGEN,NT(1),NT(2),WBAR(1),WBAR(2)
30		CONTINUE
C ******* Start evolution 
		IP=0
		IPTP=0
		IPRC=0
		IPX=IPRT/DT
		IG=0
		IGPT=IGEN/IPRT
		NT(1)=XN1
		NT(2)=XN2
		NTT=NT(1)+NT(2)
		DO 500 II=1, IGEN/DT
		IP=IP+1
		IPTP=IPTP+1
		IPRC=IPRC+1

C	Mutation loop - decides if a mutant is generated at II=II+DT
		DO 400 JN=1,2
		DEN=1/NT(JN)
		PMUT=NT(JN)*XMU*DT
		RAN=RANDOM()
		IF (RAN .LE. PMUT) GO TO 110
		GO TO 400
C 	CHOOSE CLONE FOR MUTATION
110		CP=0
		RAN=RANDOM()
		DO 120 I=1,3
		DO 120 J=1,3
		DO 120 K=1,3
		DO 120 L=1,3
		DO 120 M=1,3
		CP=CP+XPX(I,J,K,L,M,JN)
		IF (RAN .LE. CP) GO TO 130

120		CONTINUE

C	CHOOSE LOCUS FOR MUTATION
130		CM=0
		JX=0
c		WRITE (1,*) 'MUTANT CLONE', I,J,K,L,M,JN 		
		RAN=RANDOM()
		DO 140 JJ=1,5
		CM=CM+0.2
		JX=JX+1
		IF (RAN .LE.CM) GO TO 150
140 	CONTINUE
C	Chose allele of mutation
150		CONTINUE
	
		JY=0
		RAN=RANDOM()
		IF (RAN .LT. XM1) JY=1
		IF (RAN .GT. XM1 .AND. RAN .LT. XM1+XM2) JY=2
		IF (RAN .GE. XM1+XM2) JY=3

C	INCREMENT MUTANT GENOTYPE BY 1/NT
		IF (JX .EQ. 1)  XPX(JY,J,K,L,M,JN) = XPX(JY,J,K,L,M,JN)+1/NT(JN)
		IF (JX .EQ. 2 ) XPX(I,JY,K,L,M,JN) = XPX(I,JY,K,L,M,JN)+1/NT(JN)
		IF (JX .EQ. 3 ) XPX(I,J,JY,L,M,JN) = XPX(I,J,JY,L,M,JN)+1/NT(JN)
		IF (JX .EQ. 4 ) XPX(I,J,K,JY,M,JN) = XPX(I,J,K,JY,M,JN)+1/NT(JN)
		IF (JX .EQ. 5 ) XPX(I,J,K,L,JY,JN) = XPX(I,J,K,L,JY,JN)+1/NT(JN)
		
C		REDUCE THE FREQUENCY OF THE CHOSEN RECIPIENT GENOTYPE FREQUENCY BY 1/N
		XPX(I,J,K,L,M,JN) = XPX(I,J,K,L,M,JN)-1/NT(JN)
		TX=XPX(I,J,K,L,M,JN)
		IF (TX .LE.0 ) XPX(I,J,K,L,M,JN)=0

C 		*** MUTATION TEST LOOP	
c		WRITE (1,*) 'MUT LOCUS AND ALLELE', JX,JY
c		WRITE (1,50) I,J,K,L,M,XPX(I,J,K,L,M,1),XPX(I,J,K,L,M,2) 
		
400		CONTINUE
C		WRITE (1,*) 'MUTANT CLONE JX', I,J,K,L,M,JN,JX


C	******* Recombination Loop ***********
		DO 800 JQ=1,2
		RCK=NT(JQ)*NT(JQ)*X(JQ)*DT
		IF (RCK .GT.1.0) GO TO 705 
C      Random Determinaion if recombination occurs
	
		PREC=NT(JQ)*NT(JQ)*X(JQ)*DT		
		RAN=RANDOM()
		IF (RAN .GT. PREC) GO TO 800
C 	CHOOSE RECIPIENT CLONE
705 	IREC=RCK
		IF (IREC .LE. 1) IREC=1
		DO 707 IRC=1,IREC

709		CONTINUE
		CP=0
		RAN=RANDOM()
		DO 720 I=1,3
		DO 720 J=1,3
		DO 720 K=1,3
		DO 720 L=1,3
		DO 720 M=1,3
		CP=CP+XPX(I,J,K,L,M,JQ)
		IF (RAN .LE. CP) GO TO 730
720		CONTINUE
730		CONTINUE
C		WRITE(1,*) 'RECIPIENT      ', I,J,K,L,M,JQ

C	Choose Donor Clone				
		CP=0
		RAN=RANDOM()
		DO 740 ID=1,3
		DO 740 JD=1,3
		DO 740 KD=1,3
		DO 740 LD=1,3
		DO 740 MD=1,3
		CP=CP+XPX(ID,JD,KD,LD,MD,JQ)
		IF (RAN .LE. CP) GO TO 750
740		CONTINUE	

C		Chose Locus for recombination	
750		CM=0
C	    WRITE (1,*) 'DONOR          ' , ID,JD,KD,LD,MD,JQ
		JX=0
		RAN=RANDOM()
		DO 755 JJ=1,5
		CM=CM+0.2
		JX=JX+1
		IF (RAN .LE.CM) GO TO 760
755		CONTINUE
760 	CONTINUE
		
c		REPLACE ALLELE
		IF (JX .EQ. 1) XPX(ID,J,K,L,M,JQ)=XPX(ID,J,K,L,M,JQ)+ 1/NT(JQ)
		IF (JX .EQ. 2) XPX(I,JD,K,L,M,JQ)=XPX(I,JD,K,L,M,JQ)+ 1/NT(JQ)
		IF (JX .EQ. 3) XPX(I,J,KD,L,M,JQ)=XPX(I,J,KD,L,M,JQ)+ 1/NT(JQ)				
		IF (JX .EQ. 4) XPX(I,J,K,LD,M,JQ)=XPX(I,J,K,LD,M,JQ)+ 1/NT(JQ)
		IF (JX .EQ. 5) XPX(I,J,K,L,MD,JQ)=XPX(I,J,K,L,MD,JQ)+ 1/NT(JQ) 

C		Reduce recipient clone frequency
		XPX(I,J,K,L,M,JQ) = XPX(I,J,K,L,M,JQ)-1/NT(JQ)
		IF (XPX(I,J,K,L,M,JQ) .LE.0) XPX(I,J,K,L,M,JQ)=0
** End of recombination loop
C		WRITE(1,*) 'RECOMBINANT    ', I,J,K,L,M,JQ
707		CONTINUE		
800 	CONTINUE
c 	CALCULATE MEAN FITNESS
		WBAR(1)=0
		WBAR(2)=0
		DO 300 I=1, 3
		DO 300 J=1, 3
		DO 300 K=1, 3
		DO 300 L=1, 3
		DO 300 M=1, 3
		DO 300 N=1,2
		WBAR(N) = WBAR(N) + XPX(I,J,K,L,M,N)*F(I,J,K,L,M,N)
300		CONTINUE

C		CALACULATE POST SELECTION WITHIN-POPULATION FREQUENCES T+DT
	
		DO 310 I=1, 3
		DO 310 J=1, 3
		DO 310 K=1, 3
		DO 310 L=1, 3
		DO 310 M=1, 3
		DO 310 N=1,2
		P=XPX(I,J,K,L,M,N)
		XPX(I,J,K,L,M,N)=P+P*DT*(F(I,J,K,L,M,N)-WBAR(N))/WBAR(N)
310		CONTINUE
C ******** CALACULATE THE NEW VALUES FOR THE NTs
		P1=NT(1)/NTT
		P2=NT(2)/NTT
		WBCLONES= (P1*WBAR(1)+P2*WBAR(2))
		P1P=P1+P1*DT*(WBAR(1)-WBCLONES)/WBCLONES
		NT(1)=NTT*P1P
		NT(2)=NTT-NT(1)	
C		TEMPORARY PRINT LOOP
		IF (IPTP .LT. IPX) GO TO 842
		JGEN=II*DT
c		WRITE (1,*) 'GEN',JGEN,NT(1),NT(2),WBAR(1),WBAR(2),WBCLONES
		IPTP=0
842		CONTINUE
350		CONTINUE	

C  PRINT LOOP
		IF (IP .LT. IPX) GO TO 190
		IGX=IGX+1
		WB(1,IRN,IGX)=WBAR(1)
		WB(2,IRN,IGX)=WBAR(2)
		NTP(1,IRN,IGX)=NT(1)
		NTP(2,IRN,IGX)=NT(2)
190		CONTINUE
		IF (NT(1).LT. NMIN) GO TO 505
		IF (NT(2) .LT. NMIN) GO TO 505
500		CONTINUE
505		JGEN = II*DT
		WRITE (1,*) JGEN,NT(1),NT(2),WBAR(1),WBAR(2)		
C		Test print loop
c		DO 40 I=1,3
c		DO 40 J=1,3
c		DO 40 K=1,3
c		DO 40 L=1,3
c		DO 40 M=1,3
c		WRITE (1,50) I,J,K,L,M,XPX(I,J,K,L,M,1),XPX(I,J,K,L,M,2)
50		FORMAT (5I3,2X,F8.6,2X,F8.6)
40 		CONTINUE

550		CONTINUE

		PRINT *,'  MY JOB IS DONE, NOW DO YOURS'
	 	END 
	
C******************************************************************************
C     MODULE:   RANDOM
C     VERSION:  1.0
C     AUTHOR:   NEAL J. BOGDANOVICH
C     CREATED:  JANUARY 1987
C-----COPYRIGHT (C) 1987 BY NEAL J. BOGDANOVICH, AMHERST, MASSACHUSETTS
	FUNCTION RANDOM ()
C     RANDOM GENERATES A PSEUDO-RANDOM NUMBER.  THE PROCEDURE ASSUMES THAT
C     INTEGER ARITHMETIC UP TO 125 * 2796203 = 348525375 < 2 ** 29 IS
C     AVAILABLE.  THE START VALUE MUST BE WITHIN THE LIMITS 1 TO 2796202
C     INCLUSIVE, AND MUST NOT BE CHANGED BETWEEN SUCCESSIVE CALLS.
C
C     PIKE, M. C., AND HILL, I. D.  ALGORITHM 266, PSEUDO-RANDOM NUMBERS
C     COMM. ACM 8 (OCT. 1965), 605
C     HANSSON, L.  REMARK ON ALGORITHM 266
C
C     Z(I+1) = (AZ(I) + B)(MOD C)     I=0,1,2
C     R(I+1) = Z(I+1) / C
 
C     ARGUMENT DECLARATIONS
      REAL RANDOM
      INTEGER Y
      COMMON /RANDSUB/ Y
 
      Y = 125 * Y
      Y = Y - INT(Y / 2796203) * 2796203
      RANDOM = Y / 2796203.0
 
      RETURN
      END