subroutine WriteResults( unit, EquilOrProd, Nblocks, StartConf, EndConf, &
						 Nsp, BAVG_RHO, BAVG_U, BAVG_NMOL, BAVG_MOLE_FRAC, &
						 BAVG_P, AVG_RHO, AVG_U, AVG_NMOL, &
						 AVG_MOLE_FRAC, AVG_P, ERR_RHO, ERR_U, &
						 ERR_NMOL, ERR_MOLE_FRAC, ERR_P,  MaxSp, &
						 NATT, NSUC, NAMEsp, CDV, MOVETIME, tot_time, Ensemble )

implicit none

integer, parameter									:: Nphases = 2

integer, intent(in)									:: unit

integer, intent(in)									:: EquilOrProd

integer, intent(in)									:: Nblocks

integer, intent(in)									:: StartConf
integer, intent(in)									:: EndConf

integer, intent(in)									:: Nsp

real, dimension(Nblocks, Nphases), intent(in)		:: BAVG_RHO
real, dimension(Nblocks, Nphases), intent(in)		:: BAVG_U
real, dimension(Nblocks, Nphases), intent(in)		:: BAVG_NMOL
real, dimension(Nblocks, Nsp, Nphases), intent(in)	:: BAVG_MOLE_FRAC
real, dimension(Nblocks, Nphases), intent(in)		:: BAVG_P

real, dimension(Nphases), intent(in)				:: AVG_RHO
real, dimension(Nphases), intent(in)				:: AVG_U
real, dimension(Nphases), intent(in)				:: AVG_NMOL
real, dimension(Nsp, Nphases), intent(in)			:: AVG_MOLE_FRAC
real, dimension(Nphases), intent(in)				:: AVG_P

real, dimension(Nphases), intent(in)				:: ERR_RHO
real, dimension(Nphases), intent(in)				:: ERR_U
real, dimension(Nphases), intent(in)				:: ERR_NMOL
real, dimension(Nsp, Nphases), intent(in)			:: ERR_MOLE_FRAC
real, dimension(Nphases), intent(in)				:: ERR_P

integer, intent(in)									:: MaxSp	

integer, dimension(5,MaxSp,Nphases), intent(in)		:: NATT
integer, dimension(5,MaxSp,Nphases), intent(in)		:: NSUC

character*15, dimension(Nsp), intent(in)			:: NAMEsp

real, dimension(Nphases), intent(in)				:: CDV

real, dimension(5)									:: MOVETIME

real												:: tot_time					

character*5, intent(in)								:: Ensemble

! Local Stuff

integer									:: i, j, k

real									:: ratio, total

character*10, dimension(Nsp*Nphases)	:: LABEL



write(unit,"(80A)") ('=', j=1, 80)

if( EquilOrProd == 1 ) then

	write(unit,*) ' Equilibration period averages:  Steps ', StartConf, &
					' to ', EndConf

else

	write(unit,*) ' Production period averages:  Steps ', StartConf, &
					' to ', EndConf

	if( Ensemble == 'NVT' ) then

		write(unit,"(A, F6.2, A)") '  Constant Volume Change ', CDV(1), ' A^3 '

	else

		write(unit,"(A, 2F6.2, A)") '  Constant Volume Change ', CDV(1:2), ' A^3 '

	end if

end if

write(unit,"(80A)") ('=', j=1, 80)

if( Ensemble == 'NVT' .AND. EquilOrProd == 2 ) then

	write(unit, "(1x,A6, 35(2X,A10))") ' ', 'Rho 1      ', 'Rho 2      ', &
				' U 1      ', ' U 2      ',   ' P 1      ', ' P 2      '
	write(unit, "(1x,A6, 35(A12))") 'Block ', '  kg/m^3    ', '  kg/m^3    ', &
		'  kJ/kg     ', '  kJ/kg     ', '   bar      ', '   bar      '
	write(unit, "(1x,7A)") '----- ', ('----------- ', j=1,6)
	
	do j = 1, Nblocks
	
		write(unit, "(1x,I3, 2X, 35(1X,G11.4))") j, BAVG_RHO(j,1:2), &
										BAVG_U(j,1:2), BAVG_P(j,1:2)
			
	end do

	write(unit, "(1x,7A)") '----- ', ('----------- ', j=1,6)

	write(unit, "(1x,A5, 35(1X,G11.4))") 'Aver ', AVG_RHO(1:2), &
									AVG_U(1:2), AVG_P(1:2)

	write(unit, "(1x,A5, 35(1X,G11.4))") '  +- ', ERR_RHO(1:2), &
									ERR_U(1:2), ERR_P(1:2)

