C
C This is the FORTRAN version of program 3.5 from page 94 of
C "Modeling Infectious Disease in humans and animals"
C by Keeling & Rohani.
C 
C It is the S(E)IR model with multiple stages to create gamma-distributed exposed and infectious periods
C
C n is the total number of infection classes; classes 1 to m are exposed, classes m+1 to n are infectious.
C
C Birth and death rates are assumed equal.
C
C This code is written to be simple, transparent and readily compiled.
C Far more elegant and efficient code can be written.
C
C This code can be compiled using the intel fortran compiler:
C ifort -Vaxlib  -o Program_3_5  Program_3_5.f  
C 


C Main program starts here.

      program main

      IMPLICIT NONE


      INTEGER  MaxGroups, j
      PARAMETER(MaxGroups=200)

      REAL*8 beta
      REAL*8 gamma
      REAL*8 mu
      INTEGER n,m

      REAL*8 S0
      REAL*8 I0

      REAL*8 S, I(MaxGroups)

      REAL*8 MaxTime
      REAL*8 EVERY, step, t

      INTEGER GivesName

      CHARACTER*2000 str, FileName

      COMMON /parameters/ beta, gamma, mu
      COMMON /basics/ n, m
      COMMON /variables/ S, I


      GivesName=iargc()

      if (GivesName.eq.0) then
         beta=17.0/5.0
         gamma=1.0/13.0
         mu=1.0/(55.0*365.0)
         n=13
         m=8
         S0=0.05
         I0=1.0d-6
         MaxTime=30*365
      else
 
c     
c     READ IN ALL THE VARIABLES
c

         call getarg(1,FileName)

         open(1,file=FileName,STATUS='OLD',ACCESS='SEQUENTIAL')
  
         
         read(1,*) str
         read(1,*) n

         if (n.gt.MaxGroups) then
            write(*,*) "ERROR: Number of risk groups ",n," exceeds the 
     .           preset maximum"
            STOP
         endif
         
         read(1,*) str
         read(1,*) m

         read(1,*) str
         read(1,*) beta

         read(1,*) str
         read(1,*) gamma

         read(1,*) str
         read(1,*) mu
         
         read(1,*) str
         read(1,*) S0
         
         read(1,*) str
         read(1,*) I0
         
         read(1,*) str
         read(1,*) MaxTime 
         
         close(1)
      endif
C
C     Check all variables are OK & set up intitial conditions */
C     
      
      if ( S0.le.0) then
         write(*,*) "ERROR: Initial level of susceptibles (",S0,") is 
     .        less than or equal to zero."
         STOP
      endif

      if ( I0.le.0) then
         write(*,*) "ERROR: Initial level of infectious (",I0,") is 
     .        less than or equal to zero."
         STOP
      endif

      if ( beta.le.0) then
         write(*,*) "ERROR: Transmission rate beta (",beta,") is 
     .        less than or equal to zero."
         STOP
      endif

      if ( gamma.le.0) then
         write(*,*) "ERROR: Recovery rate gamma (",gamma,") is 
     .        less than or equal to zero."
         STOP
      endif

      if ( mu.le.0) then
         write(*,*) "ERROR: Recovery rate mu (",mu,") is 
     .        less than or equal to zero."
         STOP
      endif

      if ( MaxTime.le.0) then
         write(*,*) "ERROR: Maximum run time (",MaxTime,") is 
     .        less than or equal to zero."
         STOP
      endif

     
      if (n.le.0) then
         write(*,*) "ERROR: Number of infected stages (",n,") is less
     .        than or equal to zero."
         STOP
      endif

        if (m.ge.n) then
         write(*,*) "ERROR: Number of infectious stages (",m-n,") is less
     .        than or equal to zero."
         write(*,*) "m is probably too big."
         STOP
      endif    

      S=S0
      do j=1,n
         I(j)=I0/n
      enddo
          
