Getting NaN for a write output in FORTRAN

[antony1103]

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!

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)?: "
read *, diff_const

write (*,'(a,$)') "Macroscopic cross section (1/cm)?: "
read *, macro_cross_sec

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



! 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

[D H]

First off, please use the [ code ] ... [ /code ] tags, sans the spaces, when you post code. I edited your post for you.

As far as NaNs are concerned:
NaNs can arise from some computations such as 0/0. They can also result from the compiler. I suspect you are using a compiler that initializes uninitialized floating point variables with NaNs. In particular, you did not initialize the array flux2 or the variable e. Since e is not set, the then part of the if statement never executes, and flux2 keeps its (uninitialized) initial values.

The best way to find these kinds of errors is to use your debugger. Learn to use it.