else

	write(unit, "(1x,A6, 35(2X,A10))") ' ', 'Rho 1      ', 'Rho 2      ', &
				' U 1      ', ' U 2      '
	write(unit, "(1x,A6, 35(A12))") 'Block ', '  kg/m^3    ', '  kg/m^3    ', &
		'  kJ/kg     ', '  kJ/kg     '
	
	write(unit, "(1x,7A)") '----- ', ('----------- ', j=1,4)
	
	do j = 1, Nblocks
	
		write(unit, "(1x,I3, 2X, 35(1X,G11.4))") j, BAVG_RHO(j,1:2), BAVG_U(j,1:2)
			
	end do

	if( EquilOrProd == 2 ) then

		write(unit, "(1x,7A)") '----- ', ('----------- ', j=1,6)

		write(unit, "(1x,A5, 35(1X,G11.4))") 'Aver ', AVG_RHO(1:2), AVG_U(1:2)

		write(unit, "(1x,A5, 35(1X,G11.4))") '  +- ', ERR_RHO(1:2), ERR_U(1:2)

	end if

end if

if( Nsp > 1 ) then

	k = 0

	do i = 1, Nphases

		do j = 1, Nsp

			k = k + 1

			LABEL(k) = '  x_'//char(48+j)//' '//char(48+i)

		end do

	end do

	write(unit, "(/,1x,A5, 35(2X,A10))") ' ',' Ntot 1   ', ' Ntot 2   ', LABEL(1:Nsp*Nphases)
	write(unit, "(1x,A5, 35(A12))") 'Block', ' Molecules ', ' Molecules ', &
									(' Mol %   ', j=1,Nsp*Nphases)

	write(unit, "(1x,35A)") '----- ', ('----------- ', j=1,(Nsp+1)*Nphases)

	do j = 1, Nblocks

		write(unit, "(1x,I3, 2X, 35(1X,G11.4))") j, BAVG_Nmol(j,1:2), &
											BAVG_MOLE_FRAC(j,1:Nsp,1:2)

	end do

	if( EquilOrProd == 2 ) then

		write(unit, "(1x,35A)") '----- ', ('----------- ', j=1,(Nsp+1)*Nphases)

		write(unit, "(1x,A5, 35(1X,G11.4))") 'Aver ', AVG_Nmol(1:2), AVG_MOLE_FRAC(1:Nsp,1:2)

		write(unit, "(1x,A5, 35(1X,G11.4))") '  +- ', ERR_Nmol(1:2), ERR_MOLE_FRAC(1:Nsp,1:2)

	end if

else

	write(unit, "(/,1x,A5, 35(2X,A10))") ' ',' Ntot 1   ', ' Ntot 2   '
	write(unit, "(1x,A5, 35(A12))") 'Block', '  Molecules ', '  Molecules '

	write(unit, "(1x,35A)") '----- ', ('----------- ', j=1,Nphases)

	do j = 1, Nblocks

		write(unit, "(1x,I3, 2X, 35(1X,G11.4))") j, BAVG_Nmol(j,1:2)

	end do

	if( EquilOrProd == 2 ) then

		write(unit, "(1x,35A)") '----- ', ('----------- ', j=1,Nphases)

		write(unit, "(1x,A5, 35(1X,G11.4))") 'Aver ', AVG_Nmol(1:2)

		write(unit, "(1x,A5, 35(1X,G11.4))") '  +- ', ERR_Nmol(1:2)

	end if

end if

write(unit,*)


write(unit, "(A/A/A)") '  Monte Carlo Move Statistics ', &
  '  Move Type   Phase/Donor   Species        Attempted   Successful     % Acceptance  ', &
  ' ----------- ------------- ---------      ----------- ------------   ------------- '
	 
if( Ensemble == 'NPT' ) then

	do i = 1, Nphases

		ratio = 0.0

		if( NATT(3,1,i) /= 0 ) ratio = real( NSUC(3,1,i) ) / real( NATT(3,1,i) ) * 100.0

		write(unit, "(2X, A, T19, I1, T44, I8, T56, I8, T72, F6.2)") &
			'Volume', i, NATT(3,1,i), NSUC(3,1,i), ratio

	end do

else if( Ensemble == 'NVT' ) then

	ratio = 0.0

	if( NATT(3,1,1) /= 0 ) ratio = real( NSUC(3,1,1) ) / real( NATT(3,1,1) ) * 100.0

	write(unit, "(2X, A, T44, I8, T56, I8, T72, F6.2)") &
		'Volume', NATT(3,1,1), NSUC(3,1,1), ratio

end if

do i = 1, Nphases

	do j = 1, Nsp

		ratio = 0.0

		if( NATT(1,j,i) /= 0 ) ratio = real( NSUC(1,j,i) ) / real( NATT(1,j,i) ) * 100.0

		write(unit, "(2X, A, T19, I1, T29, A, T44, I8, T56, I8, T72, F6.2)") &
			'Displace', i, NAMEsp(j), NATT(1,j,i), NSUC(1,j,i), ratio

	end do

