一段代码放在subroutine中结果为nan，主程序中运算正常

石子

如题。由于计算需要，我将一段代码放入subtoutine计算，方便后续omp并行。这部分代码可以正常运行并且输出正确的结果，当我把他们放入子程序中然后通过call subroutine调用来更新某些数值时，关键的量总是NAN （经逐步输出发现是无穷大）。
详细信息：subroutine中i,j,Er,...E2i均为数组，在subroutine中不改变其值，并且在整个main程序中的初始值设置也是0。然而，经过若干次的subroutine中的迭代后会变为nan。在subroutine中（这个subroutine本质上是吸收边界，初始该区域并无源项，因此该区域应该始终为0），很多量来源于这几个数组的梯度、拉普拉斯算子等，而这些量Hr,Hi以及数组Er,Ei...会随着迭代次数t的增加而发散。

由于该段代码一直正常且放到subroutine外之后才出现该问题，因此大概率是程序问题，但我还是不太理解原因是什么。求助。

[Fortran] 纯文本查看 复制代码

call cal_cor1(i,j,Er,Ei,E2r,E2i,V,uper,upei,uper2,upei2)

[Fortran] 纯文本查看 复制代码

    subroutine cal_cor1(i,j,vEr,vEi,vE2r,vE2i,V,uper,upei,uper2,upei2)
integer,intent(in)::i,j
    double precision,dimension(1:Nx,1:Ny),intent(in)::vEr,vEi,vE2r,vE2i,V
    double precision,intent(inout)::uper,upei,uper2,upei2
    double precision::tempuper,tempupei,tempuper2,tempupei2
    double precision::coe1,coe2,coe3,coe4,coe11,coe22,coe33,coe44,temp,ytemp
    double precision::Hr,Hi,Hr2,Hi2,yHr,yHi,yHr2,yHi2
        temp=sigma*((i*dx-xleft)/depth)**index
        coe1=-(index*temp*(depth/(i*dx-xleft))*(1+6*temp*(1+temp)))/(depth*(1+2*temp*(1+temp))**3)
        coe2=(index*temp*depth/(i*dx-xleft)*(-1+6*temp**2+4*temp**3))/(depth*(1+2*temp+2*temp**2)**3)
        coe3=(1+2*temp)/(1+2*temp+2*temp**2)**2
        coe4=-temp*(1+temp)/(1+2*temp+2*temp**2)**2
        Hr=-0.5*coe2*((vEr(i+1,j)-vEr(i-1,j))/2/dx)-coe4*(vEr(i+1,j)+vEr(i-1,j)-2*vEr(i,j))/dx**2
        Hi=(-0.5*coe1*((vEi(i+1,j)-vEi(i-1,j))/2/dx)-0.5*coe3*(vEi(i+1,j)+vEi(i-1,j)-2*vEi(i,j))/dx**2)-V(i,j)*vEi(i,j)
        Hr2=(-0.5*coe2*((vE2r(i+1,j)-vE2r(i-1,j))/2/dx)-coe4*(vE2r(i+1,j)+vE2r(i-1,j)-2*vE2r(i,j))/dx**2)
        Hi2=(-0.5*coe1*((vE2i(i+1,j)-vE2i(i-1,j))/2/dx)-0.5*coe3*(vE2i(i+1,j)+vE2i(i-1,j)-2*vE2i(i,j))/dx**2)&
            -V(i,j)*vE2i(i,j)
        ytemp=ysigma*((j*dy-ybottom)/ydepth)**index
        coe11=-(index*ytemp*(ydepth/(j*dy-ybottom))*(1+6*ytemp*(1+ytemp)))/(ydepth*(1+2*ytemp*(1+ytemp))**3)
        coe22=(index*ytemp*ydepth/(j*dy-ybottom)*(-1+6*ytemp**2+4*ytemp**3))/(ydepth*(1+2*ytemp+2*ytemp**2)**3)
        coe33=(1+2*ytemp)/(1+2*ytemp+2*ytemp**2)**2
        coe44=-ytemp*(1+ytemp)/(1+2*ytemp+2*ytemp**2)**2
        yHr=(-0.5*coe22*((vEr(i,j+1)-vEr(i,j-1))/2/dy)-coe44*(vEr(i,j+1)+vEr(i,j-1)-2*vEr(i,j))/dy**2)
        yHi=(-0.5*coe11*((vEi(i,j+1)-vEi(i,j-1))/2/dy)-0.5*coe33*(vEi(i,j+1)+vEi(i,j-1)-2*vEi(i,j))/dy**2)
        yHr2=(-0.5*coe22*((vE2r(i,j+1)-vE2r(i,j-1))/2/dy)-coe44*(vE2r(i,j+1)+vE2r(i,j-1)-2*vE2r(i,j))/dy**2)
        yHi2=(-0.5*coe11*((vE2i(i,j+1)-vE2i(i,j-1))/2/dy)-0.5*coe33*(vE2i(i,j+1)+vE2i(i,j-1)-2*vE2i(i,j))/dy**2)
        tempuper=(Hr+Hi+yHr+yHi)
        tempuper2=(Hr2+Hi2+yHr2+yHi2)
            ! imaginary part
。。。
tempupei=Hr+Hi+yHr+yHi
        tempupei2=Hr2+Hi2+yHr2+yHi2
        if (.not. isnan(tempuper) .and. .not. isnan(tempupei) .and. .not. isnan(tempuper2) .and. .not. isnan(tempupei2)) then
            uper=vEr(i,j)+dt*tempuper
            upei=vEi(i,j)+dt*tempupei
            uper2=vE2r(i,j)+dt*tempuper2
            upei2=vE2i(i,j)+dt*tempupei2
        else
            ! print*, "NaN detected in tempuper at (", i, ",", j, ")"
        end if
    end subroutine cal_cor1