C
C This is the FORTRAN version of program 3.3 from page 79 of
C "Modeling Infectious Disease in humans and animals"
C by Keeling & Rohani.
C 
C It is the SIR model with 2 age groups (children and adults).
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_3  Program_3_3.f  
C 


C Main program starts here.

      program main

      REAL*8 beta(2,2)
      REAL*8 gamma
      REAL*8 lC
      REAL*8 mu(2)
      REAL*8 S0(2),I0(2)
      REAL*8 nu , n(2)

      REAL*8 S(2),I(2)

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

      INTEGER GivesName, j, k

      CHARACTER*2000 str, FileName

      COMMON /parameters/ beta, gamma, lC, mu, nu, n
      COMMON /variables/ S, I

      GivesName=iargc()

      if (GivesName.eq.0) then
         beta(1,1)=100
         beta(1,2)=10
         beta(2,1)=10
         beta(2,2)=20
         gamma=10
         lC=0.0666667
         mu(1:2)=(/ 0.0, 0.0166667 /)
         S0(1:2)=(/ 0.1, 0.1 /)
         I0(1:2)=(/ 0.0001, 0.0001 /)
         MaxTime=100
      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,*) beta(1,1:2)
         read(1,*) beta(2,1:2)
         
         read(1,*) str
         read(1,*) gamma
         
         read(1,*) str
         read(1,*) lC
         
         read(1,*) str
         read(1,*) mu(1)

         read(1,*) str
         read(1,*) mu(2)
         
         read(1,*) str
         read(1,*) S0(1)
         
         read(1,*) str
         read(1,*) I0(1)
         
         read(1,*) str
         read(1,*) S0(2)
         
         read(1,*) str
         read(1,*) I0(2)
       
         read(1,*) str
         read(1,*) MaxTime 
         
         close(1)
      endif
C
C     Check all variables are OK & set up intitial conditions */
C     
      
      n(1)=mu(2)/(lC+mu(2))
      n(2)=1-n(1)
      nu=(lC+mu(2))*n(1);
      S(1)=S0(1)
      S(2)=S0(2)
      I(1)=I0(1)
      I(2)=I0(2)

      do j=1,2
         do k=1,2
            if ( beta(j,k).lt.0 ) then
               write(*,*) "ERROR: Transmission rate beta(",j,",",k,")=
     .              ",beta(j,k)," is less than zero."
               STOP
            endif
         enddo
         
         if ( I0(j).lt.0) then
            write(*,*) "ERROR: Initial level of infection in age 
     .           group ",j," (",I0(j),") is less than zero."
            STOP
         endif
         
         if ( S0(j).lt.0) then
            write(*,*) "ERROR: Initial level of susceptibles in age 
     .           group ",j," (",S0(j),") is less than zero."
            STOP
         endif
         
         if ( I0(j)+S0(j).gt.n(j)) then
            write(*,*) "ERROR: Initial level of susceptibles and 
     .           infected in age group ",j," (",S0(j)+I0(j),") is 
     .           greater than the calculate group size (",
     .           n(j),")."
            STOP
         endif
         
         if ( mu(j).lt.0) then
            write(*,*) "ERROR: Death rate in age 
     .           group ",j," (",mu(j),") is less than zero."
            STOP
         endif
         
      enddo
      
      if ( gamma.le.0) then
         write(*,*) "ERROR: Recovery rate gamma (",gamma,") is 
     .        less than or equal to zero."
         STOP
      endif

      if ( lC.le.0) then
         write(*,*) "ERROR: Maturing rate l_C (",lC,") 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

      
C     
C     Find a suitable time-scale for outputs 
C     
      step=0.01/((beta(1,1)+beta(1,2)+beta(2,1)+beta(2,2)+gamma+nu))

      Every=1.0/((beta(1,1)+beta(1,2)+beta(2,1)+beta(2,2)+gamma+nu))
      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(1),I(1),S(2),I(2)
         endif
      
      ENDDO
      
      END
  


      SUBROUTINE Runge_Kutta(step)

      REAL*8 InitialPop(4), tmpPop(4)
      REAL*8 dPop1(4), dPop2(4), dPop3(4), dPop4(4)

      REAL*8 S(2),I(2),step
      COMMON /variables/ S,I

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(1)
      InitialPop(2)=I(1)
      InitialPop(3)=S(2)
      InitialPop(4)=I(2)
 
      CALL Diff(InitialPop,dPop1)
      do k=1,4
         tmpPop(k)=InitialPop(k)+step*dPop1(k)/2.0
      ENDDO 

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

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

      CALL Diff(tmpPop,dPop4)

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

      S(1)=tmpPop(1)
      I(1)=tmpPop(2)
      S(2)=tmpPop(3)
      I(2)=tmpPop(4)
      
      RETURN
      END


 
C The Main Differential Equations.

      SUBROUTINE Diff(Pop, dPop)

      REAL*8 Pop(4), dPop(4)
      REAL*8 beta(2,2)
      REAL*8 gamma
      REAL*8 lC
      REAL*8 mu(2)
      REAL*8 nu , n(2)

      COMMON /parameters/ beta, gamma, lC, mu, nu, n


C     Set up temporary variables to make the equations look neater

      REAL*8 tmpSC, tmpSA, tmpIC, tmpIA

      tmpSC=Pop(1)
      tmpIC=Pop(2)
      tmpSA=Pop(3)
      tmpIA=Pop(4)
      
C
C     The differential equations
C

C     dSC/dt =
      dPop(1) = nu - (beta(1,1)*tmpIC + beta(1,2)*tmpIA)*tmpSC
     .     - mu(1)*tmpSC - lC*tmpSC

C     dIC/dt =
      dPop(2) = (beta(1,1)*tmpIC + beta(1,2)*tmpIA)*tmpSC
     .     - gamma*tmpIC - mu(1)*tmpIC - lC*tmpIC

C     dSA/dt =
      dPop(3) = lC*tmpSC - (beta(2,1)*tmpIC + beta(2,2)*tmpIA)*tmpSA
     .     - mu(2)*tmpSA

C     dIA/dt =
      dPop(4) = lC*tmpIC + (beta(2,1)*tmpIC + beta(2,2)*tmpIA)*tmpSA
     .     - gamma*tmpIA - mu(2)*tmpIA


      RETURN
      END