C     
C     Find a suitable time-scale for outputs 
C     
      step=0.01/((beta+gamma*n+mu))

      Every=1.0/((beta+gamma*n+mu))
      if (Every.gt.1) then
         Every=10.0**INT(log10(Every))
      else
         Every=10.0**INT(log10(Every)-1)
      endif
      DO WHILE (MaxTime/Every.gt.10000)
         Every=Every*10.0  
      ENDDO

      open(1,recl=3000,file='Output',ACCESS='SEQUENTIAL')
C     for F77 use
C      open(1,file='Output_Risk',ACCESS='SEQUENTIAL')
C

C
C     The main iteration routine
C     

      t=0

      DO WHILE (t.lt.MaxTime)
         

         CALL Runge_Kutta(step)
         
         t=t+step
         
C     If time has moved on sufficiently, output the current data
      
         if( INT(t/Every).gt.INT((t-step)/Every) ) then
            write(1,*) t,S,I(1:n)
         endif
      
      ENDDO
      
      END
  


      SUBROUTINE Runge_Kutta(step)

      IMPLICIT NONE


      INTEGER MaxGroups
      PARAMETER(MaxGroups=200)

      REAL*8 step

      REAL*8 InitialPop(MaxGroups+1), tmpPop(MaxGroups+1)
      REAL*8 dPop1(MaxGroups+1), dPop2(MaxGroups+1)
      REAL*8 dPop3(MaxGroups+1), dPop4(MaxGroups+1)

      REAL*8 S,I(MaxGroups)
      INTEGER n,m,k

      COMMON /variables/ S, I
      COMMON /basics/ n, m

C
C     Integrates the equations one step, using Runge-Kutta 4 
C     Note: we work with arrays rather than variables to make the
C     coding easier
C
      InitialPop(1)=S
      do k=1,n
         InitialPop(1+k)=I(k)
      enddo

      CALL Diff(InitialPop,dPop1)
      do k=1,n+1
         tmpPop(k)=InitialPop(k)+step*dPop1(k)/2.0
      ENDDO 

      CALL Diff(tmpPop,dPop2)
      do k=1,n+1
         tmpPop(k)=InitialPop(k)+step*dPop2(k)/2.0
      ENDDO

      CALL Diff(tmpPop,dPop3)
      do k=1,n+1
         tmpPop(k)=InitialPop(k)+step*dPop2(k)
      ENDDO

      CALL Diff(tmpPop,dPop4)

      do k=1,n+1
         tmpPop(k)=InitialPop(k)+step*(dPop1(k)/6 + dPop2(k)/3 +
     .        dPop3(k)/3 + dPop4(k)/6)
      ENDDO

      S=tmpPop(1)
      do k=1,n
         I(k)=tmpPop(1+k)
      enddo

      RETURN
      END


 
C The Main Differential Equations.

      SUBROUTINE Diff(Pop, dPop)

      IMPLICIT NONE


      INTEGER  MaxGroups, k
      PARAMETER(MaxGroups=200)

      REAL*8 Pop(MaxGroups+1), dPop(MaxGroups+1)

      REAL*8 beta
      REAL*8 gamma
      REAL*8 mu
      REAL*8 Inf
      INTEGER n,m

      COMMON  /parameters/ beta, gamma, mu
      COMMON /basics/ n, m

C
C     Calculate proportion infectious
C
      Inf=0.0
      do k=(m+1),n
         Inf=Inf+Pop(k+1)
      enddo

C
C     The differential equations
C

C     dS/dt =
      dPop(1) = mu - beta*Inf*Pop(1) - mu*Pop(1)

C     dI_1/dt =
      dPop(2) = beta*Inf*Pop(1) - gamma*n*Pop(2) - mu*Pop(2)

C     dI_k/dt =
      do k=2,n
         dPop(k+1) =  gamma*n*Pop(k) - gamma*n*Pop(k+1) - mu*Pop(k+1)
      enddo

      RETURN
      END