******************************************************************************
*                                                                            *
      PROGRAM projec
*                                                                            *
* Fortran 77 program to calculate prevalence and incidence projections of a  *
*    chronic disease in the United States, such as Alzheimer's disease,      *
*    which affects the elderly. Projections are made for the years 1995-2050 *
*    and the potential impact of public health interventions on the          *
*    projections can be evalualated.                                         *
*                                                                            *
* See Brookmeyer and Gray (1998) for details.                                *
*                                                                            *
* Johns Hopkins University                                                   *
* School of Hygiene and Public Health                                        *
* Baltimore, Maryland                                                        *
*                                                                            *
* Last Modified: February 11, 1999                                           *
*                                                                            *
******************************************************************************
c  
c 
*****************************
* Definition of variables.  *
*****************************
c 
      INTEGER i,j,yr,yearf,yearm,reqyear
      PARAMETER (n=59,m1=114,m2=56,n2=60)
      DOUBLE PRECISION ratem(n,m1),ratef(n,m1),poplm(n,m2),poplf(n,m2)
      DOUBLE PRECISION poplm2(n2,m2),poplf2(n2,m2)
      DOUBLE PRECISION year(m2),age(n),incidm(n),incidf(n),lamm,lamf
      DOUBLE PRECISION prevm(n),prevf(n),sumpm(m2),sumpf(m2)
      DOUBLE PRECISION keeppm(n),keeppf(n)
      DOUBLE PRECISION sumim(m2),sumif(m2),prevm2(n+1),prevf2(n+1)
      DOUBLE PRECISION prev(n),summ,sumf,summ2,sumf2
      DOUBLE PRECISION thetam,thetaf,gammam,gammaf
      DOUBLE PRECISION incidm2(n,m1),incidf2(n,m1)
      REAL resp,sev
      CHARACTER*30 filem
      CHARACTER*30 filef
      CHARACTER*6 yrunits 
      INTEGER tmppm,tmppf,tmpim,tmpif
c 
c 
************************************************************************
* Read in the smoothed and projected death rate data for 1937-2050 for *
* both males and females.                                              *
************************************************************************
c 
      open(12,file='dthm.dat',status='old')
      open(13,file='dthf.dat',status='old')
      do 10 i=1,n
	read(12,*) (ratem(i,j),j=1,m1)
	read(13,*) (ratef(i,j),j=1,m1)
   10 continue
c 
c 
***********************************************************************
* Read in the "middle series" population data for 1995-2050 for both  *
* males and females (smoothed for ages 100+).                         *
***********************************************************************
c 
      open(14,file='midpopm.dat',status='old')
      open(15,file='midpopf.dat',status='old')
      do 15 i=1,n2
	read(14,*) (poplm2(i,j),j=1,56)
	read(15,*) (poplf2(i,j),j=1,56)
   15 continue
c 
      do 18 i=1,n
        do 16 j=1,m2
	  poplm(i,j)=poplm2(i+1,j)
	  poplf(i,j)=poplf2(i+1,j)
   16   continue
   18 continue
