This file contains the codes used to generate the results described in:

 The effects of weak selection on neutral diversity at linked sites (B. Charlesworth)

These are all written in GFortran for the MacIntosh (https://gnuc.org.wiki/GFortranBinariesMacOS)

     Program weaksweep6
!    program for weak sweep effects- fixations
!    diversity summed over course of fixation
!    set of gamma values for fixed recombination rate
!    locus B is neutral; A has a mutant allele A2

    CHARACTER*20 FOUT
    integer :: N,NH,nrep,ngen1,NX,int
    real :: theta,AN,r,r1,q,q1,ak11,ak12,ak22,u2
    write (*,*) 'Output file?'
    read(*,*)  FOUT
    OPEN (1,FILE=FOUT)

    write (*,*) 'population size?'
    read (*,*) N
    write (*,*) 'number of reps'
    read (*,*) nrep
    write (*,*) 'Initial gamma and h?'
    read (*,*) gam1,h2
    write(*,*) 'No. of gamma values and increment?'
    read (*,*) ngam,delgam
    write (*,*) 'Scaled mutation rate?'
    read (*,*) theta
    write (*,*) 'Rho?'
    read (*,*) rho
!    scaled selection, recombination and mutation parameters

    anrep=nrep
    NH=2*N
    AN=NH

    r1=rho/AN
    u2=theta/AN
    u=0.5*u2
    s2=gam1
!    unscaled parameters

    write(1,*) 'Program for effect of fixation of weakly selected mutations on expected diversity'
    write(1,*) 'Locus A is subject to selection'
    write(1,*) 'Diversities for genotypes at A summed over paths to fixation'
    write(1,*) ''
    write(1,*) 'Scaled parameters'
    write(1,*) 'Initial gamma= ',gam1,' h= ',h2
    write(1,*) 'No. of gamma values= ',ngam,'increment= ',delgam
    write(1,*) 'Rho= ',rho
    write(1,*) 'theta= ',theta
    write(1,*) 'u= ',u
    write(1,*) 'Population size= ',N
    write(1,*) 'Number of replicate fixations= ',nrep
    write(1,*) ''

    write(*,*) 'Population size= ',N
    write(*,*) 'Initial gamma= ',gam1,' h= ',h2
    write(*,*) 'No. of gamma values= ',ngam,'increment= ',delgam
    write(*,*) 'rho= ',rho
    write(*,*) 'theta= ',theta
    write(*,*) ''
    write(*,*) 'Number of replicates= ',nrep
    write(*,*) ''

    q0=1/AN
    write (*,*) 'Initial A2 allele frequency= ',q0
    write (1,*) 'Initial A2 allele frequency= ',q0
    write (1,*) ''

    a=h2
    b=1.0-2*h2
    r1=rho/AN
    
    call init_random_seed()

    do 350 igam=1,ngam

    gam=gam1+(igam-1)*delgam
    write(1,*) ''
    write(1,*) 'gamma= ',gam
    write(*,*) ''
    write(*,*) 'gamma= ',gam
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!    starts simulations
    write(1,*) ''
    ros=rho/gam
    s2=gam/AN
    write(1,*) 'r/s= ',ros
    write(1,*) ''
    write(*,*) ''
    write(*,*) 'r/s= ',ros
    write(*,*) ''

    i=0
    ak11gs=0.0
    ak12gs=0.0
    ak22gs=0.0
    akgms=0.0

    ak11gss=0.0
    ak12gss=0.0
    ak22gss=0.0
    akgmss=0.0

    ak22s=0.0
    ak22ss=0.0

    delpis=0.0
    delpiss=0.0

    stf=0.0
    sstf=0.0
    sH=0.0
    ssH=0.0
    nrun=0

20    i=i+1
    nrun=nrun+1
    ngen1=0

10    if(ngen1.eq.0) then
    q=q0
    p=1-q
    nlA2=0
    ak11=theta
    ak12=theta
    ak22=0.0
!    akij is neutral diversity between alleles Ai and Aj

    ak11g=0.0
    ak12g=0.0
    ak22g=0.0
    akgm=0.0
!    sums of neutral diversities within A locus genotypes over path of A2 to fixation
        end if

15    q1=q
!    stores previous generation frequency

    delsq=p*q*s2*(q*(1-2*h2)+h2)
    delsq=delsq/(1-2*p*q*s2*h2-(q**2)*s2)
    q=q+delsq
    p=1-q
!    deterministic allele frequency change
    ngen1=ngen1+1

!    write(*,*)  'Generation ',ngen1,' q= ',q,'delq= ',delsq
    N11=0
  do 120 j1=1,NH
    
    call random_number(r)
    z=r
    if(z.le.p) then
    N11=N11+1
    go to 120
    end if
120     continue
!    binomial sampling of A1 alleles

    p=N11/AN
    q=1-p
!    new allele frequencies
    H=H+2*p*q
!    cumulative diversity at selected locus over run

    if(N11.eq.NH) then
!    write(*,*) ''
!    write(*,*) 'A2 is lost'
!    write(*,*) ''
    i=i-1
    go to 20
    end if
!    A2 is lost
!    run is terminated and i is reset

    call swdiv(AN,r1,q1,ak11,ak12,ak22,u2)

    ak11a=ak11/theta
    ak12a=ak12/theta
    ak22a=ak22/theta

    ak11g=ak11g+ak11a
    ak12g=ak12g+ak12a
    ak22g=ak22g+ak22a
    akgm=akgm+(p**2)*ak11a+(2*p*q)*ak12a+(q**2)*ak22a
!    sum of mean neutral diversity (relative to theta)

    if(N11.eq.0) then
!    write(*,*) ''
!    write(*,*) 'Replicate= ',i
!    write(*,*) 'ak11= ',ak11a,' ak12= ',ak12a,' ak22= ',ak22a
!    write(*,*) ''
!    pause
!    A2 is fixed and run is completed

    tfix=ngen1

    stf=stf+tfix
    sstf=sstf+tfix**2
!    sum and sum of squares of time to fixation

    ak11gs=ak11gs+ak11g
    ak12gs=ak12gs+ak12g
    ak22gs=ak22gs+ak22g
    akgms=akgms+akgm

    ak11gss=ak11gss+ak11g**2
    ak12gss=ak12gss+ak12g**2
    ak22gss=ak22gss+ak22g**2
    akgmss=akgmss+akgm**2
!    sums and sums of squares of mean diversities over paths to fixation

    ak22s=ak22s+ak22a
    ak22ss=ak22ss+ak22a**2
!    sum and sum of squares of diversities of A2/A2 population at time of fixation of A2 (relative to theta)

    delpi=akgm-tfix+(AN*(ak22a-1))
!    measure of potential change in diversity due to recurrent sweeps

    delpis=delpis+delpi
    delpiss=delpiss+delpi**2

  if(i.eq.nrep) then
   go to 200
      else
  go to 20
  end if
    end if

    go to 15
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

200    anrun=nrun
    nloss=nrun-nrep
    pfix=nrep/anrun
    ploss=1.0-pfix

    tfm=stf/anrep
    vtf=(sstf-(stf**2)/anrep)/(anrep-1.0)
    setf=sqrt(vtf/anrep)

    ak11gm=ak11gs/stf
    ak12gm=ak12gs/stf
    ak22gm=ak22gs/stf
    akgmm=akgms/stf
    vak11g=(ak11gss-(ak11gs**2)/anrep)/(anrep-1.0)
    vak11g=vak11g/(tfm**2)
    vak12g=(ak12gss-(ak12gs**2)/anrep)/(anrep-1.0)
    vak12g=vak12g/(tfm**2)
    vak22g=(ak22gss-(ak22gs**2)/anrep)/(anrep-1.0)
    vak22g=vak22g/(tfm**2)
    vakgm=(akgmss-(akgms**2)/anrep)/(anrep-1.0)
    vakgm=vakgm/(tfm**2)
    seak11g=sqrt(vak11g/anrep)
    seak12g=sqrt(vak12g/anrep)
    seak22g=sqrt(vak22g/anrep)
    seakgm=sqrt(vakgm/anrep)
    
    ak22m=ak22s/anrep
    vak22=(ak22ss-(ak22s**2/anrep))/(anrep-1.0)
    seak22=sqrt(vak22/anrep)

    tfm=tfm/AN
    setf=setf/AN
!    time in 2N generation units

    delpim=delpis/(AN*anrep)
    vdelpi=(delpiss-(delpis**2)/anrep)/(anrep-1.0)
    sedelpi=sqrt(vdelpi/anrep)/AN

    write(*,*) 'Frequency of fixations of A2= ',pfix
    write(*,*) 'Mean final relative diversity= ',ak22m,' s.e.= ',seak22
    delak=1.0-ak22m
    sedef=seak22
    write(*,*) 'Mean final diversity reduction= ',delak,' s.e.= ',sedef
    
    write(1,*) 'Total number of runs= ',nrun
    write(1,*) 'Frequency of fixations of A2= ',pfix

    write (1,*) 'Mean time to fixation (2N generation units)= ',tfm,' s.e.= ',setf
    timp=2*(1-0.05556*gam*b-0.05556*(gam**2)*(a*(a+b)+0.19667*b**2))
    write(1,*) 'Predicted approximate mean time to fixation= ',timp
    write(1,*) ''
    write(1,*) 'Mean weighted relative diversities over paths to fixation'
    write(1,*) 'A1A1= ',ak11gm,'s.e.= ',seak11g
    write(1,*) 'A1A2= ',ak12gm,'s.e.= ',seak12g
    write(1,*) 'A2A2= ',ak22gm,'s.e.= ',seak22g
    write(1,*) 'Mean= ',akgmm,'s.e.= ',seakgm
    write(1,*) ''
    write(1,*) 'Mean final relative diversity= ',ak22m,' s.e.= ',seak22
    write(1,*) 'Mean final diversity reduction= ',delak,' s.e.= ',sedef
    write(1,*) ''
    write(1,*) 'Weighted measure of potential recurrent sweep effect= ',delpim,' s.e.= ',sedelpi
    write(1,*) ''
    delpim1=tfm*(akgmm-1)-delak
    write(1,*) 'Alternative weighted measure of potential recurrent sweep effect= ',delpim1
    write(1,*) ''

350    continue

    end program weaksweep6

    subroutine swdiv(AN,r,x,ak11,ak12,ak22,u2)
    real AN,r,x,ak11,ak12,ak22,u2
!    recursions for pairwise diversities
    
    y=1-x
    c1=1.0/(AN*y)
    c2=1.0/(AN*x)
    r1=2.0*r*x
    r2=2.0*r*y

    ak11n=ak11+u2+r1*ak12-(r1+(u2*(1.0-c1)+c1))*ak11
    ak12n=(1-r)*ak12+0.5*r2*ak11+0.5*r1*ak22+u2
    ak22n=ak22+u2+r2*ak12-(r2+(u2*(1.0-c2)+c2))*ak22
!    write(9,*) 'x=',x,' c1= ',c1,' c2= ',c2
!    write(9,*) 'r1= ',r1,' r2= ',r2
!    write(9,*) 'Old ak values',ak11,ak12,ak22
!    write(9,*) 'New ak values',ak11n,ak12n,ak22n
!    write(9,*) ''
!    pause

    ak11=ak11n
    ak12=ak12n
    ak22=ak22n

    end subroutine swdiv

    subroutine init_random_seed()
    use iso_fortran_env, only: int64
    implicit none
    integer, allocatable :: seed(:)
    integer :: i, n, un, istat, dt(8), pid
    integer(int64) :: t

    call random_seed(size = n)
    allocate(seed(n))
! First try if the OS provides a random number generator
    open(newunit=un, file="/dev/urandom", access="stream", &
     form="unformatted", action="read", status="old", iostat=istat)
    if (istat == 0) then
   read(un) seed
   close(un)
    else
   ! Fallback to XOR:ing the current time and pid. The PID is
   ! useful in case one launches multiple instances of the same
   ! program in parallel.
   call system_clock(t)
   if (t == 0) then
      call date_and_time(values=dt)
      t = (dt(1) - 1970) * 365_int64 * 24 * 60 * 60 * 1000 &
           + dt(2) * 31_int64 * 24 * 60 * 60 * 1000 &
           + dt(3) * 24_int64 * 60 * 60 * 1000 &
           + dt(5) * 60 * 60 * 1000 &
           + dt(6) * 60 * 1000 + dt(7) * 1000 &
           + dt(8)
   end if
   pid = getpid()
   t = ieor(t, int(pid, kind(t)))
   do i = 1, n
      seed(i) = lcg(t)
   end do
    end if
    call random_seed(put=seed)
    contains
! This simple PRNG might not be good enough for real work, but is
! sufficient for seeding a better PRNG.
    function lcg(s)
  integer :: lcg
  integer(int64) :: s
  if (s == 0) then
     s = 104729
  else
     s = mod(s, 4294967296_int64)
  end if
  s = mod(s * 279470273_int64, 4294967291_int64)
  lcg = int(mod(s, int(huge(0), int64)), kind(0))
    end function lcg
    end subroutine init_random_seed

    Program weaksweep8
!    program for weak sweep effects- losses
!    diversity summed over course of loss
!    set of gamma values for fixed recombination rate
!    locus B is neutral; A has a mutant allele A2

    CHARACTER*20 FOUT
    CHARACTER*10 CONT
    integer :: N,NH,nrep,ngen1,NX,int
    real :: theta,AN,r,r1,q,q1,ak11,ak12,ak22,u2
    write (*,*) 'Output file?'
    read(*,*)  FOUT
    OPEN (1,FILE=FOUT)

    write (*,*) 'population size?'
    read (*,*) N
    write (*,*) 'number of reps'
    read (*,*) nrep
    write (*,*) 'Initial gamma and h?'
    read (*,*) gam1,h2
    write(*,*) 'No. of gamma values and increment?'
    read (*,*) ngam,delgam
    write (*,*) 'Scaled mutation rate?'
    read (*,*) theta
    write (*,*) 'Rho?'
    read (*,*) rho
!    scaled selection, recombination and mutation parameters

    anrep=nrep
    NH=2*N
    AN=NH

    r1=rho/AN
    u2=theta/AN
    u=0.5*u2
    s2=gam1
!    unscaled parameters

    write(1,*) 'Program for effect of loss of weakly selected mutations on expected diversity'
    write(1,*) 'Locus A is subject to selection'
    write(1,*) 'Diversities for genotypes at A summed over paths to loss’
    write(1,*) ''
    write(1,*) 'Scaled parameters'
    write(1,*) 'Initial gamma= ',gam1,' h= ',h2
    write(1,*) 'No. of gamma values= ',ngam,'increment= ',delgam
    write(1,*) 'Rho= ',rho
    write(1,*) 'theta= ',theta
    write(1,*) 'u= ',u
    write(1,*) 'Population size= ',N
    write(1,*) 'Number of replicate fixations= ',nrep
    write(1,*) ''

    write(*,*) 'Population size= ',N
    write(*,*) 'Initial gamma= ',gam1,' h= ',h2
    write(*,*) 'No. of gamma values= ',ngam,'increment= ',delgam
    write(*,*) 'rho= ',rho
    write(*,*) 'theta= ',theta
    write(*,*) ''
    write(*,*) 'Number of replicates= ',nrep
    write(*,*) ''

    q0=1/AN
    write (*,*) 'Initial A2 allele frequency= ',q0
    write (1,*) 'Initial A2 allele frequency= ',q0
    write (1,*) ''
    write(*,*) 'Continue?'
    read (*,*) CONT
    a=h2
    b=1.0-2*h2
    r1=rho/AN
    
    call init_random_seed()

    do 350 igam=1,ngam

    gam=gam1+(igam-1)*delgam
    write(1,*) ''
    write(1,*) 'gamma= ',gam
    write(*,*) ''
    write(*,*) 'gamma= ',gam
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!    starts simulations
    write(1,*) ''
    ros=rho/gam
    s2=gam/AN
    write(1,*) 'r/s= ',ros
    write(1,*) ''
    write(*,*) ''
    write(*,*) 'r/s= ',ros
    write(*,*) ''

    i=0
    ak11gs=0.0
    ak12gs=0.0
    ak22gs=0.0
    akgms=0.0

    ak11gss=0.0
    ak12gss=0.0
    ak22gss=0.0
    akgmss=0.0

    ak11s=0.0
    ak11ss=0.0

    delpis=0.0
    delpiss=0.0

    stl=0.0
    sstl=0.0
    sH=0.0
    ssH=0.0
    nrun=0

20    i=i+1
    nrun=nrun+1
    ngen1=0

10    if(ngen1.eq.0) then
    q=q0
    p=1-q
    nlA2=0
    ak11=theta
    ak12=theta
    ak22=0.0
!    akij is neutral diversity between alleles Ai and Aj

    H=0.0
    ak11g=0.0
    ak12g=0.0
    ak22g=0.0
    akgm=0.0
!    sums of neutral diversities within A locus genotypes over path of A2 to fixation
        end if

15    q1=q
!    stores previous generation frequency

    delsq=p*q*s2*(q*(1-2*h2)+h2)
    delsq=delsq/(1-2*p*q*s2*h2-(q**2)*s2)
    q=q+delsq
    p=1-q
!    deterministic allele frequency change
    ngen1=ngen1+1

!    write(*,*)  'Generation ',ngen1,' q= ',q,'delq= ',delsq
    N11=0
  do 120 j1=1,NH
    
    call random_number(r)
    z=r
    if(z.le.p) then
    N11=N11+1
    go to 120
    end if
120     continue
!    binomial sampling of A1 alleles

    p=N11/AN
    q=1-p
!    new allele frequencies
    H=H+2*p*q
!    cumulative diversity at selected locus over run

    if(N11.eq.0) then
    i=i-1
    go to 20
    end if
!    A2 is fixed
!    run is terminated and i is reset

    call swdiv(AN,r1,q1,ak11,ak12,ak22,u2)

    ak11a=ak11/theta
    ak12a=ak12/theta
    ak22a=ak22/theta

    ak11g=ak11g+ak11a
    ak12g=ak12g+ak12a
    ak22g=ak22g+ak22a
    akgm=akgm+(p**2)*ak11a+(2*p*q)*ak12a+(q**2)*ak22a
!    sum of mean neutral diversity (relative to theta)

    if(N11.eq.NH) then
!    write(*,*) ''
!    write(*,*) 'Replicate= ',i
!    write(*,*) 'ak11= ',ak11a,' ak12= ',ak12a,' ak22= ',ak22a
!    write(*,*) ''
!    pause
!    A2 is fixed and run is completed

    tloss=ngen1
    stl=stl+tloss
    sstl=sstl+tloss**2
!    sum and sum of squares of time to loss
!    sum and sum of squares of time to loss

    ak11gs=ak11gs+ak11g
    ak12gs=ak12gs+ak12g
    ak22gs=ak22gs+ak22g
    akgms=akgms+akgm

    ak11gss=ak11gss+ak11g**2
    ak12gss=ak12gss+ak12g**2
    ak22gss=ak22gss+ak22g**2
    akgmss=akgmss+akgm**2
!    sums and sums of squares of mean diversities over paths to fixation

    ak11s=ak11s+ak11a
    ak11ss=ak11ss+ak11a**2
!    sum and sum of squares of diversities of A2/A2 population at time of fixation of A2 (relative to theta)

    delpi=akgm-tloss+(AN*(ak11a-1))
!    measure of potential change in diversity due to recurrent sweeps

    delpis=delpis+delpi
    delpiss=delpiss+delpi**2

  if(i.eq.nrep) then
   go to 200
      else
  go to 20
  end if
    end if
!    A2 is lost and run is completed

    go to 15
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

200    anrun=nrun
    nloss=nrun-nrep
    ploss=nrep/anrun
    pfix=1.0-ploss

    tlm=stl/anrep
    vtl=(sstl-(stl**2)/anrep)/(anrep-1.0)
    setl=sqrt(vtl/anrep)

    ak11gm=ak11gs/stl
    ak12gm=ak12gs/stl
    ak22gm=ak22gs/stl
    akgmm=akgms/stl
    vak11g=(ak11gss-(ak11gs**2)/anrep)/(anrep-1.0)
    vak11g=vak11g/(tlm**2)
    vak12g=(ak12gss-(ak12gs**2)/anrep)/(anrep-1.0)
    vak12g=vak12g/(tlm**2)
    vak22g=(ak22gss-(ak22gs**2)/anrep)/(anrep-1.0)
    vak22g=vak22g/(tlm**2)
    vakgm=(akgmss-(akgms**2)/anrep)/(anrep-1.0)
    vakgm=vakgm/(tlm**2)
    seak11g=sqrt(vak11g/anrep)
    seak12g=sqrt(vak12g/anrep)
    seak22g=sqrt(vak22g/anrep)
    seakgm=sqrt(vakgm/anrep)
    
    ak11m=ak11s/anrep
    vak11=(ak11ss-(ak11s**2/anrep))/(anrep-1.0)
    seak11=sqrt(vak11/anrep)

    tlm=tlm/AN
    setl=setl/AN
!    time in 2N generation units

    delpim=delpis/(AN*anrep)
    vdelpi=(delpiss-(delpis**2)/anrep)/(anrep-1.0)
    sedelpi=sqrt(vdelpi/anrep)/AN

    write(*,*) 'Frequency of losses of A2= ',ploss
    write(*,*) 'Mean final relative diversity= ',ak11m,' s.e.= ',seak11
    delak=1.0-ak11m
    sedef=seak11
    write(*,*) 'Mean final diversity reduction= ',delak,' s.e.= ',sedef
    
    write(1,*) 'Total number of runs= ',nrun
    write(1,*) 'Frequency of losses of A2= ',ploss
    write (1,*) 'Mean time to loss (2N generation units)= ',tlm,' s.e.= ',setl
    timp=2*(log(AN)-0.4228-0.2778*gam*b-0.05556*(gam**2)*(a*(5*a+3.5*b)+0.17*b**2))
    timp=timp/AN
    write(1,*) 'Predicted approximate mean time to loss= ',timp
    write(1,*) ''
    write(1,*) 'Mean weighted relative diversities over paths to loss'
    write(1,*) 'A1A1= ',ak11gm,'s.e.= ',seak11g
    write(1,*) 'A1A2= ',ak12gm,'s.e.= ',seak12g
    write(1,*) 'A2A2= ',ak22gm,'s.e.= ',seak22g
    write(1,*) 'Mean= ',akgmm,'s.e.= ',seakgm
    write(1,*) ''
    write(1,*) 'Mean final relative diversity= ',ak11m,' s.e.= ',seak11
    write(1,*) 'Mean final diversity reduction= ',delak,' s.e.= ',sedef
    write(1,*) ''
    write(1,*) 'Weighted measure of potential recurrent sweep effect= ',delpim,' s.e.= ',sedelpi
    write(1,*) ''
    delpim1=tlm*(akgmm-1)-delak
    write(1,*) 'Alternative weighted measure of potential recurrent sweep effect= ',delpim1
    write(1,*) ''

350    continue

    end program weaksweep8

    subroutine swdiv(AN,r,x,ak11,ak12,ak22,u2)
    real AN,r,x,ak11,ak12,ak22,u2
!    recursions for pairwise diversities
    
    y=1-x
    c1=1.0/(AN*y)
    c2=1.0/(AN*x)
    r1=2.0*r*x
    r2=2.0*r*y

    ak11n=ak11+u2+r1*ak12-(r1+(u2*(1.0-c1)+c1))*ak11
    ak12n=(1-r)*ak12+0.5*r2*ak11+0.5*r1*ak22+u2
    ak22n=ak22+u2+r2*ak12-(r2+(u2*(1.0-c2)+c2))*ak22
!    write(9,*) 'x=',x,' c1= ',c1,' c2= ',c2
!    write(9,*) 'r1= ',r1,' r2= ',r2
!    write(9,*) 'Old ak values',ak11,ak12,ak22
!    write(9,*) 'New ak values',ak11n,ak12n,ak22n
!    write(9,*) ''
!    pause

    ak11=ak11n
    ak12=ak12n
    ak22=ak22n

    end subroutine swdiv

    subroutine init_random_seed()
    use iso_fortran_env, only: int64
    implicit none
    integer, allocatable :: seed(:)
    integer :: i, n, un, istat, dt(8), pid
    integer(int64) :: t

    call random_seed(size = n)
    allocate(seed(n))
! First try if the OS provides a random number generator
    open(newunit=un, file="/dev/urandom", access="stream", &
     form="unformatted", action="read", status="old", iostat=istat)
    if (istat == 0) then
   read(un) seed
   close(un)
    else
   ! Fallback to XOR:ing the current time and pid. The PID is
   ! useful in case one launches multiple instances of the same
   ! program in parallel.
   call system_clock(t)
   if (t == 0) then
      call date_and_time(values=dt)
      t = (dt(1) - 1970) * 365_int64 * 24 * 60 * 60 * 1000 &
           + dt(2) * 31_int64 * 24 * 60 * 60 * 1000 &
           + dt(3) * 24_int64 * 60 * 60 * 1000 &
           + dt(5) * 60 * 60 * 1000 &
           + dt(6) * 60 * 1000 + dt(7) * 1000 &
           + dt(8)
   end if
   pid = getpid()
   t = ieor(t, int(pid, kind(t)))
   do i = 1, n
      seed(i) = lcg(t)
   end do
    end if
    call random_seed(put=seed)
    contains
! This simple PRNG might not be good enough for real work, but is
! sufficient for seeding a better PRNG.
    function lcg(s)
  integer :: lcg
  integer(int64) :: s
  if (s == 0) then
     s = 104729
  else
     s = mod(s, 4294967296_int64)
  end if
  s = mod(s * 279470273_int64, 4294967291_int64)
  lcg = int(mod(s, int(huge(0), int64)), kind(0))
    end function lcc
    end subroutine init_random_seed

    program weaksweep9
!    program for autosomal weak sweep effect approximations- losses and fixation
!    uses series representation for set of h values

    CHARACTER*20 FOUT
    CHARACTER*10 CONT
    integer :: N,NH,NMAX,nh1
    real :: gam,AN,ttot1,ttot2,Htot1,Htot2,P01,P11,a,b,tnet,Hnet
    write(*,*) 'Output file?'
    read(*,*)  FOUT
    OPEN (1,FILE=FOUT)

    write (*,*) 'Population size?'
    read (*,*) N
    write (*,*) 'Maximum size of series?'
    read (*,*) NMAX
    write (*,*) 'Gamma and initial h?'
    read (*,*) gam,h1
!   gam is the scaled selection parameter, h1 is the initial dominance coefficient
    write(*,*) 'No. of h values and increment?'
    read (*,*) nh1,delh

    NH=2*N
    AN=NH

    write(1,*) 'Program for effects of fixations and losses of weakly selected mutations'
    write(1,*) 'Autosomal locus with random mating'
    write(1,*) ''
    write(1,*) 'Population size= ',N
    write(1,*) 'Maximum size of series= ',NMAX
    write(*,*) 'Maximum size of series= ',NMAX
    write(1,*) 'Initial h= ',h1,' gamma= ',gam
    write(1,*) 'No. of h values= ',nh1,'increment= ',delh
    write(*,*) 'No. of h values= ',nh1,'increment= ',delh
    write(*,*) ''
    write(*,*) 'Continue?'
    read (*,*) CONT
!  allows program to be cancelled if desired

    write(1,*) ''
    write(1,*) 'Times are in units of coalescent time(2Ne generations)'
    write(1,*) 'Diversities over fixation or loss are relative to 2Ne'
    write(*,*) ''

    do 10 ih=1,nh1

    h=h1+(ih-1)*delh
    a=h
    b=1.0-2*h
    c=a+b

    write(1,*) ''
    write(1,*) 'h= ',h
    write(*,*) ''
    write(*,*) 'h= ',h
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    timf=2.0*(1-0.05556*gam*b-0.05556*(gam**2)*(a*(a+b)+0.1967*b**2))
    write(1,*) 'Approximate mean time to fixation= ',timf
    write(*,*) 'Approximate mean time to fixation= ',timf
    write(1,*) ''
    hetf=0.6667*(1-0.06667*gam*b-0.06667*(gam**2)*(a*(a+b)+0.1905*b**2))
    hetf=hetf
    write(1,*) 'Approximate mean relative diversity over fixation= ',hetf
    write(1,*) ''
    timl=2*(log(AN)-0.4228-0.2778*gam*b-0.05556*(gam**2)*(a*(5*a+3.5*b)+0.17*b**2))
    timl=timl/AN
    write(1,*) 'Approximate mean time to loss= ',timl
    write(1,*) ''
    hetl=1.3333*(1-0.1333*gam*b-(gam**2)*(0.1333*(a**2)+0.1*a*b+0.00476*b**2))
    hetl=hetl/AN
    write(1,*) 'Approximate mean relative diversity over loss= ',hetl
    write(1,*) ''
    tsoj=2*(log(AN)+0.5772+gam*a-0.02778*b*(gam**2)*(5.0*a+b))/AN
    write(1,*) 'Approximate mean sojourn time= ',tsoj
    hetsoj=2*(1+0.3333*gam*a-0.05556*b*(gam**2)*(a+0.2*b))
    write(1,*) ''
    write(1,*) 'Approximate mean absolute diversity over sojourn= ',hetsoj
!   write(*,*) 'Approximate mean absolute diversity over sojourn= ',hetsoj
    write(1,*) ''

    b1=abs(b)
    gam2=abs(gam)
    
!   write(*,*) 'gam2= ',gam2

    if(b1.le.0.001) then
    call swadd(AN,NH,gam,ttot1,ttot2,Htot1,Htot2,tnet,Hnet,P01,P11)
!   results for semi-dominance or neutrality
    go to 5
    end if
    
    if(gam2.ge.0.001) then
    call swintegs(AN,N,NMAX,gam,a,b,ttot1,ttot2,Htot1,Htot2,tnet,Hnet,P01,P11)
    else
    call swadd(AN,NH,gam,ttot1,ttot2,Htot1,Htot2,tnet,Hnet,P01,P11)
    end if

5   write(1,*) ''
    write(1,*) 'Results of integrations'
    write(1,*) ''
    write(1,*) 'Fixation probability= ',P11
    write(1,*) 'Mean time to fixation= ',ttot1
    Htot1=Htot1
    write(1,*) 'Mean diversity over fixation= ',Htot1
    write(1,*) ''
    write(1,*) 'Loss probability= ',P01
    write(1,*) 'Mean time to loss= ',ttot2
    write(1,*) 'Mean diversity over loss= ',Htot2
    write(1,*) ''
    write(1,*) 'Mean sojourn time= ',tnet
    write(1,*) 'Mean absolute diversity over sojourn= ',Hnet
    write(1,*) ''

10   continue

    end program weaksweep9

    subroutine swintegs(AN,N,NMAX,gam,a,b,ttot1,ttot2,Htot1,Htot2,tnet,Hnet,P01,P11)
    real :: AN,a,b,y,S0,aifac(0:2000),ip(2000),al2(2000),P01,P11,gam
    real :: psif(2000),t1(2000),t2(2000),H1(2000),H2(2000),S2(2000),S3
    real :: ttot1,ttot2,Htot1,Htot2,tnet,Hnet
    integer :: N,N1,NMAX
!   sojourn time integral calculations
!   b and gamma must be non-zero for these calculations

    N1=2*N
    a2=2*a

    ap1=gam*b
    ap=abs(ap1)
    ap=sqrt(ap)
!   factor for transforming y to u

    const=exp(gam*(a**2)/b)
!   factor for converting psi function from y to z=y+a/b
!   this means that we can use a Taylor series in z to replace psi(y)=exp-y(2a+by)

    umin=(a/b)*ap
!   lower bound for u=z/sqrt(ap)
!   note that umin < 0 if a/b < 0

    ttot1=0.0
    ttot2=0.0
    Htot1=0.0
    Htot2=0.0

    aifac(0)=0
    do 5 i=1,NMAX
    ai=i
    ip(i)=2*i+1
    al2(i)=log(2.0*ai+1.0)
    aifac(i)=aifac(i-1)+log(ai)
5   continue
!   natural logs of 2i+1 and factorial i
!   logs are used to avoid large numbers
    
    do 100 k=1,N1
    
    y=k/AN
    umax=(y+a/b)*ap
!   allele frequency and u transform for integral upper bound
    
    S21=umax
    if(umax.le.0) then
    umaxab=abs(umax)
    end if

    ifac=1
!   initialises series for y > 0

!   write (*,*) 'k= ',k
    do 40 i1=1,NMAX
    ifac=0-ifac
    
    S2(k)=S21
    if(umax.lt.0) then
    sterm=ip(i1)*log(umaxab)-aifac(i1)-al2(i1)
!   log of ith term in series for y=0
    Ser=((0.0-1.0)**ip(i1))*exp(sterm)
    else
    sterm=ip(i1)*log(umax)-aifac(i1)-al2(i1)
    Ser=exp(sterm)
    end if

    if(ap1.gt.0) then
    Ser=Ser*ifac
    end if
    
    S2(k)=S2(k)+Ser
    
!   write(*,*) 'i1=',i1,' S2= ',S2(k)
    del=(S2(k)-S21)/S21
   if(abs(del).le.0.000001) go to 130

    if(i1.eq.NMAX) then
    write(*,*) 'Series for S2 (integral from 1/(2N) to y)does not converge'
    write(1,*) 'Series for S2 (integral from 1/(2N) to y) does not converge'
    go to 130
    end if

    S21=S2(k)
40  continue

130 S22=S2(k)
!   series for integral of transformed psi function from 1/(2N) to y
    write(*,*) 'k= ',k,' Final term index= ',i1,' y= ',y,' u= ',umax,'S2= ',S22
    write(*,*) ''
100 continue

    S2MAX=S22
    
    S01=umin
    uminab=abs(umin)
!   write(*,*) 'uminab =',uminab

    if(uminab.le.0.001) then
    S0=0
    go to 235
    end if
    
    ifac=1
!   initialises series for y=0
    write(*,*) 'umin= ',umin
    
    do 50 i2=1,NMAX
    S0=S01
    ifac=0-ifac

    if(umin.lt.0) then
    sterm=ip(i2)*log(uminab)-aifac(i2)-al2(i2)
!   log of ith term in series for y=0
    Ser=((0.0-1.0)**ip(i2))*exp(sterm)
    else
    sterm=ip(i2)*log(umin)-aifac(i2)-al2(i2)
    Ser=exp(sterm)
    end if

    if(ap1.gt.0) then
    Ser=Ser*ifac
    end if

    S0=S0+Ser
    write(*,*) 'i2=',i2,'S0=',S0
    del=(S0-S01)/S01
   if(abs(del).le.0.000001) go to 230

    if(i2.eq.NMAX) then
    write(*,*) 'Series for S0 does not converge'
    write(1,*) 'Series for S0 does not converge'
    go to 230
    end if

    S01=S0
50  continue
!   series for lower bound of transformed psi integral

230 write(1,*) 'Results of integrations'

235 S1=S2MAX-S0
    write(*,*) 'S0=',S0,'S2MAX= ',S2MAX,'S1= ',S1
    S1=S1/ap
!   proportional to psi integral from 0 to 1
    S4=S2MAX-S2(1)
!   write(*,*)'S4= ',S4
    S4=S4/ap
!   proportional to probability of loss of new mutation
    S5=S2(1)-S0
!   write(*,*)'S5= ',S5
    S5=S5/ap
!   proportional to probability of fixation of new mutation
!   write(*,*) 'S0= ',S0,' S1= ',S1,'S4= ',S4,' S5 =',S5
!   write(*,*) ''
    aS1=log(S1)
    aS4=log(S4)
    aS5=log(S5)
   
    N1=N1-1
    do 110 i=1,N1
    y=i/AN
    x=1.0-y
    S2(i)=(S2(i)-S0)/ap
!   proportional to psi integral from 0 to i
    
    z=y+a/b
    psif(i)=ap1*z**2
    fk=psif(i)-log(y)-log(x)
!   proportional to log of inverse of psi x y(1-y)

    S3=S1-S2(i)
!   proportional to psi integral from y to 1
    aS3=log(S3)

    psi1y=S3*const
    psy=exp(-psif(i))*const
    pfixy=S2(i)/S1
!   write(*,*) ''
!   write(*,*) 'i= ',i,' y= ',y,' psy= ',psy,'psi1y= ',psi1y
!   write(*,*) 'S2= ',S2(i),' S3= ',S3
!   write(*,*) 'const= ',const,' psif= ',psif(i),' P(y)= ',pfixy,' z= ',z
    at1=log(S2(i))-log(S1)+log(S3)+fk
!   at1=log(pfixy)+log(psi1y)-log(psy)-log(y)-log(x)
    t1(i)=2.0*exp(at1)
!   sojourn time density at k at on way to fixation
    
    at2=2.0*aS3+aS5-aS1-aS4+fk
    t2(i)=2.0*exp(at2)
!   sojourn time density at frequency k at on way to loss

    H1(i)=2.0*y*(1-y)*t1(i)
    H2(i)=2.0*y*(1-y)*t2(i)
!   diversity densities at frequency y
    
!   write(*,*) 'i= ',i,' y= ',y
!   write(*,*) 't1= ',t1(i),'H1= ',H1(i)
!   write(*,*) 't2= ',t2(i),'H2= ',H2(i)

    ttot1=ttot1+t1(i)
    ttot2=ttot2+t2(i)
!   net sojourn time densities
    Htot1=Htot1+H1(i)
    Htot2=Htot2+H2(i)
!   net diversity densities

110 continue

    ttot1=ttot1/AN
    ttot2=ttot2/AN
!   sojourn time integrals
    
    P11=S5/S1
    P01=1.0-P11
!   probabilities of fixation and loss

    Hnet=P11*Htot1+P01*Htot2
    Htot1=Htot1/AN
    Htot2=Htot2/AN
    tnet=P11*ttot1+P01*ttot2
!   net sojourn time and diversity

    end subroutine swintegs

    subroutine swadd(AN,NH,gam,ttot1,ttot2,Htot1,Htot2,tnet,Hnet,P01,P11)
    real :: AN,a,y,S0,aifac(0:2000),ip(2000),al2(2000),P01,P11,gam
    real :: psif(2000),t1(2000),t2(2000),H1(2000),H2(2000),S2(2000),S3
    real :: ttot1,ttot2,Htot1,Htot2,tnet,Hnet
    integer :: NH
!   calculations for semidominant and neutral cases

    gam2=abs(gam)

!   write(*,*) 'swadd entered; gam2= ',gam2
    P11=1.0/AN
    P01=1.0-P11
    f0=exp(gam)-1.0
    f1=1.0-exp(0.0-gam)
    f2=1.0-exp(0.0-gam*P11)
    f3=1.0-exp(0.0-gam*P01)
    f4=exp(gam*P01)-1.0
    
    if(gam2.ge.0.001) then
    P11=f2/f1
    P01=1.0-P11
    end if
!   probabilities of fixation and loss

    ttot1=0.0
    ttot2=0.0
    Htot1=0.0
    Htot2=0.0

    N1=NH-1
!   upper limit to summation over segregating sites
    do 100 k=1,N1
    
    ak=k
    y=k/AN
    x=1.0-y

    if(gam2.le.0.001) then

    t1(k)=2.0
!   sojourn time density at k at on way to fixation
    t2(k)=2.0*x/ak
!   sojourn time density at frequency k at on way to loss
!   write(*,*) 'k= ',k,'t1= ',t1(k),' t2= ',t2(k)
    else
    f1y=(exp(gam*y)-1.0)/y
    f2x=exp(gam*x)-1.0
    f2y=f2x/x
    f3y=exp(-gam*x)
    f4y=f1*f4*x*y

    t1(k)=2.0*(f1y*f2y)/(gam*f0)
    t2(k)=2.0*f2*f3y*(f2x**2)
    t2(k)=t2(k)/(gam*f4y)
!   case of semidominance
    end if

    H1(k)=2.0*y*x*t1(k)
    H2(k)=2.0*y*x*t2(k)
!   diversity densities at frequency y

!   write(*,*) 'i= ',i,' y= ',y
!   write(*,*) 't1= ',t1(i),'H1= ',H1(i)
!   write(*,*) 't2= ',t2(i),'H2= ',H2(i)

    ttot1=ttot1+t1(k)
    ttot2=ttot2+t2(k)
!   net sojourn time densities
    Htot1=Htot1+H1(k)
    Htot2=Htot2+H2(k)
!   net diversity densities

100 continue
    
    ttot1=ttot1/AN
    ttot2=ttot2/AN
!   sojourn time integrals

!   write(*,*) 't1= ',tot1,' t2= ',t2

    Hnet=P11*Htot1+P01*Htot2
    Htot1=Htot1/AN
    Htot2=Htot2/AN
    tnet=P11*ttot1+P01*ttot2
!   net sojourn time and diversity

    end subroutine swadd

    program weaksweep10
!   program for autosomal weak sweep effect approximations- losses and fixation
!   uses summations over allele frequencies and set of h values
!   should be valid for any h or relatively small gamma value

    CHARACTER*20 FOUT
    CHARACTER*10 CONT
    integer :: N,NH,NMAX,nh1
    real :: gam,AN,ttot1,ttot2,Htot1,Htot2,P01,P11,a,b,tnet,Hnet
    write(*,*) 'Output file?'
    read(*,*)  FOUT
    OPEN (1,FILE=FOUT)

    write (*,*) 'Population size?'
    read (*,*) N
    write (*,*) 'Gamma and initial h?'
    read (*,*) gam,h1
!   gam is the scaled selection parameter, h1 is the initial dominance coefficient
    write(*,*) 'No. of h values and increment?'
    read (*,*) nh1,delh

    NH=2*N
    AN=NH

    write(1,*) 'Program for effects of fixations and losses of weakly selected mutations'
    write(1,*) 'Autosomal locus with random mating'
    write(1,*) ''
    write(1,*) 'Population size= ',N
    write(1,*) 'Initial h= ',h1,' gamma= ',gam
    write(1,*) 'No. of h values= ',nh1,'increment= ',delh
    write(*,*) 'No. of h values= ',nh1,'increment= ',delh
    write(*,*) ''
    write(*,*) 'Continue?'
    read (*,*) CONT
!  allows program to be cancelled if desired

    write(1,*) ''
    write(1,*) 'Times are in units of coalescent time(2Ne generations)'
    write(1,*) 'Diversities over fixation or loss are relative to 2Ne'
    write(*,*) ''

    do 10 ih=1,nh1

    h=h1+(ih-1)*delh
    a=h
    b=1.0-2*h
    c=a+b

    write(1,*) ''
    write(1,*) 'h= ',h
    write(*,*) ''
    write(*,*) 'h= ',h
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    timf=2.0*(1-0.05556*gam*b-0.05556*(gam**2)*(a*(a+b)+0.1967*b**2))
    write(1,*) 'Approximate mean time to fixation= ',timf
    write(*,*) 'Approximate mean time to fixation= ',timf
    write(1,*) ''
    hetf=0.6667*(1-0.06667*gam*b-0.06667*(gam**2)*(a*(a+b)+0.1905*b**2))
    hetf=hetf
    write(1,*) 'Approximate mean relative diversity over fixation= ',hetf
    write(1,*) ''
    timl=2*(log(AN)-0.4228-0.2778*gam*b-0.05556*(gam**2)*(a*(5*a+3.5*b)+0.17*b**2))
    timl=timl/AN
    write(1,*) 'Approximate mean time to loss= ',timl
    write(1,*) ''
    hetl=1.3333*(1-0.1333*gam*b-(gam**2)*(0.1333*(a**2)+0.1*a*b+0.00476*b**2))
    hetl=hetl/AN
    write(1,*) 'Approximate mean relative diversity over loss= ',hetl
    write(1,*) ''
    tsoj=2*(log(AN)+0.5772+gam*a-0.02778*b*(gam**2)*(5.0*a+b))/AN
    write(1,*) 'Approximate mean sojourn time= ',tsoj
    hetsoj=2*(1+0.3333*gam*a-0.05556*b*(gam**2)*(a+0.2*b))
    write(1,*) ''
    write(1,*) 'Approximate mean absolute diversity over sojourn= ',hetsoj
!   write(*,*) 'Approximate mean absolute diversity over sojourn= ',hetsoj
    write(1,*) ''

    b1=abs(b)
    gam2=abs(gam)
    
!   write(*,*) 'gam2= ',gam2

    if(b1.le.0.001) then
    call swadd(AN,NH,gam,ttot1,ttot2,Htot1,Htot2,tnet,Hnet,P01,P11)
!   results for semi-dominance or neutrality
    go to 5
    end if
    
    if(gam2.ge.0.001) then
    call swintegs(AN,N,NMAX,gam,a,b,ttot1,ttot2,Htot1,Htot2,tnet,Hnet,P01,P11)
    else
    call swadd(AN,NH,gam,ttot1,ttot2,Htot1,Htot2,tnet,Hnet,P01,P11)
    end if

5   write(1,*) ''
    write(1,*) 'Results of integrations'
    write(1,*) ''
    write(1,*) 'Fixation probability= ',P11
    write(1,*) 'Mean time to fixation= ',ttot1
    Htot1=Htot1
    write(1,*) 'Mean diversity over fixation= ',Htot1
    write(1,*) ''
    write(1,*) 'Loss probability= ',P01
    write(1,*) 'Mean time to loss= ',ttot2
    write(1,*) 'Mean diversity over loss= ',Htot2
    write(1,*) ''
    write(1,*) 'Mean sojourn time= ',tnet
    write(1,*) 'Mean absolute diversity over sojourn= ',Hnet
    write(1,*) ''

10   continue

    end program weaksweep10

    subroutine swintegs(AN,N,NMAX,gam,a,b,ttot1,ttot2,Htot1,Htot2,tnet,Hnet,P01,P11)
    real :: AN,a,b,y,S0,P01,P11,gam,u(0:2001),fu(0:2001),fum(0:2001),S11
    real :: psif(2000),t1(2000),t2(2000),H1(2000),H2(2000),S1k(0:2000)
    real :: ttot1,ttot2,Htot1,Htot2,tnet,Hnet
    integer :: N,N1,NMAX
!   sojourn time integral calculations

    N1=2*N
    a2=2*a

!   ap1=gam*b
!   ap=abs(ap1)
!   ap=sqrt(ap)
!   a1=a*ap/b
!   factors for transforming y to u

    ttot1=0.0
    ttot2=0.0
    Htot1=0.0
    Htot2=0.0

    do 100 k=0,N1
    
    y=k/AN
    
    u(k)=0.0-gam*y*(a2+b*y)
!   allele frequency and u value
    fu(k)=exp(u(k))
!   psi function for u value
    
!   write(*,*) 'k= ',k,'y= ',y,' u= ',u(k),' fu= ',fu(k)
!   write(*,*) ''
100 continue

    S1k(0)=0.0
    do 110 k1=1,N1
    k2=k1-1
    ub=0.5*(u(k2)+u(k1))
!   mean of u values between k1 and k1-1

    fum(k1)=exp(ub)
!   psi function for mean of u values
    
    S1k(k1)=S1k(k2)+fum(k1)
!   proportional to integral of psi from 0 to k1
110 continue
    
    S11=S1k(N1)
!   proportional to psi integral from 0 to 1
    S4=S1k(N1)-S1k(1)
!   proportional to probability of loss of new mutation
    S5=S1k(1)
!   proportional to probability of fixation of new mutation
!   write(1,*) 'Results of integrations'
    write(*,*) 'Results of integrations'
    write(*,*) 'S11= ',S11,'S4= ',S4,' S5 =',S5
    write(*,*) ''
    
    aS1=log(S11)
    aS4=log(S4)
    aS5=log(S5)

    dely=1.0/AN
    adely=log(dely)
!   change in y between successive values 

    y=0
    N2=N1-1
    do 120 i=1,N2
    y=y+dely
    x=1.0-y
    psif(i)=0.0-u(i)
!   proportional to -log of psi function

    fk=psif(i)-log(y)-log(x)
!   proportional to log of inverse of psi x y(1-y)

    S3=S1k(N1)-S1k(i)
!   proportional to psi integral from y to 1
    aS3=log(S3)

    psy=exp(-psif(i))
    pfixy=S1k(i)/S11
!   write(*,*) ''
!   write(*,*) 'i= ',i,' y= ',y,' psy= ',psy
!   write(*,*) 'S3= ',S3
!    write(*,*) ' psif= ',psif(i),' P(y)= ',pfixy,' y= ',y

    at1=log(S1k(i))-aS1+aS3+fk+adely
    t1(i)=2.0*exp(at1)
!   sojourn time density at k at on way to fixation
    
    at2=2.0*aS3+aS5-aS1-aS4+fk+adely
    t2(i)=2.0*exp(at2)
!   sojourn time density at frequency k at on way to loss

    H1(i)=2.0*y*(1-y)*t1(i)
    H2(i)=2.0*y*(1-y)*t2(i)
!   diversity densities at frequency y
    
!   write(*,*) 'i= ',i,' y= ',y
!   write(*,*) 't1= ',t1(i),'H1= ',H1(i)
!   write(*,*) 't2= ',t2(i),'H2= ',H2(i)

    ttot1=ttot1+t1(i)
    ttot2=ttot2+t2(i)
!   net sojourn time densities
    Htot1=Htot1+H1(i)
    Htot2=Htot2+H2(i)
!   net diversity densities

120 continue

    ttot1=ttot1/AN
    ttot2=ttot2/AN
!   sojourn time integrals
    
    P11=S5/S11
    P01=1.0-P11
!   probabilities of fixation and loss

    Hnet=P11*Htot1+P01*Htot2
    Htot1=Htot1/AN
    Htot2=Htot2/AN
    tnet=P11*ttot1+P01*ttot2
!   net sojourn time and diversity

    end subroutine swintegs

    subroutine swadd(AN,NH,gam,ttot1,ttot2,Htot1,Htot2,tnet,Hnet,P01,P11)
    real :: AN,a,y,S0,aifac(0:2000),ip(2000),al2(2000),P01,P11,gam
    real :: psif(2000),t1(2000),t2(2000),H1(2000),H2(2000),S2(2000),S3
    real :: ttot1,ttot2,Htot1,Htot2,tnet,Hnet
    integer :: NH
!   calculations for semidominant and neutral cases

    gam2=abs(gam)

!   write(*,*) 'swadd entered; gam2= ',gam2
    P11=1.0/AN
    P01=1.0-P11
    f0=exp(gam)-1.0
    f1=1.0-exp(0.0-gam)
    f2=1.0-exp(0.0-gam*P11)
    f3=1.0-exp(0.0-gam*P01)
    f4=exp(gam*P01)-1.0
    
    if(gam2.ge.0.001) then
    P11=f2/f1
    P01=1.0-P11
    end if
!   probabilities of fixation and loss

    ttot1=0.0
    ttot2=0.0
    Htot1=0.0
    Htot2=0.0

    N1=NH-1
!   upper limit to summation over segregating sites
    do 100 k=1,N1
    
    ak=k
    y=k/AN
    x=1.0-y

    if(gam2.le.0.001) then

    t1(k)=2.0
!   sojourn time density at k at on way to fixation
    t2(k)=2.0*x/ak
!   sojourn time density at frequency k at on way to loss
!   write(*,*) 'k= ',k,'t1= ',t1(k),' t2= ',t2(k)
    else
    f1y=(exp(gam*y)-1.0)/y
    f2x=exp(gam*x)-1.0
    f2y=f2x/x
    f3y=exp(-gam*x)
    f4y=f1*f4*x*y

    t1(k)=2.0*(f1y*f2y)/(gam*f0)
    t2(k)=2.0*f2*f3y*(f2x**2)
    t2(k)=t2(k)/(gam*f4y)
!   case of semidominance
    end if

    H1(k)=2.0*y*x*t1(k)
    H2(k)=2.0*y*x*t2(k)
!   diversity densities at frequency y

!   write(*,*) 'i= ',i,' y= ',y
!   write(*,*) 't1= ',t1(i),'H1= ',H1(i)
!   write(*,*) 't2= ',t2(i),'H2= ',H2(i)

    ttot1=ttot1+t1(k)
    ttot2=ttot2+t2(k)
!   net sojourn time densities
    Htot1=Htot1+H1(k)
    Htot2=Htot2+H2(k)
!   net diversity densities

100 continue
    
    ttot1=ttot1/AN
    ttot2=ttot2/AN
!   sojourn time integrals

!   write(*,*) 't1= ',tot1,' t2= ',t2

    Hnet=P11*Htot1+P01*Htot2
    Htot1=Htot1/AN
    Htot2=Htot2/AN
    tnet=P11*ttot1+P01*ttot2
!   net sojourn time and diversity

    end subroutine swadd

    program weaksweep11
!   program for autosomal weak sweep effect appoximations- losses and fixation
!   uses summations over allele frequencies and set of h values
!   site frequency spectra relative to neutrality are displayed
!   should be valid for any h or relatively small gamma value

    CHARACTER*20 FOUT
    CHARACTER*10 CONT
    integer :: N,NH,NMAX,nh1,inprint
    real :: gam,AN,ttot1,ttot2,Htot1,Htot2,P01,P11,a,b,tnet,Hnet
    write(*,*) 'Output file?'
    read(*,*)  FOUT
    OPEN (1,FILE=FOUT)

    write (*,*) 'Population size?'
    read (*,*) N
    write (*,*) 'Gamma and initial h?'
    read (*,*) gam,h1
!   gam is the scaled selection parameter, h1 is the initial dominance coefficient
    write(*,*) 'No. of h values and increment?'
    read (*,*) nh1,delh

    NH=2*N
    AN=NH

    write(1,*) 'Program for effects of fixations and losses of weakly selected mutations'
    write(1,*) 'Autosomal locus with random mating'
    write(1,*) 'Site frequency spectra relative to neutrality are displayed'
    write(1,*) ''
    write(1,*) 'Population size= ',N
    write(1,*) 'Initial h= ',h1,' gamma= ',gam
    write(1,*) 'No. of h values= ',nh1,'increment= ',delh
    write(*,*) 'No. of h values= ',nh1,'increment= ',delh
    write(*,*) ''
    write(*,*) 'Print interval for frequency spectrum?'
    read(*,*)  iprint
    write(1,*) 'Print interval for frequency spectrum= ',iprint
    write(*,*) 'Continue?'
    read (*,*) CONT
!  allows program to be cancelled if desired

    write(1,*) ''
    write(1,*) 'Times are in units of coalescent time(2Ne generations)'
    write(1,*) 'Diversities over fixation or loss are relative to 2Ne'
    write(*,*) ''

    do 10 ih=1,nh1

    h=h1+(ih-1)*delh
    a=h
    b=1.0-2*h
    c=a+b

    write(1,*) ''
    write(1,*) 'h= ',h
    write(*,*) ''
    write(*,*) 'h= ',h
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    timf=2.0*(1-0.05556*gam*b-0.05556*(gam**2)*(a*(a+b)+0.1967*b**2))
    write(1,*) 'Approximate mean time to fixation= ',timf
    write(*,*) 'Approximate mean time to fixation= ',timf
    write(1,*) ''
    hetf=0.6667*(1-0.06667*gam*b-0.06667*(gam**2)*(a*(a+b)+0.1905*b**2))
    hetf=hetf
    write(1,*) 'Approximate mean relative diversity over fixation= ',hetf
    write(1,*) ''
    timl=2*(log(AN)-0.4228-0.2778*gam*b-0.05556*(gam**2)*(a*(5*a+3.5*b)+0.17*b**2))
    timl=timl/AN
    write(1,*) 'Approximate mean time to loss= ',timl
    write(1,*) ''
    hetl=1.3333*(1-0.1333*gam*b-(gam**2)*(0.1333*(a**2)+0.1*a*b+0.00476*b**2))
    hetl=hetl/AN
    write(1,*) 'Approximate mean relative diversity over loss= ',hetl
    write(1,*) ''
    tsoj=2*(log(AN)+0.5772+gam*a-0.02778*b*(gam**2)*(5.0*a+b))/AN
    write(1,*) 'Approximate mean sojourn time= ',tsoj
    hetsoj=2*(1+0.3333*gam*a-0.05556*b*(gam**2)*(a+0.2*b))
    write(1,*) ''
    write(1,*) 'Approximate mean absolute diversity over sojourn= ',hetsoj
!   write(*,*) 'Approximate mean absolute diversity over sojourn= ',hetsoj
    write(1,*) ''

    b1=abs(b)
    gam2=abs(gam)
    
!   write(*,*) 'gam2= ',gam2

    if(b1.le.0.001) then
    call swadd(AN,NH,gam,a,b,ttot1,ttot2,Htot1,Htot2,tnet,Hnet,P01,P11,iprint)
!   results for semi-dominance or neutrality
    go to 5
    end if
    
    if(gam2.ge.0.001) then
    call swintegs(AN,N,NMAX,gam,a,b,ttot1,ttot2,Htot1,Htot2,tnet,Hnet,P01,P11,iprint)
    else
    call swadd(AN,NH,gam,a,b,ttot1,ttot2,Htot1,Htot2,tnet,Hnet,P01,P11,iprint)
    end if

5   write(1,*) ''
    write(1,*) 'Results of integrations'
    write(1,*) ''
    write(1,*) 'Fixation probability= ',P11
    write(1,*) 'Mean time to fixation= ',ttot1
    Htot1=Htot1
    write(1,*) 'Mean diversity over fixation= ',Htot1
    write(1,*) ''
    write(1,*) 'Loss probability= ',P01
    write(1,*) 'Mean time to loss= ',ttot2
    write(1,*) 'Mean diversity over loss= ',Htot2
    write(1,*) ''
    write(1,*) 'Mean sojourn time= ',tnet
    write(1,*) 'Mean absolute diversity over sojourn= ',Hnet
    write(1,*) ''

10   continue

    end program weaksweep11

    subroutine swintegs(AN,N,NMAX,gam,a,b,ttot1,ttot2,Htot1,Htot2,tnet,Hnet,P01,P11,iprint)
    real :: AN,a,b,y,S0,P01,P11,gam,u(0:2001),fu(0:2001),fum(0:2001),S11
    real :: psif(2000),t1(2000),t2(2000),H1(2000),H2(2000),S1k(0:2000)
    real :: ttot1,ttot2,Htot1,Htot2,tnet,Hnet,t1a(2000),t2a(2000),ttot1a,ttot1b
    integer :: N,N1,NMAX,inprint
!   sojourn time integral calculations

    N1=2*N
    a2=2*a
!   ap1=gam*b
!   ap=abs(ap1)
!   ap=sqrt(ap)
!   a1=a*ap/b
!   factors for transforming y to u

    ttot1=0.0
    ttot2=0.0
    ttot1a=0.0
    ttot2a=0.0
    Htot1=0.0
    Htot2=0.0
    
    do 100 k=0,N1
    y=k/AN
    
    u(k)=0.0-gam*y*(a2+b*y)
!   allele frequency and u value
    fu(k)=exp(u(k))
!   psi function for u value
    
!   write(*,*) 'k= ',k,'y= ',y,' u= ',u(k),' fu= ',fu(k)
!   write(*,*) ''
100 continue

    S1k(0)=0.0
    do 110 k1=1,N1
    k2=k1-1
    ub=0.5*(u(k2)+u(k1))
!   mean of u values between k1 and k1-1

    fum(k1)=exp(ub)
!   psi function for mean of u values
    
    S1k(k1)=S1k(k2)+fum(k1)
!   proportional to integral of psi from 0 to k1
110 continue
    
    S11=S1k(N1)
!   proportional to psi integral from 0 to 1
    S4=S1k(N1)-S1k(1)
!   proportional to probability of loss of new mutation
    S5=S1k(1)
!   proportional to probability of fixation of new mutation
!   write(1,*) 'Results of integrations'
!   write(*,*) 'Results of integrations'
!   write(*,*) 'S11= ',S11,'S4= ',S4,' S5 =',S5
!   write(*,*) ''
    
    aS1=log(S11)
    aS4=log(S4)
    aS5=log(S5)

    dely=1.0/AN
    adely=log(dely)
!   change in y between successive values 

    y=0
    N2=N1-1
    do 120 i=1,N2
    y=y+dely
    x=1.0-y
    psif(i)=0.0-u(i)
!   proportional to -log of psi function

    fk=psif(i)-log(y)-log(x)
!   proportional to log of inverse of psi x y(1-y)

    S3=S1k(N1)-S1k(i)
!   proportional to psi integral from y to 1
    aS3=log(S3)

    psy=exp(-psif(i))
    pfixy=S1k(i)/S11
!   write(*,*) ''
!   write(*,*) 'i= ',i,' y= ',y,' psy= ',psy
!   write(*,*) 'S3= ',S3
!    write(*,*) ' psif= ',psif(i),' P(y)= ',pfixy,' y= ',y

    at1=log(S1k(i))-aS1+aS3+fk+adely
    t1(i)=2.0*exp(at1)
!   sojourn time density at y at on way to fixation

    yxp=y*x
    t1a(i)=1-0.3333*gam*b*yxp
    t1a(i)=t1a(i)-0.3333*(gam**2)*yxp*((a**2)+0.5*a*b*(1.0+2.0*y)+0.03333*(b**2)*(1+y*(2.0+13.0*y)))
!   approximate relative sojourn time density at y at on way to fixation
    
    at2=2.0*aS3+aS5-aS1-aS4+fk+adely
    t2(i)=2.0*exp(at2)
!   sojourn time density at y at on way to loss

    t2a(i)=1.0-0.3333*gam*b*y*(2.0-y)
    fy1=(a**2)*(2.0-y)+a*b*(1+yxp)
    fy2=0.03333*(b**2)*(y+(2.0-3.0*y)*(1.0-3.0*y)+13.0*(y**2)*x)
    t2a(i)=t2a(i)-0.3333*y*(gam**2)*(fy1+fy2)
!   approximate relative sojourn time density at y at on way to loss

    H1(i)=2.0*y*(1-y)*t1(i)
    H2(i)=2.0*y*(1-y)*t2(i)
!   diversity densities at frequency y

    ttot1=ttot1+t1(i)
    ttot2=ttot2+t2(i)
    ttot1a=ttot1a+t1a(i)
    ttot2a=ttot2a+t2a(i)
!   net sojourn time densities
    Htot1=Htot1+H1(i)
    Htot2=Htot2+H2(i)
!   net diversity densities

120 continue

    icount=0
    write(1,*) ''
    write(1,*) 'Frequency spectra relative to neutrality for fixations and losses'

    do 130 i1=1,N2
    ai=i1
    y=ai/AN
    x=1.0-y
    icount=icount+1
 
    t1(i1)=t1(i1)/2.0
    t2(i1)=AN*y*t2(i1)/(2.0*x)
    t1a(i1)=t1a(i1)
    t2a(i1)=t2a(i1)
!   gives site frequency spectra relative to neutraliy

    if(icount.eq.iprint) then
    write(1,*) ''
    write(1,*) 'i= ',i1,' y= ',y
    write(1,*) 't1= ',t1(i1),' t2= ',t2(i1)
    write(1,*) 't1approx= ',t1a(i1),' t2approx= ',t2a(i1)
    icount=0
    end if

130 continue
    
    ttot1=ttot1/AN
    ttot2=ttot2/AN
!   sojourn time integrals

    P11=S5/S11
    P01=1.0-P11
!   probabilities of fixation and loss

    Hnet=P11*Htot1+P01*Htot2
    Htot1=Htot1/AN
    Htot2=Htot2/AN
    tnet=P11*ttot1+P01*ttot2
!   net sojourn time and diversity

    end subroutine swintegs

    subroutine swadd(AN,NH,gam,a,b,ttot1,ttot2,Htot1,Htot2,tnet,Hnet,P01,P11,iprint)
    real :: AN,a,y,S0,P01,P11,gam,b
    real :: psif(2000),t1(2000),t2(2000),H1(2000),H2(2000),S2(2000),S3
    real :: ttot1,ttot2,Htot1,Htot2,tnet,Hnet,t1a(2000),t2a(2000),ttot1a,ttot1b
    integer :: NH,iprint
!   calculations for semidominant and neutral cases

    gam2=abs(gam)

!   write(*,*) 'swadd entered; gam2= ',gam2
    P11=1.0/AN
    P01=1.0-P11
    f0=exp(gam)-1.0
    f1=1.0-exp(0.0-gam)
    f2=1.0-exp(0.0-gam*P11)
    f3=1.0-exp(0.0-gam*P01)
    f4=exp(gam*P01)-1.0
    
    if(gam2.ge.0.001) then
    P11=f2/f1
    P01=1.0-P11
    end if
!   probabilities of fixation and loss

    ttot1=0.0
    ttot2=0.0
    Htot1=0.0
    Htot2=0.0

    N1=NH-1
!   upper limit to summation over segregating sites
    do 100 k=1,N1
    if(k.gt.0) then
    icount=icount+1
    end if
    
    ak=k
    y=k/AN
    x=1.0-y

    if(gam2.le.0.001) then
    t1(k)=2.0
!   sojourn time density at k at on way to fixation
    t2(k)=2.0*x/ak
!   sojourn time density at frequency k at on way to loss
!   write(*,*) 'k= ',k,'t1= ',t1(k),' t2= ',t2(k)
    else
    f1y=(exp(gam*y)-1.0)/y
    f2x=exp(gam*x)-1.0
    f2y=f2x/x
    f3y=exp(-gam*x)
    f4y=f1*f4*x*y

    t1(k)=2.0*(f1y*f2y)/(gam*f0)
    t2(k)=2.0*f2*f3y*(f2x**2)
    t2(k)=t2(k)/(gam*f4y)
!   case of semidominance
    end if

    yxp=y*x
    t1a(k)=1-0.3333*gam*b*yxp
    t1a(k)=t1a(k)-0.3333*(gam**2)*yxp*((a**2)+0.5*a*b*(1.0+2.0*y)+0.03333*(b**2)*(1+y*(2.0+13.0*y)))
!   approximate sojourn time density relative to neutrality at y at on way to fixation

    t2a(k)=1-0.3333*gam*b*y*(2.0-y)
    fy1=(a**2)*(2.0-y)+a*b*(1+yxp)
    fy2=0.03333*(b**2)*(y+(2.0-3.0*y)*(1.0-3.0*y)+13.0*(y**2)*x)
    t2a(k)=t2a(k)-0.3333*y*(gam**2)*(fy1+fy2)
!   approximate sojourn time density relative to neutrality at y at on way to loss

    H1(k)=2.0*y*x*t1(k)
    H2(k)=2.0*y*x*t2(k)
!   diversity densities at frequency y

    ttot1=ttot1+t1(k)
    ttot2=ttot2+t2(k)
    ttot1a=ttot1a+t1a(k)
    ttot2a=ttot2a+t2(k)
!   net sojourn time densities
    Htot1=Htot1+H1(k)
    Htot2=Htot2+H2(k)
!   net diversity densities

100 continue
    
    icount=0
    write(1,*) ''
    write(1,*) 'Frequency spectra for fixations and losses relative to neutrality'

    do 130 i1=1,N1
    ai=i1
    y=ai/AN
    x=1.0-y
    icount=icount+1

    t1(i1)=t1(i1)/2.0
    t2(i1)=AN*y*t2(i1)/(2.0*x)
    t1a(i1)=t1a(i1)
    t2a(i1)=t2a(i1)
!   gives site frequency spectra relative to neutraliy

    if(icount.eq.iprint) then
    write(1,*) ''
    write(1,*) 'i= ',i1,' y= ',y
    write(1,*) 't1= ',t1(i1),' t2= ',t2(i1)
    write(1,*) 't1approx= ',t1a(i1),' t2approx= ',t2a(i1)
    icount=0
    end if

130 continue

    ttot1=ttot1/AN
    ttot2=ttot2/AN
!   sojourn time integrals

    Hnet=P11*Htot1+P01*Htot2
    Htot1=Htot1/AN
    Htot2=Htot2/AN
    tnet=P11*ttot1+P01*ttot2
!   net sojourn time and diversity

    end subroutine swadd

    Program weaksweep12
!   program for weak sweep effects- stationary effect of recurrent fixations and losses
!   set of gamma and h values for fixed recombination rate and N
!   locus B is neutral; A has a mutant allele A2

    CHARACTER*20 FOUT
    CHARACTER*20 FIN

    integer :: i,N,NH,ndat
    real :: theta,AN,r,rho,u,v,akap,gam1,h2

    write (*,*) 'Input file?'
    read(*,*)  FIN
    OPEN (2,file=FIN)

    write (*,*) 'Output file?'
    read(*,*)  FOUT
    OPEN (1,file=FOUT)

!   read (*,*) 'Population size?'
    read (2,*) N
!   write (*,*) 'Scaled mutation rate?'
    read (2,*) theta
!   write (*,*) 'Mutational bias?'
    read (2,*) akap
!   write (*,*) 'Rho?'
    read (2,*) rho
!   scaled selection, recombination and mutation parameters
    
!   write (*,*) 'Number of data items?'
    read (2,*) ndat

    write(1,*) 'Program for stationary effect of fixations and losses of weakly selected mutations on expected      	diversity'
    write(1,*) 'A1 is favourable and A2 is deleterious'
    write(1,*) ''
    write(1,*) 'Rho= ',rho
    write(1,*) 'theta= ',theta
    write(1,*) 'u= ',u
    write(1,*) 'Mutational bias= ',akap
    write(1,*) 'Population size= ',N
    write(1,*) ''
    write(1,*) 'No. of data items= ',ndat

    NH=2*N
    AN=NH
    
    do 10 i=1,ndat
    read(2,*) gam1,h2
!   write(1,*) gam1,h2
    read(2,*) P10
!   loss probability for A1>A2
    P11=1.0-P10
    read(2,*) delpw10,delpw11
!   net changes in relative diversity for A1>A2 loss and fixation events
    read(2,*) se10,se11
!   their standard errors
    read(2,*) P20
!   loss probability for A2>A1
    P21=1.0-P20
    read(2,*) delpw20,delpw21
!   net changes in relative diversity after A2>A1 loss and and fixation events
    read(2,*) se20,se21
!   their standard errors

!   write(*,*) ' Data file is read'

    write(1,*) 'gamma= ',gam1,' h= ',h2
    write(1,*) 'P10= ',P10,' P11= ',P11
    write(1,*) 'delpm10= ',delpw10,' delpm11= ',delpw11
    write(1,*) 'Standard errors of delpm for 1-0 and 1-1 events'
    write(1,*) se10,se11
    var10=se10**2
    var11=se11**2
    write(1,*) ''
    write(1,*) 'P20= ',P20,' P21= ',P21
    write(1,*) 'delpm20= ',delpw20,' delpw21= ',delpw21
    write(1,*) 'Standard errors of delpm for 2-0 and 2-1 events'
    write(1,*) se20,se21
    write(1,*) ''
    var20=se20**2
    var21=se21**2

    alam10=P10/P11
!   number of losses of A1>A2 before A1>A2 fixation
    alam20=P20/P21
!   number of losses of A2>A1 before A2>A1 fixation

    delp1=alam10*delpw10
    delp1=delp1+delpw11
    vdelp1=alam10*var10
    vdelp1=vdelp1+var11
    
    delp2=alam20*delpw20
    delp2=delp2+delpw21
    vdelp2=alam20*var20
    vdelp2=vdelp2+var21

    alam11=0.5*AN*theta*P11
!   substitution rate A1>A2
    alam21=0.5*AN*theta*akap*P21
!   substitution rate A2>A1

    TS=(1.0/alam11)+(1.0/alam21)
!   time between successive A1>A2 substitutions
    delp11=(delp1+delp2)/TS
!   net stationary state change in diversity

    write(1,*) 'delp1= ',delp1,' delp2= ',delp2
    delp=delp1+delp2
    vdelp=vdelp1+vdelp2
    sedelp=sqrt(vdelp)
    write(1,*) 'Net change in diversity per cycle= ',delp, 's.e.= ',sedelp
    write(1,*)  ''
    write(1,*) 'alam10= ',alam10,' alam20= ',alam20
    write(1,*) 'alam11= ',alam11,' alam21= ',alam21
    sedelp=sqrt(vdelp)/TS
    write(1,*) 'Time over whole cycle= ',TS
    write(1,*)  'Stationary change in diversity= ',delp11,'s.e.= ',sedelp
    write(1,*)  ''
    write(1,*)  ''
10  continue

    end program weaksweep12

    Program weaksweep13
!   Solution for critical h from Eq. 18 of Zhao & Charlesworth 2016

    CHARACTER*20 FOUT
    CHARACTER*10 CONT
    integer :: N,NH,NMAX,nh1
    real :: gam,AN,ttot1,ttot2,Htot1,Htot2,P01,P11,a,b,tnet,Hnet
    write(*,*) 'Output file?'
    read(*,*)  FOUT
    OPEN (1,FILE=FOUT)

    write(1,*) 'Solution for critical h from Eq. 18 of Zhao & Charlesworth 2016'
    write(1,*) 'Case of no recombination'
    write(1,*) 'Autosomal locus with random mating'
    write(1,*) ''
    write (*,*) 'initial gamma?'
    read (*,*) gam1
!   gam1 is the scaled selection parameter
    write(1,*) 'Initial gamma= ',gam1
    write(*,*) 'No. of gamma values and increment?'
    read(*,*) ng1,delg
    write(1,*) 'No. of gamma values= ',ng1,'increment= ',delg
    write(*,*) ''
    write(*,*) 'Maximum number of iterations? '
    read (*,*) NMAX
    write(1,*) 'Maximum number of iterations= ',NMAX
    write(*,*) 'Continue?'
    read (*,*) CONT
!  allows program to be cancelled if desired
    
    do 10 ig=1,ng1
    nit=0
    write(1,*) ''
    gam=gam1+(ig-1)*delg
    h0=0.5+0.1*gam
!   initial h for NR iteration
    h=h0
    a3=2.0*(gam**2)
    a2=0.0-gam*(2.0-3.0*gam)
    a1=0.0-2.0*(6.0-gam*(1.0-gam))
    a0=6.0+gam
!   coefficients for cubic equation for h

5   Fh=a3*(h**3)+a2*(h**2)+a1*h+a0
!   cubic equation for h
    
    b2=3.0*a3
    b1=2.0*a2
    b0=a1
!   coefficients for derivative of Fh
    Gh=b2*(h**2)+b1*h+b0
!   derivative of Fh
    write(*,*) 'h= ',h,' F= ',Fh,' G= ',Gh
    h2=h-0.1*Fh/Gh
    del=(h2-h)/h
    nit=nit+1
    h=h2
    if(abs(del).le.0.000001) then
    go to 20
    else
        if(nit.ge.NMAX) then
        write(1,*) 'Iteration does not converge'
        go to 20
        end if
   end if
    go to 5
    
20  write(1,*) 'gamma= ',gam,' h= ',h
    write(1,*) 'F(h)= ',Fh,' G(h)= ',Gh
    write(1,*) 'Number of iterations= ',nit
    write(*,*) ''
    write(*,*) 'Number of iterations= ',nit
    write(*,*) 'h= ',h

10  continue
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
    end program weaksweep13