! creates a block with a specified map length, trait value, number of loci, and genic variance per unit map length. program genome integer(KIND=8) logpieces,i,j,pieces real(kind=8) y0,z0,V0,temp1,temp2,x3,x4,delta real(KIND=8),allocatable,dimension(:):: eff1,eff2 REAL(KIND=8), PARAMETER :: Pi = 3.1415926535897932 y0=0.4D+0 ! map length of introduced block z0=0.0D+0 ! total trait value of introduced block V0=0.1d+0 ! genic variance per unit map length logpieces=11 pieces=2**logpieces ! total number of loci on block allocate(eff1(pieces),eff2(pieces)) call init_random_seed ! choose contributions of loci such that they sum to z0. eff1(1)=z0 do i=1,logpieces print*, i temp2=sqrt(V0*0.25D+0*(y0/(2**(i-1)))) do j=1,2**(i-1) print*,j temp1=eff1(j)*0.5D+0 call random_number(x3) call random_number(x4) eff2(2*j-1)=(sqrt(-2.0D+0*log(x3))*cos(2.0D+0*Pi*x4)*temp2)+temp1 eff2(2*j)=eff1(j)-eff2(2*j-1) enddo do j=1,2**i eff1(j)=eff2(j) enddo enddo ! output file: column 1--> locus number i; column 2--> effect of locus i; column 3--> cumulative effect of all loci from [0 to i]. open(4,file='genome_L2048_z0_V0y0_0.04.txt') temp1=eff1(1) do i=1,pieces write(4,93) i,eff1(i),temp1 93 Format(I8,1X,F12.9,1X,F15.9) temp1=temp1+eff1(i+1) enddo close(4) end SUBROUTINE init_random_seed() INTEGER :: z5, m, clock INTEGER, DIMENSION(:), ALLOCATABLE :: seed, seedold CALL RANDOM_SEED(size = m) ALLOCATE(seed(m), seedold(m)) CALL SYSTEM_CLOCK(COUNT=clock) seed = clock + 37 * (/ (z5 - 1, z5 = 1, m) /) CALL RANDOM_SEED(PUT = seed) CALL RANDOM_SEED(GET = seedold) DEALLOCATE(seed) END SUBROUTINE init_random_seed