Getting NaN for a write output in FORTRAN

• Fortran

Main Question or Discussion Point

Hi all,

I'm writing a program in FORTRAN to calculate neutron flux in a reactor, but when I try to get a numerical output it simply says NaN. Anyone know what this means or how to fix it? Code is below, and this is happening on the flux2(N) output on the last write statement. Thanks!

Code:
      program slabflux

implicit none

! VARIABLES

integer, parameter :: N=30		! Number of slab increments

real :: thickness         		! slab thickness in meters
real :: diff_const        		! diffusion constant of the slab
real :: macro_cross_sec   		! macroscopic cross section of the slab
real :: source            		! neutron source rate
real :: distance          		! distance of neutron flux calculation in the slab
real, dimension(-1:N+1) :: flux1       ! flux calculation at a certain distance in the slab
real, dimension(0:N) :: flux2		! used for iterative flux calculation
real, dimension(0:N) :: S		! source array
real :: del				! increment length
real :: a				! flux coefficient
real :: b				! flux coefficient
real, dimension(-1:N) :: prev_flux	! previous flux summation for loop
real, dimension(0:N) :: prev_iter	! flux summation from previous iteration
integer :: i				! do loop counter
integer :: j				! do loop counter
integer :: k				! do loop counter
real :: e				! iteration error

! INPUT

write (*,'(a,$)') "Slab thickness (cm)?: " read *, thickness write (*,'(a,$)') "Diffusion constant (cm)?: "

write (*,'(a,$)') "Macroscopic cross section (1/cm)?: " read *, macro_cross_sec write (*,'(a,$)') "Neutron source rate (neutrons/cm^3*s)?: "

! CALCULATION

del=thickness/N					! increment length

S(N/2)=source				         	! source centered in slab

flux1(1:29)=1						! initial flux
prev_flux(-1:N)=0
prev_flux(1:29)=0

a=-(diff_const/(del**2))				! flux coefficients

b=(diff_const/(del**2))+macro_cross_sec

if (0.997<e .AND. e<1.003) then
flux2(0:N)=flux1(0:N)
do i=0,N,1						! first iteration
prev_flux(i)=flux2(i-1)+prev_flux(i-1)
do j=i,N,1
prev_iter(i)=prev_flux(i)+flux1(i+1)
end do
flux2(i)=(1/b)*(S(j)-(a*prev_flux(i))-(a*prev_iter(i)))
end do
e=abs(flux2(15)/flux1(15))
else
do k=0,N,1
distance=del*N
write (*,'(a,f6.2,3x,a,f6.2)') "Distance (cm): ",
+                                       distance,
+                                       "Flux (neutrons/cm^2*s): ",
+                                       flux2(N)
end do
end if

end program slabflux

Last edited by a moderator:

Related Programming and Computer Science News on Phys.org
D H
Staff Emeritus