- #1

s_hy

- 61

- 0

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)

Code:

```
!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
```