[Fortran 90] Output is NaN in variables

[s_hy]

Hi all,

For this problem, I have checked all possible factor that may cause NaN in variables such as
- division by zero --> not found
- undefined variables - set initialization for variables - ep, hr,ht
- parameter setting --> double precision problems checked, no issue

but still couldn't find the cause of the NaN in variable ep, ht and hr as well. need an advice from professionals ( i am an amateur anyway)

!this is code from berenger 2002 for vlf propagation in earth-ionosphere waveguide.
!equations 7a-7c
!to date, ionosphere is assumed as free space -29/10/2013

!**************************************************************************************

subroutine spherical
implicit none

integer                                              :: N, L, i, j, t, k
integer                                              :: ir !, rr
double precision                                :: c, f, lamda, dx, dt, sigma, eps0,u0
double precision                                :: epsr,sigmar,omegap
double precision                                :: dr, dth,rmh,rph
double precision,dimension(300,300)   :: ep, ht, hr, ga, gb
double precision,dimension(300)         :: r
double precision, parameter               :: pi = 3.14159265d0
character(len=20)                              :: filename


!parameter
 k = 201
 L = 100

 f = 3.e9
 u0 = 4*pi*1e-7
 eps0 = 8.8e-12
 c = 1D0/sqrt(u0*eps0)
 lamda = c/f
! print *,'c=',c
 
 dx = lamda/20 
 dt = dx / (2*c)
 dr = 0.05*lamda
 dth = 1*pi/180          !theta  = pi
 
 !initialize r
 do j = 1,k+1
   r(j) = 0.0
   !print *, 'j=',j,'r(j)=',r(j)
 end do


!do r = 1,dr*(k-1),dr
 do ir = 1, (k-1)
    r(ir) = (ir - 1)*dr
 end do !IR

 !material parameter , assumed to be free space
  epsr = 1d0
  sigma = 1d0
  sigmar = 1d0
         

!initialization for updated equations
 do i = 1,k+1
  do j =1 ,k+1
    ep(i,j) = 0
    ht(i,j) = 0
    hr(i,j) = 0
  end do
 end do


!material properties
do i = 1, k                                  
    do j = 1,k
        ga(i,j) = exp(-(sigmar*dt)/(eps0*epsr)) 
        gb(i,j) = (1-exp(-(sigmar*dt)/(eps0*epsr)))/sigmar              
        !print *,'ga(i,j)=',ga(i,j),'gb(i,j)=',gb(i,j)
    end do
end do

!**********************************************************************
!**********************************************************************
 ! update equations
 do t=1,L
 
      write (filename, "('data',I3.3,'.dat')") t
      open (unit=130,file=filename)
            
      ! source
      ep(1,1) = sin(2*pi*f*t*dt)
      !print *, 'ep(1,1)=',ep(1,1)
     
    !  update Ez field
    do i = 2,k
     do j = 2,k
       rph = r(j)+dr/2
       rmh = r(j)-dr/2
       ep(i,j) = ga(i,j)*ep(i,j) + gb(i,j) * (((1/r(j)*dth)*(hr(i,j) - hr(i-1,j)) - (1/r(j)*dr)*(rph*ht(i,j) + rmh*ht(i,j-1))))
       !print*, 'i=',i,'j=',j,'r(j) = ',r(j)
       write (130,*) i,j,ep(i,j)
       if (j == 200) write (130,*) ' '
     end do
    end do
     
     !update Ht (theta) field =Hx
     do  i = 1,k-1
      do j = 1,k-1
         rph = r(j)+dr
         rmh = r(j)
         ht(i,j) = ht(i,j) - (dt/u0/r(j)/dr)*(rph*ep(i,j+1) - rmh*ep(i,j))
         hr(i,j) = hr(i,j) + (dt/u0/r(j)/dth/sin(dth))*(sin(dth)*ep(i+1,j) - sin(dth)*ep(i,j))
         print*,'i=',i,'j=',j,'hr(i,j)  =',hr(i,j) 
       end do
      end do
   
   
    close (unit=130)
           
 end do

end subroutine

 

[D H]

You are dividing by zero (infinity) and then multiplying by zero, which results in NaN.

!initialize r
 do j = 1,k+1
   r(j) = 0.0
   !print *, 'j=',j,'r(j)=',r(j)
 end do

!do r = 1,dr*(k-1),dr
 do ir = 1, (k-1)
    r(ir) = (ir - 1)*dr
 end do !IR

The first loop sets the first k+1 elements of r to zero. The second loop sets the first k-1 elements of r to (ir-1)*dr. Elements k and k+1? They're still zero. (Aside: r(1) is also zero.)

    !  update Ez field
    do i = 2,k
     do j = 2,k
       rph = r(j)+dr/2
       rmh = r(j)-dr/2
       ep(i,j) = ga(i,j)*ep(i,j) + gb(i,j) * (((1/r(j)*dth)*(hr(i,j) - hr(i-1,j)) - (1/r(j)*dr)*(rph*ht(i,j) + rmh*ht(i,j-1))))
       !print*, 'i=',i,'j=',j,'r(j) = ',r(j)
       write (130,*) i,j,ep(i,j)
       if (j == 200) write (130,*) ' '
     end do
    end do

When j=k that 1/r(j) is 1.0/0.0., or +infinity. On the first pass, hr(i,j) - hr(i-1,j) is 0, and infinity*0 is NaN.

 

[s_hy]

solved! thank you