end do

do i = 1, Nphases

	do j = 1, Nsp

		ratio = 0.0

		if( NATT(2,j,i) /= 0 ) ratio = real( NSUC(2,j,i) ) / real( NATT(2,j,i) ) * 100.0

		write(unit, "(2X, A, T19, I1, T29, A, T44, I8, T56, I8, T72, F6.2)") &
			'Rotate', i, NAMEsp(j), NATT(2,j,i), NSUC(2,j,i), ratio

	end do

end do

do i = 1, Nphases
	
	do j = 1, Nsp

		ratio = 0.0

		if( NATT(4,j,i) /= 0 ) ratio = real( NSUC(4,j,i) ) / real( NATT(4,j,i) ) * 100.0

		write(unit, "(2X, A, T19, I1, T29, A, T44, I8, T56, I8, T72, F6.2)") &
			'Transfer', i, NAMEsp(j), NATT(4,j,i), NSUC(4,j,i), ratio

	end do

end do

do i = 1, Nphases
	
	do j = 1, Nsp

		ratio = 0.0

		if( NATT(5,j,i) /= 0 ) ratio = real( NSUC(5,j,i) ) / real( NATT(5,j,i) ) * 100.0

		write(unit, "(2X, A, T19, I1, T29, A, T44, I8, T56, I8, T72, F6.2)") &
			'Regrowth', i, NAMEsp(j), NATT(5,j,i), NSUC(5,j,i), ratio

	end do

end do

total = sum( MOVETIME )

write(unit,*)

write(unit, "(A/A/A)") '  Time Distribution', &
					   '  Move Type        Move Time       % of total ', &
					   ' -----------    ---------------   ------------ '

if ( MOVETIME(3) < 60.0 ) then
	
	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Volume', MOVETIME(3), ' sec.', MOVETIME(3) / total	* 100.0

else if ( MOVETIME(3) < 3600.0 ) then

	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Volume', MOVETIME(3)/60.0, ' min.', MOVETIME(3) / total	* 100.0

else

	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Volume', MOVETIME(3)/3600.0, ' hrs.', MOVETIME(3) / total	* 100.0

endif

if ( MOVETIME(1) < 60.0 ) then
	
	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Displace', MOVETIME(1), ' sec.', MOVETIME(1) / total	* 100.0

else if ( MOVETIME(1) < 3600.0 ) then

	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Displace', MOVETIME(1)/60.0, ' min.', MOVETIME(1) / total	* 100.0

else

	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Displace', MOVETIME(1)/3600.0, ' hrs.', MOVETIME(1) / total	* 100.0

endif

if ( MOVETIME(2) < 60.0 ) then
	
	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Rotate', MOVETIME(2), ' sec.', MOVETIME(2) / total	* 100.0

else if ( MOVETIME(2) < 3600.0 ) then

	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Rotate', MOVETIME(2)/60.0, ' min.', MOVETIME(2) / total	* 100.0

else

	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Rotate', MOVETIME(2)/3600.0, ' hrs.', MOVETIME(2) / total	* 100.0

endif

if ( MOVETIME(4) < 60.0 ) then
	
	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Transfer', MOVETIME(4), ' sec.', MOVETIME(4) / total	* 100.0

else if ( MOVETIME(4) < 3600.0 ) then

	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Transfer', MOVETIME(4)/60.0, ' min.', MOVETIME(4) / total	* 100.0

else

	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Transfer', MOVETIME(4)/3600.0, ' hrs.', MOVETIME(4) / total	* 100.0

endif

if ( MOVETIME(5) < 60.0 ) then
	
	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Regrow', MOVETIME(5), ' sec.', MOVETIME(5) / total	* 100.0

else if ( MOVETIME(5) < 3600.0 ) then

	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Regrow', MOVETIME(5)/60.0, ' min.', MOVETIME(5) / total	* 100.0

else

	write(unit, "(1X, A, T17, F8.2, A, T37, F6.2)")  &
		' Regrow', MOVETIME(5)/3600.0, ' hrs.', MOVETIME(5) / total	* 100.0

endif

write(unit,*)

if ( tot_time < 60.0 ) then
	
	write(unit, "(1X, A, F8.2, A)")  &
		' Total Time = ', tot_time, ' sec.'

else if ( tot_time < 3600.0 ) then

	write(unit, "(1X, A, F8.2, A)")  &
		' Total Time = ', tot_time/60.0, ' min.'

else

	write(unit, "(1X, A, F8.2, A)")  &
		' Total Time = ', tot_time/3600.0, ' hrs.'

endif

return

end subroutine WriteResults