c
c 
**************************************************************
* Read in the AD incidence rates for both males and females. *
**************************************************************
c 
      print*, ' '
      print*, ' '
      print*, 'Names of incidence data files required.'
      write(*,21)
   21 format(4x,'Enter filename containing incidence data (/100) for MAL
     +ES in quotes'/,5x,'("q" to quit): ',$)
      read*, filem
      if (filem.EQ.'q') then 
        STOP
      end if
      open(16,file=filem,status='old',iostat=ios)
   22 if (ios.NE.0) then
        print*, '   Error: Unable to open the specified file.'
        write(*,23)
   23   format(4x,'Please re-enter filename with incidence data for MALE
     +S ("q" to quit): '/,5x,'',$)
        read*, filem
        if (filem.EQ.'q') then 
          STOP
        end if
        open(16,file=filem,status='old',iostat=ios)
        go to 22
      end if
c
      write(*,26)
   26 format(4x,'Enter filename containing incidence data (/100) for FEM
     +ALES in quotes'/,5x,'("q" to quit): ',$)
      read*, filef
      if (filef.EQ.'q') then 
        STOP
      end if
      open(17,file=filef,status='old',iostat=ios)
   28 if (ios.NE.0) then
        print*,'   Error: Unable to open the specified file.'
        write(*,29)
   29   format(4x,'Please re-enter filename with incidence data for FEMA
     +LES ("q" to quit): '/,5x,'',$)
        read*, filef
        if (filef.EQ.'q') then 
          STOP
        end if
        open(17,file=filef,status='old',iostat=ios)
        go to 28
      end if
      do 30 i=1,n
	read(16,*) incidm(i)
	read(17,*) incidf(i)
   30 continue
c 
c 
**********************************************************************
* Prompt user for relative risk of death, diseased vs. non-diseased. *
**********************************************************************
c 
      print*, ' '
      print*, ' '
      print*,'Relative risk of death (diseased vs. non diseased) require 
     +d. (e.g. A RR'
      print*,'of 1.5 implies that individuals with disease have mortalit
     +y rates 50%'
      print*,'higher than individuals of the same age without disease.)'
      write(*,40)
   40 format(4x,'Enter RR of death for MALES (-99 to quit): ',$)
      read*, lamm
      if (lamm.EQ.-99) then
        STOP
      end if
   41 if (lamm.LE.0) then
	print*,'   Error: Relative risk must be positive.'
        write(*,42)
   42   format(4x,'Please re-enter RR of death for MALES (-99 to quit): 
     +',$)
        read*, lamm
        if (lamm.EQ.-99) then
          STOP
        end if
        go to 41
      end if
c
      write(*,50)
   50 format(4x,'Enter RR of death for FEMALES (-99 to quit): ',$)
      read*, lamf
      if (lamf.EQ.-99) then
        STOP
      end if
   51 if (lamf.LE.0) then
	print*,'   Error: Relative risk must be positive.'
        write(*,52)
   52   format(4x,'Please re-enter RR of death for FEMALES (-99 to quit 
     +): ',$)
        read*, lamf
        if (lamf.EQ.-99) then
          STOP
        end if
        go to 51
      end if
c 
c 
*********************************************
* Prompt user for year intervention begins. *
*********************************************
c 
      print*, ' '
      print*, ' '
      print*,'Calendar year of intervention required.'
      write(*,60)
   60 format(4x,'Enter year intervention begins for MALES (-99 to quit):
     + ',$)
      read*, yearm
      if (yearm.EQ.-99) then
        STOP
      end if
   61 if (yearm.LE.1900) then
	print*,'   Error: Intervention year must be after 1900.'
        write(*,62)
   62   format(4x,'Please re-enter intervention year for MALES (-99 to q
     +uit): ',$)
        read*,yearm 
        if (yearm.EQ.-99) then
          STOP
        end if
        go to 61
      end if
      write(*,70)
c
   70 format(4x,'Enter year intervention begins for FEMALES (-99 to quit
     +): ',$)
      read*, yearf
      if (yearf.EQ.-99) then
        STOP
      end if
   71 if (yearf.LE.1900) then
	print*,'   Error: Intervention year must be after 1900.'
        write(*,72)
   72   format(4x,'Please re-enter intervention year for FEMALES (-99 to
     + quit): ',$)
        read*,yearf 
        if (yearf.EQ.-99) then
          STOP
        end if
        go to 71
      end if
c 
c 
*********************************************************************
* Prompt user for intervention value - effect on hazard of disease. *
*********************************************************************
c 
      print*,' '
      print*,' '
      print*,'Relative risk of disease (treated vs. untreated) required.
     + (e.g. A RR of'
      print*,'0.90 implies intervention reduces age-specific incidence r
     +ates by 10%.)'
      write(*,80)
   80 format(4x,'Enter RR of disease associated with the intervention fo
     +r MALES (-99'/,5x,'to quit): ',$)
      read*, thetam
      if (thetam.EQ.-99) then
        STOP
      end if
   81 if (thetam.LE.0) then
	print*,'   Error: RR must be positive.'
        write(*,82)
   82   format(4x,'Please re-enter RR of disease for MALES (-99 to quit)
     +: ', $)
        read*,thetam
        if (thetam.EQ.-99) then
          STOP
        end if
        go to 81
      end if
      write(*,83)
c
   83 format(4x,'Enter RR of disease associated with the intervention fo
     +r FEMALES (-99'/,5x,'to quit): ',$)
      read*, thetaf
      if (thetaf.EQ.-99) then
        STOP
      end if
   84 if (thetaf.LE.0) then
	print*,'   Error: RR must be positive.'
        write(*,85)
   85   format(4x,'Please re-enter RR of disease for FEMALES (-99 to qui
     +t): ',$)
        read*,thetaf
        if (thetaf.EQ.-99) then
          STOP
        end if
        go to 84
      end if
c 
c 
*************************************************************
* Prompt user for intervention value - effect on mortality. *
*************************************************************
c 
      print*,' '
      print*,' '
      print*,'Relative risk of death among those with disease (treated v
     +s. untreated)'
      print*,'required. (e.g. A RR of 0.80 implies intervention reduces 
     +disease'
      print*,'mortality by 20%.)'
      write(*,87)
   87 format(4x,'Enter RR of death among those with disease associated wi
     +th the'/,5x,'intervention for MALES (-99 to quit): ',$)
      read*, gammam
      if (gammam.EQ.-99) then
        STOP
      end if
   88 if (gammam.LE.0) then
	print*,'   Error: RR must be positive.'
        write(*,89)
   89   format(4x,'Please re-enter RR of death for MALES (-99 to quit): 
     +',$)
        read*,gammam
        if (gammam.EQ.-99) then
          STOP
        end if
        go to 88
      end if
c
      write(*,91)
   91 format(4x,'Enter RR of death among those with disease associated wi
     +th the'/,5x,'intervention for FEMALES (-99 to quit): ',$)
      read*, gammaf
      if (gammaf.EQ.-99) then
        STOP
      end if
   92 if (gammaf.LE.0) then
	print*,'   Error: RR must be positive.'
        write(*,93)
   93   format(4x,'Please re-enter RR of death for FEMALES (-99 to quit)
     +: ',$)
        read*,gammaf
        if (gammaf.EQ.-99) then
          STOP
        end if
        go to 92
      end if
c 
c 
*****************************************
* Prompt user for duration of illness.  *
*****************************************
c
      print*,' '
      print*,' '
      print*,'Duration of illness for prevalence estimates required. (e.
     +g. If you want'
      print*,'prevalence estimates of all people with disease, enter 0. 
     +If you want'
      print*,'prevalence estimates of people with disease at least 2 yea
     +rs, enter 2.)'
      write(*,95)
   95 format(4x,'Enter duration of illness (-99 to quit): ',$)
      read*,sev
      if (sev.EQ.-99) then
        STOP
      end if
   96 if (sev.LT.0.OR.sev.NE.INT(sev)) then
	print*,'   Error: Duration must be a non-negative integer.'
        write(*,97)
   97   format(4x,'Please re-enter duration of illness (-99 to quit): ',
     +$)
        read*,sev
        if (sev.EQ.-99) then
          STOP
        end if
        go to 96
      end if
c
c 
*************************************************************
* Ask user if would like to request age-specific prevalence *
* rates for a given year. If yes, request year.             *
*************************************************************
c 
      print*,' '
      print*,' '
      write(*,104)
  104 format(1x,'Would you like age-specific prevalence rates printed fo
     +r a specific'/,1x,'year? (1=yes, 0=no, -99 to quit): ',$)
      read*,resp 
      if (resp.EQ.-99) then
        STOP
      end if
  105 if (resp.NE.1.AND.resp.NE.0) then
	print*,'   Error: Invalid response.'
        write(*,106)
  106   format(4x,'Please re-enter response (1=yes, 0=no, -99 to quit): 
     +',$)  
        read*,resp 
        if (resp.EQ.-99) then
          STOP
        end if
        go to 105
      end if
      if (resp.EQ.1) then
        write(*,107)
  107   format(4x,'Enter calendar year for which you would like age-spec
     +ific prevalence'/,5x,'rates (-99 to quit): ',$)
        read*,reqyear
        if (reqyear.EQ.-99) then
          STOP
        end if
  108   if (reqyear.LE.1997.OR.reqyear.GE.2051) then
          print*,'   Error: Requested year must be between 1998 and 2050
     +.'
          write(*,109)
  109     format(4x,'Please re-enter year (-99 to quit): ',$)
          read*,reqyear
          if (reqyear.EQ.-99) then
            STOP
          end if
          go to 108
        end if
      end if
c
c
***********************
* Begin calculations. *
***********************
c 
      print*, ' '
      print*, ' '
      print*, 'Please wait...'
      print*, ' '
      print*, ' '
      do 110 i=1,59
	age(i)=59+i
  110 continue
      do 113 i=1,56
	year(i)=1994+i
  113 continue
c 
      do 120 i=1,n
	do 115 j=1,m1
	  if (j.LT.(yearm-1936)) then
	    incidm2(i,j)=incidm(i)/100
	    ratem(i,j)=ratem(i,j)/100000
	  else 
	    incidm2(i,j)=thetam*incidm(i)/100
	    ratem(i,j)=ratem(i,j)/100000
	  end if
	  if (j.LT.(yearf-1936)) then
	    incidf2(i,j)=incidf(i)/100
	    ratef(i,j)=ratef(i,j)/100000
	  else 
	    incidf2(i,j)=thetaf*incidf(i)/100
	    ratef(i,j)=ratef(i,j)/100000
	  end if
  115   continue
  120 continue
c
      do 160 yr=95,150
        if (sev.EQ.0) then
          call calcp(prev,incidm2,ratem,lamm,gammam,yr,yearm,n,m1)
        else 
          call calcpsev(prev,incidm2,ratem,lamm,gammam,yr,yearm,n,
     +                    m1,INT(sev))
        end if
c
	do 130 i=1,n
	  prevm(i)=prev(i)
  130   continue
c
        if (sev.EQ.0) then
          call calcp(prev,incidf2,ratef,lamf,gammaf,yr,yearf,n,m1)
        else
          call calcpsev(prev,incidf2,ratef,lamf,gammaf,yr,yearf,n,
     +                    m1,INT(sev))
        end if
c
	do 140 i=1,n
	  prevf(i)=prev(i)
  140   continue
c
        if (yr.EQ.(reqyear-1900)) then
	  do 145 i=1,n
	    keeppm(i)=prevm(i)
            keeppf(i)=prevf(i)
  145     continue
        end if
        prevm2(1)=0
        prevf2(1)=0
        do 147 i=1,n
          prevm2(i+1)=prevm(i)
          prevf2(i+1)=prevf(i)
  147   continue
	summ=0
	sumf=0
	summ2=0
	sumf2=0
        do 150 i=1,n
     	  summ=summ+prevm(i)*poplm(i,yr-94)
     	  sumf=sumf+prevf(i)*poplf(i,yr-94)
          if (yr.NE.95) then
            summ2=summ2+(poplm2(i,yr-94-1)-prevm2(i)*poplm2(i,yr-94-1))
     +            *incidm2(i,yr-36)
            sumf2=sumf2+(poplf2(i,yr-94-1)-prevf2(i)*poplf2(i,yr-94-1))
     +            *incidf2(i,yr-36)
	  end if
  150   continue
        sumpm(yr-94)=summ
        sumpf(yr-94)=sumf
        if (yr.NE.95) then
          sumim(yr-94)=summ2
          sumif(yr-94)=sumf2
        else
          sumim(1)=0
          sumif(1)=0
        end if
  160 continue
c 
c 
******************
* Create output. *
******************
c 
      print*, 'Output written to file projec.out.'
      print*, ' '
      open(18,file='projec.out',status='unknown')
      write(18,*), ' '
      write(18,*), 'Incidence and prevalence projections based on progra
     +m projec.f.'
      write(18,*),'Brookmeyer and Gray, 1998'
      write(18,*),'Johns Hopkins University'
      write(18,*),'School of Hygiene and Public Health'
      write(18,*),'Baltimore, Maryland'
      write(18,*),' '
      write(18,*),' '
      write(18,*),'Calculations are based on age- and gender-specific U.
     +S. mortality rates'
      write(18,*),'and U.S. Census Bureau population projections ("middl
     +e series"). In'
      write(18,*),'addition, the following user-specified values were us
     +ed:'
      write(18,*),' '
      write(18,170),filem,filef
  170 format(3x,'Incidence data file (male): ',a/,21x,'(female): ',a)
      write(18,*),' '
      write(18,190),lamm,lamf
  190 format(3x,'Relative risk of death, diseased vs. non diseased (male
     +): ',F7.4/,51x,'(female): ',F7.4)
      write(18,*),' '
      write(18,210),yearm,yearf
  210 format(3x,'Intervention year (male): ',I4/,19x,'(female): ',I4)
      write(18,*),' '
      write(18,230),thetam,thetaf
  230 format(3x,'Relative risk of disease, treated vs. untreated (male):
     + ',F7.4/,49x,'(female): ',F7.4)
      write(18,*),' '
      write(18,241),gammam,gammaf
  241 format(3x,'Relative risk of death among diseased, treated vs. untr
     +eated'/,54x,' (male): ',F7.4/,53x,'(female): ',F7.4)
      write(18,*),' '
      write(18,245),INT(sev)
  245 format(3x,'Duration of illness (years) for prevalence estimates: '
     +,I3)
      write(18,*),' '
      if (resp.EQ.1) then
        write(18,246),reqyear
  246   format(3x,'Age-specific prevalence rates requested for year: ',I
     +4)
      end if
      write(18,*),' '
      write(18,*),' '
      write(18,*),' '
      write(18,*),'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     +~~~~~~~~~~~~~~~~~~~'
      write(18,*),'              INCIDENCE PROJECTIONS            PREVAL
     +ENCE PROJECTIONS'
      if (INT(sev).EQ.1) then
        yrunits=' YEAR'
      else
        yrunits=' YEARS'
      end if
      write(18,247),INT(sev),yrunits
  247 format(48x,'FOR THOSE WITH DISEASE '/,50x,'AT LEAST ',I3,a)
      write(18,*),' Year        Male     Female    Total         Male   
     +  Female    Total'
      write(18,*),'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     +~~~~~~~~~~~~~~~~~~~'
      do 260 i=4,m2
        tmppm=NINT(sumpm(i)/10000)*10000
        tmppf=NINT(sumpf(i)/10000)*10000
        tmpim=NINT(sumim(i)/10000)*10000
        tmpif=NINT(sumif(i)/10000)*10000
        write(18,250),INT(year(i)),tmpim,tmpif,tmpim+tmpif,tmppm,tmppf,
     +tmppm+tmppf
  250 format(2x,I4,3x,3I10,3x,3I10)
      if (year(i)/10.EQ.INT(year(i)/10)) then
        write(18,*),' '
      end if
  260 continue
      write(18,*),'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     +~~~~~~~~~~~~~~~~~~~'
      write(18,*),' '
      write(18,*),' '
      write(18,*),' '
      write(18,*),' '
      if (resp.EQ.1) then
        write(18,*),'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     +~~~~~~~~~~'
        write(18,263),INT(sev),yrunits,INT(reqyear)
  263   format(2x,'ESTIMATED AGE-SPECIFIC PREVALENCE RATES (/100 PY) FOR
     + THOSE'/,7x,'WITH DISEASE AT LEAST ',I3,a,' IN THE YEAR ',I4)
        write(18,*),'               Age          Male       Female'
        write(18,*),'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     +~~~~~~~~~~'
        do 270 i=1,36
          write(18,265),INT(age(i)),100*keeppm(i),100*keeppf(i)
  265     format(16x,I3,8x,F6.2,6x,F6.2)
          if ((age(i)+1)/10.EQ.INT((age(i)+1)/10)) then
            write(18,*),' '
          end if
  270   continue
        write(18,*),'~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     +~~~~~~~~~~'
        write(18,*),' '
      end if
c
      END
******************************************************************************
c 
c 
****************************************
*      SUBROUTINE calcp                *
****************************************
c
	SUBROUTINE calcp(prev,prob,dr,lam,gam,y,iyr,rows,cols)
c
c	received parameters
	INTEGER rows,cols,y,iyr
	DOUBLE PRECISION lam,gam,prev(59),prob(59,114),dr(rows,cols)
c
c	local parameters
	INTEGER a,i,j,k,t1,t2
	DOUBLE PRECISION piece1,piece2,temp1,temp2,mult,tmp0,tmp1
c
	  DO 40 a=60,118
	    piece1=0
	    piece2=1
	      DO 30 i=60,a
   	        temp1=1 
	        temp2=1
                if (i.ne.60) then
		  DO 10 j=60,(i-1)
		    t1=j-59
		    t2=y-a+j-36
                    temp1=temp1*(1-prob(t1,t2))*(1-dr(t1,t2))
   10		  continue
		end if
		DO 20 k=i,a
		  t1=k-59
		  t2=y-a+k-36
                  if (t2.LT.(iyr-1936)) then
		    mult=lam
		  else
		    mult=gam*lam
		  end if
                  tmp0=0.0001
                  tmp1=1-mult*dr(t1,t2)
                  temp2=temp2*(max(tmp0,tmp1))
   20		continue
		t1=i-59
		t2=y-a+i-36
	        piece1=piece1+prob(t1,t2)*temp1*temp2
                piece2=piece2*(1-prob(t1,t2))*(1-dr(t1,t2))
   30	      continue 
              prev(a-59)=piece1/(piece1+piece2)
   40	  continue
	return
	END
******************************************************************************
c 
c 
****************************************
*      SUBROUTINE calcpsev             *
****************************************
c
	SUBROUTINE calcpsev(prev,prob,dr,lam,gam,y,iyr,rows,cols,sev)
c
c	received paramters
	INTEGER rows,cols,y,iyr,sev
	DOUBLE PRECISION lam,gam,prev(59),prob(59,114),dr(rows,cols)
c
c	local parameters
	INTEGER a,i,j,k,t1,t2
	DOUBLE PRECISION piece1,piece2,piece3,temp1,temp2
	DOUBLE PRECISION mult,tmp0,tmp1
c
	  DO  5 a=60,(60+sev-1)
	    prev(a-59)=0
    5	  continue
c
	  DO 100 a=(60+sev),118
	    piece1=0
	    piece2=1
	      DO 30 i=60,a
   	        temp1=1 
	        temp2=1
                if (i.ne.60) then
		  DO 10 j=60,(i-1)
		    t1=j-59
		    t2=y-a+j-36
                    temp1=temp1*(1-prob(t1,t2))*(1-dr(t1,t2))
   10		  continue
		end if
		DO 20 k=i,a
		  t1=k-59
		  t2=y-a+k-36
                  if (t2.LT.(iyr-1936)) then
                    mult=lam
                  else
                    mult=gam*lam
                  end if
                  tmp0=0.0001
                  tmp1=1-mult*dr(t1,t2)
                  temp2=temp2*(max(tmp0,tmp1))
   20		continue
		t1=i-59
		t2=y-a+i-36
	        piece1=piece1+prob(t1,t2)*temp1*temp2
                piece2=piece2*(1-prob(t1,t2))*(1-dr(t1,t2))
   30	      continue 
	    piece3=0
	      DO 90 i=60,(a-sev)
   	        temp1=1 
	        temp2=1
                if (i.ne.60) then
		  DO 40 j=60,(i-1)
		    t1=j-59
		    t2=y-a+j-36
                    temp1=temp1*(1-prob(t1,t2))*(1-dr(t1,t2))
   40		  continue
		end if
		DO 50 k=i,a
		  t1=k-59
		  t2=y-a+k-36
                  if (t2.LT.(iyr-1936)) then
                    mult=lam
                  else
                    mult=gam*lam
                  end if
                  tmp0=0.0001
                  tmp1=1-mult*dr(t1,t2)
                  temp2=temp2*(max(tmp0,tmp1))
   50		continue
		t1=i-59
		t2=y-a+i-36
	        piece3=piece3+prob(t1,t2)*temp1*temp2
   90	      continue 
              prev(a-59)=piece3/(piece1+piece2)
  100	  continue
	return
	END
******************************************************************************