92 行
3.7 KiB
Fortran
92 行
3.7 KiB
Fortran
!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 GET_COORDINATES
|
|
!> @brief This subroutine calculates the coordinates of each grid node in the 3D domain,
|
|
!>including the Yee grid nodes and the source-centered coordinate system.
|
|
USE CONSTANTPARAMETERS
|
|
IMPLICIT NONE
|
|
|
|
INTEGER(KIND=4) :: ii,jj,kk
|
|
!================================================
|
|
Coordiz(NZS) = -Cdelz(NZS)/2
|
|
!>In the Coordiz coordinate system, the coordinates are on the grid edge
|
|
!! Coordinates are defined at cell centers relative to source location
|
|
do kk=NZS-1,1,-1
|
|
Coordiz(kk)=Coordiz(kk+1)-(Cdelz(kk+1)+Cdelz(kk))/2
|
|
end do
|
|
do kk=NZS+1,NZ,1
|
|
Coordiz(kk)=Coordiz(kk-1)+(Cdelz(kk-1)+Cdelz(kk))/2
|
|
end do
|
|
!> Initialize X and Y coordinates based on source grid length being odd or even
|
|
IF(Logi_Sourcelenth) THEN
|
|
!> For odd source grid length, origin is at the grid center
|
|
Coordix(NXS)=0
|
|
do ii=NXS-1,1,-1
|
|
Coordix(ii)=Coordix(ii+1)-(Cdelx(ii)+Cdelx(ii+1))/2
|
|
end do
|
|
do ii=NXS+1,NX,1
|
|
Coordix(ii)=Coordix(ii-1)+(Cdelx(ii-1)+Cdelx(ii))/2
|
|
end do
|
|
|
|
Coordiy(NYS)=0
|
|
do jj=NYS-1,1,-1
|
|
Coordiy(jj)=Coordiy(jj+1)-(Cdely(jj)+Cdely(jj+1))/2
|
|
end do
|
|
do jj=NYS+1,NY,1
|
|
Coordiy(jj)=Coordiy(jj-1)+(Cdely(jj)+Cdely(jj-1))/2
|
|
end do
|
|
ELSE
|
|
!> For even source grid length, origin spans between two grid centers
|
|
Coordix(NXS)=-GridSize/2.0
|
|
Coordix(NXS+1)=GridSize/2.0
|
|
do ii=NXS-1,1,-1
|
|
Coordix(ii)=Coordix(ii+1)-(Cdelx(ii)+Cdelx(ii+1))/2
|
|
end do
|
|
do ii=NXS+2,nx,1
|
|
Coordix(ii)=Coordix(ii-1)+(Cdelx(ii-1)+Cdelx(ii))/2
|
|
end do
|
|
|
|
Coordiy(NYS)=-GridSize/2.0
|
|
Coordiy(NYS+1)=GridSize/2.0
|
|
do jj=NYS-1,1,-1
|
|
Coordiy(jj)=Coordiy(jj+1)-(Cdely(jj)+Cdely(jj+1))/2
|
|
end do
|
|
do jj=NYS+2,ny,1
|
|
Coordiy(jj)=Coordiy(jj-1)+(Cdely(jj)+Cdely(jj-1))/2
|
|
end do
|
|
ENDIF
|
|
|
|
!===============================Calculate the Yee node coordinates========================================
|
|
!> Allocate and calculate Yee grid node coordinates
|
|
ALLOCATE(coordinates_x(NXB),coordinates_y(NYB),coordinates_z(NZB))
|
|
|
|
DO ii=1,NX
|
|
coordinates_x(ii)=Coordix(ii)-Cdelx(ii)/2.0
|
|
ENDDO
|
|
coordinates_x(NXB)=Coordix(NX)+Cdelx(NX)/2.0
|
|
|
|
DO jj=1,NY
|
|
coordinates_y(jj)=Coordiy(jj)-Cdely(jj)/2.0
|
|
ENDDO
|
|
coordinates_y(NYB)=Coordiy(NY)+Cdely(NY)/2.0
|
|
|
|
DO kk=1,NZ
|
|
coordinates_z(kk)=Coordiz(kk)-Cdelz(kk)/2
|
|
ENDDO
|
|
coordinates_z(NZB)=coordinates_z(NZ)+Cdelz(NZ)
|
|
|
|
!---------------------Create a global coordinate system about HZ-----------------------!
|
|
! The HZ grid planes are located at the Yee nodes, whose source-centered
|
|
! coordinates are stored in coordinates_x/y/z (origin at the loop source /
|
|
! ground surface). The receiver coordinates read from input.dat are also
|
|
! source-centered, so Coord_HZ_* must use the SAME origin, otherwise the
|
|
! receiver-to-grid search in Get_Receiver_Gridlabel yields index 0 (or an
|
|
! uninitialized value) and the observer interpolation in Iteration.f90 reads
|
|
! out-of-bounds EX/EY/CDELX/CDELY entries -> NaN in the dBzdt output files.
|
|
ALLOCATE(Coord_HZ_X(Nx),Coord_HZ_Y(NY),Coord_HZ_Z(NZB))
|
|
Coord_HZ_X(1:NX) = coordinates_x(1:NX) !Record the HZ coordinate information in the x direction.
|
|
Coord_HZ_Y(1:NY) = coordinates_y(1:NY) !Record the HZ coordinate information in the y direction.
|
|
Coord_HZ_Z(1:NZB) = coordinates_z(1:NZB) !Record the HZ coordinate information in the z direction.
|
|
END SUBROUTINE |