program ipums_2000_totpop

!	Tetyana Dubovyk
!	July 17, 2003
!	Revised on January 12, 2004
!
!     INPUT:    census data are organized in the following data files:
!				td-1950-totpop.dat
!				td-1960-totpop.dat
!				td-1970-totpop.dat
!				td-1980-totpop.dat
!				td-1990-totpop.dat
!				td-2000-totpop.dat
!				The layout of each of the above files is the same:
!				There are 192 rows and columns are as follows:
!                     Field      Type   Size   Definition
!                     -----      ----   ----   ----------
!                     Sex        C      5      0=Total population;
!											   1 = Male; 
!											   2 = Female
!                     Agegrp     C      5      0 = Total (15 yrs and over)
!                                              1 = 15 to 24 years 
!                                              2 = 25 to 34 years 
!                                              3 = 35 to 44 years 
!                                              4 = 45 to 54 years 
!                                              5 = 55 to 64 years 
!                                              6 = 65 to 74 years 
!                                              7 = 75 years and over
!                     Mstat      C      5      0 = Total
!                                              1 = Single
!                                              2 = Married,spouse present
!                                              3 = Married,s.p.,youngest<6
!                                              4 = Married,s.p.,youngest 6-17
!                                              5 = Married,spouse absent
!                                              6 = Widowed
!                                              7 = Divorced
!                     Totpop     N      9      Total population             (1)
!                     Empl       N      9      Number employed              (2)
!                     Totpt      N      9      Total employed part-time     (3)
!                     Hrs1_14    N      9      Employed 1  to 14 hrs/week   (4)
!                     Hrs15_29   N      9      Employed 15 to 29 hrs/week   (5)
!                     Hrs30_34   N      9      Employed 30 to 34 hrs/week   (6)
!                     Totft      N      9      Total employed full-time     (7)
!                     Hrs35_39   N      9      Employed 35 to 39 hrs/week   (8)
!                     Hrs40      N      9      Employed 40 hours            (9)
!                     Hrs41_48   N      9      Employed 41 to 48 hrs/week  (10)
!                     Hrs49_59   N      9      Employed 49 to 59 hrs/week  (11)
!                     Hrs60pl    N      9      Employed 60+ hours          (12)
!
!				To incorporate scaling of hours for 2000, the following 
!				input files are used:
!					matrix_h80.dat
!					matrix_h90.dat
!					matrix_h00.dat
!				The layout of each of the above files is the same:
!				There are 192 rows and columns are as follows:
!					1 = population in group i
!					2 = total hours for group i calculated with hrswork2
!					3 = average hours for group i based on hrswork2
!					  = column 2/column 1
!					4 = total hours for group i calculated with uhrswork
!					5 = average hours for group i based on uhrswork
!					  = column 4/column 1
!
!     OUTPUT:         ipums_5per_totpop.out:  Tables summarizing hours, employment 
!                                for various characteristics 
!                                (e.g., age, sex, marital status).
!					  ipums-extrap_totpop.out: The same tables as in ipums_5per_totpop.out
!						but with extrapolated values
!					  paper_tables.out: Tables that are included into published
!						version of the update of QRw98
!					  ipums_extrap_totpop.m: Matlab file for plotting life-cycle
!						profiles for different groups
!					  figure1sr317.m: Matlab file to plot figure 1 from
!						SR 317 (paper by McGrattan, Jones, Manuelli)
!					  comparison_TS.m: Matlab file to plot the comparison of
!						time-series for average hours calculated based on
!						hrswork2 for 1950-1990 and uhrswork for 1980-2000					  
!
!     REFERENCES:     Integrated Public Use Microdata Series , Bureau of the Census
!						www.ipums.org
!
!                     U S. Department of Commerce, Bureau of the Census,
!                     Subject Reports, Employment Status and Work 
!                     Experience 1950,1960,1970,1990,2000
!
!     CONTACTS:       Matt Sobek 
!                     Minnesota Population Center -- University of Minnesota 
!                     (612) 627-4509
!
!                     Wendy Treadwell 
!                     Machine Readable Data Center -- Univ. Minnesota
!                     (612) 624-4389
!     Key matrices:
!
!     ryr(i,j) = real-valued characteristics for group i, year yr
!		  j=1:	 hours (hours = reported hours x #empl./#reporting hours)	
!         j=2:   hours per person  
!         j=3:   hours per worker
!         j=4:   employed/population in group i in %
!         j=5:   cell population/total population in group i, in %
!
!     hyr(i,j) = real-valued characteristics for group i, year yr
!		  j=1:	 hours calculated using hrswork2
!				(hours = reported hours from hrswork2 x #empl./#reporting hours)	
!		  j=2:   average hours, hrswork2
!		  j=3:	 hours calculated using uhrswork
!				(hours = reported hours from uhrswork x #empl./#reporting hours)	
!		  j=4:   average hours, uhrswork
!         j=5:   scaling coefficient found using 1990
!				 = h90(i,1)/h90(i,3)
!				 for 2000 this column = 0
!         j=6:   for 1900 this column = 0
!				 for 2000 this column = scaled hours from method 4
!									  = h00(i,3)*h90(i,5)
! 
!     iyr(i,j) = integer-valued characteristics for group i, year yr
!         j=1:   total population
!         j=2:   number employed
!         j=3:   number employed  1-14 hours
!         j=4:   number employed 15-29 hours
!         j=5:   number employed 30-34 hours
!         j=6:   number employed 35-39 hours
!         j=7:   number employed  40   hours
!         j=8:   number employed 41-48 hours
!         j=9:   number employed 49-59 hours
!         j=10:  number employed  60+  hours
!         j=11:  number reporting hours
!
!     imar(i,j,b,k) = integer-valued characteristics for different
!				groups of marital status, year k
!		  i=1:	 Married S.P.
!		  i=2:   Married S.A.
!		  i=3:   Single
!		  i=4:   Widowed
!		  i=5:   Divorced
!
!         j=1:   population in the group
!         j=2:   number employed
!         j=3:   number employed  1-14 hours
!         j=4:   number employed 15-29 hours
!         j=5:   number employed 30-34 hours
!         j=6:   number employed 35-39 hours
!         j=7:   number employed  40   hours
!         j=8:   number employed 41-48 hours
!         j=9:   number employed 49-59 hours
!         j=10:  number employed  60+  hours
!         j=11:  number reporting hours
!
!		  b=1:   total population
!		  b=2:   males
!		  b=3:   females
!
!		  k=1:	 year 1950 and so on
!         k=6:   year 2000
!
!     rmar(i,j,b,k) = real-valued characteristics for different
!				groups of marital status, year k
!		  i=1:	 Married S.P.
!		  i=2:   Married S.A.
!		  i=3:   Single
!		  i=4:   Widowed
!		  i=5:   Divorced
!
!		  j=1:	 hours (hours = reported hours x #empl./#reporting hours)	
!         j=2:   hours per person  
!         j=3:   hours per worker
!         j=4:   employed/population in group i in %
!         j=5:   cell population/total population in group i, in %
!
!		  b=1:   total population
!		  b=2:   males
!		  b=3:   females
!
!		  k=1:	 year 1950 and so on
!         k=6:   year 2000
!
!	cat=category:
!	   =1 if we print Table 2: average hours by age
!	   =2 if we print Table 3: hours per worker by age
!	   =3 if we print Table 4: employment/population by age
!	   =4 if we print Table 5: cell population/total population by age


	implicit integer(i,j,k) 
	integer								:: cat,dog,b,rat,count
    real, parameter						:: epsilon=10e-8
	real(kind = 8), dimension(192,5)	:: r50, r60, r70, r80, r90, r00
	real(kind = 8), dimension(192,5,6)	:: r
	real(kind = 8), dimension(192,6)	:: h80,h90, h00
	real(kind = 8), dimension(5,5,3,6)	:: rmar
	real(kind = 8)						:: tab(6,7),tab1(6,3),tab2(6,5),&
										   h(5,3,3), tab3(5),tab4(3),&
										   tabfig(6,4,4)
	real(kind = 8)						:: tpop50,tpop60,tpop70,tpop80,&
										   tpop90,tpop00,pop50,pop60,pop70,&
										   pop80,pop90,pop00
	real(kind = 8), dimension(6)		:: tpop
	integer, dimension(192,11)			:: i50, i60, i70, i80, i90, i00
	integer, dimension(5,11,3,6)		:: imar
	integer, dimension(192)				:: s, a, m
	character*70						:: title,label,thead,note
	character*3							:: name
! real(kind = 8) refers to real variable of double precision.
! tab(6,7) contains numbers printed into output tables,
!		7 = number of age groups, 15 yrs and up,
!		6 = number of census years.
! tab1(6,3) contains numbers printed into table with aggregate statistics,
!		3 = number of aggregate statistics:
!				average hours,
!				hours/worker,
!				employment/population,
!		6 = number of census years.
! s = sex, a = age, m = marital status

	! output files
	open(unit=2, file='ipums_extrap_totpop.m')
	open(unit=3, file='ipums-extrap_totpop.out')
	open(unit=4, file='ipums_5per_totpop.out')
	open(unit=5, file='comparison_TS.m')
	open(unit=30, file='matrix_r70.dat')
	open(unit=31, file='matrix_r00.dat')
	open(unit=32, file='paper_tables.out')
	open(unit=33, file='figure1sr317.m')
	open(unit=34, file='figure1sr317_4graphs.m')

	! input files
	open(unit=7, file='td-1950-totpop.dat')
	open(unit=8, file='td-1960-totpop.dat')
	open(unit=9, file='td-1970-totpop.dat')
	open(unit=10, file='td-1980-totpop.dat')
	open(unit=11, file='td-1990-totpop.dat')
	open(unit=12, file='td-2000-totpop.dat')
	open(unit=20, file='matrix_h80.dat')
	open(unit=21, file='matrix_h90.dat')
	open(unit=22, file='matrix_h00.dat')



! Read in IPUMS data files

	do i = 1,192
		read(7,'(I6, 2I5, 12I15)')  s(i),a(i),m(i),(i50(i,j),j=1,2),i4, &
                  (i50(i,j),j=3,5),i5,(i50(i,j),j=6,10)
        read(8,'(I6, 2I5, 12I15)')  i1,i2,i3,      (i60(i,j),j=1,2),i4, &
                  (i60(i,j),j=3,5),i5,(i60(i,j),j=6,10)
		read(9,'(I6, 2I5, 12I15)')  i1,i2,i3,      (i70(i,j),j=1,2),i4, &
                  (i70(i,j),j=3,5),i5,(i70(i,j),j=6,10)
        read(10,'(I6, 2I5, 12I15)')  i1,i2,i3,      (i80(i,j),j=1,2),i4, &
                  (i80(i,j),j=3,5),i5,(i80(i,j),j=6,10)
        read(11,'(i6, 2I5, 12I15)')  i1,i2,i3,      (i90(i,j),j=1,2),i4, &
                  (i90(i,j),j=3,5),i5,(i90(i,j),j=6,10)
        read(12,'(i6, 2I5, 12I15)')  i1,i2,i3,      (i00(i,j),j=1,2),i4, &
                  (i00(i,j),j=3,5),i5,(i00(i,j),j=6,10)
		read(20,'(I12,4F17.3)') i1,(h80(i,j),j=1,4) ! hours for methods 2 and 4
		read(21,'(I12,4F17.3)') i1,(h90(i,j),j=1,4) ! hours for methods 2 and 4
		read(22,'(I12,4F17.3)') i1,(h00(i,j),j=1,4) ! hours for methods 2 and 4
	
		! Calculate number of people reporting hours
		i50(i,11) = i50(i,3) + i50(i,4) + i50(i,5) + i50(i,6) + i50(i,7) + &
					i50(i,8) + i50(i,9) + i50(i,10)
		i60(i,11) = i60(i,3) + i60(i,4) + i60(i,5) + i60(i,6) + i60(i,7) + &
					i60(i,8) + i60(i,9) + i60(i,10)
		i70(i,11) = i70(i,3) + i70(i,4) + i70(i,5) + i70(i,6) + i70(i,7) + &
					i70(i,8) + i70(i,9) + i70(i,10)
		i80(i,11) = i80(i,3) + i80(i,4) + i80(i,5) + i80(i,6) + i80(i,7) + &
					i80(i,8) + i80(i,9) + i80(i,10)
		i90(i,11) = i90(i,3) + i90(i,4) + i90(i,5) + i90(i,6) + i90(i,7) + &
					i90(i,8) + i90(i,9) + i90(i,10)		 
		i00(i,11) = i00(i,3) + i00(i,4) + i00(i,5) + i00(i,6) + i00(i,7) + &
					i00(i,8) + i00(i,9) + i00(i,10)		 
	end do

! Calculate columns of matrices h90 and h00
	do i=1,192
		h90(i,5) = h90(i,1)/h90(i,3) ! hrswork2/uhrswork
		h00(i,6) = h00(i,3)*h90(i,5)
	end do

	print *, "h90(1,5)=",h90(1,5), "    h90(9,5)=",h90(9,5)

! For each year calculate total population = population of males + pop of females

	tpop50 = dfloat(i50(1,1)); tpop(1) = tpop50
	tpop60 = dfloat(i60(1,1)); tpop(2) = tpop60
	tpop70 = dfloat(i70(1,1)); tpop(3) = tpop70
	tpop80 = dfloat(i80(1,1)); tpop(4) = tpop80
	tpop90 = dfloat(i90(1,1)); tpop(5) = tpop90
	tpop00 = dfloat(i00(1,1)); tpop(6) = tpop00
	! dfloat(x) converts integer x into Real(8) (real of double precision)

!	print *, "i50(1,11) = ", i50(1,11), "  tpop50 = ", tpop50

! Calculate variables of matrix ryr
! Hours for group i (i=1,128) are defined as:
!            hours(i) = (7.5*[#emp 1-14](i) + ...+62.5*[#emp 60+](i))
!                       x  #employed(i)/#reporting hours(i) 

      do  i=1,192
		! Calculate hours for each group
        r50(i,1)  = (7.5d0*dfloat(i50(i,3)) & 
                   +22.0d0*dfloat(i50(i,4)) &
                   +32.0d0*dfloat(i50(i,5)) &
                   +37.0d0*dfloat(i50(i,6)) &
                   +40.0d0*dfloat(i50(i,7)) &
                   +44.5d0*dfloat(i50(i,8)) &
                   +54.0d0*dfloat(i50(i,9)) &
                   +62.5d0*dfloat(i50(i,10))) &
                   * (dfloat(i50(i,2))+epsilon)/(dfloat(i50(i,11))+epsilon)

        r60(i,1)  = (7.5d0*dfloat(i60(i,3)) & 
                   +22.0d0*dfloat(i60(i,4)) &
                   +32.0d0*dfloat(i60(i,5)) &
                   +37.0d0*dfloat(i60(i,6)) &
                   +40.0d0*dfloat(i60(i,7)) &
                   +44.5d0*dfloat(i60(i,8)) &
                   +54.0d0*dfloat(i60(i,9)) &
                   +62.5d0*dfloat(i60(i,10))) &
                   * dfloat(i60(i,2))/(dfloat(i60(i,11))+epsilon)

        r70(i,1)  = (7.5d0*dfloat(i70(i,3)) & 
                   +22.0d0*dfloat(i70(i,4)) &
                   +32.0d0*dfloat(i70(i,5)) &
                   +37.0d0*dfloat(i70(i,6)) &
                   +40.0d0*dfloat(i70(i,7)) &
                   +44.5d0*dfloat(i70(i,8)) &
                   +54.0d0*dfloat(i70(i,9)) &
                   +62.5d0*dfloat(i70(i,10))) &
                   * dfloat(i70(i,2))/(dfloat(i70(i,11))+epsilon)

        r80(i,1)  = (7.5d0*dfloat(i80(i,3)) & 
                   +22.0d0*dfloat(i80(i,4)) &
                   +32.0d0*dfloat(i80(i,5)) &
                   +37.0d0*dfloat(i80(i,6)) &
                   +40.0d0*dfloat(i80(i,7)) &
                   +44.5d0*dfloat(i80(i,8)) &
                   +54.0d0*dfloat(i80(i,9)) &
                   +62.5d0*dfloat(i80(i,10))) &
                   * dfloat(i80(i,2))/(dfloat(i80(i,11))+epsilon)

        r90(i,1)  = (7.5d0*dfloat(i90(i,3)) & 
                   +22.0d0*dfloat(i90(i,4)) &
                   +32.0d0*dfloat(i90(i,5)) &
                   +37.0d0*dfloat(i90(i,6)) &
                   +40.0d0*dfloat(i90(i,7)) &
                   +44.5d0*dfloat(i90(i,8)) &
                   +54.0d0*dfloat(i90(i,9)) &
                   +62.5d0*dfloat(i90(i,10))) &
                   * dfloat(i90(i,2))/(dfloat(i90(i,11))+epsilon)

        r00(i,1)  = h00(i,6)
        
		! Calculate average hours = hours/ population in the group
		r50(i,2)  = r50(i,1)/(dfloat(i50(i,1))+epsilon)
        r60(i,2)  = r60(i,1)/(dfloat(i60(i,1))+epsilon)
        r70(i,2)  = r70(i,1)/(dfloat(i70(i,1))+epsilon)
        r80(i,2)  = r80(i,1)/(dfloat(i80(i,1))+epsilon)
        r90(i,2)  = r90(i,1)/(dfloat(i90(i,1))+epsilon)
        r00(i,2)  = r00(i,1)/(dfloat(i00(i,1))+epsilon)
        
		! Calculate hours per worker = hours/ # employed
		r50(i,3)  = r50(i,1)/(dfloat(i50(i,2))+epsilon)
        r60(i,3)  = r60(i,1)/(dfloat(i60(i,2))+epsilon)
        r70(i,3)  = r70(i,1)/(dfloat(i70(i,2))+epsilon)
        r80(i,3)  = r80(i,1)/(dfloat(i80(i,2))+epsilon)
        r90(i,3)  = r90(i,1)/(dfloat(i90(i,2))+epsilon)
        r00(i,3)  = r00(i,1)/(dfloat(i00(i,2))+epsilon)
        
		! Calculate employment/population ratio for each group in %
		r50(i,4)  = (dfloat(i50(i,2))/(dfloat(i50(i,1))+epsilon))*100
        r60(i,4)  = (dfloat(i60(i,2))/(dfloat(i60(i,1))+epsilon))*100
        r70(i,4)  = (dfloat(i70(i,2))/(dfloat(i70(i,1))+epsilon))*100
        r80(i,4)  = (dfloat(i80(i,2))/(dfloat(i80(i,1))+epsilon))*100
        r90(i,4)  = (dfloat(i90(i,2))/(dfloat(i90(i,1))+epsilon))*100
        r00(i,4)  = (dfloat(i00(i,2))/(dfloat(i00(i,1))+epsilon))*100

		! Calculate cell population/total population ratio for each group
		r50(i,5)  = (dfloat(i50(i,1))/tpop50)*100
        r60(i,5)  = (dfloat(i60(i,1))/tpop60)*100
        r70(i,5)  = (dfloat(i70(i,1))/tpop70)*100
        r80(i,5)  = (dfloat(i80(i,1))/tpop80)*100
        r90(i,5)  = (dfloat(i90(i,1))/tpop90)*100
        r00(i,5)  = (dfloat(i00(i,1))/tpop00)*100
	end do

	! Matrix r for all years
	do i=1,192
		do j=1,5
			r(i,j,1) = r50(i,j)
			r(i,j,2) = r60(i,j)
			r(i,j,3) = r70(i,j)
			r(i,j,4) = r80(i,j)
			r(i,j,5) = r90(i,j)
			r(i,j,6) = r00(i,j)
		end do
	end do

	!Check matrices r70 and r00
	do i=1,192
	write(30,'(3I2,5F17.4)') s(i),a(i),m(i),(r50(i,j),j=1,5)
	write(31,'(5F17.4)') (r00(i,j),j=1,5)
	end do

! Print tables with aggregate statistics for different groups into ipums.out.
	title='TABLE 1: AGGREGATE STATISTICS '
	write(4,'(A42)') title
	write(4,*) 
	write(4,*) 
	label='a. TOTAL '

	call table2(label,r(1,2:4,1),r(1,2:4,2),r(1,2:4,3), &
			r(1,2:4,4), r(1,2:4,5),r(1,2:4,6))

	! Print output file with tables that go into published paper
	title='TABLE 1 '
	write(32,'(A42)') title
	title='A LOOK BEHIND AN AGGREGATE FACT '
	write(32,'(A42)') title
	title='In the United States, 1950-2000 '
	write(32,'(A42)') title
	write(32,*) 
	call table2p1(r(1,2:4,1),r(1,2:4,2),r(1,2:4,3), &
			r(1,2:4,4), r(1,2:4,5),r(1,2:4,6))



	do i = 1,192
		IF ((s(i).eq.1) .and. (a(i).eq.0) .and. (m(i).eq.0)) THEN
			label='b. ALL MALES '
			call table2(label,r50(i,2:4),r60(i,2:4),r70(i,2:4),r80(i,2:4), &
							r90(i,2:4),r00(i,2:4))
		elseif ((s(i).eq.2) .and. (a(i).eq.0) .and. (m(i).eq.0)) then
			label='c. ALL FEMALES '
			call table2(label,r50(i,2:4),r60(i,2:4),r70(i,2:4),r80(i,2:4), &
							r90(i,2:4),r00(i,2:4))
		endif
	end do
	do i=1,192
		if ((s(i).eq.1) .and. (a(i).eq.0) .and. (m(i).eq.1)) then
			label='d. SINGLE MALES '
			call table2(label,r50(i,2:4),r60(i,2:4),r70(i,2:4),r80(i,2:4), &
							r90(i,2:4),r00(i,2:4))
 		elseif ((s(i).eq.1) .and. (a(i).eq.0) .and. (m(i).eq.2)) then
			label='e. MARRIED MALES, SPOUSE PRESENT '
			call table2(label,r50(i,2:4),r60(i,2:4),r70(i,2:4),r80(i,2:4), &
							r90(i,2:4),r00(i,2:4))
 		elseif ((s(i).eq.1) .and. (a(i).eq.0) .and. (m(i).eq.5)) then
			label='f. MARRIED MALES, SPOUSE ABSENT '
			call table2(label,r50(i,2:4),r60(i,2:4),r70(i,2:4),r80(i,2:4), &
							r90(i,2:4),r00(i,2:4))
 		elseif ((s(i).eq.1) .and. (a(i).eq.0) .and. (m(i).eq.6)) then
			label='g. WIDOWED MALES '
			call table2(label,r50(i,2:4),r60(i,2:4),r70(i,2:4),r80(i,2:4), &
							r90(i,2:4),r00(i,2:4))
 		elseif ((s(i).eq.1 ) .and. (a(i).eq.0) .and. (m(i).eq.7)) then
			label='h. DIVORCED MALES '
			call table2(label,r50(i,2:4),r60(i,2:4),r70(i,2:4),r80(i,2:4), &
							r90(i,2:4),r00(i,2:4))
		ENDIF
	end do
	do i = 1,192
 		IF ((s(i).eq.2) .and. (a(i).eq.0) .and. (m(i).eq.1)) then
			label='i. SINGLE FEMALES '
			call table2(label,r50(i,2:4),r60(i,2:4),r70(i,2:4),r80(i,2:4), &
							r90(i,2:4),r00(i,2:4))
 		elseif ((s(i).eq.2) .and. (a(i).eq.0) .and. (m(i).eq.2)) then
			label='j. MARRIED FEMALES, SPOUSE PRESENT '
			call table2(label,r50(i,2:4),r60(i,2:4),r70(i,2:4),r80(i,2:4), &
							r90(i,2:4),r00(i,2:4))
 		elseif ((s(i).eq.2) .and. (a(i).eq.0) .and. (m(i).eq.3)) then
			label='k. MARRIED FEMALES, S.P., YOUNGEST < 6 '
			call table2(label,r50(i,2:4),r60(i,2:4),r70(i,2:4),r80(i,2:4), &
							r90(i,2:4),r00(i,2:4))
 		elseif ((s(i).eq.2) .and. (a(i).eq.0) .and. (m(i).eq.4)) then
			label='l. MARRIED FEMALES, S.P., YOUNGEST 6-17 '
			call table2(label,r50(i,2:4),r60(i,2:4),r70(i,2:4),r80(i,2:4), &
							r90(i,2:4),r00(i,2:4))
 		elseif ((s(i).eq.2) .and. (a(i).eq.0) .and. (m(i).eq.5)) then
			label='m. MARRIED FEMALES, SPOUSE ABSENT '
			call table2(label,r50(i,2:4),r60(i,2:4),r70(i,2:4),r80(i,2:4), &
							r90(i,2:4),r00(i,2:4))
  		elseif ((s(i).eq.2) .and. (a(i).eq.0) .and. (m(i).eq.6)) then
			label='n. WIDOWED FEMALES '
			call table2(label,r50(i,2:4),r60(i,2:4),r70(i,2:4),r80(i,2:4), &
							r90(i,2:4),r00(i,2:4))
   		elseif ((s(i).eq.2) .and. (a(i).eq.0) .and. (m(i).eq.7)) then
			label='o. DIVORCED FEMALES '
			call table2(label,r50(i,2:4),r60(i,2:4),r70(i,2:4),r80(i,2:4), &
							r90(i,2:4),r00(i,2:4))
		END IF
	end do

	! writing into file ipums-extrap.m
		write(2,'(''a=20:10:80;'')')

		! writing into file paper_tables.out
		title='TABLES 8-10'
		write(32,'(A42)') title
		title='EXTRAPOLATED LIFE-CYCLE PROFILES OF HOURS WORKED'
		write(32,'(A70)') title
		title='U.S. Census Data Extrapolated as Explained in Appendix D'
		write(32,'(A70)') title
		write(32,*) 


!	cat=category:
!	   =1 if we print Table 2: average hours by age
!	   =2 if we print Table 3: hours per worker by age
!	   =3 if we print Table 4: employment/population by age
!	   =4 if we print Table 5: cell population/total population by age
do cat = 1,4
	if (cat.eq.1) then
		title='TABLE 2: AVERAGE HOURS '
		thead='          Average Hours by Age'
		name ='ave'
	elseif (cat.eq.2) then
		title='TABLE 3: HOURS PER WORKER '
		thead='          Hours/Worker by Age'
		name ='hou'
	elseif (cat.eq.3) then
		title='TABLE 4: EMPLOYMENT/POPULATION '
		thead='     Employment/Population by Age, in %'
		name ='emp'
	elseif (cat.eq.4) then
		title='TABLE 5: CELL POPULATION/TOTAL POPULATION '
		thead='          % Population by Age'
		name ='cel'
	endif

	! writing into file ipums.out
	write(4,'(A42)') title
	write(4,*) 
	write(4,*)
	
	! writing into file ipums-extrap.out
		write(3,'(A42)') title
		write(3,*) 
		write(3,*)

		
	label='a. TOTAL '
		do i=1,7
		tab(1,i) = r(i+1,cat+1,6) ! year 2000
		tab(2,i) = r(i+1,cat+1,5)
		tab(3,i) = r(i+1,cat+1,4)
		tab(4,i) = r(i+1,cat+1,3)
		tab(5,i) = r(i+1,cat+1,2)
		tab(6,i) = r(i+1,cat+1,1) ! year 1950
		end do	
	k = 9 ! total table
	i = 9 ! total table
	count=15*(cat-1)+1 ! variable is used to number figures in ipums_extrap.m
	call table1(tab,label,thead)
		call extrapolation(tab,label,thead,name,k,i,cat,count)

	do k=1,2  ! correspond to sex
		if (k.eq.1) then 
			label='b. ALL MALES '
		elseif (k.eq.2) then 
			label='c. ALL FEMALES '
		endif
			do i=1,192
				if ((s(i).eq.k) .and. (m(i).eq.0) .and. (a(i).ne.0)) then
					tab(1,a(i)) = r00(i,cat+1)
					tab(2,a(i)) = r90(i,cat+1)
					tab(3,a(i)) = r80(i,cat+1)
					tab(4,a(i)) = r70(i,cat+1)
					tab(5,a(i)) = r60(i,cat+1)
					tab(6,a(i)) = r50(i,cat+1)
				endif
			end do	
			call table1(tab,label,thead)
				j=0 ! correspond to age group=all
				count=count+1
				call extrapolation(tab,label,thead,name,k,j,cat,count)
	end do

	do k=1,2  ! correspond to sex 
		do j=1,7  ! correspond to marital status
			! insert labels
			if ((k.eq.1).and.(j.eq.1)) then
				label='d. SINGLE MALES '
			elseif ((k.eq.1).and.(j.eq.2)) then
				label='e. MARRIED MALES, SPOUSE PRESENT '
			elseif ((k.eq.1).and.(j.eq.3)) then
				goto 70
			elseif ((k.eq.1).and.(j.eq.4)) then
				goto 70
			elseif ((k.eq.1).and.(j.eq.5)) then
				label='f. MARRIED MALES, SPOUSE ABSENT '
			elseif ((k.eq.1).and.(j.eq.6)) then
				label='g. WIDOWED MALES '
			elseif ((k.eq.1).and.(j.eq.7)) then
				label='h. DIVORCED MALES '
			elseif ((k.eq.2).and.(j.eq.1)) then
				label='i. SINGLE FEMALES '
			elseif ((k.eq.2).and.(j.eq.2)) then
				label='j. MARRIED FEMALES, SPOUSE PRESENT '
			elseif ((k.eq.2).and.(j.eq.3)) then
				label='k. MARRIED FEMALES, S.P., YOUNGEST < 6 '
			elseif ((k.eq.2).and.(j.eq.4)) then
				label='l. MARRIED FEMALES, S.P., YOUNGEST 6-17 '
			elseif ((k.eq.2).and.(j.eq.5)) then
				label='m. MARRIED FEMALES, SPOUSE ABSENT '
			elseif ((k.eq.2).and.(j.eq.6)) then
				label='n. WIDOWED FEMALES '
			elseif ((k.eq.2).and.(j.eq.7)) then
				label='o. DIVORCED FEMALES '
			endif
			
			do i=1,192
				if ((s(i).eq.k) .and. (m(i).eq.j) .and. (a(i).ne.0)) then
					tab(1,a(i)) = r00(i,cat+1)
					tab(2,a(i)) = r90(i,cat+1)
					tab(3,a(i)) = r80(i,cat+1)
					tab(4,a(i)) = r70(i,cat+1)
					tab(5,a(i)) = r60(i,cat+1)
					tab(6,a(i)) = r50(i,cat+1)
				endif
			end do
			call table1(tab,label,thead)
				count=count+1
				call extrapolation(tab,label,thead,name,k,j,cat,count)
70		continue
		end do
	end do
end do ! corresponds to cat=1,4

! Tables 2,3,4 in QRw98
! cat=category
!    =1: average hours
!    =2: hours per worker
!    =3: employment/population ratio

! Table 2 in QRw98: A distribution of cat By Sex
	title='TABLE 6: DISTRIBUTION BY SEX '
	write(4,'(A42)') title
	title='This is table 2 in QRw98 '
	write(4,'(A42)') title
	write(4,*) 
	write(4,*) 
	do cat=1,3
		if (cat.eq.1) then
			label='a. Average Hours '
			thead='Average Hours'
		elseif (cat.eq.2) then
			label='b. Hours per Worker '
			thead='Hours per Worker'
		elseif (cat.eq.3) then
			label='c. Employment-Population ratio '
			thead='Employment/Population, in %'
		endif
		do k=1,6
			tab1(k,1) = r(1,cat+1,k) ! total population column
			tab1(k,2) = r(65,cat+1,k) ! males column
			tab1(k,3) = r(129,cat+1,k) ! females column
		end do
		call table6(label,thead,tab1)		
	end do ! corresponds to cat=1,3

! Table 2 in update of QRw98: A distribution of cat By Sex
	title='TABLE 2 '
	write(32,'(A42)') title
	title='DISAGGREGATING FIRST BY SEX '
	write(32,'(A42)') title
	title='Average Weekly Hours Worked per Person '
	write(32,'(A50)') title
	title='and per Worker and the employment-to-population ratio '
	write(32,'(A70)') title
	write(32,*) 
	do cat=1,3
		if (cat.eq.1) then
			thead='a. Average Hours'
		elseif (cat.eq.2) then
			thead='b. Hours per Worker'
		elseif (cat.eq.3) then
			thead='c. Employment-to-Population Ratio, in %'
		endif
		do k=1,6
			tab1(k,1) = r(1,cat+1,k) ! total population column
			tab1(k,2) = r(65,cat+1,k) ! males column
			tab1(k,3) = r(129,cat+1,k) ! females column
		end do
		call table6p1(thead,tab1)		
	end do ! corresponds to cat=1,3


! Table 3 in QRw98: A distribution of cat By Age
	title='TABLE 7: DISTRIBUTION BY AGE '
	write(4,'(A42)') title
	title='This is table 3 in QRw98 '
	write(4,'(A42)') title
	write(4,*) 

do dog=1,3
!  dog=1: total population
!     =2: males
!     =3: females
	if (dog.eq.1) then
		title='Table 7.1: Distribution by age for TOTAL POPULATION '
		write(4,'(A60)') title
		write(4,*) 
	elseif (dog.eq.2) then
		title='Table 7.2: Distribution by age for MALES '
		write(4,'(A60)') title
		write(4,*) 
	elseif (dog.eq.3) then
		title='Table 7.3: Distribution by age for FEMALES '
		write(4,'(A60)') title
		write(4,*) 
	endif
	do cat=1,3
		if (cat.eq.1) then
			label='a. Average Hours '
			thead='Average Hours'
		elseif (cat.eq.2) then
			label='b. Hours per Worker '
			thead='Hours per Worker'
		elseif (cat.eq.3) then
			label='c. Employment-Population ratio '
			thead='Employment/Population, in %'
		endif
		if (dog.eq.1) then
			do k=1,6
				do i=1,7
				tab(k,i) = r(i+1,cat+1,k) ! matrix for total population 
				end do
			end do
		elseif (dog.eq.2) then
			do k=1,6
				do i=1,7
				tab(k,i) = r(i+65,cat+1,k) ! matrix for males 
				end do
			end do
		elseif (dog.eq.3) then
			do k=1,6
				do i=1,7
				tab(k,i) = r(i+129,cat+1,k) ! matrix for females 
				end do
			end do
		endif
		call table7(label,thead,tab)		
	end do ! corresponds to cat=1,3
end do ! corresponds to dog=1,3

! Table 3-5 in update of QRw98: A distribution of cat By Age
	title='TABLES 3-5 '
	write(32,'(A42)') title
	title='AND THEN BY AGE '
	write(32,'(A42)') title
	write(32,*) 

do dog=1,3
!  dog=1: total population
!     =2: males
!     =3: females
	if (dog.eq.1) then
		title='Table 3: TOTAL  '
		write(32,'(A60)') title
		write(4,*) 
	elseif (dog.eq.2) then
		title='Table 4: MALES '
		write(32,'(A60)') title
		write(32,*) 
	elseif (dog.eq.3) then
		title='Table 5: FEMALES '
		write(32,'(A60)') title
		write(32,*) 
	endif
	do cat=1,3
		if (cat.eq.1) then
			thead='          a. Weekly Hours Worked per Person by Age (in Years) '
		elseif (cat.eq.2) then
			thead='          b. Weekly Hours Worked per Worker by Age (in Years) '
		elseif (cat.eq.3) then
			thead='          c. Employment-to-Population Ratio by Age, in % (in Years) '
		endif
		if (dog.eq.1) then
			do k=1,6
				do i=1,7
				tab(k,i) = r(i+1,cat+1,k) ! matrix for total population 
				end do
			end do
		elseif (dog.eq.2) then
			do k=1,6
				do i=1,7
				tab(k,i) = r(i+65,cat+1,k) ! matrix for males 
				end do
			end do
		elseif (dog.eq.3) then
			do k=1,6
				do i=1,7
				tab(k,i) = r(i+129,cat+1,k) ! matrix for females 
				end do
			end do
		endif
		call table7p1(thead,tab)		
	end do ! corresponds to cat=1,3
end do ! corresponds to dog=1,3

! Insert printing of tables 6-10 into paper_tables.out here
	title='TABLES 6-7 '
	write(32,'(A42)') title
	title='A MORE COMPREHENSIVE DISTRIBUTION OF HOURS WORKED PER PERSON * '
	write(32,'(A70)') title
	title='Average Weekly Hours Worked per Person '
	write(32,'(A50)') title
	title='for Sets of Demographic Categories in the United States, 1950-2000 '
	write(32,'(A70)') title
	write(32,*)
	title='TABLE 6: MARRIED... '
	write(32,'(A70)') title
	write(32,*)
	thead='          Weekly Hours Worked per Person by Age (in Years)'

	do i=1,6
		do j=1,7
			tab(i,j) = 0.d0
		end do
	end do
	do j=2,4  ! correspond to MARRIED SPOUSE PRESENT
		do k=1,2  ! correspond to sex 
			! insert labels
			if ((k.eq.1).and.(j.eq.2)) then
				label='a. MARRIED MALES, SPOUSE PRESENT '
				write(32,'(A70)') label
			elseif ((k.eq.1).and.(j.eq.3)) then
				goto 80
			elseif ((k.eq.1).and.(j.eq.4)) then
				goto 80
			elseif ((k.eq.2).and.(j.eq.2)) then
				label='b. MARRIED FEMALES, SPOUSE PRESENT '
				write(32,'(A70)') label
			elseif ((k.eq.2).and.(j.eq.3)) then
				label='c. MARRIED FEMALES, S.P., YOUNGEST < 6 '
				write(32,'(A70)') label
			elseif ((k.eq.2).and.(j.eq.4)) then
				label='d. MARRIED FEMALES, S.P., YOUNGEST 6-17 '
				write(32,'(A70)') label
			endif
			
			do i=1,192
				if ((s(i).eq.k) .and. (m(i).eq.j) .and. (a(i).ne.0)) then
					do rat=1,6
						tab(rat,a(i)) = r(i,2,rat)
					end do
				endif
			end do
		call table7p1(thead,tab)
80		continue
		end do
	end do
	! correspond to MARRIED SPOUSE ABSENT
		do k=1,2  ! correspond to sex 
			! insert labels
			if (k.eq.1) then
				label='e. MARRIED MALES, S.A. '
				write(32,'(A70)') label
			elseif (k.eq.2) then
				label='f. MARRIED FEMALES, S.A. '
				write(32,'(A70)') label
			endif
			
			do i=1,192
				if ((s(i).eq.k) .and. (m(i).eq.5) .and. (a(i).ne.0)) then
					do rat=1,6
						tab(rat,a(i)) = r(i,2,rat)
					end do
				endif
			end do
			call table7p1(thead,tab)
		end do
	title='* See website for a decomposition of Hours Worked per Person '
	write(32,'(A70)') title
	title='into per Worker and Employment-to-Population Ratio'
	write(32,'(A70)') title
	write(32,*)
	write(32,*)

	title='TABLE 7: ...AND NOT MARRIED '
	write(32,'(A70)') title
	write(32,*)
	thead='          Weekly Hours Worked per Person by Age (in Years)  '

	do i=1,6
		do j=1,7
			tab(i,j) = 0.d0
		end do
	end do
	! correspond to SINGLE
		do k=1,2  ! correspond to sex 
			! insert labels
			if (k.eq.1) then
				label='a. SINGLE MALES '
				write(32,'(A70)') label
			elseif (k.eq.2) then
				label='b. SINGLE FEMALES '
				write(32,'(A70)') label
			endif
			
			do i=1,192
				if ((s(i).eq.k) .and. (m(i).eq.1) .and. (a(i).ne.0)) then
					do rat=1,6
						tab(rat,a(i)) = r(i,2,rat)
					end do
				endif
			end do
			call table7p1(thead,tab)
		end do
	do j=6,7  ! correspond to WIDOWED AND DIVORCED
		do k=1,2  ! correspond to sex 
			! insert labels
			if ((k.eq.1).and.(j.eq.6)) then
				label='c. WIDOWED MALES '
				write(32,'(A70)') label
			elseif ((k.eq.2).and.(j.eq.6)) then
				label='d. WIDOWED FEMALES '
				write(32,'(A70)') label
			elseif ((k.eq.1).and.(j.eq.7)) then
				label='e. DIVORCED MALES '
				write(32,'(A70)') label
			elseif ((k.eq.2).and.(j.eq.7)) then
				label='F. DIVORCED FEMALES  '
				write(32,'(A70)') label
			endif
			
			do i=1,192
				if ((s(i).eq.k) .and. (m(i).eq.j) .and. (a(i).ne.0)) then
					do rat=1,6
						tab(rat,a(i)) = r(i,2,rat)
					end do
				endif
			end do
		call table7p1(thead,tab)
		end do
	end do



! Table 4 in QRw98: A distribution of cat By Marital Status
! Calculate matrices imar and rmar (for description of matrices, see 
!	the beginning of the program)
! NOTE: for table 4, we consider age 25+
do rat=1,2
! rat=1: age group 15+
!    =2: age group 25+
	do i=1,5
	do j=1,11
		do l=1,3
		do k=1,6
			imar(i,j,l,k)=0
		end do
		end do
	end do
	end do
	! h(i,j,k) = real-valued matrix of scaling coefficients for 2000
	! the matrix is constructed from the matrix h90
	!	i= marital status dimension: 5 groups of marital status
	!   j=1: total population
	!	j=2: males
	!	j=3: females
	!	k=1: hours from method 2
	!	k=2: hours from method 4
	!	k=3: scaling coefficients
	!		 h(:,:,3) = h(:,:,1)/h(:,:,2)
	do i=1,5
		do j=1,3
		do k=1,3
			h(i,j,k)=0.d0
		end do
		end do
	end do
	do i=1,64 ! this for layer=1: total population
		if ((m(i).eq.1).and.(a(i).ge.rat)) then ! Single
			do j=1,11
				imar(3,j,1,1)=imar(3,j,1,1) + i50(i,j)
				imar(3,j,1,2)=imar(3,j,1,2) + i60(i,j)
				imar(3,j,1,3)=imar(3,j,1,3) + i70(i,j)
				imar(3,j,1,4)=imar(3,j,1,4) + i80(i,j)
				imar(3,j,1,5)=imar(3,j,1,5) + i90(i,j)
				imar(3,j,1,6)=imar(3,j,1,6) + i00(i,j)
				h(3,1,1) = h(3,1,1) + h90(i,1) ! hours for method 2
				h(3,1,2) = h(3,1,2) + h90(i,3) ! hours for method 4
			end do
		elseif ((m(i).eq.2).and.(a(i).ge.rat)) then ! Married S.P.
			do j=1,11
				imar(1,j,1,1)=imar(1,j,1,1) + i50(i,j)
				imar(1,j,1,2)=imar(1,j,1,2) + i60(i,j)
				imar(1,j,1,3)=imar(1,j,1,3) + i70(i,j)
				imar(1,j,1,4)=imar(1,j,1,4) + i80(i,j)
				imar(1,j,1,5)=imar(1,j,1,5) + i90(i,j)
				imar(1,j,1,6)=imar(1,j,1,6) + i00(i,j)
				h(1,1,1) = h(1,1,1) + h90(i,1)
				h(1,1,2) = h(1,1,2) + h90(i,3)
			end do
		elseif ((m(i).eq.5).and.(a(i).ge.rat)) then ! Married S.A.
			do j=1,11
				imar(2,j,1,1)=imar(2,j,1,1) + i50(i,j)
				imar(2,j,1,2)=imar(2,j,1,2) + i60(i,j)
				imar(2,j,1,3)=imar(2,j,1,3) + i70(i,j)
				imar(2,j,1,4)=imar(2,j,1,4) + i80(i,j)
				imar(2,j,1,5)=imar(2,j,1,5) + i90(i,j)
				imar(2,j,1,6)=imar(2,j,1,6) + i00(i,j)
				h(2,1,1) = h(2,1,1) + h90(i,1)
				h(2,1,2) = h(2,1,2) + h90(i,3)
			end do
		elseif ((m(i).eq.6).and.(a(i).ge.rat)) then ! Widowed
			do j=1,11
				imar(4,j,1,1)=imar(4,j,1,1) + i50(i,j)
				imar(4,j,1,2)=imar(4,j,1,2) + i60(i,j)
				imar(4,j,1,3)=imar(4,j,1,3) + i70(i,j)
				imar(4,j,1,4)=imar(4,j,1,4) + i80(i,j)
				imar(4,j,1,5)=imar(4,j,1,5) + i90(i,j)
				imar(4,j,1,6)=imar(4,j,1,6) + i00(i,j)
				h(4,1,1) = h(4,1,1) + h90(i,1)
				h(4,1,2) = h(4,1,2) + h90(i,3)
			end do
		elseif ((m(i).eq.7).and.(a(i).ge.rat)) then ! Divorced
			do j=1,11
				imar(5,j,1,1)=imar(5,j,1,1) + i50(i,j)
				imar(5,j,1,2)=imar(5,j,1,2) + i60(i,j)
				imar(5,j,1,3)=imar(5,j,1,3) + i70(i,j)
				imar(5,j,1,4)=imar(5,j,1,4) + i80(i,j)
				imar(5,j,1,5)=imar(5,j,1,5) + i90(i,j)
				imar(5,j,1,6)=imar(5,j,1,6) + i00(i,j)
				h(5,1,1) = h(5,1,1) + h90(i,1)
				h(5,1,2) = h(5,1,2) + h90(i,3)
			end do
		endif
	end do ! i=1,64

do k=1,2 ! sex=Male, Female
	     ! this is for layers 2,3: males, females respectively
	do i=1,192 ! this for layer=1: total population
		if ((s(i).eq.k).and.(m(i).eq.1).and.(a(i).ge.rat)) then ! Single
			do j=1,11
				imar(3,j,1+k,1)=imar(3,j,1+k,1) + i50(i,j)
				imar(3,j,1+k,2)=imar(3,j,1+k,2) + i60(i,j)
				imar(3,j,1+k,3)=imar(3,j,1+k,3) + i70(i,j)
				imar(3,j,1+k,4)=imar(3,j,1+k,4) + i80(i,j)
				imar(3,j,1+k,5)=imar(3,j,1+k,5) + i90(i,j)
				imar(3,j,1+k,6)=imar(3,j,1+k,6) + i00(i,j)
				h(3,1+k,1) = h(3,1+k,1) + h90(i,1)
				h(3,1+k,2) = h(3,1+k,2) + h90(i,3)
			end do
		elseif ((s(i).eq.k).and.(m(i).eq.2).and.(a(i).ge.rat)) then ! Married S.P.
			do j=1,11
				imar(1,j,1+k,1)=imar(1,j,1+k,1) + i50(i,j)
				imar(1,j,1+k,2)=imar(1,j,1+k,2) + i60(i,j)
				imar(1,j,1+k,3)=imar(1,j,1+k,3) + i70(i,j)
				imar(1,j,1+k,4)=imar(1,j,1+k,4) + i80(i,j)
				imar(1,j,1+k,5)=imar(1,j,1+k,5) + i90(i,j)
				imar(1,j,1+k,6)=imar(1,j,1+k,6) + i00(i,j)
				h(1,1+k,1) = h(1,1+k,1) + h90(i,1)
				h(1,1+k,2) = h(1,1+k,2) + h90(i,3)
			end do
		elseif ((s(i).eq.k).and.(m(i).eq.5).and.(a(i).ge.rat)) then ! Married S.A.
			do j=1,11
				imar(2,j,1+k,1)=imar(2,j,1+k,1) + i50(i,j)
				imar(2,j,1+k,2)=imar(2,j,1+k,2) + i60(i,j)
				imar(2,j,1+k,3)=imar(2,j,1+k,3) + i70(i,j)
				imar(2,j,1+k,4)=imar(2,j,1+k,4) + i80(i,j)
				imar(2,j,1+k,5)=imar(2,j,1+k,5) + i90(i,j)
				imar(2,j,1+k,6)=imar(2,j,1+k,6) + i00(i,j)
				h(2,1+k,1) = h(2,1+k,1) + h90(i,1)
				h(2,1+k,2) = h(2,1+k,2) + h90(i,3)
			end do
		elseif ((s(i).eq.k).and.(m(i).eq.6).and.(a(i).ge.rat)) then ! Widowed
			do j=1,11
				imar(4,j,1+k,1)=imar(4,j,1+k,1) + i50(i,j)
				imar(4,j,1+k,2)=imar(4,j,1+k,2) + i60(i,j)
				imar(4,j,1+k,3)=imar(4,j,1+k,3) + i70(i,j)
				imar(4,j,1+k,4)=imar(4,j,1+k,4) + i80(i,j)
				imar(4,j,1+k,5)=imar(4,j,1+k,5) + i90(i,j)
				imar(4,j,1+k,6)=imar(4,j,1+k,6) + i00(i,j)
				h(4,1+k,1) = h(4,1+k,1) + h90(i,1)
				h(4,1+k,2) = h(4,1+k,2) + h90(i,3)
			end do
		elseif ((s(i).eq.k).and.(m(i).eq.7).and.(a(i).ge.rat)) then ! Divorced
			do j=1,11
				imar(5,j,1+k,1)=imar(5,j,1+k,1) + i50(i,j)
				imar(5,j,1+k,2)=imar(5,j,1+k,2) + i60(i,j)
				imar(5,j,1+k,3)=imar(5,j,1+k,3) + i70(i,j)
				imar(5,j,1+k,4)=imar(5,j,1+k,4) + i80(i,j)
				imar(5,j,1+k,5)=imar(5,j,1+k,5) + i90(i,j)
				imar(5,j,1+k,6)=imar(5,j,1+k,6) + i00(i,j)
				h(5,1+k,1) = h(5,1+k,1) + h90(i,1)
				h(5,1+k,2) = h(5,1+k,2) + h90(i,3)
			end do
		endif
	end do ! i=1,192
end do ! k=1,2: sex

	!Calculate scaling coefficients used for 2000
	h(:,:,3) = h(:,:,1)/h(:,:,2)
	print *, "scaling coefficients"
	print *, "h(1,1,3)=", h(1,1,3)
	print *, "h(1,2,3)=", h(1,2,3)
	print *, "h(1,3,3)=", h(1,3,3)

	do i=1,5
	do j=1,5
		do l=1,3
		do k=1,6
			rmar(i,j,l,k)=0.d0
		end do
		end do
	end do
	end do
	do k=1,6 ! year-dimension
		do i=1,5 ! marital status-dimension
        rmar(i,1,:,k)  = (7.5d0*dfloat(imar(i,3,:,k)) & 
                   +22.0d0*dfloat(imar(i,4,:,k)) &
                   +32.0d0*dfloat(imar(i,5,:,k)) &
                   +37.0d0*dfloat(imar(i,6,:,k)) &
                   +40.0d0*dfloat(imar(i,7,:,k)) &
                   +44.5d0*dfloat(imar(i,8,:,k)) &
                   +54.0d0*dfloat(imar(i,9,:,k)) &
                   +62.5d0*dfloat(imar(i,10,:,k))) &
                   * dfloat(imar(i,2,:,k))/dfloat(imar(i,11,:,k))
		! Scaling of hours for 2000, I use average for scaling coefficients
		! for males and females
		if (k.eq.6) then
			rmar(i,1,:,6) = rmar(i,1,:,6)*h(i,:,3)
		endif
		! Calculate average hours = hours/ population in the group
		rmar(i,2,:,k)  = rmar(i,1,:,k)/dfloat(imar(i,1,:,k))
		! Calculate hours per worker = hours/ # employed
		rmar(i,3,:,k)  = rmar(i,1,:,k)/dfloat(imar(i,2,:,k))
		! Calculate employment/population ratio for each group in %
		rmar(i,4,:,k)  = (dfloat(imar(i,2,:,k))/dfloat(imar(i,1,:,k)))*100
		! Calculate cell population/total population ratio for each group
		rmar(i,5,:,k)  = (dfloat(imar(i,1,:,k))/tpop(k))*100
		end do ! correspond to i: marital status-dimension
	end do ! correspond to k: year-dimension



	title='TABLE 8: DISTRIBUTION BY MARITAL STATUS '
	write(4,'(A42)') title
	title='This is table 4 in QRw98 '
	write(4,'(A42)') title
	if (rat.eq.1) then
		note='Note: This includes age group 15 plus'
		write(4,'(A42)') note
	elseif (rat.eq.2) then
		note='Note: This includes age group 25 plus'
		write(4,'(A42)') note
	endif
	write(4,*) 
	write(4,*) 

do dog=1,3
!  dog=1: total population
!     =2: males
!     =3: females
	if (dog.eq.1) then
		title='Table 8.1: Distr. by marital status for TOTAL POPULATION '
		write(4,'(A60)') title
		write(4,*) 
	elseif (dog.eq.2) then
		title='Table 8.2: Distribution by marital status for MALES '
		write(4,'(A60)') title
		write(4,*) 
	elseif (dog.eq.3) then
		title='Table 8.3: Distribution by marital status for FEMALES '
		write(4,'(A60)') title
		write(4,*) 
	endif
	do cat=1,4
	! cat=1: average hours
	!    =2: hours per worker
	!    =3: employment-population ratio
	!    =4: cel population-total population ratio
		if (cat.eq.1) then
			label='a. Average Hours '
			thead='Average Hours'
		elseif (cat.eq.2) then
			label='b. Hours per Worker '
			thead='Hours per Worker'
		elseif (cat.eq.3) then
			label='c. Employment-Population ratio '
			thead='Employment/Population, in %'
		elseif (cat.eq.4) then
			label='d. Cell Population-Total Population ratio '
			thead='Cell Pop/Total Population, in %'
		endif
		if (dog.eq.1) then
			do k=1,6
				do i=1,5
				tab2(k,i) = rmar(i,cat+1,1,k) ! matrix for total population 
				end do
			end do
		elseif (dog.eq.2) then
			do k=1,6
				do i=1,5
				tab2(k,i) = rmar(i,cat+1,2,k) ! matrix for males 
				end do
			end do
		elseif (dog.eq.3) then
			do k=1,6
				do i=1,5
				tab2(k,i) = rmar(i,cat+1,3,k) ! matrix for females 
				end do
			end do
		endif
		call table8(label,thead,note,tab2)		
	end do ! corresponds to cat=1,4
end do ! corresponds to dog=1,3
end do ! corresponds to rat=1,2



! Table for Appendix B of QRw98
	title='TABLE 9: TABLE OF APPENDIX B '
	write(4,'(A42)') title
	title='This is table of Appendix B in QRw98 '
	write(4,'(A42)') title
	write(4,*)
	write(4,*)

	do i=1,6
		do j=1,3
			tab1(i,j)=0.d0
		end do
	end do
	! tab1(i,j)=real-valued matrix
	!	i=1: 1950 and so on
	!	i=6: 2000
	!	j=1: average hours given in Table 1 of QRw98
	!	j=2: average hours calculated with the distribution of population
	!			across categories fixed at 1950
	!	j=3: average hours calculated with the distribution of hours
	!			across categories fixed at 1950

	! 1st column of tab1 is average hours given in Table 1 of QRw98
	do i=1,6
		tab1(i,1)=r(1,2,i) 
		! column 1 = average hours for total population in year i
	end do
	do i=65,192
        if (((m(i).lt.3).or.(m(i).gt.4)).and.(a(i).gt.0) &
                                        .and.(m(i).gt.0)) then
		  ! 2nd column of tab1
          tab1(1,2)  = tab1(1,2) + r50(i,2)*dfloat(i50(i,1))/tpop50
          tab1(2,2)  = tab1(2,2) + r60(i,2)*dfloat(i50(i,1))/tpop50
          tab1(3,2)  = tab1(3,2) + r70(i,2)*dfloat(i50(i,1))/tpop50
          tab1(4,2)  = tab1(4,2) + r80(i,2)*dfloat(i50(i,1))/tpop50
          tab1(5,2)  = tab1(5,2) + r90(i,2)*dfloat(i50(i,1))/tpop50
          tab1(6,2)  = tab1(6,2) + r00(i,2)*dfloat(i50(i,1))/tpop50
		  ! 3rd column of tab1
          tab1(1,3)  = tab1(1,3) + r50(i,2)*dfloat(i50(i,1))/tpop50
          tab1(2,3)  = tab1(2,3) + r50(i,2)*dfloat(i60(i,1))/tpop60
          tab1(3,3)  = tab1(3,3) + r50(i,2)*dfloat(i70(i,1))/tpop70
          tab1(4,3)  = tab1(4,3) + r50(i,2)*dfloat(i80(i,1))/tpop80
          tab1(5,3)  = tab1(5,3) + r50(i,2)*dfloat(i90(i,1))/tpop90
          tab1(6,3)  = tab1(6,3) + r50(i,2)*dfloat(i00(i,1))/tpop00
        endif
	end do
	call table9(tab1)

! Table for Appendix B of update of QRw98
	title='APPENDIX B '
	write(32,'(A42)') title
	title='A DECOMPOSITION OF AVERAGE WEEKLY HOURS '
	write(32,'(A42)') title
	title='WORKED PER PERSON '
	write(32,'(A42)') title
	write(32,*)

	call table9p1(tab1)


! Make graph for comparison of average hours calculated using 
!	method 2 (based on hrswork2) and method 4 (based on uhrswork).
	title='Chart B1: Two Methods of Calculating Average Hours: Total Population'
	count=1
	do i=1,5
		tab3(i) = 0.d0
	end do
	do i=1,3
		tab4(i) = 0.d0
	end do
	! Method 2
	do i=1,5
		tab3(i) = r(1,2,i) ! average hours for total population
	end do
	! Method 4
		tab4(1) = h80(1,4)
		tab4(2) = h90(1,4)
		tab4(3) = h00(1,4)
	call comparison(title,count,tab3,tab4)

	title='Chart B2: Two Methods of Calculating Average Hours: Males'
	count=2
	do i=1,5
		tab3(i) = 0.d0
	end do
	do i=1,3
		tab4(i) = 0.d0
	end do
	! Method 2
	do i=1,5
		tab3(i) = r(65,2,i) ! average hours for males
	end do
	! Method 4
		tab4(1) = h80(65,4)
		tab4(2) = h90(65,4)
		tab4(3) = h00(65,4)
	call comparison(title,count,tab3,tab4)

	title='Chart B3: Two Methods of Calculating Average Hours: Females'
	count=3
	do i=1,5
		tab3(i) = 0.d0
	end do
	do i=1,3
		tab4(i) = 0.d0
	end do
	! Method 2
	do i=1,5
		tab3(i) = r(129,2,i) ! average hours for females
	end do
	! Method 4
		tab4(1) = h80(129,4)
		tab4(2) = h90(129,4)
		tab4(3) = h00(129,4)
	call comparison(title,count,tab3,tab4)

! This part of the program produces figure 1 from paper by McGrattan,
!	Jones, Mannuelli at
! http://www.minneapolisfed.org/research/sr/sr317.html
! Construct matrix tabfig with dimensions (6,4,4):
!	6 rows correspond to census years: 1950, and so on
!	4 columns are
!		1 = Single Males
!		2 = Married Males, S.P.
!		3 = Single Females
!		4 = Married Females, S.P.
!	4 layers are
!		1 = average hours
!		2 = hours per worker
!		3 = employment-to-population ratio
!		4 = percentage of population in the group
! Age group is 25-64 years old. 

count = 0
do cat=1,4 ! cat is third dimension of tabfig
	count = count + 1
		do i=1,6
			do Col=1,4
				do j=1,4
				tabfig(i,Col,j) = 0.d0
				end do
			end do
		end do
	do k = 1,2  ! correspond to sex (Male, Female)
		do j = 1,2  ! correspond to marst (Single, Married,S.P.)
		pop50 = 0.d0
		pop60 = 0.d0
		pop70 = 0.d0
		pop80 = 0.d0
		pop90 = 0.d0
		pop00 = 0.d0
			Col = 2*(k-1) + j ! Column # for matrix tabfig
			do i = 1,192
			if ((s(i).eq.k).and.(m(i).eq.j).and.(a(i).ge.2) &
					.and.(a(i).le.5)) then
				! calculate total hours for each group
				tabfig(1,Col,cat) = tabfig(1,Col,cat) + r50(i,cat+1)*dfloat(i50(i,1))
				tabfig(2,Col,cat) = tabfig(2,Col,cat) + r60(i,cat+1)*dfloat(i60(i,1))
				tabfig(3,Col,cat) = tabfig(3,Col,cat) + r70(i,cat+1)*dfloat(i70(i,1))
				tabfig(4,Col,cat) = tabfig(4,Col,cat) + r80(i,cat+1)*dfloat(i80(i,1))
				tabfig(5,Col,cat) = tabfig(5,Col,cat) + r90(i,cat+1)*dfloat(i90(i,1))
				tabfig(6,Col,cat) = tabfig(6,Col,cat) + r00(i,cat+1)*dfloat(i00(i,1))
				! calculate population in the group
				pop50 = pop50 + dfloat(i50(i,1))
				pop60 = pop60 + dfloat(i60(i,1))
				pop70 = pop70 + dfloat(i70(i,1))
				pop80 = pop80 + dfloat(i80(i,1))
				pop90 = pop90 + dfloat(i90(i,1))
				pop00 = pop00 + dfloat(i00(i,1))

			endif				
			end do
			! calculate average hours
			tabfig(1,Col,cat) = tabfig(1,Col,cat)/pop50
			tabfig(2,Col,cat) = tabfig(2,Col,cat)/pop60
			tabfig(3,Col,cat) = tabfig(3,Col,cat)/pop70
			tabfig(4,Col,cat) = tabfig(4,Col,cat)/pop80
			tabfig(5,Col,cat) = tabfig(5,Col,cat)/pop90
			tabfig(6,Col,cat) = tabfig(6,Col,cat)/pop00
		end do
	end do
	if (cat.eq.1) then
	call figure1(tabfig(:,:,1)) ! graph for average hours
	endif
	if (cat.eq.1) then
		title = 'U.S. Average Hours since 1950, age group 25-64'
		thead = 'average hours'
		call figure1a(tabfig(:,:,1),count,title,thead)
	elseif (cat.eq.2) then
		title = 'U.S. Hours per Worker since 1950, age group 25-64'
		thead = 'hours per worker'
		call figure1a(tabfig(:,:,2),count,title,thead)
	elseif (cat.eq.3) then
		title = 'U.S. Employment-to-Population Ratio since 1950, age group 25-64'
		thead = 'employment-to-population ratio'
		call figure1a(tabfig(:,:,3),count,title,thead)
	elseif (cat.eq.4) then
		title = 'U.S. Percentage of Population in the group since 1950, age group 25-64'
		thead = '% of population in the group'
		call figure1a(tabfig(:,:,4),count,title,thead)
	endif

end do ! corresponds to cat=1,4

end program ipums_2000_totpop

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!	SUBROUTINES

subroutine table2(label,tab50,tab60,tab70,tab80,tab90,tab00)
! subroutine is used to print tables with aggregate statistics into ipums.out
! input variables are:
!	label of the table
!	ryr(3) correspond to columns 2 to 4 of matrix ryr
	
	real(kind = 8), intent(in) 	:: tab50(3),tab60(3),tab70(3),tab80(3),tab90(3),tab00(3) 
	character*56				:: line1, head1, head2
	character*42, intent(in)	:: label
	line1='-------------------------------------------'
	head1='       Average    Hours/     Employment/'
	head2='Year    Hours     Worker     Population, %'

	write(4,'(A42)') label
	write(4,'(A56)') line1
	write(4,'(A56)') head1
	write(4,'(A56)') head2
	write(4,'(A56)') line1
	write(4,'("1950 |", F7.4, 2F12.4)') (tab50(i), i=1,3)
	write(4,'("1960 |", F7.4, 2F12.4)') (tab60(i), i=1,3)
	write(4,'("1970 |", F7.4, 2F12.4)') (tab70(i), i=1,3)
	write(4,'("1980 |", F7.4, 2F12.4)') (tab80(i), i=1,3)
	write(4,'("1990 |", F7.4, 2F12.4)') (tab90(i), i=1,3)
	write(4,'("2000 |", F7.4, 2F12.4)') (tab00(i), i=1,3)
	write(4,'(A56)') line1
	write(4,*) 
	write(4,*) 
	return

end subroutine table2

subroutine table2p1(tab50,tab60,tab70,tab80,tab90,tab00)
! subroutine is used to print tables with aggregate statistics into 
! paper_tables.out, table that goes into published paper
! input variables are:
!	ryr(3) correspond to columns 2 to 4 of matrix ryr
	
	real(kind = 8), intent(in) 	:: tab50(3),tab60(3),tab70(3),tab80(3),tab90(3),tab00(3) 
	real(kind = 8)				:: change(3)
	character*56				:: line1, head1, head2
	line1='-----------------------------------------------'
	head1='           Average    Hours/     Employment/'
	head2='Year        Hours     Worker     Population, %'
	do i=1,3
		change(i) = ((tab00(i)-tab50(i))/tab50(i))*100
	end do
	
	write(32,'(A56)') line1
	write(32,'(A56)') head1
	write(32,'(A56)') head2
	write(32,'(A56)') line1
	write(32,'("1950     |", F7.2, 2F12.2)') (tab50(i), i=1,3)
	write(32,'("1960     |", F7.2, 2F12.2)') (tab60(i), i=1,3)
	write(32,'("1970     |", F7.2, 2F12.2)') (tab70(i), i=1,3)
	write(32,'("1980     |", F7.2, 2F12.2)') (tab80(i), i=1,3)
	write(32,'("1990     |", F7.2, 2F12.2)') (tab90(i), i=1,3)
	write(32,'("2000     |", F7.2, 2F12.2)') (tab00(i), i=1,3)
	write(32,'(A56)') line1
	write(32,'("% Change |", F7.2, 2F12.2)') (change(i), i=1,3)
	write(32,'("1950-2000|")') 
	write(32,'(A56)') line1
	write(32,*) 
	write(32,*) 
	return

end subroutine table2p1


subroutine table1(tab,label,thead)
! subroutine is used to print tables with aggregate statistics into ipums.out
	
	real(kind = 8), intent(in)	:: tab(6,7)
	character*70				:: line1, head1, head2
	character*42, intent(in)	:: label, thead
	line1='---------------------------------------------------------------'
	head1=' Age     ------------------------------------------------------'
	head2='in 50     15-24   25-34   35-44   45-54   55-64   65-74   75-99'
	write(4,'(A42)') label
	write(4,'(A70)') line1
	write(4,'(A56)') thead	
	write(4,'(A70)') head1
	write(4,'(A70)') head2
	write(4,'(A70)') line1
 	
	write(4,'(''-35:-26'', 6F8.4)')     tab(1,1)
    write(4,'(''-25:-16'', 6F8.4)')     tab(2,1),tab(1,2)
    write(4,'(''-15: -6'', 6F8.4)')     tab(3,1),tab(2,2),tab(1,3)
    write(4,'('' -5:  4'', 6F8.4)')     tab(4,1),tab(3,2),tab(2,3),tab(1,4)
    write(4,'(''  5: 14'', 6F8.4)')     tab(5,1),tab(4,2),tab(3,3), &
                                          tab(2,4),tab(1,5)
    write(4,'('' 15: 24'', 6F8.4)')  tab(6,1),tab(5,2),tab(4,3), &
                                          tab(3,4),tab(2,5),tab(1,6)
    write(4,'('' 25: 34'',8X,6F8.4)')   tab(6,2),tab(5,3),tab(4,4), &
                                          tab(3,5),tab(2,6),tab(1,7)
    write(4,'('' 35: 44'',16X,5F8.4)')  tab(6,3),tab(5,4),tab(4,5),	&
                                          tab(3,6),tab(2,7)
    write(4,'('' 45: 54'',24X,5F8.4)')  tab(6,4),tab(5,5),tab(4,6),tab(3,7)
    write(4,'('' 55: 64'',32X,5F8.4)')  tab(6,5),tab(5,6),tab(4,7)
    write(4,'('' 65: 74'',40X,5F8.4)')  tab(6,6),tab(5,7)
	write(4,'('' 75: 99'',48X,5F8.4)')  tab(6,7)	
	
	write(4,'(A70)') line1
	write(4,*) 
	write(4,*) 
	return

end subroutine table1

subroutine table6(label,thead,tab)
! subroutine is used to print tables with aggregate statistics into ipums.out
! this is to reproduce table 2 in QRw98: Distribution by sex
	
	real(kind = 8), intent(in) 	:: tab(6,3) 
	character*56				:: line1,line2, head1, head2
	character*42, intent(in)	:: label,thead
	line1='-------------------------------------------'
	line2='      -------------------------------------'
	head1='        Total              Sex             '
	head2='Year    Population    Males     Females'

	write(4,'(A42)') label
	write(4,'(A56)') line1
	write(4,'(A60)') thead	
	write(4,'(A56)') line2
	write(4,'(A56)') head1
	write(4,'(A56)') head2
	write(4,'(A56)') line1
	write(4,'("1950 |", F9.4, 2F12.4)') (tab(1,i), i=1,3)
	write(4,'("1960 |", F9.4, 2F12.4)') (tab(2,i), i=1,3)
	write(4,'("1970 |", F9.4, 2F12.4)') (tab(3,i), i=1,3)
	write(4,'("1980 |", F9.4, 2F12.4)') (tab(4,i), i=1,3)
	write(4,'("1990 |", F9.4, 2F12.4)') (tab(5,i), i=1,3)
	write(4,'("2000 |", F9.4, 2F12.4)') (tab(6,i), i=1,3)
	write(4,'(A56)') line1
	write(4,*) 
	write(4,*) 
	return

end subroutine table6

subroutine table6p1(thead,tab)
! subroutine is used to print tables with aggregate statistics into ipums.out
! this is to reproduce table 2 in QRw98: Distribution by sex
	
	real(kind = 8), intent(in) 	:: tab(6,3) 
	real(kind = 8)				:: change(3)
	character*56				:: line1,line2, head1, head2
	character*42, intent(in)	:: thead
	line1='-----------------------------------------------'
	line2='          -------------------------------------'
	head1='            Total              Sex             '
	head2='Year        Population    Males     Females'
	do i=1,3
		change(i) = ((tab(6,i)-tab(1,i))/tab(1,i))*100
	end do

	write(32,'(A56)') line1
	write(32,'(A60)') thead	
	write(32,'(A56)') line2
	write(32,'(A56)') head1
	write(32,'(A56)') head2
	write(32,'(A56)') line1
	write(32,'("1950     |", F9.2, 2F12.2)') (tab(1,i), i=1,3)
	write(32,'("1960     |", F9.2, 2F12.2)') (tab(2,i), i=1,3)
	write(32,'("1970     |", F9.2, 2F12.2)') (tab(3,i), i=1,3)
	write(32,'("1980     |", F9.2, 2F12.2)') (tab(4,i), i=1,3)
	write(32,'("1990     |", F9.2, 2F12.2)') (tab(5,i), i=1,3)
	write(32,'("2000     |", F9.2, 2F12.2)') (tab(6,i), i=1,3)
	write(32,'(A56)') line1
	write(32,'("% Change |", F9.2, 2F12.2)') (change(i), i=1,3)
	write(32,'("1950-2000|")') 
	write(32,'(A56)') line1
	write(32,*) 
	write(32,*) 
	return

end subroutine table6p1


subroutine table7(label,thead,tab)
! subroutine is used to print tables with aggregate statistics into ipums.out
! this is to reproduce table 3 in QRw98: Distribution by age
	
	real(kind = 8), intent(in) 	:: tab(6,7) 
	character*70				:: line1, head1, head2
	character*42, intent(in)	:: label,thead
	line1='-------------------------------------------------------'
	head1='        -----------------------------------------------'
	head2=' Year   15-24  25-34  35-44  45-54  55-64  65-74  75-99'

	write(4,'(A42)') label
	write(4,'(A70)') line1
	write(4,'(A60)') thead	
	write(4,'(A70)') head1
	write(4,'(A70)') head2
	write(4,'(A70)') line1
	write(4,'("1950 |", 7F7.2)') (tab(1,i), i=1,7)
	write(4,'("1960 |", 7F7.2)') (tab(2,i), i=1,7)
	write(4,'("1970 |", 7F7.2)') (tab(3,i), i=1,7)
	write(4,'("1980 |", 7F7.2)') (tab(4,i), i=1,7)
	write(4,'("1990 |", 7F7.2)') (tab(5,i), i=1,7)
	write(4,'("2000 |", 7F7.2)') (tab(6,i), i=1,7)
	write(4,'(A70)') line1
	write(4,*) 
	write(4,*) 
	return

end subroutine table7

subroutine table7p1(thead,tab)
! subroutine is used to print tables with aggregate statistics into ipums.out
! this is to reproduce table 3 in QRw98: Distribution by age
	
	real(kind = 8), intent(in) 	:: tab(6,7) 
	real(kind = 8)				:: change(7)
	character*80				:: line1, head1, head2
	character*80, intent(in)	:: thead
	line1='-------------------------------------------------------------------'
	head1='          ---------------------------------------------------------'
	head2=' Year        15-24   25-34   35-44   45-54   55-64   65-74   75-99'
	do i=1,7
		change(i) = ((tab(6,i)-tab(1,i))/tab(1,i))*100
	end do

	write(32,'(A70)') line1
	write(32,'(A70)') thead	
	write(32,'(A70)') head1
	write(32,'(A70)') head2
	write(32,'(A70)') line1
	write(32,'("1950     |", 7F8.2)') (tab(1,i), i=1,7)
	write(32,'("1960     |", 7F8.2)') (tab(2,i), i=1,7)
	write(32,'("1970     |", 7F8.2)') (tab(3,i), i=1,7)
	write(32,'("1980     |", 7F8.2)') (tab(4,i), i=1,7)
	write(32,'("1990     |", 7F8.2)') (tab(5,i), i=1,7)
	write(32,'("2000     |", 7F8.2)') (tab(6,i), i=1,7)
	write(32,'(A70)') line1
	write(32,'("% Change |", 7F8.2)') (change(i), i=1,7)
	write(32,'("1950-2000|")') 
	write(32,'(A70)') line1
	write(32,*) 
	write(32,*) 
	return

end subroutine table7p1


subroutine table8(label,thead,note,tab)
! subroutine is used to print tables with aggregate statistics into ipums.out
! this is to reproduce table 4 in QRw98: Distribution by marital status
	
	real(kind = 8), intent(in) 	:: tab(6,5) 
	character*70				:: line1, head1, head2
	character*42, intent(in)	:: label,thead,note
	line1='---------------------------------------------------------'
	head1='        -------------------------------------------------'
	head2=' Year   Mar S.P.   Mar S.A.   Single   Widowed   Divorced'

	write(4,'(A42)') label
	write(4,'(A70)') line1
	write(4,'(A60)') thead	
	write(4,'(A70)') head1
	write(4,'(A70)') head2
	write(4,'(A70)') line1
	write(4,'("1950 |", 5F10.2)') (tab(1,i), i=1,5)
	write(4,'("1960 |", 5F10.2)') (tab(2,i), i=1,5)
	write(4,'("1970 |", 5F10.2)') (tab(3,i), i=1,5)
	write(4,'("1980 |", 5F10.2)') (tab(4,i), i=1,5)
	write(4,'("1990 |", 5F10.2)') (tab(5,i), i=1,5)
	write(4,'("2000 |", 5F10.2)') (tab(6,i), i=1,5)
	write(4,'(A70)') line1
	write(4,'(A42)') note
	write(4,*) 
	write(4,*) 
	return

end subroutine table8

subroutine table9(tab1)
! subroutine is used to print table 9 into ipums.out
! this table corresponds to table of Appendix B of QRw98
	
	real(kind = 8), intent(in) 	:: tab1(6,3) 
	character*56				:: line1, head1, head2
	line1='-------------------------------------------'
	head1='       Actual Hours     1950         1950'
	head2='Year    per Person     Weight       Hours'

	write(4,'(A56)') line1
	write(4,'(A56)') head1
	write(4,'(A56)') head2
	write(4,'(A56)') line1
	write(4,'("1950 |", F12.4, 2F12.4)') (tab1(1,i), i=1,3)
	write(4,'("1960 |", F12.4, 2F12.4)') (tab1(2,i), i=1,3)
	write(4,'("1970 |", F12.4, 2F12.4)') (tab1(3,i), i=1,3)
	write(4,'("1980 |", F12.4, 2F12.4)') (tab1(4,i), i=1,3)
	write(4,'("1990 |", F12.4, 2F12.4)') (tab1(5,i), i=1,3)
	write(4,'("2000 |", F12.4, 2F12.4)') (tab1(6,i), i=1,3)
	write(4,'(A56)') line1
	write(4,*) 
	write(4,*) 
	return

end subroutine table9

subroutine table9p1(tab1)
! subroutine is used to print table 9 into ipums.out
! this table corresponds to table of Appendix B of QRw98
	
	real(kind = 8), intent(in) 	:: tab1(6,3) 
	real(kind = 8)				:: change(3)
	character*70				:: line1, head1, head2, &
								   line2, head3
	line1='-----------------------------------------------------------'
	line2='                         ----------------------------------'
	head3='                         Hours per Person Recalculated With'
	head1='           Actual Hours           1950           1950'
	head2='Year        per Person           Weights         Hours'
	do i=1,3
		change(i) = ((tab1(6,i)-tab1(1,i))/tab1(1,i))*100
	end do

	write(32,'(A70)') line1
	write(32,'(A70)') head3
	write(32,'(A70)') line2
	write(32,'(A70)') head1
	write(32,'(A70)') head2
	write(32,'(A70)') line1
	write(32,'("1950     |", F12.2, 2F16.2)') (tab1(1,i), i=1,3)
	write(32,'("1960     |", F12.2, 2F16.2)') (tab1(2,i), i=1,3)
	write(32,'("1970     |", F12.2, 2F16.2)') (tab1(3,i), i=1,3)
	write(32,'("1980     |", F12.2, 2F16.2)') (tab1(4,i), i=1,3)
	write(32,'("1990     |", F12.2, 2F16.2)') (tab1(5,i), i=1,3)
	write(32,'("2000     |", F12.2, 2F16.2)') (tab1(6,i), i=1,3)
	write(32,'(A70)') line1
	write(32,'("% Change |", F12.2, 2F16.2)') (change(i), i=1,3)
	write(32,'("1950-2000|")') 
	write(32,'(A70)') line1
	write(32,*) 
	write(32,*) 
	return

end subroutine table9p1


subroutine extrapolation(tab,label,thead,name,k,n,cat,count)
! Subroutine does (i) extrapolation of values in Tables 2-5,
!	(ii) writing info into ipums-extrap.m to print out graphs with 
!	extrapolated series
! Variables:	
!	tab = values from main program to put into table
!	label = label of the table
!	thead = head of the table
!	name = name of table used to call matrices in ipums-extract.m
!	k = sex
!	n = marital status
	real(kind = 8), intent(in)		:: tab(6,7)
	character*42, intent(in)		:: label, thead
	character*3, intent(in)			:: name
	integer, intent(in)				:: k,n,cat,count
	character*70					:: line1, head1, head2,title,head3
	real(kind = 8), dimension(12,7)	:: M
	real(kind = 8), parameter		:: w=0.3333
	integer							:: i,j
	
	! dimension of M: 12=# of groups of 'Year born'
	!				   7=# of age groups
	! initialize matrix M with zeros
	do i=1,12
		do j=1,7
			M(i,j) = 0
		end do
	end do

	! insert values for some rows of M from tab(6,7)
	! correspondance between rows of M and 'Year born' is
	!	row  Year born
	!	1		1976-85
	!	2		1966-75
	!	and so on
	!	11		1876-85
	!	12		1866-75
	M(1,1)   = tab(1,1)
    M(2,1:2) = (/tab(2,1),tab(1,2)/)
    M(3,1:3) = (/tab(3,1),tab(2,2),tab(1,3)/)
    M(4,1:4) = (/tab(4,1),tab(3,2),tab(2,3),tab(1,4)/)
    M(5,1:5) = (/tab(5,1),tab(4,2),tab(3,3),tab(2,4),tab(1,5)/)
    M(6,1:6) = (/tab(6,1),tab(5,2),tab(4,3),tab(3,4),tab(2,5),tab(1,6)/)
    M(7,2:7) = (/tab(6,2),tab(5,3),tab(4,4),tab(3,5),tab(2,6),tab(1,7)/)
    M(8,3:7) = (/tab(6,3),tab(5,4),tab(4,5),tab(3,6),tab(2,7)/)
    M(9,4:7) = (/tab(6,4),tab(5,5),tab(4,6),tab(3,7)/)
    M(10,5:7) = (/tab(6,5),tab(5,6),tab(4,7)/)
    M(11,6:7) = (/tab(6,6),tab(5,7)/)
	M(12,7)   = tab(6,7)
	
	! extrapolation of values, explaned in QRw98, appendix C
	do i=6,1,-1
		do j=i+1,7
			M(i,j) = M(i+1,j) + w*(M(i+1,j) - M(i+2,j)) &
							  + w*(M(i+2,j) - M(i+3,j)) &
							  + w*(M(i+3,j) - M(i+4,j)) 
		end do
	end do
	do i=7,12
		do j=1,i-6
			M(i,j) = M(i-1,j) + w*(M(i-1,j) - M(i-2,j)) &
							  + w*(M(i-2,j) - M(i-3,j)) &
							  + w*(M(i-3,j) - M(i-4,j)) 
		end do	
	end do	

	! get rid of negative values in matrix M
	do i=1,12
		do j=1,7
			if (M(i,j).le.0) then
				M(i,j) = 0
			endif
		end do
	end do
	! get rid of bogus values in matrix M, 
	!	bogus values are present to avoid divicion by zero
	do i=1,12
		do j=1,7
			if (M(i,j) .eq. 299.50) then
				M(i,j) = 0
			endif
		end do
	end do

! making table
	line1='----------------------------------------------------------'
	head1=' Year      -----------------------------------------------'
	head2=' born      15-24  25-34  35-44  45-54  55-64  65-74  75-99'
	write(3,'(A42)') label
	write(3,'(A60)') line1
	write(3,'(A60)') thead	
	write(3,'(A60)') head1
	write(3,'(A60)') head2
	write(3,'(A60)') line1
	write(3,'(''1866-75  '', 7F7.2)')   M(12,:)
	write(3,'(''1876-85  '', 7F7.2)')   M(11,:)
	write(3,'(''1886-95  '', 7F7.2)')   M(10,:)
	write(3,'(''1896-1905'', 7F7.2)')   M(9,:)
	write(3,'(''1906-15  '', 7F7.2)')   M(8,:)
	write(3,'(''1916-25  '', 7F7.2)')   M(7,:)
	write(3,'(''1926-35  '', 7F7.2)')   M(6,:)
	write(3,'(''1936-45  '', 7F7.2)')   M(5,:)
	write(3,'(''1946-55  '', 7F7.2)')   M(4,:)
	write(3,'(''1956-65  '', 7F7.2)')   M(3,:)
	write(3,'(''1966-75  '', 7F7.2)')   M(2,:)
	write(3,'(''1976-85  '', 7F7.2)')   M(1,:)
	write(3,'(A60)') line1
	write(3,*) 
	write(3,*) 

! Writing into paper_tables.out
if (cat.eq.1) then
	if ((k.eq.9).and.(n.eq.9)) then
		title='TABLE 8: TOTAL'
		write(32,'(A60)') title
		write(32,*)
	elseif ((k.eq.1).and.(n.eq.0)) then
		title='TABLE 9: MALES'
		write(32,'(A60)') title
		write(32,*)
	elseif ((k.eq.2).and.(n.eq.0)) then
		title='TABLE 10: FEMALES'
		write(32,'(A60)') title
		write(32,*)
	endif

	if (((k.eq.9).and.(n.eq.9)).or. &
	    ((k.eq.1).and.(n.eq.0)).or. &
		((k.eq.2).and.(n.eq.0))) then
	head3='Average Weekly Hours Worked per Person by Age (in Years)'
	line1='-----------------------------------------------------------------'
	head1=' Year      ------------------------------------------------------'
	head2=' born       15-24   25-34   35-44   45-54   55-64   65-74   75-99'
	write(32,'(A65)') line1
	write(32,'(A80)') head3	
	write(32,'(A65)') head1
	write(32,'(A65)') head2
	write(32,'(A65)') line1
	write(32,'(''1866-75  '', 7F8.2)')   M(12,:)
	write(32,'(''1876-85  '', 7F8.2)')   M(11,:)
	write(32,'(''1886-95  '', 7F8.2)')   M(10,:)
	write(32,'(''1896-1905'', 7F8.2)')   M(9,:)
	write(32,'(''1906-15  '', 7F8.2)')   M(8,:)
	write(32,'(''1916-25  '', 7F8.2)')   M(7,:)
	write(32,'(''1926-35  '', 7F8.2)')   M(6,:)
	write(32,'(''1936-45  '', 7F8.2)')   M(5,:)
	write(32,'(''1946-55  '', 7F8.2)')   M(4,:)
	write(32,'(''1956-65  '', 7F8.2)')   M(3,:)
	write(32,'(''1966-75  '', 7F8.2)')   M(2,:)
	write(32,'(''1976-85  '', 7F8.2)')   M(1,:)
	write(32,'(A65)') line1
	write(32,*) 
	write(32,*) 

	endif

endif ! correspond to cat=1

! Writing into ipums-extrap.m
	write(2,'(A3,I1,I1,''=[ ... '')') name,k,n
    write(2,'( 7F7.2)')   M(12,:)
	write(2,'( 7F7.2)')   M(11,:)
	write(2,'( 7F7.2)')   M(10,:)
	write(2,'( 7F7.2)')   M(9,:)
	write(2,'( 7F7.2)')   M(8,:)
	write(2,'( 7F7.2)')   M(7,:)
	write(2,'( 7F7.2)')   M(6,:)
	write(2,'( 7F7.2)')   M(5,:)
	write(2,'( 7F7.2)')   M(4,:)
	write(2,'( 7F7.2)')   M(3,:)
	write(2,'( 7F7.2)')   M(2,:)
	write(2,'( 7F7.2)')   M(1,:)
    write(2,'(''];'')')

	write(2,'("figure(",I2,")")') count
    write(2,'(''title('',A1,A42,A1,'')'')') "'",label,"'"
    write(2,'(''xlabel('',A5,'')'')') "'age'"
    write(2,'(''ylabel('',A1,A42,A1,'')'')') "'",thead,"'"

	write(2,'("line(a,",A3,I1,I1,"(1,:),",A11,",2,",A7, &
				",[0,0,0],",A11,",",A3,")")') name,k,n, &
				"'LineWidth'","'Color'","'LineStyle'","':'"
	write(2,'("line(a,",A3,I1,I1,"(2,:),",A11,",2,",A7, &
				",[0,0,0])")') name,k,n, &
				"'LineWidth'","'Color'"
	write(2,'("line(a,",A3,I1,I1,"(3,:),",A11,",2,",A7, &
				",[0.5,0.5,0.5])")') name,k,n, &
				"'LineWidth'","'Color'"
	write(2,'("line(a,",A3,I1,I1,"(4,:),",A11,",2,",A7, &
				",[0,0,0],",A11,",",A4,")")') name,k,n, &
				"'LineWidth'","'Color'","'LineStyle'","'--'"
	write(2,'("line(a,",A3,I1,I1,"(5,:),",A11,",1,",A7, &
				",[0,0,0])")') name,k,n, &
				"'LineWidth'","'Color'"
	write(2,'("line(a,",A3,I1,I1,"(6,:),",A11,",1,",A7, &
				",[0,0,0],",A11,",",A4,")")') name,k,n, &
				"'LineWidth'","'Color'","'LineStyle'","'--'"
 	write(2,'("line(a,",A3,I1,I1,"(7,:),",A11,",1,",A7, &
				",[0,0,1],",A11,",",A4,")")') name,k,n, &
				"'LineWidth'","'Color'","'LineStyle'","'--'"
	write(2,'("line(a,",A3,I1,I1,"(8,:),",A11,",1,",A7, &
				",[0,0,1])")') name,k,n, &
				"'LineWidth'","'Color'"
	write(2,'("line(a,",A3,I1,I1,"(9,:),",A11,",2,",A7, &
				",[0,0,1],",A11,",",A4,")")') name,k,n, &
				"'LineWidth'","'Color'","'LineStyle'","'--'"
	write(2,'("line(a,",A3,I1,I1,"(10,:),",A11,",2,",A7, &
				",[1,0,0],",A11,",",A4,")")') name,k,n, &
				"'LineWidth'","'Color'","'LineStyle'","'--'"
	write(2,'("line(a,",A3,I1,I1,"(11,:),",A11,",2,",A7, &
				",[1,0,0])")') name,k,n, &
				"'LineWidth'","'Color'"
 	write(2,'("line(a,",A3,I1,I1,"(12,:),",A11,",2,",A7, &
				",[1,0,0],",A11,",",A3,")")') name,k,n, &
				"'LineWidth'","'Color'","'LineStyle'","':'"

	write(2,'("legend(",A70)') "'1866-75','1876-85', &
				'1886-95','1896-1905','1906-15','1916-25',..."
	write(2,'(A70,")")') "'1926-35','1936-45', &
				'1946-55','1956-65','1966-75','1976-85'"

    if (cat .eq. 1) then
 		write(2,'(''axis([15,85,0,50]);'')')
    elseif (cat .eq. 2) then
!		write(2,'(''axis([10,80,0,50]);'')')
    elseif (cat .eq. 3) then
		write(2,'(''axis([15,85,0,100]);'')')
    elseif (cat .eq. 4) then
!		write(2,'(''axis([10,80,0,30]);'')')
	endif
    write(2,'(''grid'')')
    write(2,'(''  print -dpsc fig_'',A3,I1,I1)') name,k,n
!    write(2,'(''pause(5)'')')
    write(2,*)
	return
end subroutine extrapolation


SUBROUTINE comparison(title,count,tab1,tab2)
! The subroutine writes into matlab file used to produce figure 1
! from QRw98, time series for average hours

real(kind=8),intent(in) :: tab1(5),tab2(3)
integer, intent(in)		:: count
character*70, intent(in):: title

write(5,'("yr1=1950:10:1990;")')
write(5,'("tab1=[ ... ")') 
do i=1,5
	write(5,'(2F20.4)') tab1(i)
end do
write(5,'(''];'')')
write(5,*)
write(5,'("yr2=1980:10:2000;")')
write(5,'("tab2=[ ... ")') 
do i=1,3
	write(5,'(4F20.4)') tab2(i)
end do
write(5,'(''];'')')
write(5,*)
	! plot the graph
write(5,'("figure(",I1,")")') count


	write(5,'("set(gca,",A7,",[1950 1960 1970 1980  &
			1990 2000])")') "'XTick'"
    write(5,'(''title('',A1,A70,A2,A14,'')'')') "'",title,"',",&
					"'FontSize',16"
    write(5,'(''xlabel('',A8,A14,'')'')') "'year',",&
					"'FontSize',14"
    write(5,'(''ylabel('',A17,A14,'')'')') "'average hours',",&
					"'FontSize',14"
	write(5,'("line(yr1,tab1,",A11,",2,",A7, &
				",[0,0,1]",A50,")")')  &
				"'LineWidth'","'Color'", &
				",'Marker','s','MarkerFaceColor',[0,0,1]"
	write(5,'("line(yr2,tab2,",A11,",2,",A7, &
				",[1,0,0]",A50,")")')  &
				"'LineWidth'","'Color'", &
				",'Marker','s','MarkerFaceColor',[1,0,0]"
	write(5,'("legend(",A70,",1)")') "'based on hrswork2',&
				'based on uhrswork'"
	write(5,'(''axis([1950,2000,10,35]);'')')
    write(5,'(''grid'')')

    write(5,'(''  print -dpsc comparison_'',I1)') count
	write(5,*)
END SUBROUTINE comparison

SUBROUTINE figure1(tab)
! The subroutine writes into matlab file used to produce figure 1
! from SR317

real(kind=8),dimension(6,4),intent(in) :: tab

write(33,'("yr=1950:10:2000;")')
write(33,'("tab=[ ... ")') 
do i=1,6
	write(33,'(4F20.4)') tab(i,:)
end do
write(33,'(''];'')')
write(33,*)
	! plot the graph
	write(33,'("set(gca,",A7,",[1950 1960 1970 1980  &
			1990 2000])")') "'XTick'"
    write(33,'(''title('',A50,A14,'')'')') "'U.S. Average Hours since 1950, age group 25-64',",&
					"'FontSize',16"
    write(33,'(''xlabel('',A6,'')'')') "'year'"
    write(33,'(''ylabel('',A15,'')'')') "'average hours'"
	write(33,'("line(yr,tab(:,1),",A11,",2,",A7, &
				",[0,0,1],",A11,",",A4,A40,")")')  &
				"'LineWidth'","'Color'","'LineStyle'","'--'",&
				",'Marker','s','MarkerFaceColor',[0,0,1]"
	write(33,'("line(yr,tab(:,2),",A11,",2,",A7, &
				",[0,0,1]",A40,")")')  &
				"'LineWidth'","'Color'", &
				",'Marker','s','MarkerFaceColor',[0,0,1]"
	write(33,'("line(yr,tab(:,3),",A11,",2,",A7, &
				",[1,0,0],",A11,",",A4,A40,")")')  &
				"'LineWidth'","'Color'","'LineStyle'","'--'",&
				",'Marker','s','MarkerFaceColor',[1,0,0]"
	write(33,'("line(yr,tab(:,4),",A11,",2,",A7, &
				",[1,0,0]",A40,")")')  &
				"'LineWidth'","'Color'", &
				",'Marker','s','MarkerFaceColor',[1,0,0]"
	write(33,'("text(",A10,",[1952 34],",A22,")")') "'Position'", &
				"'String','Single Men'"
	write(33,'("text(",A10,",[1952 43],",A22,")")') "'Position'", &
				"'String','Married Men'"
	write(33,'("text(",A10,",[1952 27],",A23,")")') "'Position'", &
				"'String','Single Women'"
	write(33,'("text(",A10,",[1952 7],",A25,")")') "'Position'", &
				"'String','Married Women'"
	write(33,'(''axis([1950,2000,0,50]);'')')
    write(33,'(''  print -dpsc fig1sr317'')') 
	write(33,*)

END SUBROUTINE figure1


SUBROUTINE figure1a(tab,count,title,thead)
! The matlab file will include 4 graphs
!		1 = average hours
!		2 = hours per worker
!		3 = employment-to-population ratio
!		4 = percentage of population in the group
! The subroutine writes into matlab file used to produce figure 1
! from SR317

real(kind=8),dimension(6,4),intent(in) :: tab
integer, intent(in)		:: count
character*70, intent(in):: title,thead

write(34,'("yr=1950:10:2000;")')
write(34,'("tab=[ ... ")') 
do i=1,6
	write(34,'(4F20.4)') tab(i,:)
end do
write(34,'(''];'')')
write(34,*)
	! plot the graph
	write(34,'("figure(",I1,")")') count
	write(34,'("set(gca,",A7,",[1950 1960 1970 1980  &
			1990 2000])")') "'XTick'"
    write(34,'(''title('',A1,A70,A16,'')'')') "'",title,&
					"','FontSize',14"
    write(34,'(''xlabel('',A6,'')'')') "'year'"
    write(34,'(''ylabel('',A1,A30,A1,'')'')') "'",thead,"'"
	write(34,'("line(yr,tab(:,1),",A11,",2,",A7, &
				",[0,0,1],",A11,",",A4,A40,")")')  &
				"'LineWidth'","'Color'","'LineStyle'","'--'",&
				",'Marker','s','MarkerFaceColor',[0,0,1]"
	write(34,'("line(yr,tab(:,2),",A11,",2,",A7, &
				",[0,0,1]",A40,")")')  &
				"'LineWidth'","'Color'", &
				",'Marker','s','MarkerFaceColor',[0,0,1]"
	write(34,'("line(yr,tab(:,3),",A11,",2,",A7, &
				",[1,0,0],",A11,",",A4,A40,")")')  &
				"'LineWidth'","'Color'","'LineStyle'","'--'",&
				",'Marker','s','MarkerFaceColor',[1,0,0]"
	write(34,'("line(yr,tab(:,4),",A11,",2,",A7, &
				",[1,0,0]",A40,")")')  &
				"'LineWidth'","'Color'", &
				",'Marker','s','MarkerFaceColor',[1,0,0]"
	write(34,'("legend(",A70,",-1)")') "'Single Men','Married Men',&
				'Single Women','Married Women'"
!	write(34,'(''axis([1950,2000,0,50]);'')')
    write(34,'(''grid'')')
    write(34,'(''  print -dpsc fig1sr317_'',I1)') count
	write(34,*)

END SUBROUTINE figure1a