This is for temperature: 298

This is for reac- 1

program check
implicit none


! Setting up parameters and necessary conversion factors
real,parameter::h=4.135667516e-15;            		    ! planck's constant in eVs
real,parameter::c=3.0e10;                           	            ! velocity of light in cm/s
real,parameter::kb=8.6173324e-5;                            ! boltzmann's constant in eV/K

							
integer, parameter:: dp = selected_real_kind(15, 307)
real(dp),allocatable::zper(:)					! zero point energy array 
real(dp),allocatable::zpcr(:)				! zero point corrected energy
real(dp),allocatable::qvibr(:)
real(dp)::zpcrtot=0.0					! total zero point corrected energy 

integer::i,j
integer::r
character(len=20)::filename	
integer::ierror,alloc_err,nlines,nless,array_t
real(dp)::scf,vibmulti,freqtot
real(dp),allocatable::freq(:),vib(:)
real(dp)::T(29) =(/(array_t,array_t=273,973,25)/)



write(*,*)"How many reactants do you have?"
read(*,*)r


!allocate proper dimension to all zero point arrays
allocate(zper(r))
allocate(zpcr(r))
allocate(qvibr(r))

! Loop through reactant files

do j= 1,29 !temp loop

write(*,*)"This is for temperature: ",T(j)
    
  do i=1,r	! file loop
  write(filename,100)i
  100 format('reac-',I1)											! get filenames

	write(*,*)"This is for reac-",i
     
	freqtot=0														! initializing total frequency value
	vibmulti=1														!
    ! Count no of lines	
	nlines=0
	OPEN(unit=10,file=filename,status='old',action='read',&
	 iostat=ierror)
			do	! count loop
            read(10,*,iostat=ierror)
             if (ierror /= 0) exit
  			nlines=nlines+1
            end do !end count loop
	CLOSE(unit=10)
	
	! Count no of frequency elements
	nless=nlines-1
	
	!allocate a dimension to frquency array
	allocate(freq(nless))
	allocate(vib(nless))
	
	OPEN(unit=10,file=filename,status='old',action='read',&
	 iostat=ierror)

			read(10,*,iostat=ierror) scf
			
			do  	! freq loop
            read(10,*,iostat=ierror) freq(nless)
            if (ierror /= 0) exit
          	if (freq(nless) < 100) then
              	freq(nless) = 100
            end if            
            
			freqtot=freqtot+freq(nless)
            			           
			vib(nless)= 1/(1-exp(-h*c*freq(nless)/kb/T(j)))
            
            vibmulti=vib(nless)*vibmulti
			
	     	end do	!end freq loop
     
     deallocate (freq,stat=alloc_err) 
     deallocate (vib,stat=alloc_err)      		
     CLOSE(10)


zper(i)= 0.5*h*c*freqtot												! zero point energy
zpcr(i)= scf + zper(i)													! zero point correction
zpcrtot= zpcrtot+zpcr(i)												! total sum of zero point corrected energies of reactants												

write(*,*)"ZPER", zper(i)
write(*,*)"ZPCR", zpcr(i)
write(*,*)"ZPCRTOT", zpcrtot
write(*,*)"Vib.P.F for reac",i,vibmulti  
qvibr(i)=vibmulti
write(*,*)"This is test qvib: ", qvibr(i) 

	end do	! end file loop
    
write(*,*) "Try product: ",PRODUCT(qvibr)

deallocate (qvibr,stat=alloc_err)     
deallocate (zper,stat=alloc_err) 
deallocate (zpcr,stat=alloc_err)

end do  ! end temp loop

end program check