[其他行业算法] MT一维正演，用SOR法不收敛

program main
implicit none
integer,parameter::n=1000                       !地上n层 地下n层
complex*16,parameter::i=(0d0,1d0)               
real*16,parameter::PI=3.1415926535897932384626  
real*16,parameter::miu=4d0*PI*1E-7              
real*16,parameter::w=1.5                        !SOR迭代因子
integer::j,k
real*16::o(-n-1:n),o1(-n-1:n),dz(-n-1:n+1),s,f,E14,ro,fa
!o地层电导率，o1分界面电导率，dz层厚度，f频率，E14四分之一厚度处电场，ro视电阻率，fa相位
complex*16::E(-n-1:n+1),H(-n-1:n),E1(-n-1:n+1),Hs,a,b
!E电场，H磁场，Hs地表磁场
open(123,file='a.dat')
dz=1d0                    !单层厚度1m
dz(-n-1)=1.0E6            !最上层厚度1000km，视作无限远
dz(n)=1.0E6               !最底层厚度1000km，视作无限远
forall (k=-n-1:-1:1) o(k)=1E-10      !空气中电导率1E-10，视作无限大
forall(k=0:300:1) o(k)=0.1d0         !0到300米，电导率0.1
forall(k=301:600:1) o(k)=0.005d0     !301到600米，电导率0.005
forall(k=601:n:1) o(k)=0.1d0         !601到无穷远，电导率0.1
do k=-n,n
    o1(k)=(o(k)*dz(k)+o(k-1)*dz(k-1))/(dz(k)+dz(k-1))  !计算层界面电导率
enddo
do j=20,2000,20                     !j为厚度
    f=(500**2)*100/j                  !趋附深度求取频率
    E=(1d0,1d0)                     !迭代初值
    E(-n-1)=(1.0E-6,0d0)            !地上无穷远处电场1E-6v/m
    E(n+1)=(0d0,0d0)                !地下无穷远处电场为0
    do                              !SOR迭代法
        E1=E
        s=0
        do k=-n,n
            a=2d0*(E(k+1)/(dz(k)*dz(k)+dz(k-1)*dz(k))+E(k-1)/(dz(k-1)*dz(k)+dz(k-1)*dz(k-1)))
            b=2d0/(dz(k)*dz(k)+dz(k-1)*dz(k))+2d0/(dz(k-1)*dz(k)+dz(k-1)*dz(k-1))+i*f*miu*o1(k)
            E(k)=(1-w)*E(k)+w*a/b   
            s=s+abs(E(k)-E1(k))
        enddo
        if(s<1E-20)exit             !迭代结束条件
    enddo
    do k=-n-1,n
        H(k)=(E(k+1)-E(k))/(dz(k)*i*f*miu)      !计算各地层中点磁场
    enddo
    E14=0.75*E(0)+0.25d0*E(1)
    Hs=E14*o1(0)*dz(0)/2d0+H(0)                 !地表处磁场
    ro=(abs(E(0)/Hs))**2/f/miu                  !卡尼亚视电阻率
    write(*,*)f,ro
    write(123,*)f,ro
enddo
pause
endprogram

program main
implicit none
integer,parameter::n=1000                       !地上n层 地下n层
complex*16,parameter::i=(0d0,1d0)               
real*8,parameter::PI=3.1415926535897932384626d0  
real*8,parameter::miu=4.0d0*PI*1.0d-7              
real*8,parameter::w=0.8                        !SOR迭代因子
integer::j,k
real*8::o(-n-1:n),o1(-n:n),dz(-n-1:n+1),s,f,E14,ro,fa
!o地层电导率，o1分界面电导率，dz层厚度，f频率，E14四分之一厚度处电场，ro视电阻率，fa相位
complex*16::E(-n-1:n+1),H(-n-1:n),E1(-n-1:n+1),Hs,a,b
!E电场，H磁场，Hs地表磁场
open(123,file='a.dat')
dz=1d0                    !单层厚度1m
dz(-n-1)=1.0E6            !最上层厚度1000km，视作无限远
dz(n)=1.0E6               !最底层厚度1000km，视作无限远
forall (k=-n-1:-1:1) o(k)=1E-10      !空气中电导率1E-10，视作无限大
forall(k=0:300:1) o(k)=0.1d0         !0到300米，电导率0.1
forall(k=301:600:1) o(k)=0.005d0     !301到600米，电导率0.005
forall(k=601:n:1) o(k)=0.1d0         !601到无穷远，电导率0.1
do k=-n,n
    o1(k)=(o(k)*dz(k)+o(k-1)*dz(k-1))/(dz(k)+dz(k-1))  !计算层界面电导率
enddo
do j=20,2000,20                     !j为厚度
    f=2*PI*(500**2)*100/j           !趋附深度求取频率
    E=(1d0,1d0)                     !迭代初值
    !E(-n-1)=(1.0E-6,0d0)           !地上无穷远处电场1E-6v/m
    E(n+1)=(0d0,0d0)                !地下无穷远处电场为0
    do                              !SOR迭代法
        E1=E
        s=0
        do k=-n,n
            a=2d0*(E(k+1)/(dz(k)*dz(k)+dz(k-1)*dz(k))+E(k-1)/(dz(k-1)*dz(k)+dz(k-1)*dz(k-1)))
            b=2d0/(dz(k)*dz(k)+dz(k-1)*dz(k))+2d0/(dz(k-1)*dz(k)+dz(k-1)*dz(k-1))+i*f*miu*o1(k)
            E(k)=(1-w)*E(k)+w*a/b   
            s=s+abs(E(k)-E1(k))**2
        enddo
        s = dsqrt( s/n/2 ) 
        if(s<1E-10)exit             !迭代结束条件
    enddo
    do k=-n-1,n
        H(k)=(E(k+1)-E(k))/(dz(k)*i*f*miu)      !计算各地层中点磁场
    enddo
    E14=0.75*E(0)+0.25d0*E(1)
    Hs=E14*o1(0)*dz(0)/2d0+H(0)                 !地表处磁场
    ro=(abs(E(0)/Hs))**2/f/miu                  !卡尼亚视电阻率
    write(*,*)f,ro
    write(123,*)f,ro
enddo
pause
endprogram