openSource/tem3dfdtd-open
4
0
派生 0
代码 项目 活动
文件
a1bfcdd26bae133c936d11298d48f7571b83576a
tem3dfdtd-open/tem3dfdtd/lib/Iteration.f90
T
sunhf 89ff871b0b 修正代码中的注释错误
2026-08-14 20:51:58 +08:00

706 行
42 KiB
Fortran
原始文件 Blame 文件历史

!Copyright (c) 2013 by https://git.em3d.cn/ under guide of Xiu Li(lixiu@chd.edu.cn)
!written by Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn) and Xushan Lu(luxushan@gmail.com)
!Code distribution @ https://git.em3d.cn/
! --------------------------------Subroutine part---------------------------------------------!
subroutine Iteration_cpml
use constantparameters
USE CONSTANTPARAMETERS
USE ELECTROMAGNETIC_VARIABLES
USE RES_MODEL_PARAMETER
USE TIME_PARAMETER
USE PML_PARAMETER
implicit none
real::t1,t2,t,t_start,t_end,t_total !t1 denotes original cpu time at the beginning of each computation fraction, t2 denotes the end cpu time and t=t2-t1
REAL*8 CA,CB,DELX1,DELY1,DELZ1 !ca, cb, delx1, dely1, delz1 are all middle variables used in the computation of EM field
REAL*8 TEMP_SIG,temp_cacb,data_rec(point_num) !Temp_sig and temp_cacb are middle variables used in the computation of EM field
REAL*8 DELY2,DELZ2,delx2 !They are all middle variables as above ones.
integer num,i,j,k,ii,iii,jj,kk,idx_write,x_pos_observer(8),y_pos_observer(8),z_pos_observer(8) !num is the number of computation fraction
integer :: N_hight=0
real*8,allocatable::Meps_r(:),Mdelt(:),Msource(:),Mcq(:) !They are local substitution of eps_r, delt and cq
REAL*8 hz_observer(8)
CHARACTER*20::string,str_num
WRITE(*,*)'[Iteration_cpml] Boundary condition: CPML absorbing boundary (PML = unbounded absorbing layer, uniform grid required)'
WRITE(*,*)'[Iteration_cpml] Iteration starts .. .. .. ..'
!Create output files
idx_start=12000
do iii=1,point_num
idx_write=idx_start+iii
IF(iii<10) THEN
write(str_num,"(I1)")iii
ELSEif(iii<100) THEN
write(str_num,"(I2)")iii
ELSEif(iii<1000) THEN
write(str_num,"(I3)")iii
ENDIF
string='dBzdt'//"_"//trim(str_num)//'.txt'
open(idx_write,file=string)
write(idx_write,*)"point_"//trim(str_num)
write(idx_write,*)Points_observer(iii)%local_coord_to_source%coord_x,Points_observer(iii)%local_coord_to_source%coord_y,&
Points_observer(iii)%local_coord_to_source%coord_z
enddo
call cpu_time(t_start)
!OPEN(20250220,file='dBzdt.txt')
do num=1,num_fra_com,1 !The outer loop which begins from the first fraction ends at the last fraction
call cpu_time(t1) !Record the cpu time at the beginning of each computing fraction
allocate(mdelt(0:mstop(num)),meps_r(mstop(num)),mcq(mstop(num)),msource(mstop(num)))
! The memory of mdelt, meps_r, mcq and msource are allocated at the begining of fraction
do ii=mstart(num),mstart(num)+mstop(num)-1,1
mdelt(ii-mstart(num)+1)=delt(ii)
meps_r(ii-mstart(num)+1)=eps_r(ii)
mcq(ii-mstart(num)+1)=cq(ii)
msource(ii-mstart(num)+1)=source(ii) !Link the local value of mdelt, meps_r, mcq and msorce to the global value of delt, eps_r, cq and source array.
end do
print*,'Now computing fraction:',num
mdelt(0)=mdelt(1)
do loop=1,mstop(num),1
! --------------------------------CPML coefficients b/c of the current time step-------------------------------!
! b = exp(-(sig/kappa+alpha)*delt/eps0), c = sig*(b-1)/(sig+kappa*alpha)/kappa
! They are recomputed at every step because MDELT changes between the raise, wave, ramp and off phases.
! The sigma/alpha/kappa profiles are built once by Get_pml_parameters.
! When Logic_PML=0 this whole block is skipped and the original Dirichlet boundary scheme is used.
DO i=1,PML_X1
b_e_x1(i)=DEXP(-(sig_PML_e_x1(i)/kappa_PML_e_x1(i)+alpha_PML_e_x1(i))*MDELT(LOOP-1)/EPS0)
IF(sig_PML_e_x1(i)==0.0 .AND. alpha_PML_e_x1(i)==0.0 .AND. i==PML_X1)THEN
c_e_x1(i)=0.0
ELSE
c_e_x1(i)=sig_PML_e_x1(i)*(b_e_x1(i)-1.0)/(sig_PML_e_x1(i)+kappa_PML_e_x1(i)*alpha_PML_e_x1(i))/kappa_PML_e_x1(i)
ENDIF
ENDDO
DO ii=1,PML_X1-1
b_h_x1(ii)=DEXP(-(sig_PML_h_x1(ii)/kappa_PML_h_x1(ii)+alpha_PML_h_x1(ii))*MDELT(LOOP-1)/EPS0)
c_h_x1(ii)=sig_PML_h_x1(ii)*(b_h_x1(ii)-1.0)/(sig_PML_h_x1(ii)+kappa_PML_h_x1(ii)*alpha_PML_h_x1(ii))/kappa_PML_h_x1(ii)
ENDDO
DO i=1,PML_X2
b_e_x2(i)=DEXP(-(sig_PML_e_x2(i)/kappa_PML_e_x2(i)+alpha_PML_e_x2(i))*MDELT(LOOP-1)/EPS0)
IF(sig_PML_e_x2(i)==0.0 .AND. alpha_PML_e_x2(i)==0.0 .AND. i==PML_X2)THEN
c_e_x2(i)=0.0
ELSE
c_e_x2(i)=sig_PML_e_x2(i)*(b_e_x2(i)-1.0)/(sig_PML_e_x2(i)+kappa_PML_e_x2(i)*alpha_PML_e_x2(i))/kappa_PML_e_x2(i)
ENDIF
ENDDO
DO ii=1,PML_X2-1
b_h_x2(ii)=DEXP(-(sig_PML_h_x2(ii)/kappa_PML_h_x2(ii)+alpha_PML_h_x2(ii))*MDELT(LOOP-1)/EPS0)
c_h_x2(ii)=sig_PML_h_x2(ii)*(b_h_x2(ii)-1.0)/(sig_PML_h_x2(ii)+kappa_PML_h_x2(ii)*alpha_PML_h_x2(ii))/kappa_PML_h_x2(ii)
ENDDO
DO j=1,PML_Y1
b_e_y1(j)=DEXP(-(sig_PML_e_y1(j)/kappa_PML_e_y1(j)+alpha_PML_e_y1(j))*MDELT(LOOP-1)/EPS0)
IF(sig_PML_e_y1(j)==0.0 .AND. alpha_PML_e_y1(j)==0.0 .AND. j==PML_Y1)THEN
c_e_y1(j)=0.0
ELSE
c_e_y1(j)=sig_PML_e_y1(j)*(b_e_y1(j)-1.0)/(sig_PML_e_y1(j)+kappa_PML_e_y1(j)*alpha_PML_e_y1(j))/kappa_PML_e_y1(j)
ENDIF
ENDDO
DO jj=1,PML_Y1-1
b_h_y1(jj)=DEXP(-(sig_PML_h_y1(jj)/kappa_PML_h_y1(jj)+alpha_PML_h_y1(jj))*MDELT(LOOP-1)/EPS0)
c_h_y1(jj)=sig_PML_h_y1(jj)*(b_h_y1(jj)-1.0)/(sig_PML_h_y1(jj)+kappa_PML_h_y1(jj)*alpha_PML_h_y1(jj))/kappa_PML_h_y1(jj)
ENDDO
DO j=1,PML_Y2
b_e_y2(j)=DEXP(-(sig_PML_e_y2(j)/kappa_PML_e_y2(j)+alpha_PML_e_y2(j))*MDELT(LOOP-1)/EPS0)
IF(sig_PML_e_y2(j)==0.0 .AND. alpha_PML_e_y2(j)==0.0 .AND. j==PML_Y2)THEN
c_e_y2(j)=0.0
ELSE
c_e_y2(j)=sig_PML_e_y2(j)*(b_e_y2(j)-1.0)/(sig_PML_e_y2(j)+kappa_PML_e_y2(j)*alpha_PML_e_y2(j))/kappa_PML_e_y2(j)
ENDIF
ENDDO
DO jj=1,PML_Y2-1
b_h_y2(jj)=DEXP(-(sig_PML_h_y2(jj)/kappa_PML_h_y2(jj)+alpha_PML_h_y2(jj))*MDELT(LOOP-1)/EPS0)
c_h_y2(jj)=sig_PML_h_y2(jj)*(b_h_y2(jj)-1.0)/(sig_PML_h_y2(jj)+kappa_PML_h_y2(jj)*alpha_PML_h_y2(jj))/kappa_PML_h_y2(jj)
ENDDO
DO k=1,PML_Z1
b_e_z1(k)=DEXP(-(sig_PML_e_z1(k)/kappa_PML_e_z1(k)+alpha_PML_e_z1(k))*MDELT(LOOP-1)/EPS0)
IF(sig_PML_e_z1(k)==0.0 .AND. alpha_PML_e_z1(k)==0.0 .AND. k==PML_Z1)THEN
c_e_z1(k)=0.0
ELSE
c_e_z1(k)=sig_PML_e_z1(k)*(b_e_z1(k)-1.0)/(sig_PML_e_z1(k)+kappa_PML_e_z1(k)*alpha_PML_e_z1(k))/kappa_PML_e_z1(k)
ENDIF
ENDDO
DO kk=1,PML_Z1-1
b_h_z1(kk)=DEXP(-(sig_PML_h_z1(kk)/kappa_PML_h_z1(kk)+alpha_PML_h_z1(kk))*MDELT(LOOP-1)/EPS0)
c_h_z1(kk)=sig_PML_h_z1(kk)*(b_h_z1(kk)-1.0)/(sig_PML_h_z1(kk)+kappa_PML_h_z1(kk)*alpha_PML_h_z1(kk))/kappa_PML_h_z1(kk)
ENDDO
DO k=1,PML_Z2
b_e_z2(k)=DEXP(-(sig_PML_e_z2(k)/kappa_PML_e_z2(k)+alpha_PML_e_z2(k))*MDELT(LOOP-1)/EPS0)
IF(sig_PML_e_z2(k)==0.0 .AND. alpha_PML_e_z2(k)==0.0 .AND. k==PML_Z2)THEN
c_e_z2(k)=0.0
ELSE
c_e_z2(k)=sig_PML_e_z2(k)*(b_e_z2(k)-1.0)/(sig_PML_e_z2(k)+kappa_PML_e_z2(k)*alpha_PML_e_z2(k))/kappa_PML_e_z2(k)
ENDIF
ENDDO
DO kk=1,PML_Z2-1
b_h_z2(kk)=DEXP(-(sig_PML_h_z2(kk)/kappa_PML_h_z2(kk)+alpha_PML_h_z2(kk))*MDELT(LOOP-1)/EPS0)
c_h_z2(kk)=sig_PML_h_z2(kk)*(b_h_z2(kk)-1.0)/(sig_PML_h_z2(kk)+kappa_PML_h_z2(kk)*alpha_PML_h_z2(kk))/kappa_PML_h_z2(kk)
ENDDO
!Assemble c_h_zz used by the Hz recursion in the z direction.
DO k=1,PML_Z1-1
c_h_zz(k)=c_h_z1(k)
ENDDO
DO k=NZ+2-PML_Z2,NZ
c_h_zz(k)=c_h_z2(NZ+1-k)
ENDDO
!Precompute the inverse denominator of the Hz recursion once per step
!(bit-for-bit neutral when Logic_PML=0: inv_hz_den stays 1.0 from ZERO).
DO k=1,NZ
inv_hz_den(k)=1.0D0/(den_hz(k)+c_h_zz(k))
ENDDO
! --------------------------------update the value of Ex ---------------------------------------!
! 忠实移植自参考版 tem3dfdtd_第二版（孙师兄版）：psi 内嵌在场更新循环内，
! 循环结构 DO I / DO K / DO J；源项仅在源平面 K=NZS+1 施加（当前项目 2D 掩码，
! 等价参考版 3D 掩码在非源平面层为 0）。
DO I=1,NX
DO K=2,NZB-1
DO J=2,NYB-1
DELY1=(CDELY(J-1)+CDELY(J))/2.0D0
DELZ1=(CDELZ(K-1)+CDELZ(K))/2.0D0
CA=(2.0D0*Meps_r(loop)-Mdelt(LOOP-1)*CCSIGX(I,J,K))/(2.0D0*Meps_r(loop)+Mdelt(LOOP-1)*CCSIGX(I,J,K))
CB=(2.0D0*MDELT(LOOP-1))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGX(I,J,K))
EX(I,J,K)=CA*EX(I,J,K)+CB*((HZ(I,J,K)-HZ(I,J-1,K))*den_ey(J)/DELY1&
&-(HY(I,J,K)-HY(I,J,K-1))*den_ez(K)/DELZ1)
IF(K==NZS+1)THEN
EX(I,J,K)=EX(I,J,K)-CB*Msource(loop)*is_ex_in_source(I,J)
ENDIF
! PML for Ex, y-direction
IF(J<=PML_Y1)THEN
psi_Exy_1(I,J,K)=b_e_y1(J)*psi_Exy_1(I,J,K)+c_e_y1(J)*(HZ(I,J,K)-HZ(I,J-1,K))/DELY1
EX(I,J,K)=EX(I,J,K)+CB*psi_Exy_1(I,J,K)
ELSEIF(J>=NY+2-PML_Y2)THEN
psi_Exy_2(I,NY+2-J,K)=b_e_y2(NY+2-J)*psi_Exy_2(I,NY+2-J,K)+c_e_y2(NY+2-J)*(HZ(I,J,K)-HZ(I,J-1,K))/DELY1
EX(I,J,K)=EX(I,J,K)+CB*psi_Exy_2(I,NY+2-J,K)
ENDIF
! PML for Ex, z-direction
IF(K<=PML_Z1)THEN
psi_Exz_1(I,J,K)=b_e_z1(K)*psi_Exz_1(I,J,K)+c_e_z1(K)*(HY(I,J,K)-HY(I,J,K-1))/DELZ1
EX(I,J,K)=EX(I,J,K)-CB*psi_Exz_1(I,J,K)
ELSEIF(K>=NZ+2-PML_Z2)THEN
psi_Exz_2(I,J,NZ+2-K)=b_e_z2(NZ+2-K)*psi_Exz_2(I,J,NZ+2-K)+c_e_z2(NZ+2-K)*(HY(I,J,K)-HY(I,J,K-1))/DELZ1
EX(I,J,K)=EX(I,J,K)-CB*psi_Exz_2(I,J,NZ+2-K)
ENDIF
ENDDO
ENDDO
ENDDO
!===============end of updating Ex=========================!
! --------------------------------update the value of Ey ---------------------------------------!
DO J=1,NY
DO K=2,NZB-1
DO I=2,NXB-1
DELX1=(CDELX(I-1)+CDELX(I))/2.0D0
DELZ1=(CDELZ(K-1)+CDELZ(K))/2.0D0
CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*CCSIGY(I,J,K))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGY(I,J,K))
CB=(2.0D0*MDELT(LOOP-1))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGY(I,J,K))
EY(I,J,K)=CA*EY(I,J,K)+CB*((HX(I,J,K)-HX(I,J,K-1))*den_ez(K)/DELZ1&
&-(HZ(I,J,K)-HZ(I-1,J,K))*den_ex(I)/DELX1)
IF(K==NZS+1)THEN
EY(I,J,K)=EY(I,J,K)-CB*Msource(loop)*is_ey_in_source(I,J)
ENDIF
! PML for Ey, z-direction
IF(K<=PML_Z1)THEN
psi_Eyz_1(I,J,K)=b_e_z1(K)*psi_Eyz_1(I,J,K)+c_e_z1(K)*(HX(I,J,K)-HX(I,J,K-1))/DELZ1
EY(I,J,K)=EY(I,J,K)+CB*psi_Eyz_1(I,J,K)
ELSEIF(K>=NZ+2-PML_Z2)THEN
psi_Eyz_2(I,J,NZ+2-K)=b_e_z2(NZ+2-K)*psi_Eyz_2(I,J,NZ+2-K)+c_e_z2(NZ+2-K)*(HX(I,J,K)-HX(I,J,K-1))/DELZ1
EY(I,J,K)=EY(I,J,K)+CB*psi_Eyz_2(I,J,NZ+2-K)
ENDIF
! PML for Ey, x-direction
IF(I<=PML_X1)THEN
psi_Eyx_1(I,J,K)=b_e_x1(I)*psi_Eyx_1(I,J,K)+c_e_x1(I)*(HZ(I,J,K)-HZ(I-1,J,K))/DELX1
EY(I,J,K)=EY(I,J,K)-CB*psi_Eyx_1(I,J,K)
ELSEIF(I>=NX+2-PML_X2)THEN
psi_Eyx_2(NX+2-I,J,K)=b_e_x2(NX+2-I)*psi_Eyx_2(NX+2-I,J,K)+c_e_x2(NX+2-I)*(HZ(I,J,K)-HZ(I-1,J,K))/DELX1
EY(I,J,K)=EY(I,J,K)-CB*psi_Eyx_2(NX+2-I,J,K)
ENDIF
ENDDO
ENDDO
ENDDO
!===============end of updating Ey===================!
! -------------------------------------update the value of Ez--------------------------------------!
DO K=1,NZ
DO J=2,NYB-1
DO I=2,NXB-1
DELX1=(CDELX(I-1)+CDELX(I))/2.0D0
DELY1=(CDELY(J-1)+CDELY(J))/2.0D0
TEMP_CACB=2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGZ(I,J,K)
CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*CCSIGZ(I,J,K))/TEMP_CACB
CB=(2.0D0*MDELT(LOOP-1))/TEMP_CACB
EZ(I,J,K)=CA*EZ(I,J,K)+CB*((HY(I,J,K)-HY(I-1,J,K))*den_ex(I)/DELX1&
&-(HX(I,J,K)-HX(I,J-1,K))*den_ey(J)/DELY1)
! PML for Ez, x-direction
IF(I<=PML_X1)THEN
psi_Ezx_1(I,J,K)=b_e_x1(I)*psi_Ezx_1(I,J,K)+c_e_x1(I)*(HY(I,J,K)-HY(I-1,J,K))/DELX1
EZ(I,J,K)=EZ(I,J,K)+CB*psi_Ezx_1(I,J,K)
ELSEIF(I>=NX+2-PML_X2)THEN
psi_Ezx_2(NX+2-I,J,K)=b_e_x2(NX+2-I)*psi_Ezx_2(NX+2-I,J,K)+c_e_x2(NX+2-I)*(HY(I,J,K)-HY(I-1,J,K))/DELX1
EZ(I,J,K)=EZ(I,J,K)+CB*psi_Ezx_2(NX+2-I,J,K)
ENDIF
! PML for Ez, y-direction
IF(J<=PML_Y1)THEN
psi_Ezy_1(I,J,K)=b_e_y1(J)*psi_Ezy_1(I,J,K)+c_e_y1(J)*(HX(I,J,K)-HX(I,J-1,K))/DELY1
EZ(I,J,K)=EZ(I,J,K)-CB*psi_Ezy_1(I,J,K)
ELSEIF(J>=NY+2-PML_Y2)THEN
psi_Ezy_2(I,NY+2-J,K)=b_e_y2(NY+2-J)*psi_Ezy_2(I,NY+2-J,K)+c_e_y2(NY+2-J)*(HX(I,J,K)-HX(I,J-1,K))/DELY1
EZ(I,J,K)=EZ(I,J,K)-CB*psi_Ezy_2(I,NY+2-J,K)
ENDIF
ENDDO
ENDDO
ENDDO
!===============end of updating Ez=========================!
! ------------------------------------update the value of Hx-----------------------------------------------!
DO I=1,NXB
DO K=1,NZ
DO J=1,NY
HX(I,J,K)=HX(I,J,K)-MCQ(LOOP)*((EZ(I,J+1,K)-EZ(I,J,K))*den_hy(J))/CDELY(J)&
&+MCQ(LOOP)*((EY(I,J,K+1)-EY(I,J,K))*den_hz(K))/CDELZ(K)
! PML for Hx, y-direction
IF(J<=PML_Y1-1)THEN
psi_Hxy_1(I,J,K)=b_h_y1(J)*psi_Hxy_1(I,J,K)+c_h_y1(J)*(EZ(I,J+1,K)-EZ(I,J,K))/CDELY(J)
HX(I,J,K)=HX(I,J,K)-MCQ(LOOP)*psi_Hxy_1(I,J,K)
ELSEIF(J>=NY+2-PML_Y2)THEN
psi_Hxy_2(I,NY+1-J,K)=b_h_y2(NY+1-J)*psi_Hxy_2(I,NY+1-J,K)+c_h_y2(NY+1-J)*(EZ(I,J+1,K)-EZ(I,J,K))/CDELY(J)
HX(I,J,K)=HX(I,J,K)-MCQ(LOOP)*psi_Hxy_2(I,NY+1-J,K)
ENDIF
! PML for Hx, z-direction
IF(K<=PML_Z1-1)THEN
psi_Hxz_1(I,J,K)=b_h_z1(K)*psi_Hxz_1(I,J,K)+c_h_z1(K)*(EY(I,J,K+1)-EY(I,J,K))/CDELZ(K)
HX(I,J,K)=HX(I,J,K)+MCQ(LOOP)*psi_Hxz_1(I,J,K)
ELSEIF(K>=NZ+2-PML_Z2)THEN
psi_Hxz_2(I,J,NZ+1-K)=b_h_z2(NZ+1-K)*psi_Hxz_2(I,J,NZ+1-K)+c_h_z2(NZ+1-K)*(EY(I,J,K+1)-EY(I,J,K))/CDELZ(K)
HX(I,J,K)=HX(I,J,K)+MCQ(LOOP)*psi_Hxz_2(I,J,NZ+1-K)
ENDIF
ENDDO
ENDDO
ENDDO
!================end of updating Hx=======================!
! -------------------------------------update the value of Hy---------------------------------------------!
DO J=1,NYB
DO K=1,NZ
DO I=1,NX
HY(I,J,K)=HY(I,J,K)-MCQ(LOOP)*((EX(I,J,K+1)-EX(I,J,K))*den_hz(K))/CDELZ(K)&
&+MCQ(LOOP)*((EZ(I+1,J,K)-EZ(I,J,K))*den_hx(I))/CDELX(I)
! PML for Hy, x-direction
IF(I<=PML_X1-1)THEN
psi_Hyx_1(I,J,K)=b_h_x1(I)*psi_Hyx_1(I,J,K)+c_h_x1(I)*(EZ(I+1,J,K)-EZ(I,J,K))/CDELX(I)
HY(I,J,K)=HY(I,J,K)+MCQ(LOOP)*psi_Hyx_1(I,J,K)
ELSEIF(I>=NX+2-PML_X2)THEN
psi_Hyx_2(NX+1-I,J,K)=b_h_x2(NX+1-I)*psi_Hyx_2(NX+1-I,J,K)+c_h_x2(NX+1-I)*(EZ(I+1,J,K)-EZ(I,J,K))/CDELX(I)
HY(I,J,K)=HY(I,J,K)+MCQ(LOOP)*psi_Hyx_2(NX+1-I,J,K)
ENDIF
! PML for Hy, z-direction
IF(K<=PML_Z1-1)THEN
psi_Hyz_1(I,J,K)=b_h_z1(K)*psi_Hyz_1(I,J,K)+c_h_z1(K)*(EX(I,J,K+1)-EX(I,J,K))/CDELZ(K)
HY(I,J,K)=HY(I,J,K)-MCQ(LOOP)*psi_Hyz_1(I,J,K)
ELSEIF(K>=NZ+2-PML_Z2)THEN
psi_Hyz_2(I,J,NZ+1-K)=b_h_z2(NZ+1-K)*psi_Hyz_2(I,J,NZ+1-K)+c_h_z2(NZ+1-K)*(EX(I,J,K+1)-EX(I,J,K))/CDELZ(K)
HY(I,J,K)=HY(I,J,K)-MCQ(LOOP)*psi_Hyz_2(I,J,NZ+1-K)
ENDIF
ENDDO
ENDDO
ENDDO
!===============end of updating Hy========================!
!-------------------------------------update the value of Hz----------------------------------------------!
! 上扫 k=1..NZS-1（参考版原样）：HZ(I,J,K+1) 由 HZ(I,J,K) 推出，psi 修正内嵌。
DO K=1,NZs-1
DO I=1,NX
DO J=1,NY
HZ(I,J,K+1)=HZ(I,J,K)-((CDELZ(K)*den_hx(I))/(den_hz(K)+c_h_zz(K)))*((HX(I+1,J,K)-HX(I,J,K))/CDELX(I))&
&-((CDELZ(K)*den_hy(J))/(den_hz(K)+c_h_zz(K)))*((HY(I,J+1,K)-HY(I,J,K))/CDELY(J))
! PML for Hz up, x-direction
IF(I<=PML_X1-1)THEN
psi_Hzx_1(I,J,K)=b_h_x1(I)*psi_Hzx_1(I,J,K)+c_h_x1(I)*((HX(I+1,J,K)-HX(I,J,K))/CDELX(I))
HZ(I,J,K+1)=HZ(I,J,K+1)-((CDELZ(K))/(den_hz(K)+c_h_zz(K)))*psi_Hzx_1(I,J,K)
ELSEIF(I>=NX+2-PML_X2)THEN
psi_Hzx_2(NX+1-I,J,K)=b_h_x2(NX+1-I)*psi_Hzx_2(NX+1-I,J,K)+c_h_x2(NX+1-I)*((HX(I+1,J,K)-HX(I,J,K))/CDELX(I))
HZ(I,J,K+1)=HZ(I,J,K+1)-((CDELZ(K))/(den_hz(K)+c_h_zz(K)))*psi_Hzx_2(NX+1-I,J,K)
ENDIF
! PML for Hz up, y-direction
IF(J<=PML_Y1-1)THEN
psi_Hzy_1(I,J,K)=b_h_y1(J)*psi_Hzy_1(I,J,K)+c_h_y1(J)*((HY(I,J+1,K)-HY(I,J,K))/CDELY(J))
HZ(I,J,K+1)=HZ(I,J,K+1)-((CDELZ(K))/(den_hz(K)+c_h_zz(K)))*psi_Hzy_1(I,J,K)
ELSEIF(J>=NY+2-PML_Y2)THEN
psi_Hzy_2(I,NY+1-J,K)=b_h_y2(NY+1-J)*psi_Hzy_2(I,NY+1-J,K)+c_h_y2(NY+1-J)*((HY(I,J+1,K)-HY(I,J,K))/CDELY(J))
HZ(I,J,K+1)=HZ(I,J,K+1)-((CDELZ(K))/(den_hz(K)+c_h_zz(K)))*psi_Hzy_2(I,NY+1-J,K)
ENDIF
ENDDO
ENDDO
ENDDO
! 下扫 k=NZ..NZS+1（参考版原样）：HZ(I,J,K) 由 HZ(I,J,K+1) 推出（HZ(NZ+1) 保持 0），psi 修正内嵌。
DO K=NZ,NZs+1,-1
DO I=1,NX
DO J=1,NY
HZ(I,J,K)=HZ(I,J,K+1)+((CDELZ(K)*den_hx(I))/(den_hz(K)+c_h_zz(K)))*((HX(I+1,J,K)-HX(I,J,K))/CDELX(I))&
&+((CDELZ(K)*den_hy(J))/(den_hz(K)+c_h_zz(K)))*((HY(I,J+1,K)-HY(I,J,K))/CDELY(J))
! PML for Hz down, x-direction
IF(I<=PML_X1-1)THEN
psi_Hzx_1(I,J,K)=b_h_x1(I)*psi_Hzx_1(I,J,K)+c_h_x1(I)*((HX(I+1,J,K)-HX(I,J,K))/CDELX(I))
HZ(I,J,K)=HZ(I,J,K)+((CDELZ(K))/(den_hz(K)+c_h_zz(K)))*psi_Hzx_1(I,J,K)
ELSEIF(I>=NX+2-PML_X2)THEN
psi_Hzx_2(NX+1-I,J,K)=b_h_x2(NX+1-I)*psi_Hzx_2(NX+1-I,J,K)+c_h_x2(NX+1-I)*((HX(I+1,J,K)-HX(I,J,K))/CDELX(I))
HZ(I,J,K)=HZ(I,J,K)+((CDELZ(K))/(den_hz(K)+c_h_zz(K)))*psi_Hzx_2(NX+1-I,J,K)
ENDIF
! PML for Hz down, y-direction
IF(J<=PML_Y1-1)THEN
psi_Hzy_1(I,J,K)=b_h_y1(J)*psi_Hzy_1(I,J,K)+c_h_y1(J)*((HY(I,J+1,K)-HY(I,J,K))/CDELY(J))
HZ(I,J,K)=HZ(I,J,K)+((CDELZ(K))/(den_hz(K)+c_h_zz(K)))*psi_Hzy_1(I,J,K)
ELSEIF(J>=NY+2-PML_Y2)THEN
psi_Hzy_2(I,NY+1-J,K)=b_h_y2(NY+1-J)*psi_Hzy_2(I,NY+1-J,K)+c_h_y2(NY+1-J)*((HY(I,J+1,K)-HY(I,J,K))/CDELY(J))
HZ(I,J,K)=HZ(I,J,K)+((CDELZ(K))/(den_hz(K)+c_h_zz(K)))*psi_Hzy_2(I,NY+1-J,K)
ENDIF
ENDDO
ENDDO
ENDDO
!===================end of updating Hz==========================!
DO i=1,point_num
!=================================point1=======================================
x_pos_observer(1)=Points_observer(i)%global_coordmesh(1)%coordmesh_x
y_pos_observer(1)=Points_observer(i)%global_coordmesh(1)%coordmesh_y
z_pos_observer(1)=Points_observer(i)%global_coordmesh(1)%coordmesh_z
hz_observer(1)=(EX(x_pos_observer(1),y_pos_observer(1)+1,z_pos_observer(1))-EX(x_pos_observer(1),y_pos_observer(1),z_pos_observer(1)))/CDELY(y_pos_observer(1))-&
(EY(x_pos_observer(1)+1,y_pos_observer(1),z_pos_observer(1))-EY(x_pos_observer(1),y_pos_observer(1),z_pos_observer(1)))/CDELX(x_pos_observer(1))
!=================================point2=======================================
x_pos_observer(2)=Points_observer(i)%global_coordmesh(2)%coordmesh_x
y_pos_observer(2)=Points_observer(i)%global_coordmesh(2)%coordmesh_y
z_pos_observer(2)=Points_observer(i)%global_coordmesh(2)%coordmesh_z
hz_observer(2)=(EX(x_pos_observer(2),y_pos_observer(2)+1,z_pos_observer(2))-EX(x_pos_observer(2),y_pos_observer(2),z_pos_observer(2)))/CDELY(y_pos_observer(2))-&
(EY(x_pos_observer(2)+1,y_pos_observer(2),z_pos_observer(2))-EY(x_pos_observer(2),y_pos_observer(2),z_pos_observer(2)))/CDELX(x_pos_observer(2))
!=================================point3=======================================
x_pos_observer(3)=Points_observer(i)%global_coordmesh(3)%coordmesh_x
y_pos_observer(3)=Points_observer(i)%global_coordmesh(3)%coordmesh_y
z_pos_observer(3)=Points_observer(i)%global_coordmesh(3)%coordmesh_z
hz_observer(3)=(EX(x_pos_observer(3),y_pos_observer(3)+1,z_pos_observer(3))-EX(x_pos_observer(3),y_pos_observer(3),z_pos_observer(3)))/CDELY(y_pos_observer(3))-&
(EY(x_pos_observer(3)+1,y_pos_observer(3),z_pos_observer(3))-EY(x_pos_observer(3),y_pos_observer(3),z_pos_observer(3)))/CDELX(x_pos_observer(3))
!=================================point4=======================================
x_pos_observer(4)=Points_observer(i)%global_coordmesh(4)%coordmesh_x
y_pos_observer(4)=Points_observer(i)%global_coordmesh(4)%coordmesh_y
z_pos_observer(4)=Points_observer(i)%global_coordmesh(4)%coordmesh_z
hz_observer(4)=(EX(x_pos_observer(4),y_pos_observer(4)+1,z_pos_observer(4))-EX(x_pos_observer(4),y_pos_observer(4),z_pos_observer(4)))/CDELY(y_pos_observer(4))-&
(EY(x_pos_observer(4)+1,y_pos_observer(4),z_pos_observer(4))-EY(x_pos_observer(4),y_pos_observer(4),z_pos_observer(4)))/CDELX(x_pos_observer(4))
!=================================point5=======================================
x_pos_observer(5)=Points_observer(i)%global_coordmesh(5)%coordmesh_x
y_pos_observer(5)=Points_observer(i)%global_coordmesh(5)%coordmesh_y
z_pos_observer(5)=Points_observer(i)%global_coordmesh(5)%coordmesh_z
hz_observer(5)=(EX(x_pos_observer(5),y_pos_observer(5)+1,z_pos_observer(5))-EX(x_pos_observer(5),y_pos_observer(5),z_pos_observer(5)))/CDELY(y_pos_observer(5))-&
(EY(x_pos_observer(5)+1,y_pos_observer(5),z_pos_observer(5))-EY(x_pos_observer(5),y_pos_observer(5),z_pos_observer(5)))/CDELX(x_pos_observer(5))
!=================================point6=======================================
x_pos_observer(6)=Points_observer(i)%global_coordmesh(6)%coordmesh_x
y_pos_observer(6)=Points_observer(i)%global_coordmesh(6)%coordmesh_y
z_pos_observer(6)=Points_observer(i)%global_coordmesh(6)%coordmesh_z
hz_observer(6)=(EX(x_pos_observer(6),y_pos_observer(6)+1,z_pos_observer(6))-EX(x_pos_observer(6),y_pos_observer(6),z_pos_observer(6)))/CDELY(y_pos_observer(6))-&
(EY(x_pos_observer(6)+1,y_pos_observer(6),z_pos_observer(6))-EY(x_pos_observer(6),y_pos_observer(6),z_pos_observer(6)))/CDELX(x_pos_observer(6))
!=================================point7=======================================
x_pos_observer(7)=Points_observer(i)%global_coordmesh(7)%coordmesh_x
y_pos_observer(7)=Points_observer(i)%global_coordmesh(7)%coordmesh_y
z_pos_observer(7)=Points_observer(i)%global_coordmesh(7)%coordmesh_z
hz_observer(7)=(EX(x_pos_observer(7),y_pos_observer(7)+1,z_pos_observer(7))-EX(x_pos_observer(7),y_pos_observer(7),z_pos_observer(7)))/CDELY(y_pos_observer(7))-&
(EY(x_pos_observer(7)+1,y_pos_observer(7),z_pos_observer(7))-EY(x_pos_observer(7),y_pos_observer(7),z_pos_observer(7)))/CDELX(x_pos_observer(7))
!=================================point8=======================================
x_pos_observer(8)=Points_observer(i)%global_coordmesh(8)%coordmesh_x
y_pos_observer(8)=Points_observer(i)%global_coordmesh(8)%coordmesh_y
z_pos_observer(8)=Points_observer(i)%global_coordmesh(8)%coordmesh_z
hz_observer(8)=(EX(x_pos_observer(8),y_pos_observer(8)+1,z_pos_observer(8))-EX(x_pos_observer(8),y_pos_observer(8),z_pos_observer(8)))/CDELY(y_pos_observer(8))-&
(EY(x_pos_observer(8)+1,y_pos_observer(8),z_pos_observer(8))-EY(x_pos_observer(8),y_pos_observer(8),z_pos_observer(8)))/CDELX(x_pos_observer(8))
data_rec(i) = hz_observer(1) * Points_observer(i)%coeff(1) + hz_observer(2) * Points_observer(i)%coeff(2) +&
hz_observer(3) * Points_observer(i)%coeff(3) + hz_observer(4) * Points_observer(i)%coeff(4) +&
hz_observer(5) * Points_observer(i)%coeff(5) + hz_observer(6) * Points_observer(i)%coeff(6) +&
hz_observer(7) * Points_observer(i)%coeff(7) + hz_observer(8) * Points_observer(i)%coeff(8)
ENDDO
enddo
deallocate(meps_r,mcq,msource,mdelt)
print*,mstop(num),'steps have just finished'
IF(Ctime(mstart(num)+mstop(num)-1)>(RAISETIME+WAVE+RAMP))THEN
DO i=1,point_num
idx_write=idx_start+i
WRITE(idx_write,*)mstart(num)+mstop(num)-1,Ctime(mstart(num)+mstop(num)-1)-(RAISETIME+WAVE+RAMP),data_rec(i)
ENDDO
write(*,'(a,i8,a,i8,a,f6.2,a)') 'Progress: [', num, '/', num_fra_com, '] (',100.0*num/num_fra_com, '%)'
ENDIF
ENDDO
call cpu_time(t_end)
t_total=t_end-t_start
print*,'The computing time is:', t_total
end subroutine Iteration_cpml
!===============================================================================================!
! ITERATION (below): original Dirichlet boundary iteration (Logic_PML=0).
! ITERATION_CPML (above): CPML absorbing boundary iteration (Logic_PML=1).
! main.f90 dispatches to either one according to Logic_PML read from input.dat.
!===============================================================================================!
! --------------------------------Subroutine part---------------------------------------------!
subroutine Iteration
use constantparameters
USE CONSTANTPARAMETERS
USE ELECTROMAGNETIC_VARIABLES
USE RES_MODEL_PARAMETER
USE TIME_PARAMETER
USE PML_PARAMETER
implicit none
real::t1,t2,t,t_start,t_end,t_total !t1 denotes original cpu time at the beginning of each computation fraction, t2 denotes the end cpu time and t=t2-t1
REAL*8 CA,CB,DELX1,DELY1,DELZ1 !ca, cb, delx1, dely1, delz1 are all middle variables used in the computation of EM field
REAL*8 TEMP_SIG,temp_cacb,data_rec(point_num) !Temp_sig and temp_cacb are middle variables used in the computation of EM field
REAL*8 DELY2,DELZ2,delx2 !They are all middle variables as above ones.
integer num,i,j,k,ii,iii,jj,kk,idx_write,x_pos_observer(8),y_pos_observer(8),z_pos_observer(8) !num is the number of computation fraction
integer :: N_hight=0
real*8,allocatable::Meps_r(:),Mdelt(:),Msource(:),Mcq(:) !They are local substitution of eps_r, delt and cq
REAL*8 hz_observer(8)
CHARACTER*20::string,str_num
WRITE(*,*)'[Iteration] Boundary condition: Dirichlet (zero-field) boundary (field fixed to zero at the outer grid faces)'
WRITE(*,*)'[Iteration] Iteration starts .. .. .. ..'
!Create output files
idx_start=12000
do iii=1,point_num
idx_write=idx_start+iii
IF(iii<10) THEN
write(str_num,"(I1)")iii
ELSEif(iii<100) THEN
write(str_num,"(I2)")iii
ELSEif(iii<1000) THEN
write(str_num,"(I3)")iii
ENDIF
string='dBzdt'//"_"//trim(str_num)//'.txt'
open(idx_write,file=string)
write(idx_write,*)"point_"//trim(str_num)
write(idx_write,*)Points_observer(iii)%local_coord_to_source%coord_x,Points_observer(iii)%local_coord_to_source%coord_y,&
Points_observer(iii)%local_coord_to_source%coord_z
enddo
call cpu_time(t_start)
!OPEN(20250220,file='dBzdt.txt')
do num=1,num_fra_com,1 !The outer loop which begins from the first fraction ends at the last fraction
call cpu_time(t1) !Record the cpu time at the beginning of each computing fraction
allocate(mdelt(0:mstop(num)),meps_r(mstop(num)),mcq(mstop(num)),msource(mstop(num)))
! The memory of mdelt, meps_r, mcq and msource are allocated at the begining of fraction
do ii=mstart(num),mstart(num)+mstop(num)-1,1
mdelt(ii-mstart(num)+1)=delt(ii)
meps_r(ii-mstart(num)+1)=eps_r(ii)
mcq(ii-mstart(num)+1)=cq(ii)
msource(ii-mstart(num)+1)=source(ii) !Link the local value of mdelt, meps_r, mcq and msorce to the global value of delt, eps_r, cq and source array.
end do
print*,'Now computing fraction:',num
mdelt(0)=mdelt(1)
do loop=1,mstop(num),1
! --------------------------------update the value of Ex and Ey in source area---------------------------------------!
DO J=2,NYB-1
DO I=1,NX
K=NZS+1-N_hight
DELY1=(CDELY(J-1)+CDELY(J))/2.0D0
DELZ1=CDELZ(NZ/2+1)
CA=(2.0D0*Meps_r(loop)-Mdelt(LOOP-1)*CCSIGX(I,J,K))/(2.0*Meps_r(loop)+Mdelt(LOOP-1)*CCSIGX(I,J,K))
CB=(2.0D0*MDELT(LOOP-1))/(2.0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGX(I,J,K))
EX(I,J,K)=CA*EX(I,J,K)+CB*((HZ(I,J,K)-HZ(I,J-1,K))*den_ey(J)/DELY1&
&-(HY(I,J,K)-HY(I,J,K-1))*den_ez(K)/DELZ1)-cb*Msource(loop)*is_ex_in_source(i,j)
ENDDO
ENDDO
! end of updating Ex while k=Nzs+1
! update the value of Ey while k=Nzs+1
DO J=1,NY
DO I=2,NX
K=NZS+1-N_hight
DELX1=(CDELX(I-1)+CDELX(I))/2.0
DELZ1=CDELZ(NZ/2+1)
CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*CCSIGY(I,J,K))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGY(I,J,K))
CB=(2.0D0*MDELT(LOOP-1))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGY(I,J,K))
EY(I,J,K)=CA*EY(I,J,K)+CB*((HX(I,J,K)-HX(I,J,K-1))*den_ez(K)/DELZ1&
&-(HZ(I,J,K)-HZ(I-1,J,K))*den_ex(I)/DELX1)-cb*Msource(loop)*is_ey_in_source(i,j)
ENDDO
ENDDO
! end of uptating Ey while k=Nzs+1
! ---------------------------------------------------Ex Part-------------------------------------------------------------!
DO K=NZS+2-N_hight,NZ
DO J=2,NY
DO I=1,NX
DELY1=(CDELY(J-1)+CDELY(J))/2.0D0
DELZ1=(CDELZ(K-1)+CDELZ(K))/2.0D0
CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*CCSIGX(I,J,K))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGX(I,J,K))
CB=(2.0D0*MDELT(LOOP-1))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGX(I,J,K))
EX(I,J,K)=CA*EX(I,J,K)+CB*((HZ(I,J,K)-HZ(I,J-1,K))*den_ey(J)/DELY1&
&-(HY(I,J,K)-HY(I,J,K-1))*den_ez(K)/DELZ1)
ENDDO
ENDDO
ENDDO
DO K=2,NZS-N_hight
DO J=2,NY
DO I=1,NX
DELY1=(CDELY(J-1)+CDELY(J))/2.0D0
DELZ1=(CDELZ(K-1)+CDELZ(K))/2.0D0
CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*CCSIGX(I,J,K))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGX(I,J,K))
CB=(2.0D0*MDELT(LOOP-1))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGX(I,J,K))
EX(I,J,K)=CA*EX(I,J,K)+CB*((HZ(I,J,K)-HZ(I,J-1,K))*den_ey(J)/DELY1&
&-(HY(I,J,K)-HY(I,J,K-1))*den_ez(K)/DELZ1)
ENDDO
ENDDO
ENDDO
! ================end of updating Ex==================!
! -----------------------------------------update the value of Ey--------------------------------!
DO K=NZS+2-N_hight,NZ
DO J=1,NY
DO I=2,NX
DELX1=(CDELX(I-1)+CDELX(I))/2.0D0
DELZ1=(CDELZ(K-1)+CDELZ(K))/2.0D0
CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*CCSIGY(I,J,K))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGY(I,J,K))
CB=(2.0D0*MDELT(LOOP-1))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGY(I,J,K))
EY(I,J,K)=CA*EY(I,J,K)+CB*((HX(I,J,K)-HX(I,J,K-1))*den_ez(K)/DELZ1&
&-(HZ(I,J,K)-HZ(I-1,J,K))*den_ex(I)/DELX1)
ENDDO
ENDDO
ENDDO
DO K=2,NZS-N_hight
DO J=1,NY
DO I=2,NXB-1
DELX1=(CDELX(I-1)+CDELX(I))/2.0D0
DELZ1=(CDELZ(K-1)+CDELZ(K))/2.0D0
CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*CCSIGY(I,J,K))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGY(I,J,K))
CB=(2.0D0*MDELT(LOOP-1))/(2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGY(I,J,K))
EY(I,J,K)=CA*EY(I,J,K)+CB*((HX(I,J,K)-HX(I,J,K-1))*den_ez(K)/DELZ1&
&-(HZ(I,J,K)-HZ(I-1,J,K))*den_ex(I)/DELX1)
ENDDO
ENDDO
ENDDO
!===============end of updating Ey===================!
! -------------------------------------update the value of Ez--------------------------------------!
DO K=1,NZ
DO J=2,NYB-1
DO I=2,NXB-1
DELX1=(CDELX(I-1)+CDELX(I))/2.0D0
DELY1=(CDELY(J-1)+CDELY(J))/2.0D0
TEMP_CACB=2.0D0*Meps_r(loop)+MDELT(LOOP-1)*CCSIGZ(I,J,K)
CA=(2.0D0*Meps_r(loop)-MDELT(LOOP-1)*CCSIGZ(I,J,K))/TEMP_CACB
CB=(2.0D0*MDELT(LOOP-1))/TEMP_CACB
EZ(I,J,K)=CA*EZ(I,J,K)+CB*((HY(I,J,K)-HY(I-1,J,K))*den_ex(I)/DELX1&
&-(HX(I,J,K)-HX(I,J-1,K))*den_ey(J)/DELY1)
ENDDO
ENDDO
ENDDO
!===============end of updating Ez=========================!
! ------------------------------------update the value of Hx-----------------------------------------------!
DO K=1,NZ
DO J=1,NY
DO I=1,NXB
DELY2=CDELY(J)
DELZ2=CDELZ(K)
HX(I,J,K)=HX(I,J,K)-MCQ(LOOP)*((EZ(I,J+1,K)-EZ(I,J,K))*den_hy(J)/DELY2&
&-(EY(I,J,K+1)-EY(I,J,K))*den_hz(K)/DELZ2)
ENDDO
ENDDO
ENDDO
!================end of updating Hx=======================!
! -------------------------------------update the value of Hy---------------------------------------------!
DO K=1,NZ
DO J=1,NYB
DO I=1,NX
DELZ2=CDELZ(K)
DELX2=CDELX(I)
HY(I,J,K)=HY(I,J,K)-MCQ(LOOP)*((EX(I,J,K+1)-EX(I,J,K))*den_hz(K)/DELZ2&
&-(EZ(I+1,J,K)-EZ(I,J,K))*den_hx(I)/DELX2)
ENDDO
ENDDO
ENDDO
!===============end of updating Hy========================!
!-------------------------------------update the value of Hz----------------------------------------------!
DO J=1,NY
DO I=1,NX
DO K=NZ,NZS+1,-1
DELX2=CDELX(I)
DELY2=CDELY(J)
DELZ2=CDELZ(K)
HZ(I,J,K)=HZ(I,J,K+1)+DELZ2*((HX(I+1,J,K)-HX(I,J,K))*den_hx(I)/DELX2&
&+(HY(I,J+1,K)-HY(I,J,K))*den_hy(J)/DELY2)*inv_hz_den(K)
ENDDO
ENDDO
ENDDO
DO K=1,NZS-1
DO J=1,NY
DO I=1,NX
DELX2=CDELX(I)
DELY2=CDELY(J)
DELZ2=CDELZ(K)
HZ(I,J,K+1)=HZ(I,J,K)-DELZ2*((HX(I+1,J,K)-HX(I,J,K))*den_hx(I)/DELX2&
&+(HY(I,J+1,K)-HY(I,J,K))*den_hy(J)/DELY2)*inv_hz_den(K)
ENDDO
ENDDO
ENDDO
!===================end of updating Hz==========================!
DO i=1,point_num
!=================================point1=======================================
x_pos_observer(1)=Points_observer(i)%global_coordmesh(1)%coordmesh_x
y_pos_observer(1)=Points_observer(i)%global_coordmesh(1)%coordmesh_y
z_pos_observer(1)=Points_observer(i)%global_coordmesh(1)%coordmesh_z
hz_observer(1)=(EX(x_pos_observer(1),y_pos_observer(1)+1,z_pos_observer(1))-EX(x_pos_observer(1),y_pos_observer(1),z_pos_observer(1)))/CDELY(y_pos_observer(1))-&
(EY(x_pos_observer(1)+1,y_pos_observer(1),z_pos_observer(1))-EY(x_pos_observer(1),y_pos_observer(1),z_pos_observer(1)))/CDELX(x_pos_observer(1))
!=================================point2=======================================
x_pos_observer(2)=Points_observer(i)%global_coordmesh(2)%coordmesh_x
y_pos_observer(2)=Points_observer(i)%global_coordmesh(2)%coordmesh_y
z_pos_observer(2)=Points_observer(i)%global_coordmesh(2)%coordmesh_z
hz_observer(2)=(EX(x_pos_observer(2),y_pos_observer(2)+1,z_pos_observer(2))-EX(x_pos_observer(2),y_pos_observer(2),z_pos_observer(2)))/CDELY(y_pos_observer(2))-&
(EY(x_pos_observer(2)+1,y_pos_observer(2),z_pos_observer(2))-EY(x_pos_observer(2),y_pos_observer(2),z_pos_observer(2)))/CDELX(x_pos_observer(2))
!=================================point3=======================================
x_pos_observer(3)=Points_observer(i)%global_coordmesh(3)%coordmesh_x
y_pos_observer(3)=Points_observer(i)%global_coordmesh(3)%coordmesh_y
z_pos_observer(3)=Points_observer(i)%global_coordmesh(3)%coordmesh_z
hz_observer(3)=(EX(x_pos_observer(3),y_pos_observer(3)+1,z_pos_observer(3))-EX(x_pos_observer(3),y_pos_observer(3),z_pos_observer(3)))/CDELY(y_pos_observer(3))-&
(EY(x_pos_observer(3)+1,y_pos_observer(3),z_pos_observer(3))-EY(x_pos_observer(3),y_pos_observer(3),z_pos_observer(3)))/CDELX(x_pos_observer(3))
!=================================point4=======================================
x_pos_observer(4)=Points_observer(i)%global_coordmesh(4)%coordmesh_x
y_pos_observer(4)=Points_observer(i)%global_coordmesh(4)%coordmesh_y
z_pos_observer(4)=Points_observer(i)%global_coordmesh(4)%coordmesh_z
hz_observer(4)=(EX(x_pos_observer(4),y_pos_observer(4)+1,z_pos_observer(4))-EX(x_pos_observer(4),y_pos_observer(4),z_pos_observer(4)))/CDELY(y_pos_observer(4))-&
(EY(x_pos_observer(4)+1,y_pos_observer(4),z_pos_observer(4))-EY(x_pos_observer(4),y_pos_observer(4),z_pos_observer(4)))/CDELX(x_pos_observer(4))
!=================================point5=======================================
x_pos_observer(5)=Points_observer(i)%global_coordmesh(5)%coordmesh_x
y_pos_observer(5)=Points_observer(i)%global_coordmesh(5)%coordmesh_y
z_pos_observer(5)=Points_observer(i)%global_coordmesh(5)%coordmesh_z
hz_observer(5)=(EX(x_pos_observer(5),y_pos_observer(5)+1,z_pos_observer(5))-EX(x_pos_observer(5),y_pos_observer(5),z_pos_observer(5)))/CDELY(y_pos_observer(5))-&
(EY(x_pos_observer(5)+1,y_pos_observer(5),z_pos_observer(5))-EY(x_pos_observer(5),y_pos_observer(5),z_pos_observer(5)))/CDELX(x_pos_observer(5))
!=================================point6=======================================
x_pos_observer(6)=Points_observer(i)%global_coordmesh(6)%coordmesh_x
y_pos_observer(6)=Points_observer(i)%global_coordmesh(6)%coordmesh_y
z_pos_observer(6)=Points_observer(i)%global_coordmesh(6)%coordmesh_z
hz_observer(6)=(EX(x_pos_observer(6),y_pos_observer(6)+1,z_pos_observer(6))-EX(x_pos_observer(6),y_pos_observer(6),z_pos_observer(6)))/CDELY(y_pos_observer(6))-&
(EY(x_pos_observer(6)+1,y_pos_observer(6),z_pos_observer(6))-EY(x_pos_observer(6),y_pos_observer(6),z_pos_observer(6)))/CDELX(x_pos_observer(6))
!=================================point7=======================================
x_pos_observer(7)=Points_observer(i)%global_coordmesh(7)%coordmesh_x
y_pos_observer(7)=Points_observer(i)%global_coordmesh(7)%coordmesh_y
z_pos_observer(7)=Points_observer(i)%global_coordmesh(7)%coordmesh_z
hz_observer(7)=(EX(x_pos_observer(7),y_pos_observer(7)+1,z_pos_observer(7))-EX(x_pos_observer(7),y_pos_observer(7),z_pos_observer(7)))/CDELY(y_pos_observer(7))-&
(EY(x_pos_observer(7)+1,y_pos_observer(7),z_pos_observer(7))-EY(x_pos_observer(7),y_pos_observer(7),z_pos_observer(7)))/CDELX(x_pos_observer(7))
!=================================point8=======================================
x_pos_observer(8)=Points_observer(i)%global_coordmesh(8)%coordmesh_x
y_pos_observer(8)=Points_observer(i)%global_coordmesh(8)%coordmesh_y
z_pos_observer(8)=Points_observer(i)%global_coordmesh(8)%coordmesh_z
hz_observer(8)=(EX(x_pos_observer(8),y_pos_observer(8)+1,z_pos_observer(8))-EX(x_pos_observer(8),y_pos_observer(8),z_pos_observer(8)))/CDELY(y_pos_observer(8))-&
(EY(x_pos_observer(8)+1,y_pos_observer(8),z_pos_observer(8))-EY(x_pos_observer(8),y_pos_observer(8),z_pos_observer(8)))/CDELX(x_pos_observer(8))
data_rec(i) = hz_observer(1) * Points_observer(i)%coeff(1) + hz_observer(2) * Points_observer(i)%coeff(2) +&
hz_observer(3) * Points_observer(i)%coeff(3) + hz_observer(4) * Points_observer(i)%coeff(4) +&
hz_observer(5) * Points_observer(i)%coeff(5) + hz_observer(6) * Points_observer(i)%coeff(6) +&
hz_observer(7) * Points_observer(i)%coeff(7) + hz_observer(8) * Points_observer(i)%coeff(8)
ENDDO
enddo
deallocate(meps_r,mcq,msource,mdelt)
print*,mstop(num),'steps have just finished'
IF(Ctime(mstart(num)+mstop(num)-1)>(RAISETIME+WAVE+RAMP))THEN
DO i=1,point_num
idx_write=idx_start+i
WRITE(idx_write,*)mstart(num)+mstop(num)-1,Ctime(mstart(num)+mstop(num)-1)-(RAISETIME+WAVE+RAMP),data_rec(i)
ENDDO
write(*,'(a,i8,a,i8,a,f6.2,a)') 'Progress: [', num, '/', num_fra_com, '] (',100.0*num/num_fra_com, '%)'
ENDIF
ENDDO
call cpu_time(t_end)
t_total=t_end-t_start
print*,'The computing time is:', t_total
end subroutine Iteration
查看 Git Blame 复制永久链接