use msimsl parameter (sita=0.2) parameter (nq=300,nm=300,nw=30) parameter (nm2=int(real(nm)*sita+0.5),nm1=nm-nm2) parameter (ngener=20,nrep=10000) parameter (nlmax=1) real gfreq(0:ngener,nlmax),gfreqq(0:ngener,nlmax),gfreqm(0:ngener,nlmax) real vfreq(0:ngener,nlmax),afreq(0:ngener,nlmax) real df(ngener,nlmax),het(0:ngener,nlmax),hetfreq(0:ngener,nlmax) real ane(nlmax),anel(ngener,nlmax),ane2(nlmax) real anef(nlmax),anefl(ngener,nlmax),anef2(nlmax) real f0(nq+nm,nq+nm),f1(nq+nm,nq+nm) real f_coef(ngener),acoef(ngener),dfped(ngener),anefped(ngener) integer queen(nq,nlmax,2),male(nm,nlmax),worker(nw,nlmax,2) integer queen0(nq,nlmax,2),male0(nm,nlmax) integer male1(nm1,nlmax),male2(nm2,nlmax) integer icount(ngener,nlmax) integer fpair(nq),mpair(nm) integer f_queen(nq),m_queen(nq),f_m1(nm1),w_m2(nm2) integer f_worker(nw),m_worker(nw) write(*,*) 'seed value of random number=' read(*,*) iseed call rnset(iseed) uwm2=real(nm2)/real(nw) write(*,*) 'model= (1:modelA, 2:modelB)' read(*,*) imodel if(imodel==2) then if(nw/=nm2) then write(*,*) 'check nm2 & nw (should be nm2=nw)' stop end if end if write(*,*) 'mating system= (1:polyandry-polygyny, 2:monoandry-manogyny)' read(*,*) mtype if(mtype==2) then if(nq/=nm) then write(*,*) 'check nq & nm (should be nq=nm)' stop end if end if write(6,1000) nq,nm,nw write(6,2000) sita,nm1,nm2,uwm2 if(mtype==1) then write(6,*) 'polyandry-polygyny model' else if(mtype==2) then write(6,*) 'monoandry-monogyny model' end if write(6,*) 'iseed=',iseed 1000 format('nq=',i4,1x,'nm=',i4,1x,'nw=',i4) 2000 format('sita=',f7.4,1x,'nm1=',i4,1x,'nm2=',i4,1x,'uwm2=',f8.4) do i=1,ngener do j=1,nlmax afreq(i,j)=0.0 vfreq(i,j)=0.0 df(i,j)=0.0 hetfreq(i,j)=0.0 icount(i,j)=0 end do end do do i=1,nlmax afreq(0,i)=0.5 vfreq(0,i)=0.0 hetfreq(0,i)=0.0 end do do i=1,ngener acoef(i)=0.0 end do !---- prediction of effective population size ! if(imodel==1) then amuwm2=real(nm2)/real(nw) predne=((2.0+amuwm2*sita)*real(nq)+(2.0-sita)**2*real(nm))/ & ((3.0-sita)**2*real(nm*nq)) predne=1.0/predne else if(imodel==2) then predne=(2.0*real(nq)+(2.0-sita)**2*real(nm))/ & ((3.0-sita)**2*real(nm*nq)) predne=1.0/predne end if !********************** do irep=1,nrep !********************** do i=1,nlmax gfreq(0,i)=0.5 gfreqq(0,i)=0.5 gfreqm(0,i)=0.5 het(0,i)=0.0 end do do i=1,ngener do j=1,nlmax gfreqq(i,j)=0.0 gfreqm(i,j)=0.0 het(i,j)=0.0 end do end do do i=1,nq+nm do j=1,nq+nm f0(i,j)=0.0 end do end do do i=1,nq f0(i,i)=0.5 end do do i=nq+1,nq+nm f0(i,i)=1.0 end do do i=1,ngener f_coef(i)=0.0 end do call base(nq,nm,queen,male,gfreqq,gfreqm,het,nlmax,ngener,iseed) do i=1,nlmax hetfreq(0,i)=hetfreq(0,i)+het(0,i) end do !-------------------------- do igener=1,ngener !-------------------------- ! if(mtype==1) then ! *** polyandry-polygyny *** !------- reproduction of young queens call mate1(nq,nm,queen,male,queen0,het,f_queen,m_queen, & ngener,nlmax,iseed,igener) !------- reproduction of workers call mate2(nq,nm,nw,queen,male,worker,f_worker,m_worker, & nlmax,iseed) else if(mtype==2) then ! *** monoandry-monogyny *** !------- reproduction of young queens call mate1m(nq,nm,queen,male,queen0,fpair,mpair,het, & f_queen,m_queen,ngener,nlmax,iseed,igener) !------- reproduction of workers call mate2m(nq,nm,nw,queen,male,worker,fpair,mpair, & f_worker,m_worker,nlmax,iseed) end if !------- reproduction of young males (m1) by queens call reprom1(nq,nm1,queen,male1,f_m1,nlmax,iseed) !------- reproduction of worker-produced males (m2) if(imodel==1) then call reprom2(nw,nm2,worker,male2,w_m2,nlmax,iseed) else if(imodel==2) then call reprom2b(nw,nm2,worker,male2,w_m2,nlmax,iseed) end if do i=1,nm1 do j=1,nlmax male0(i,j)=male1(i,j) end do end do do i=1,nm2 do j=1,nlmax male0(nm1+i,j)=male2(i,j) end do end do ! !--------------------------------------------------------------------- ! !------- computation of gene frequency in queens call freqq(queen0,gfreqq,ngener,nlmax,nq,igener) !------- computation of gene frequency in males call freqm(male0,gfreqm,ngener,nlmax,nm,igener) do i=1,nlmax gfreq(igener,i)=((2.0-sita)*gfreqq(igener,i)+ & gfreqm(igener,i))/(3.0-sita) if(gfreq(igener,i)==0.0.or.gfreq(igener,i)==1.0) then icount(igener,i)=icount(igener,i)+1 end if end do !--------------------------------------------------------------------- do i=1,nq do j=1,nlmax do k=1,2 queen(i,j,k)=queen0(i,j,k) end do end do end do do i=1,nm do j=1,nlmax male(i,j)=male0(i,j) end do end do !--------------------------------------------------------------------- do i=1,nlmax afreq(igener,i)=afreq(igener,i)+gfreq(igener,i) vfreq(igener,i)=vfreq(igener,i)+gfreq(igener,i)*gfreq(igener,i) end do !--------------------------------------------------------------------- do i=1,nlmax dft=(het(igener-1,i)-het(igener,i))/het(igener-1,i) df(igener,i)=df(igener,i)+dft hetfreq(igener,i)=hetfreq(igener,i)+het(igener,i) end do !--------------------------------------------------------------------- call coanc(nq,nw,nm,nm1,nm2,f0,f1,f_queen,m_queen,f_m1, & f_worker,m_worker,w_m2) call inbr(nq,nw,nm,f0,f_queen,m_queen,f_coef,ngener,igener) ! write(*,*) igener,f_coef(igener) acoef(igener)=acoef(igener)+f_coef(igener) do i=1,nq+nm do j=1,nq+nm f0(i,j)=f1(i,j) end do end do !-------------------- end do !-------------------- !********************** end do !********************** do i=1,ngener do j=1,nlmax afreq(i,j)=afreq(i,j)/real(nrep) vfreq(i,j)=vfreq(i,j)/real(nrep)-afreq(i,j)*afreq(i,j) end do end do do i=0,ngener do j=1,nlmax hetfreq(i,j)=hetfreq(i,j)/real(nrep) end do end do do i=1,ngener do j=1,nlmax df(i,j)=(hetfreq(i-1,j)-hetfreq(i,j))/hetfreq(i-1,j) if(i==2) then df(2,j)=(hetfreq(0,j)-hetfreq(2,j))/hetfreq(0,j) end if end do end do do i=1,ngener do j=1,nlmax anel(i,j)=2.0*(vfreq(i,j)-vfreq(i-1,j))/ & (gfreq(0,j)*(1.0-gfreq(0,j)-vfreq(i-1,j))) anel(i,j)=1.0/anel(i,j) anefl(i,j)=1.0/(2.0*df(i,j)) end do end do do i=1,nlmax ane(i)=0.0 anef(i)=0.0 end do do i=1,nlmax do j=1,ngener ane(i)=ane(i)+1.0/anel(j,i) anef(i)=anef(i)+1.0/anefl(j,i) end do ane(i)=ane(i)/real(ngener) ane(i)=1.0/ane(i) anef(i)=anef(i)/real(ngener) anef(i)=1.0/anef(i) end do do i=1,nlmax ane2(i)=0.0 end do do i=1,nlmax do j=6,ngener ane2(i)=ane2(i)+1.0/anel(j,i) anef2(i)=anef2(i)+1.0/anefl(j,i) end do ane2(i)=ane2(i)/real(ngener-5) ane2(i)=1.0/ane2(i) anef2(i)=anef2(i)/real(ngener-5) anef2(i)=1.0/anef2(i) end do do i=1,ngener acoef(i)=acoef(i)/real(nrep) end do avnef_ped=0.0 do i=2,ngener dfped(i)=(acoef(i)-acoef(i-1))/(1.0-acoef(i-1)) anefped(i)=1.0/(2.0*dfped(i)) avnef_ped=avnef_ped+1.0/anefped(i) end do avnef_ped=avnef_ped/real(ngener-1) avnef_ped=1.0/avnef_ped avnef_ped2=0.0 do i=6,ngener avnef_ped2=avnef_ped2+1.0/anefped(i) end do avnef_ped2=avnef_ped2/real(ngener-5) avnef_ped2=1.0/avnef_ped2 do i=1,nlmax write(6,3000) i,ane(i),ane2(i) end do 3000 format('locus=',i,1x,'Nev=',f10.4,1x,'Nev(6-20)=',f10.4) do i=1,nlmax write(6,3500) i,anef(i),anef2(i) end do 3500 format('locus=',i,1x,'Nei=',f10.4,1x,'Nei(6-20)=',f10.4) write(6,3700) avnef_ped,avnef_ped2 3700 format('Pedigree',1x,'Nei=',f10.4,1x,'Nei(6-20)=',f10.4) write(6,4000) predne 4000 format('pred_ne=',f10.4) do i=1,ngener write(6,5000) i,anel(i,1),anefl(i,1),hetfreq(i,1),icount(i,1) end do 5000 format(i3,2x,'Nev=',f10.4,2x,'Nei=',f10.4,2x,f10.4,2x,i5) do i=2,ngener write(6,5500) i,anefped(i),acoef(i) end do 5500 format(i3,2x,'Nei=',f10.4,2x,'f=',f10.4) stop end !********************************************************************* ! Subroutine for generating base poplulation * !********************************************************************* subroutine base(nq,nm,queen,male,gfreqq,gfreqm,het,nlmax,ngener,iseed) integer queen(nq,nlmax,2),male(nm,nlmax) integer ia(1000) real gfreqq(0:ngener,nlmax),gfreqm(0:ngener,nlmax) real het(0:ngener,nlmax) do i=1,nq do j=1,nlmax do k=1,2 queen(i,j,k)=2 end do end do end do do ii=1,nlmax do i=1,nq*2 ia(i)=i end do call random(ia,nq*2,iseed) nn=int(real(nq*2)*gfreqq(0,ii)+0.5) do i=1,nn if(ia(i)<=nq) then queen(ia(i),ii,1)=1 else queen(ia(i)-nq,ii,2)=1 end if end do end do do ii=1,nlmax do i=1,nq k=queen(i,ii,1)+queen(i,ii,2) if(k==3) then het(0,ii)=het(0,ii)+1.0 end if end do het(0,ii)=het(0,ii)/real(nq) end do do i=1,nm do j=1,nlmax male(i,j)=2 end do end do do ii=1,nlmax do i=1,nm ia(i)=i end do call random(ia,nm,iseed) nn=int(real(nm)*gfreqm(0,ii)+0.5) do i=1,nn male(ia(i),ii)=1 end do end do ! do i=1,nq ! write(6,*) i,(queen(i,j,1),j=1,nlmax) ! write(6,*) i,(queen(i,j,2),j=1,nlmax) ! end do ! do i=1,nm ! write(6,*) i,(male(i,j),j=1,nlmax) ! end do return end !********************************************************************* ! Subroutine for reproducing young queens (ployandry-polygyny) * !********************************************************************* subroutine mate1(nq,nm,queen,male,queen0,het,f_queen,m_queen, & ngener,nlmax,iseed,igener) integer queen(nq,nlmax,2),male(nm,nlmax) integer queen0(nq,nlmax,2) integer f_queen(nq),m_queen(nq) real het(0:ngener,nlmax) do i=1,nq do j=1,nlmax rq=rnunf() iiq=int(rq*real(nq))+1 kk=1 if(rnunf()>0.5) then kk=2 end if queen0(i,j,1)=queen(iiq,j,kk) rm=rnunf() iim=int(rm*real(nm))+1 queen0(i,j,2)=male(iim,j) end do f_queen(i)=iiq m_queen(i)=nq+iim end do do ii=1,nlmax do i=1,nq k=queen0(i,ii,1)+queen0(i,ii,2) if(k==3) then het(igener,ii)=het(igener,ii)+1.0 end if end do het(igener,ii)=het(igener,ii)/real(nq) end do return end !********************************************************************* ! Subroutine for reproducing young queens (monoandry-monogyny) * !********************************************************************* subroutine mate1m(nq,nm,queen,male,queen0,fpair,mpair,het, & f_queen,m_queen,ngener,nlmax,iseed,igener) integer queen(nq,nlmax,2),male(nm,nlmax) integer queen0(nq,nlmax,2) integer fpair(nq),mpair(nm),ia(1000) integer f_queen(nq),m_queen(nq) real het(0:ngener,nlmax) !-------pairing of parents do i=1,nq ia(i)=i end do call random(ia,nq,iseed) do i=1,nq fpair(i)=ia(i) end do do i=1,nm ia(i)=i end do call random(ia,nm,iseed) do i=1,nm mpair(i)=ia(i) end do do i=1,nq do j=1,nlmax rq=rnunf() iiq=int(rq*real(nq))+1 jjq=fpair(iiq) jjm=mpair(iiq) kk=1 if(rnunf()>0.5) then kk=2 end if queen0(i,j,1)=queen(jjq,j,kk) queen0(i,j,2)=male(jjm,j) end do f_queen(i)=jjq m_queen(i)=nq+jjm end do do ii=1,nlmax do i=1,nq k=queen0(i,ii,1)+queen0(i,ii,2) if(k==3) then het(igener,ii)=het(igener,ii)+1.0 end if end do het(igener,ii)=het(igener,ii)/real(nq) end do return end !********************************************************************** ! Subrouting for reproducing young males (m1) by queens * !********************************************************************** subroutine reprom1(nq,nm1,queen,male1,f_m1,nlmax,iseed) integer queen(nq,nlmax,2),male1(nm1,nlmax) integer f_m1(nm1) do i=1,nm1 do j=1,nlmax rq=rnunf() iiq=int(rq*real(nq))+1 kk=1 if(rnunf()>0.5) then kk=2 end if male1(i,j)=queen(iiq,j,kk) end do f_m1(i)=iiq end do return end !********************************************************************** ! Subroutine for reproducing workers (ployandry-polygyny) * !********************************************************************** subroutine mate2(nq,nm,nw,queen,male,worker, & f_worker,m_worker,nlmax,iseed) integer queen(nq,nlmax,2),male(nm,nlmax) integer worker(nw,nlmax,2) integer f_worker(nw),m_worker(nw) do i=1,nw do j=1,nlmax rq=rnunf() iiq=int(rq*real(nq))+1 kk=1 if(rnunf()>0.5) then kk=2 end if worker(i,j,1)=queen(iiq,j,kk) rm=rnunf() iim=int(rm*real(nm))+1 worker(i,j,2)=male(iim,j) end do f_worker(i)=iiq m_worker(i)=nq+iim end do return end !********************************************************************** ! Subroutine for reproducing workers (monoandry-monogyny) * !********************************************************************** subroutine mate2m(nq,nm,nw,queen,male,worker,fpair,mpair, & f_worker,m_worker,nlmax,iseed) integer queen(nq,nlmax,2),male(nm,nlmax) integer worker(nw,nlmax,2) integer fpair(nq),mpair(nm) integer f_worker(nw),m_worker(nw) do i=1,nw do j=1,nlmax rq=rnunf() iiq=int(rq*real(nq))+1 jjq=fpair(iiq) jjm=mpair(iiq) kk=1 if(rnunf()>0.5) then kk=2 end if worker(i,j,1)=queen(jjq,j,kk) worker(i,j,2)=male(jjm,j) end do f_worker(i)=jjq m_worker(i)=nq+jjm end do return end !********************************************************************** ! Subrouting for reproducing worker-produced males (m2) * !********************************************************************** subroutine reprom2(nw,nm2,worker,male2,w_m2,nlmax,iseed) integer worker(nw,nlmax,2),male2(nm2,nlmax) integer w_m2(nm2) do i=1,nm2 do j=1,nlmax rw=rnunf() iiw=int(rw*real(nw))+1 kk=1 if(rnunf()>0.5) then kk=2 end if male2(i,j)=worker(iiw,j,kk) end do w_m2(i)=iiw end do return end !********************************************************************** ! Subrouting for reproducing worker-produced males (m2) * !********************************************************************** subroutine reprom2b(nw,nm2,worker,male2,w_m2,nlmax,iseed) integer worker(nw,nlmax,2),male2(nm2,nlmax) integer w_m2(nm2) do i=1,nm2 do j=1,nlmax iiw=i kk=1 if(rnunf()>0.5) then kk=2 end if male2(i,j)=worker(iiw,j,kk) end do w_m2(i)=iiw end do return end !********************************************************************** ! Subroutine for computating gene frequency in queens * !********************************************************************** Subroutine freqq(queen0,gfreqq,ngener,nlmax,nq,igener) real gfreqq(0:ngener,nlmax) integer queen0(nq,nlmax,2) do j=1,nlmax do i=1,nq do k=1,2 if(queen0(i,j,k)==1) then gfreqq(igener,j)=gfreqq(igener,j)+1.0 end if end do end do end do do i=1,nlmax gfreqq(igener,i)=gfreqq(igener,i)/real(2*nq) end do return end !********************************************************************** ! Subroutine for computating gene frequency in males * !********************************************************************** Subroutine freqm(male0,gfreqm,ngener,nlmax,nm,igener) real gfreqm(0:ngener,nlmax) integer male0(nm,nlmax) do j=1,nlmax do i=1,nm if(male0(i,j)==1) then gfreqm(igener,j)=gfreqm(igener,j)+1.0 end if end do end do do i=1,nlmax gfreqm(igener,i)=gfreqm(igener,i)/real(nm) end do return end !********************************************************************* subroutine random(ia,nn,iseed) integer ia(1000) real a(1000) do i=1,nn a(i)=rnunf() end do do i=1,nn-1 do j=i+1,nn if(a(i)<=a(j)) then work=a(i) iwork=ia(i) a(i)=a(j) ia(i)=ia(j) a(j)=work ia(j)=iwork end if end do end do return end !********************************************************************* ! Subrouitne for constructing coancestry matrix * !********************************************************************* subroutine coanc(nq,nw,nm,nm1,nm2,f0,f1,f_queen,m_queen,f_m1, & f_worker,m_worker,w_m2) real f0(nq+nm,nq+nm),f1(nq+nm,nq+nm) integer f_queen(nq),m_queen(nq),f_m1(nm1),w_m2(nm2) integer f_worker(nw),m_worker(nw) do i=1,nq do j=1,nq if(i==j) then f1(i,i)=0.5*(1.0+f0(f_queen(i),m_queen(i))) else f1(i,j)=0.25*(f0(f_queen(i),f_queen(j))+f0(f_queen(i),m_queen(j)) & +f0(m_queen(i),f_queen(j))+f0(m_queen(i),m_queen(j))) end if end do end do do i=nq+1,nq+nm1 k=i-nq do j=nq+1,nq+nm1 l=j-nq if(i==j) then f1(i,i)=1.0 else f1(i,j)=f0(f_m1(k),f_m1(l)) end if end do end do do i=nq+nm1+1,nq+nm1+nm2 k=i-nq-nm1 do j=nq+nm1+1,nq+nm1+nm2 l=j-nq-nm1 if(i==j) then f1(i,i)=1.0 else if(w_m2(k)==w_m2(l)) then f1(i,j)=0.5*(1.0+f0(f_worker(w_m2(k)),m_worker(w_m2(k)))) else f1(i,j)=0.25*(f0(f_worker(w_m2(k)),f_worker(w_m2(l))) & +f0(f_worker(w_m2(k)),m_worker(w_m2(l))) & +f0(m_worker(w_m2(k)),f_worker(w_m2(l))) & +f0(m_worker(w_m2(k)),m_worker(w_m2(l)))) end if end if end do end do do i=nq+1,nq+nm1 k=i-nq do j=nq+nm1+1,nq+nm1+nm2 l=j-nq-nm1 f1(i,j)=0.5*(f0(f_m1(k),f_worker(w_m2(l)))+f0(f_m1(k),m_worker(w_m2(l)))) f1(j,i)=f1(i,j) end do end do do i=1,nq do j=nq+1,nq+nm1 k=j-nq f1(i,j)=0.5*(f0(f_queen(i),f_m1(k))+f0(m_queen(i),f_m1(k))) f1(j,i)=f1(i,j) end do end do do i=1,nq do j=nq+nm1+1,nq+nm1+nm2 k=j-nq-nm1 f1(i,j)=0.25*(f0(f_queen(i),f_worker(w_m2(k))) & +f0(f_queen(i),m_worker(w_m2(k))) & +f0(m_queen(i),f_worker(w_m2(k))) & +f0(m_queen(i),m_worker(w_m2(k)))) f1(j,i)=f1(i,j) end do end do return end !******************************************************************** ! Subroutine for computing inbreeding coefficent * !******************************************************************** subroutine inbr(nq,nw,nm,f0,f_queen,m_queen,f_coef,ngener,igener) real f0(nq+nm,nq+nm),f_coef(ngener) integer f_queen(nq),m_queen(nq) integer f_worker(nw),m_worker(nw) do i=1,nq f_coef(igener)=f_coef(igener)+f0(f_queen(i),m_queen(i)) end do f_coef(igener)=f_coef(igener)/real(nq) return end