!Copyright (c) 2022 by LEEE under guide of Huaifeng Sun(sunhuaifeng@email.sdu.edu.cn)
!written by Xinyu Li(202335098@mail.sdu.edu.cn) and Qi Zhao(zhaoqi_326326@163.com)
SUBROUTINE terrain_conformal

 	USE CONSTANTPARAMETERS
	USE ELECTROMAGNETIC_VARIABLES
	USE RES_MODEL_PARAMETER
	USE TIME_PARAMETER   
    USE OMP_LIB
    IMPLICIT NONE
                                                                                      
    INTEGER :: i,ii,j,jj,k,kk,i0,j1,i1,k2,i2,k3,j3,l1,l2,l3,iii,jjj,kkk
    INTEGER :: index, recorder
    logical :: JudgmentValue,LOGICAL_1,LOGICAL_2,LOGICAL_3,LOGICAL_4,LOGICAL_5,LOGICAL_6
    CHARACTER*255 :: temp, FineNameOfTerrain
    REAL(KIND=8), DIMENSION(:), ALLOCATABLE :: Normal_temp
    CHARACTER(300) :: line	!used to read one line of the ASCII STL file
    INTEGER(KIND=4) :: ios_stl,i_face,iv_face,nv_tmp,idx_v,ip_stl
    INTEGER(KIND=4), DIMENSION(:,:), ALLOCATABLE :: tmp_face
    REAL(KIND=8), DIMENSION(:,:), ALLOCATABLE :: tmp_vert,tmp_norm
    REAL(KIND=8) :: vx,vy,vz
    REAL(KIND=8) :: ex1x,ex1y,ex1z,ex2x,ex2y,ex2z,crossx,crossy,crossz

    !>The terrain mesh file can be Complex_Terrain.dat (the original text format) or
    !!Complex_Terrain.stl (the ASCII STL format). Which one is used was decided in
    !!GETDATA, and here only the corresponding reading branch is entered. The variables
    !!filled below (Node_Label, CoordinatesX/Y/Z, Element_Label, Element_Node1/2/3,
    !!n_point, n_face) keep the same names in both formats.
    IF(Logic_TerrainDat)THEN
        !>Original .dat format:
        !!  Line 1 : "Number of Nodes and Elements:"
        !!  Line 2 : n_point (number of nodes)
        !!  Line 3 : n_face  (number of triangular elements)
        !!  Line 4 : "Nodes Coordinates:"
        !!  next n_point lines : label, Coord_X, Coord_Y, Coord_Z
        !!  then 2 title lines, then n_face lines : label, node1, node2, node3
        FineNameOfTerrain = "Complex_Terrain.dat"
        OPEN( 520,  FILE = FinenameOfTerrain )
        Read(520, *) temp
        Read(520, *) n_point !get total number of  Node
        Read(520, *) n_face !get total number of  Element
        Read(520,*) temp

        ALLOCATE( Normal(n_face,3) )
        ALLOCATE( Normal_temp(n_face) )
        ALLOCATE( Node_Label(n_point), CoordinatesX(n_point), CoordinatesY(n_point), CoordinatesZ(n_point) )
        ALLOCATE( Element_Label(n_face), Element_Node1(n_face),Element_Node2(n_face), Element_Node3(n_face) )
		ALLOCATE( vert0(3,n_face), vert1(3,n_face), vert2(3,n_face), edge1(3,n_face), edge2(3,n_face))
        ALLOCATE( Face_Triangle_NormVect(3,n_face))

        DO i =1, n_point
            Read(520,*) Node_Label(i), CoordinatesX(i), CoordinatesY(i), CoordinatesZ(i) ! get the label of Node; get the coordinates of X, Y, Z connecting the Node
        ENDDO
        Read(520,*) temp
        Read(520,*) temp
        DO i =1, n_face
            Read(520,*) Element_Label(i), Element_Node1(i),Element_Node2(i), Element_Node3(i) ! get the label of Element; get the Label of Node connecting the Node
            vert0(1,i)=CoordinatesX( Element_Node1(i) )
			vert0(2,i)=CoordinatesY( Element_Node1(i) )
			vert0(3,i)=CoordinatesZ( Element_Node1(i) )
			vert1(1,i)=CoordinatesX( Element_Node2(i) )
			vert1(2,i)=CoordinatesY( Element_Node2(i) )
			vert1(3,i)=CoordinatesZ( Element_Node2(i) )
			vert2(1,i)=CoordinatesX( Element_Node3(i) )
			vert2(2,i)=CoordinatesY( Element_Node3(i) )
			vert2(3,i)=CoordinatesZ( Element_Node3(i) )
        ENDDO
        CLOSE(520)
    ELSEIF(Logic_TerrainStl)THEN
        !>ASCII STL format:
        !!  solid <name>
        !!  facet normal nx ny nz
        !!    outer loop
        !!      vertex x y z
        !!      vertex x y z
        !!      vertex x y z
        !!    endloop
        !!  endfacet
        !!  ...
        !!  endsolid <name>
        !!In an STL file the vertices are written once per facet, so the duplicated
        !!vertices are merged into unique nodes before filling the global arrays.
        !!The vertex order of each facet is also checked against the facet normal so
        !!that the normal direction convention is the same as the .dat format.
        OPEN(520,FILE='Complex_Terrain.stl',STATUS='OLD')
        !>First pass: count the number of facets.
        n_face=0
        DO
            READ(520,'(A)',IOSTAT=ios_stl) line
            IF(ios_stl/=0) EXIT
            !>convert the line into lower case for keyword matching
            DO ii=1,LEN_TRIM(line)
                IF(line(ii:ii)>='A'.AND.line(ii:ii)<='Z') line(ii:ii)=ACHAR(IACHAR(line(ii:ii))+32)
            ENDDO
            IF(INDEX(line,'facet')>0 .AND. INDEX(line,'endfacet')==0) n_face=n_face+1
        ENDDO
        IF(n_face==0)THEN
            WRITE(*,*)'Error: no facet is found in Complex_Terrain.stl!'
            STOP
        ENDIF
        REWIND(520)
        ALLOCATE(tmp_face(3,n_face),tmp_norm(3,n_face),tmp_vert(3,3*n_face))
        nv_tmp=0
        i_face=0
        iv_face=0
        DO
            READ(520,'(A)',IOSTAT=ios_stl) line
            IF(ios_stl/=0) EXIT
            DO ii=1,LEN_TRIM(line)
                IF(line(ii:ii)>='A'.AND.line(ii:ii)<='Z') line(ii:ii)=ACHAR(IACHAR(line(ii:ii))+32)
            ENDDO
            IF(INDEX(line,'facet')>0 .AND. INDEX(line,'endfacet')==0)THEN
                !>a new facet begins
                i_face=i_face+1
                iv_face=0
                IF(INDEX(line,'normal')>0)THEN
                    READ(line(INDEX(line,'normal')+6:),*,IOSTAT=ios_stl) tmp_norm(1,i_face),tmp_norm(2,i_face),tmp_norm(3,i_face)
                    IF(ios_stl/=0)THEN
                        WRITE(*,*)'Error: failed to read the facet normal line in Complex_Terrain.stl!'
                        STOP
                    ENDIF
                ENDIF
            ELSEIF(INDEX(line,'endfacet')>0)THEN
                !>a facet is finished, check that it has exactly 3 vertices
                IF(iv_face/=3)THEN
                    WRITE(*,*)'Error: a facet with',iv_face,'vertices (instead of 3) is found in Complex_Terrain.stl!'
                    STOP
                ENDIF
            ELSEIF(INDEX(line,'vertex')>0)THEN
                iv_face=iv_face+1
                IF(iv_face>3)THEN
                    WRITE(*,*)'Error: a facet with more than 3 vertices is found in Complex_Terrain.stl!'
                    STOP
                ENDIF
                READ(line(INDEX(line,'vertex')+6:),*,IOSTAT=ios_stl) vx,vy,vz
                IF(ios_stl/=0)THEN
                    WRITE(*,*)'Error: failed to read a vertex line in Complex_Terrain.stl!'
                    STOP
                ENDIF
                !>merge the duplicated vertices
                idx_v=0
                DO ip_stl=1,nv_tmp
                    IF(ABS(tmp_vert(1,ip_stl)-vx)<eps105.AND.ABS(tmp_vert(2,ip_stl)-vy)<eps105.AND.ABS(tmp_vert(3,ip_stl)-vz)<eps105)THEN
                        idx_v=ip_stl
                        EXIT
                    ENDIF
                ENDDO
                IF(idx_v==0)THEN
                    nv_tmp=nv_tmp+1
                    tmp_vert(1,nv_tmp)=vx
                    tmp_vert(2,nv_tmp)=vy
                    tmp_vert(3,nv_tmp)=vz
                    idx_v=nv_tmp
                ENDIF
                tmp_face(iv_face,i_face)=idx_v
            ENDIF
        ENDDO
        !>Correct the vertex order of each facet: compare the cross-product normal
        !!with the facet normal stored in the STL file, swap node2/node3 if they
        !!point in opposite directions, so that the normal direction convention
        !!is the same as in the .dat format.
        DO i=1,n_face
            ex1x=tmp_vert(1,tmp_face(2,i))-tmp_vert(1,tmp_face(1,i))
            ex1y=tmp_vert(2,tmp_face(2,i))-tmp_vert(2,tmp_face(1,i))
            ex1z=tmp_vert(3,tmp_face(2,i))-tmp_vert(3,tmp_face(1,i))
            ex2x=tmp_vert(1,tmp_face(3,i))-tmp_vert(1,tmp_face(1,i))
            ex2y=tmp_vert(2,tmp_face(3,i))-tmp_vert(2,tmp_face(1,i))
            ex2z=tmp_vert(3,tmp_face(3,i))-tmp_vert(3,tmp_face(1,i))
            crossx=ex1y*ex2z-ex1z*ex2y
            crossy=ex1z*ex2x-ex1x*ex2z
            crossz=ex1x*ex2y-ex1y*ex2x
            IF(crossx*tmp_norm(1,i)+crossy*tmp_norm(2,i)+crossz*tmp_norm(3,i)<0.0D0)THEN
                idx_v=tmp_face(2,i)
                tmp_face(2,i)=tmp_face(3,i)
                tmp_face(3,i)=idx_v
            ENDIF
        ENDDO
        n_point=nv_tmp
        ALLOCATE( Normal(n_face,3) )
        ALLOCATE( Normal_temp(n_face) )
        ALLOCATE( Node_Label(n_point), CoordinatesX(n_point), CoordinatesY(n_point), CoordinatesZ(n_point) )
        ALLOCATE( Element_Label(n_face), Element_Node1(n_face),Element_Node2(n_face), Element_Node3(n_face) )
        ALLOCATE( vert0(3,n_face), vert1(3,n_face), vert2(3,n_face), edge1(3,n_face), edge2(3,n_face))
        ALLOCATE( Face_Triangle_NormVect(3,n_face))
        DO j=1,n_point
            Node_Label(j)=j
            CoordinatesX(j)=tmp_vert(1,j)
            CoordinatesY(j)=tmp_vert(2,j)
            CoordinatesZ(j)=tmp_vert(3,j)
        ENDDO
        DO i=1,n_face
            Element_Label(i)=i
            Element_Node1(i)=tmp_face(1,i)
            Element_Node2(i)=tmp_face(2,i)
            Element_Node3(i)=tmp_face(3,i)
            vert0(1,i)=CoordinatesX( Element_Node1(i) )
            vert0(2,i)=CoordinatesY( Element_Node1(i) )
            vert0(3,i)=CoordinatesZ( Element_Node1(i) )
            vert1(1,i)=CoordinatesX( Element_Node2(i) )
            vert1(2,i)=CoordinatesY( Element_Node2(i) )
            vert1(3,i)=CoordinatesZ( Element_Node2(i) )
            vert2(1,i)=CoordinatesX( Element_Node3(i) )
            vert2(2,i)=CoordinatesY( Element_Node3(i) )
            vert2(3,i)=CoordinatesZ( Element_Node3(i) )
        ENDDO
        CLOSE(520)
        DEALLOCATE(tmp_face,tmp_norm,tmp_vert)
        WRITE(*,*)'Complex_Terrain.stl read: n_point=',n_point,' n_face=',n_face
    ELSE
        WRITE(*,*)'Error: neither Complex_Terrain.dat nor Complex_Terrain.stl exists, terrain_conformal can not run!'
        RETURN
    ENDIF
    !============================================================================================================================================================
    !Find out the maximum and minimum values of the abnormal volume triangular mesh in z-axis.
    max_coord_z=maxval(CoordinatesZ)
    min_coord_z=minval(CoordinatesZ)
    DO kk=1,NZ
        IF(coordinates_z(kk)>min_coord_z)THEN
            Z_min=kk-1
            EXIT
        ENDIF
    ENDDO
    DO kk=1,NZ
        IF(coordinates_z(kk)>max_coord_z)THEN
		    Z_max=kk
		    EXIT
		ENDIF
       ENDDO
	ALLOCATE(orig_z(3,NXB*NYB),orig_y(3,NXB*NZB),orig_x(3,NXB*NZB))
	ALLOCATE(det_z(NXB*NYB),det_x(NYB*NZB),det_y(NXB*NZB))
	ALLOCATE(u_z(NXB*NYB),u_x(NYB*NZB),u_y(NXB*NZB))
	ALLOCATE(v_z(NXB*NYB),v_x(NYB*NZB),v_y(NXB*NZB))
	ALLOCATE(t_z(NXB*NYB),t_x(NYB*NZB),t_y(NXB*NZB))
	ALLOCATE(pvec_z(3,NXB*NYB),pvec_y(3,NXB*NZB),pvec_x(3,NYB*NZB))
	ALLOCATE(tvec_z(3,NXB*NYB),tvec_y(3,NXB*NZB),tvec_x(3,NYB*NZB))
	ALLOCATE(crosspoint_ZZ(50,NXB*NYB),crosspoint_YY(50,NXB*NZB),crosspoint_XX(50,NYB*NZB))
	ALLOCATE(mmz_per(NXB*NYB),mmy_per(NXB*NZB),mmx_per(NYB*NZB))
    
    !Vector of  the Triangle
	edge1 = vert1 - vert0
	edge2 = vert2 - vert0
	dir_z   = [0.D0,0.D0,1.D0]
	dir_y   = [0.D0,1.D0,0.D0]
	dir_x   = [1.D0,0.D0,0.D0]
	mmx_per=0
	mmy_per=0
	mmz_per=0
	DO i0=1,n_face
		Face_Triangle_NormVect(1,i0)=edge1(2,i0) * edge2(3,i0)-edge1(3,i0) * edge2(2,i0)
		Face_Triangle_NormVect(2,i0)=edge1(3,i0) * edge2(1,i0)-edge1(1,i0) * edge2(3,i0)
		Face_Triangle_NormVect(3,i0)=edge1(1,i0) * edge2(2,i0)-edge1(2,i0) * edge2(1,i0)
	ENDDO            
	!CALL OMP_SET_NUM_THREADS(16)
    !$OMP PARALLEL DO PRIVATE(i0,j1,i1,Rz,verts_Dotmultp,LOGICAL_1,LOGICAL_2,LOGICAL_3,LOGICAL_4,LOGICAL_5,LOGICAL_6)			
    !get the intersaction of ray and z-face and save it into "coor_z_terrain"
        Do j1=1,NYB
            Do i1=1,NXB
                Rz=(j1-1)*NXB+i1
                DO i0=1,n_face
					orig_z(1,Rz)=coordinates_x(i1)
					orig_z(2,Rz)=coordinates_y(j1)
					orig_z(3,Rz)=coordinates_z(1)
					tvec_z(1:3,Rz) = orig_z(1:3,Rz) - vert0(1:3,i0)
					pvec_z(1,Rz) = dir_z(2)*edge2(3,i0) - dir_z(3)*edge2(2,i0)
					pvec_z(2,Rz) = dir_z(3)*edge2(1,i0) - dir_z(1)*edge2(3,i0)
					pvec_z(3,Rz) = dir_z(1)*edge2(2,i0) - dir_z(2)*edge2(1,i0)
					det_z(Rz)=edge1(1,i0)*pvec_z(1,Rz)+edge1(2,i0)*pvec_z(2,Rz)+edge1(3,i0)*pvec_z(3,Rz)
					IF (abs(det_z(Rz)) < eps105) THEN
						CYCLE
					END IF
					u_z(Rz) = (tvec_z(1,Rz)*pvec_z(1,Rz)+tvec_z(2,Rz)*pvec_z(2,Rz)+tvec_z(3,Rz)*pvec_z(3,Rz))/det_z(Rz)
					IF (u_z(Rz) < 0.D0 .or. u_z(Rz) > 1.D0) THEN
						CYCLE
					END IF 
					pvec_z(1,Rz) = tvec_z(2,Rz)*edge1(3,i0) - tvec_z(3,Rz)*edge1(2,i0)
					pvec_z(2,Rz) = tvec_z(3,Rz)*edge1(1,i0) - tvec_z(1,Rz)*edge1(3,i0)
					pvec_z(3,Rz) = tvec_z(1,Rz)*edge1(2,i0) - tvec_z(2,Rz)*edge1(1,i0)
					v_z(Rz) = (dir_z(1)*pvec_z(1,Rz)+dir_z(2)*pvec_z(2,Rz)+dir_z(3)*pvec_z(3,Rz))/det_z(Rz)
					IF (v_z(Rz) <  0.D0 .or. u_z(Rz) + v_z(Rz) > 1.D0) THEN
						CYCLE
					END IF
					t_z(Rz)=(edge2(1,i0)*pvec_z(1,Rz)+edge2(2,i0)*pvec_z(2,Rz)+edge2(3,i0)*pvec_z(3,Rz))/det_z(Rz)
					mmz_per(Rz) = mmz_per(Rz) + 1 
					crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_X = orig_z(1,Rz) + t_z(Rz) * dir_z(1)
                    crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_Y = orig_z(2,Rz) + t_z(Rz) * dir_z(2)
                    crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_Z = orig_z(3,Rz) + t_z(Rz) * dir_z(3)
                    LOGICAL_1=crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_X<coordinates_x(1)
                    LOGICAL_2=crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_X>coordinates_x(NXB)
                    LOGICAL_3=crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_Y<coordinates_y(1)
                    LOGICAL_4=crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_Y>coordinates_y(NYB)
                    LOGICAL_5=crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_Z<coordinates_z(1)
                    LOGICAL_6=crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_Z>coordinates_z(NZB)
                    IF(LOGICAL_1 .OR. LOGICAL_2 .OR.LOGICAL_3 .OR.LOGICAL_4 .OR.LOGICAL_5 .OR.LOGICAL_6)THEN
                        mmz_per(Rz) = mmz_per(Rz) - 1 
						CYCLE
                    ENDIF
                    !*********************************Determine the intersection_Z attribute*******************************
                    verts_Dotmultp = dir_z(1) * Face_Triangle_NormVect(1,i0) + dir_z(2) * Face_Triangle_NormVect(2,i0) + dir_z(3) * Face_Triangle_NormVect(3,i0)
                    IF(verts_Dotmultp > 0.D0)THEN
						crosspoint_ZZ(mmz_per(Rz),Rz)%Log_In=.TRUE.
                    ELSE
						crosspoint_ZZ(mmz_per(Rz),Rz)%Log_In=.FALSE.
                    ENDIF
                ENDDO
                IF(mmz_per(Rz)>1)THEN
                DO kk=2,mmz_per(Rz)
                    DO kkk=1,kk-1
                        IF(crosspoint_ZZ(kk,Rz)%Global_Coord%Coord_X<crosspoint_ZZ(kkk,Rz)%Global_Coord%Coord_X)THEN
                            CALL SWAP(crosspoint_ZZ(kk,Rz), crosspoint_ZZ(kkk,Rz))
                        ENDIF
                    ENDDO
                ENDDO
            ENDIF
            ENDDO
        ENDDO
        !$OMP END PARALLEL DO
        !$OMP PARALLEL DO PRIVATE(i0,k2,i2,Ry,verts_Dotmultp,LOGICAL_1,LOGICAL_2,LOGICAL_3,LOGICAL_4,LOGICAL_5,LOGICAL_6)
        !get the intersaction of ray and y-face and save it into "coor_y_terrain"	
        DO k2=1,NZB
            DO i2=1,NXB
                Ry=(k2-1)*NXB+i2
                DO i0=1,n_face
					orig_y(1,Ry)=coordinates_x(i2)
					orig_y(2,Ry)=coordinates_y(1)
					orig_y(3,Ry)=coordinates_z(k2)
					tvec_y(1:3,Ry) = orig_y(1:3,Ry) - vert0(1:3,i0)
					pvec_y(1,Ry) = dir_y(2)*edge2(3,i0) - dir_y(3)*edge2(2,i0)
					pvec_y(2,Ry) = dir_y(3)*edge2(1,i0) - dir_y(1)*edge2(3,i0)
					pvec_y(3,Ry) = dir_y(1)*edge2(2,i0) - dir_y(2)*edge2(1,i0)
					det_y(Ry)=edge1(1,i0)*pvec_y(1,Ry)+edge1(2,i0)*pvec_y(2,Ry)+edge1(3,i0)*pvec_y(3,Ry)
					IF (abs(det_y(Ry)) < eps105) THEN
						CYCLE
					END IF
					u_y(Ry) = (tvec_y(1,Ry)*pvec_y(1,Ry)+tvec_y(2,Ry)*pvec_y(2,Ry)+tvec_y(3,Ry)*pvec_y(3,Ry))/det_y(Ry)
					IF (u_y(Ry) < 0.0 .or. u_y(Ry) > 1.0) THEN
						CYCLE
					END IF
					pvec_y(1,Ry) = tvec_y(2,Ry)*edge1(3,i0) - tvec_y(3,Ry)*edge1(2,i0)
					pvec_y(2,Ry) = tvec_y(3,Ry)*edge1(1,i0) - tvec_y(1,Ry)*edge1(3,i0)
					pvec_y(3,Ry) = tvec_y(1,Ry)*edge1(2,i0) - tvec_y(2,Ry)*edge1(1,i0)
					v_y(Ry) = (dir_y(1)*pvec_y(1,Ry)+dir_y(2)*pvec_y(2,Ry)+dir_y(3)*pvec_y(3,Ry))/det_y(Ry)
					IF (v_y(Ry) < 0.0 .or. u_y(Ry) + v_y(Ry) > 1.0) THEN
						CYCLE
					END IF
					t_y(Ry)=(edge2(1,i0)*pvec_y(1,Ry)+edge2(2,i0)*pvec_y(2,Ry)+edge2(3,i0)*pvec_y(3,Ry))/det_y(Ry)
					mmy_per(Ry) = mmy_per(Ry) + 1 
					crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_X=orig_y(1,Ry) + t_y(Ry) * dir_y(1)
                    crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y=orig_y(2,Ry) + t_y(Ry) * dir_y(2)
                    crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Z=orig_y(3,Ry) + t_y(Ry) * dir_y(3)
                    LOGICAL_1=crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_X<coordinates_x(1)
                    LOGICAL_2=crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_X>coordinates_x(NXB)
                    LOGICAL_3=crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y<coordinates_y(1)
                    LOGICAL_4=crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y>coordinates_y(NYB)
                    LOGICAL_5=crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Z<coordinates_z(1)
                    LOGICAL_6=crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Z>coordinates_z(NZB)
                    IF(LOGICAL_1 .OR. LOGICAL_2 .OR.LOGICAL_3 .OR.LOGICAL_4 .OR.LOGICAL_5 .OR.LOGICAL_6)THEN
                        mmy_per(Ry) = mmy_per(Ry) - 1 
						CYCLE
                    ENDIF
                    !*********************************Determine the intersection_Y attribute*******************************
                    verts_Dotmultp = dir_y(1) * Face_Triangle_NormVect(1,i0) + dir_y(2) * Face_Triangle_NormVect(2,i0) + dir_y(3) * Face_Triangle_NormVect(3,i0)
                    if(verts_Dotmultp > 0.D0)THEN
						crosspoint_YY(mmy_per(Ry),Ry)%Log_In=.TRUE.
                    ELSE
						crosspoint_YY(mmy_per(Ry),Ry)%Log_In=.FALSE.
                    ENDIF
                ENDDO
                IF(mmy_per(Ry)>1)THEN
                DO jj=2,mmy_per(Ry)
                    DO jjj=1,jj-1
                        IF(crosspoint_YY(jj,Ry)%Global_Coord%Coord_Y<crosspoint_YY(jjj,Ry)%Global_Coord%Coord_Y)THEN
                            CALL SWAP(crosspoint_YY(jj,Ry), crosspoint_YY(jjj,Ry))
                        ENDIF
                    ENDDO
                ENDDO
            ENDIF
            ENDDO
        ENDDO
        !$OMP END PARALLEL DO
        !$OMP PARALLEL DO PRIVATE(i0,Rx,verts_Dotmultp,LOGICAL_1,LOGICAL_2,LOGICAL_3,LOGICAL_4,LOGICAL_5,LOGICAL_6)
        !get the intersaction of ray and x-face and save it into "coor_x_terrain"
        Do k3=1,NZB
            Do j3=1,NYB
                Rx=(k3-1)*NYB+j3
                DO i0=1,n_face
					orig_x(1,Rx)=coordinates_x(1)
					orig_x(2,Rx)=coordinates_y(j3)
					orig_x(3,Rx)=coordinates_z(k3)
					tvec_x(1:3,Rx) = orig_x(1:3,Rx) - vert0(1:3,i0)
					pvec_x(1,Rx) = dir_x(2)*edge2(3,i0) - dir_x(3)*edge2(2,i0)
					pvec_x(2,Rx) = dir_x(3)*edge2(1,i0) - dir_x(1)*edge2(3,i0)
					pvec_x(3,Rx) = dir_x(1)*edge2(2,i0) - dir_x(2)*edge2(1,i0)
					det_x(Rx)=edge1(1,i0)*pvec_x(1,Rx)+edge1(2,i0)*pvec_x(2,Rx)+edge1(3,i0)*pvec_x(3,Rx)
					IF (abs(det_x(Rx)) < eps105) THEN
						CYCLE
					END IF
					u_x(Rx) = (tvec_x(1,Rx)*pvec_x(1,Rx)+tvec_x(2,Rx)*pvec_x(2,Rx)+tvec_x(3,Rx)*pvec_x(3,Rx))/det_x(Rx)
					IF (u_x(Rx) < 0.0 .or. u_x(Rx) > 1.0) THEN
						CYCLE
					END IF
					pvec_x(1,Rx) = tvec_x(2,Rx)*edge1(3,i0) - tvec_x(3,Rx)*edge1(2,i0)
					pvec_x(2,Rx) = tvec_x(3,Rx)*edge1(1,i0) - tvec_x(1,Rx)*edge1(3,i0)
					pvec_x(3,Rx) = tvec_x(1,Rx)*edge1(2,i0) - tvec_x(2,Rx)*edge1(1,i0)
					v_x(Rx) = (dir_x(1)*pvec_x(1,Rx)+dir_x(2)*pvec_x(2,Rx)+dir_x(3)*pvec_x(3,Rx))/det_x(Rx)
					IF (v_x(Rx) < 0.0 .or. u_x(Rx) + v_x(Rx) > 1.0) THEN
						CYCLE
					END IF
					t_x(Rx)=(edge2(1,i0)*pvec_x(1,Rx)+edge2(2,i0)*pvec_x(2,Rx)+edge2(3,i0)*pvec_x(3,Rx))/det_x(Rx)
					mmx_per(Rx) = mmx_per(Rx) + 1
					crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_X= orig_x(1,Rx) + t_x(Rx) * dir_x(1)
                    crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_Y= orig_x(2,Rx) + t_x(Rx) * dir_x(2)
                    crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_Z= orig_x(3,Rx) + t_x(Rx) * dir_x(3)
                    LOGICAL_1=crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_X<coordinates_x(1)
                    LOGICAL_2=crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_X>coordinates_x(NXB)
                    LOGICAL_3=crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_Y<coordinates_y(1)
                    LOGICAL_4=crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_Y>coordinates_y(NYB)
                    LOGICAL_5=crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_Z<coordinates_z(1)
                    LOGICAL_6=crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_Z>coordinates_z(NZB)
                    IF(LOGICAL_1 .OR. LOGICAL_2 .OR.LOGICAL_3 .OR.LOGICAL_4 .OR.LOGICAL_5 .OR.LOGICAL_6)THEN
                        mmx_per(Rx) = mmx_per(Rx) - 1 
						CYCLE
                    ENDIF
                    !*********************************Determine the intersection_X attribute*******************************
                    verts_Dotmultp = dir_x(1) * Face_Triangle_NormVect(1,i0) + dir_x(2) * Face_Triangle_NormVect(2,i0) + dir_x(3) * Face_Triangle_NormVect(3,i0)
                    if(verts_Dotmultp > 0.D0)THEN
						crosspoint_XX(mmx_per(Rx),Rx)%Log_In=.TRUE.
                    ELSE
						crosspoint_XX(mmx_per(Rx),Rx)%Log_In=.FALSE.
                    ENDIF
                ENDDO
                IF(mmx_per(Rx)>1)THEN
					DO ii=2,mmx_per(Rx)
						DO iii=1,mmx_per(Rx)-1
							IF(crosspoint_XX(ii,Rx)%Global_Coord%Coord_X<crosspoint_XX(iii,Rx)%Global_Coord%Coord_X)THEN
								CALL SWAP(crosspoint_XX(ii,Rx), crosspoint_XX(iii,Rx))
							ENDIF
						ENDDO
					ENDDO
				ENDIF
            ENDDO
		ENDDO
		!$OMP END PARALLEL DO
    PRINT*,'Ray tracing computation of terrain is complete!'
    !=========================================================================================================	
    WRITE(5141,*)"========================== This is the SUBROUTINE CONFORMALGRID =========================="
    !=========================================Terrain conformal in the x-direction==========================================
	!>The x-conductivity is being treated
    DO k=1, Z_min-1
		DO j=1, NYB
			DO i=1,NX
				CCSIGX(I,J,K)=AIR_CONDUCTIVITY
			ENDDO
		ENDDO
	ENDDO
	DO k=Z_max+1, NZB
		DO j=1, NYB
			DO i=1,NX
				CCSIGX(I,J,K)=TAR_CONDUCTIVITY(2)
			ENDDO
		ENDDO
    ENDDO
    DO k=Z_min, Z_max
        DO j=1, NYB
            Rx=(k-1)*NYB+j
            IF(mmx_per(Rx)==1)THEN  !The case with only one intersection point
                IF (crosspoint_XX(1,Rx)%Log_In) THEN !The intersection point is the penetration point, that is, the air penetrates into the stratum
                    DO i=1,NX
                        Logic_1=(crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_X > coordinates_x(i)).AND. &
                                  & (crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_X < coordinates_x(i+1))
                        IF(Logic_1) THEN
                            KIdx_1=i
                            EXIT
                        ENDIF
                    ENDDO
                    LenRatio_CCSIGX(KIdx_1,J,K) = ABS((crosspoint_XX(1,Rx)%Global_Coord%Coord_X-coordinates_x(KIdx_1)))/Cdelx(KIdx_1)
                    CCSIGX(1:(KIdx_1-1),J,K) = AIR_CONDUCTIVITY
                    CCSIGX(KIdx_1,J,K) = LenRatio_CCSIGX(KIdx_1,J,K)*AIR_CONDUCTIVITY+(1-LenRatio_CCSIGX(KIdx_1,J,K))*TAR_CONDUCTIVITY(2)
                    CCSIGX((KIdx_1+1):NX,J,K) = TAR_CONDUCTIVITY(2)
                ELSE  !The intersection point is the exit point, that is, it penetrates into the air from the stratum
                    DO i=1,NX
                        Logic_1=(crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_X > coordinates_x(i)).AND. &
                                  & (crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_X < coordinates_x(i+1))
                        IF(Logic_1) THEN
                            KIdx_1=i
                            EXIT
                        ENDIF
                    ENDDO
                    LenRatio_CCSIGX(KIdx_1,J,K) = ABS((crosspoint_XX(1,Rx)%Global_Coord%Coord_X-coordinates_x(KIdx_1)))/Cdelx(KIdx_1)
                    CCSIGX(1:(KIdx_1-1),J,K) = TAR_CONDUCTIVITY(2)
                    CCSIGX(KIdx_1,J,K) = LenRatio_CCSIGX(KIdx_1,J,K)*TAR_CONDUCTIVITY(2)+(1-LenRatio_CCSIGX(KIdx_1,J,K))*AIR_CONDUCTIVITY
                    CCSIGX((KIdx_1+1):NX,J,K) = AIR_CONDUCTIVITY
                ENDIF
            ELSEIF(mmx_per(Rx)>1)THEN  !There are intersection points and the number is greater than one
                idx_start=0;idx_end=0
                DO ii=1,mmx_per(Rx)
                    IF (ii==1)THEN      !First, determine the first intersection point
                        IF(crosspoint_XX(ii,Rx)%Log_In)THEN
                            idx_start=ii
                        ELSEIF(.NOT.crosspoint_XX(ii,Rx)%Log_In)THEN
                            idx_end=ii
                        ENDIF
                        IF (crosspoint_XX(ii,Rx)%Log_In) THEN !The intersection point is the penetration point, that is, the air penetrates into the stratum
                            DO i=1,NX
                                Logic_1=(crosspoint_XX(ii,Rx)%Global_Coord%Coord_X > coordinates_x(i)).AND. &
                                      & (crosspoint_XX(ii,Rx)%Global_Coord%Coord_X < coordinates_x(i+1))
                                IF(Logic_1) THEN
                                    KIdx_1=i
                                    EXIT
                                ENDIF
                            ENDDO
                            CCSIGX(1:(KIdx_1-1),J,K) = AIR_CONDUCTIVITY
                        ELSE  !The intersection point is the exit point, that is, it penetrates into the air from the stratum
                            DO i=1,NX
                                Logic_1=(crosspoint_XX(ii,Rx)%Global_Coord%Coord_X > coordinates_x(i)).AND. &
                                      & (crosspoint_XX(ii,Rx)%Global_Coord%Coord_X < coordinates_x(i+1))
                                IF(Logic_1) THEN
                                    KIdx_1=i
                                    EXIT
                                ENDIF
                            ENDDO
                            CCSIGX(1:(KIdx_1-1),J,K) = TAR_CONDUCTIVITY(2)
                        ENDIF
                    ENDIF
                    
                    IF(crosspoint_XX(ii,Rx)%Log_In)THEN
                        idx_start=ii
                    ELSEIF(.NOT.crosspoint_XX(ii,Rx)%Log_In)THEN
                        idx_end=ii
                    ENDIF
                    IF((idx_start>0).AND.(idx_end>0))THEN 
                        DO i=1,NX
                            Logic_1=(crosspoint_XX(idx_start,Rx)%Global_Coord%Coord_X > coordinates_x(i)).AND. &
                                  & (crosspoint_XX(idx_start,Rx)%Global_Coord%Coord_X < coordinates_x(i+1))
                            Logic_2=(crosspoint_XX(idx_end,Rx)%Global_Coord%Coord_X > coordinates_x(i)).AND. &
                                  & (crosspoint_XX(idx_end,Rx)%Global_Coord%Coord_X < coordinates_x(i+1))
                            IF(Logic_1) KIdx_1=i
                            IF(Logic_2) KIdx_2=i
                            IF(Logic_1.and.Logic_2)	EXIT
                        ENDDO
                        IF(idx_start<idx_end) THEN           !First, break through the air and enter the stratum
                            IF(KIdx_1==KIdx_2) THEN
                                LenRatio_CCSIGX(KIdx_1,J,K) = ABS((crosspoint_XX(idx_end,Rx)%Global_Coord%Coord_X-crosspoint_XX(idx_start,Rx)%Global_Coord%Coord_X))/Cdelx(KIdx_1)
                                CCSIGX(KIdx_1,J,K) =AIR_CONDUCTIVITY * (1-LenRatio_CCSIGX(KIdx_1,J,K)) + TAR_CONDUCTIVITY(2) * LenRatio_CCSIGX(KIdx_1,J,K)
                            ELSEIF(KIdx_2 > KIdx_1)THEN      
                                LenRatio_CCSIGX(KIdx_1,J,K) = ABS((crosspoint_XX(idx_start,Rx)%Global_Coord%Coord_X-coordinates_x(KIdx_1)))/Cdelx(KIdx_1)
                                LenRatio_CCSIGX(KIdx_2,J,K) = ABS((crosspoint_XX(idx_end,Rx)%Global_Coord%Coord_X-coordinates_x(KIdx_2)))/Cdelx(KIdx_2)
                                CCSIGX(KIdx_1,J,K) = AIR_CONDUCTIVITY * LenRatio_CCSIGX(KIdx_1,J,K) + TAR_CONDUCTIVITY(2) * (1-LenRatio_CCSIGX(KIdx_1,J,K))
                                CCSIGX((KIdx_1+1):(KIdx_2-1),J,K) = TAR_CONDUCTIVITY(2)
                                CCSIGX(KIdx_2,J,K) = TAR_CONDUCTIVITY(2) * LenRatio_CCSIGX(KIdx_2,J,K) + AIR_CONDUCTIVITY * (1-LenRatio_CCSIGX(KIdx_2,J,K))
                            ENDIF
                        ELSE                                 !First, break through the stratum and enter the air
                            IF(KIdx_1==KIdx_2) THEN
                                LenRatio_CCSIGX(KIdx_1,J,K) = ABS((crosspoint_XX(idx_start,Rx)%Global_Coord%Coord_X-crosspoint_XX(idx_end,Rx)%Global_Coord%Coord_X))/Cdelx(KIdx_1)
                                CCSIGX(KIdx_1,J,K) = CCSIGX(KIdx_1,J,K)+TAR_CONDUCTIVITY(2) * LenRatio_CCSIGX(KIdx_1,J,K)
                            ELSEIF(KIdx_2 < KIdx_1)THEN 
                                LenRatio_CCSIGX(KIdx_1,J,K) = ABS((coordinates_x(KIdx_1+1)-crosspoint_XX(idx_start,Rx)%Global_Coord%Coord_X))/Cdelx(KIdx_1)
                    		    LenRatio_CCSIGX(KIdx_2,J,K) = ABS((crosspoint_XX(idx_end,Rx)%Global_Coord%Coord_X-coordinates_x(KIdx_2)))/Cdelx(KIdx_2)
                                CCSIGX(KIdx_2,J,K) = TAR_CONDUCTIVITY(2) * LenRatio_CCSIGX(KIdx_2,J,K) + AIR_CONDUCTIVITY * (1-LenRatio_CCSIGX(KIdx_2,J,K))
                                CCSIGX((KIdx_2+1):(KIdx_1-1),J,K) = AIR_CONDUCTIVITY
                                CCSIGX(KIdx_1,J,K) = AIR_CONDUCTIVITY * LenRatio_CCSIGX(KIdx_1,J,K) + TAR_CONDUCTIVITY(2) * (1-LenRatio_CCSIGX(KIdx_1,J,K))
                            ENDIF
                        ENDIF
                    ENDIF
                    
                    IF(ii==mmx_per(Rx))THEN  !The last intersection point
                        IF(crosspoint_XX(ii,Rx)%Log_In)THEN
                            idx_start=ii
                        ELSEIF(.NOT.crosspoint_XX(ii,Rx)%Log_In)THEN
                            idx_end=ii
                        ENDIF
                        IF (crosspoint_XX(ii,Rx)%Log_In) THEN !The intersection point is the penetration point, that is, the air penetrates into the stratum
                            DO i=1,NX
                                Logic_1=(crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_X > coordinates_x(i)).AND. &
                                      & (crosspoint_XX(mmx_per(Rx),Rx)%Global_Coord%Coord_X < coordinates_x(i+1))
                                IF(Logic_1) THEN
                                    KIdx_1=i
                                    EXIT
                                ENDIF
                            ENDDO
                            CCSIGX((KIdx_1+1):NX,J,K) = TAR_CONDUCTIVITY(2)
                        ELSE  !The intersection point is the exit point, that is, it penetrates into the air from the stratum
                            DO i=1,NX
                                Logic_1=(crosspoint_XX(ii,Rx)%Global_Coord%Coord_X > coordinates_x(i)).AND. &
                                      & (crosspoint_XX(ii,Rx)%Global_Coord%Coord_X < coordinates_x(i+1))
                                IF(Logic_1) THEN
                                    KIdx_1=i
                                    EXIT
                                ENDIF
                            ENDDO
                            CCSIGX((KIdx_1+1):NX,J,K) = AIR_CONDUCTIVITY
                        ENDIF
                    ENDIF
                ENDDO				
            ELSEIF(mmx_per(Rx)==0)THEN  !In the absence of an intersection point, directly compare the Z-direction positional relationship between the ray and the z-intersection point on the plane
                IF(coordinates_z(k)>crosspoint_ZZ(1,(j-1)*NXB+2)%Global_Coord%Coord_Z) THEN
                    CCSIGX(:,J,K) = TAR_CONDUCTIVITY(2)
                ELSE
                    CCSIGX(:,J,K) = AIR_CONDUCTIVITY
                ENDIF
            ENDIF
        ENDDO
    ENDDO
    !=========================================Terrain conformal in the y-direction==========================================
    !>The y-conductivity is being treated
	DO k=1, Z_min-1
		DO j=1, NY
			DO i=1,NXB
				CCSIGY(I,J,K)=AIR_CONDUCTIVITY
			ENDDO
		ENDDO
	ENDDO
	DO k=Z_max+1, NZB
		DO j=1, NY
			DO i=1,NXB
				CCSIGY(I,J,K)=TAR_CONDUCTIVITY(2)
			ENDDO
		ENDDO
    ENDDO
    
	DO k=Z_min, Z_max
        DO i=1, NXB    
            Ry=(k-1)*NXB+i
            IF(mmy_per(Ry)==1)THEN  !The case with only one intersection point
                IF (crosspoint_YY(1,Ry)%Log_In) THEN !The intersection point is the penetration point, that is, the air penetrates into the stratum
                    DO j=1,NY
                        Logic_1=(crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y > coordinates_y(j)).AND. &
                                  & (crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y < coordinates_y(j+1))
                        IF(Logic_1) THEN
                            KIdx_1=j
                            EXIT
                        ENDIF
                    ENDDO
                    LenRatio_CCSIGY(i,KIdx_1,k) = ABS((crosspoint_YY(1,Ry)%Global_Coord%Coord_Y-coordinates_y(KIdx_1)))/Cdely(KIdx_1)
                    CCSIGY(i,1:(KIdx_1-1),k) = AIR_CONDUCTIVITY
                    CCSIGY(i,KIdx_1,k) = LenRatio_CCSIGY(i,KIdx_1,k)*AIR_CONDUCTIVITY+(1-LenRatio_CCSIGY(i,KIdx_1,k))*TAR_CONDUCTIVITY(2)
                    CCSIGY(i,(KIdx_1+1):NY,k) = TAR_CONDUCTIVITY(2)
                ELSE  !The intersection point is the exit point, that is, it penetrates into the air from the stratum
                    DO j=1,NY
                        Logic_1=(crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y > coordinates_y(j)).AND. &
                                  & (crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y < coordinates_y(j+1))
                        IF(Logic_1) THEN
                            KIdx_1=j
                            EXIT
                        ENDIF
                    ENDDO
                    LenRatio_CCSIGY(i,KIdx_1,k) = ABS((crosspoint_YY(1,Ry)%Global_Coord%Coord_Y-coordinates_y(KIdx_1)))/Cdely(KIdx_1)
                    CCSIGY(i,1:(KIdx_1-1),k) = TAR_CONDUCTIVITY(2)
                    CCSIGY(i,KIdx_1,K) = LenRatio_CCSIGY(i,KIdx_1,K)*TAR_CONDUCTIVITY(2)+(1-LenRatio_CCSIGY(i,KIdx_1,k))*AIR_CONDUCTIVITY
                    CCSIGY(i,(KIdx_1+1):NY,k) = AIR_CONDUCTIVITY
                ENDIF
            ELSEIF(mmy_per(Ry)>1)THEN  !There are intersection points and the number is greater than one
                idx_start=0;idx_end=0
                DO ii=1,mmy_per(Ry)
                    IF (ii==1)THEN      !First, determine the first intersection point
                        IF (crosspoint_YY(ii,Ry)%Log_In) THEN !The intersection point is the penetration point, that is, the air penetrates into the stratum
                            DO j=1,NY
                                Logic_1=(crosspoint_YY(ii,Ry)%Global_Coord%Coord_Y > coordinates_y(j)).AND. &
                                      & (crosspoint_YY(ii,Ry)%Global_Coord%Coord_Y < coordinates_y(j+1))
                                IF(Logic_1) THEN
                                    KIdx_1=j
                                    EXIT
                                ENDIF
                            ENDDO
                            CCSIGY(i,1:(KIdx_1-1),k) = AIR_CONDUCTIVITY
                        ELSE  !The intersection point is the exit point, that is, it penetrates into the air from the stratum
                            DO j=1,NY
                                Logic_1=(crosspoint_YY(ii,Ry)%Global_Coord%Coord_Y > coordinates_y(j)).AND. &
                                      & (crosspoint_YY(ii,Ry)%Global_Coord%Coord_Y < coordinates_y(j+1))
                                IF(Logic_1) THEN
                                    KIdx_1=j
                                    EXIT
                                ENDIF
                            ENDDO
                            CCSIGY(i,1:(KIdx_1-1),k) = TAR_CONDUCTIVITY(2)
                        ENDIF
                    ENDIF
                    
                    IF(crosspoint_YY(ii,Ry)%Log_In)THEN
                        idx_start=ii
                    ELSEIF(.NOT.crosspoint_YY(ii,Ry)%Log_In)THEN
                        idx_end=ii
                    ENDIF
                    IF((idx_start>0).AND.(idx_end>0))THEN 
                        DO j=1,NY
                            Logic_1=(crosspoint_YY(idx_start,Ry)%Global_Coord%Coord_Y > coordinates_y(j)).AND. &
                                  & (crosspoint_YY(idx_start,Ry)%Global_Coord%Coord_Y < coordinates_y(j+1))
                            Logic_2=(crosspoint_YY(idx_end,Ry)%Global_Coord%Coord_Y > coordinates_y(j)).AND. &
                                  & (crosspoint_YY(idx_end,Ry)%Global_Coord%Coord_Y < coordinates_y(j+1))
                            IF(Logic_1) KIdx_1=j
                            IF(Logic_2) KIdx_2=j
                            IF(Logic_1.and.Logic_2)	EXIT
                        ENDDO
                        IF(idx_start<idx_end) THEN           !The ray first penetrate the air and enter the stratum
                            IF(KIdx_1==KIdx_2) THEN
                                LenRatio_CCSIGY(i,KIdx_1,k) = ABS((crosspoint_YY(idx_end,Ry)%Global_Coord%Coord_Y-crosspoint_YY(idx_start,Ry)%Global_Coord%Coord_Y))/Cdely(KIdx_1)
                                CCSIGY(i,KIdx_1,k) = AIR_CONDUCTIVITY * (1-LenRatio_CCSIGY(i,KIdx_1,K))+TAR_CONDUCTIVITY(2) * LenRatio_CCSIGY(i,KIdx_1,k)
                            ELSEIF(KIdx_2 > KIdx_1)THEN 
                                LenRatio_CCSIGY(i,KIdx_1,k) = ABS((crosspoint_YY(idx_start,Ry)%Global_Coord%Coord_Y-coordinates_y(KIdx_1)))/Cdely(KIdx_1)
                        		LenRatio_CCSIGY(i,KIdx_2,k) = ABS((crosspoint_YY(idx_end,Ry)%Global_Coord%Coord_Y-coordinates_y(KIdx_2)))/Cdely(KIdx_2)
                                CCSIGY(i,KIdx_1,k) = AIR_CONDUCTIVITY * LenRatio_CCSIGY(i,KIdx_1,k) + TAR_CONDUCTIVITY(2) * (1-LenRatio_CCSIGY(i,KIdx_1,k))
                                CCSIGY(i,(KIdx_1+1):(KIdx_2-1),K) = TAR_CONDUCTIVITY(2)
                                CCSIGY(i,KIdx_2,k) = TAR_CONDUCTIVITY(2) * LenRatio_CCSIGY(i,KIdx_2,k) + AIR_CONDUCTIVITY * (1-LenRatio_CCSIGY(i,KIdx_2,k))
                            ENDIF
                        ELSE                                 !The ray first penetrate the stratum and enter the air
                            IF(KIdx_1==KIdx_2) THEN
                                LenRatio_CCSIGY(i,KIdx_1,k) = ABS((crosspoint_YY(idx_start,Ry)%Global_Coord%Coord_Y-crosspoint_YY(idx_end,Ry)%Global_Coord%Coord_Y))/Cdely(KIdx_1)
                                CCSIGY(i,KIdx_1,k) = TAR_CONDUCTIVITY(2) * (1-LenRatio_CCSIGY(i,KIdx_1,K))+AIR_CONDUCTIVITY * LenRatio_CCSIGY(i,KIdx_1,k)
                            ELSEIF(KIdx_2 < KIdx_1)THEN
                                LenRatio_CCSIGY(i,KIdx_1,k) = ABS((coordinates_y(KIdx_1+1)-crosspoint_YY(idx_start,Ry)%Global_Coord%Coord_Y))/Cdely(KIdx_1)
                        		LenRatio_CCSIGY(i,KIdx_2,k) = ABS((crosspoint_YY(idx_end,Ry)%Global_Coord%Coord_Y-coordinates_y(KIdx_2)))/Cdely(KIdx_2)
                                CCSIGY(i,KIdx_2,k) = TAR_CONDUCTIVITY(2) * LenRatio_CCSIGY(i,KIdx_2,K) + AIR_CONDUCTIVITY * (1-LenRatio_CCSIGY(i,KIdx_2,k))
                                CCSIGY(i,(KIdx_2+1):(KIdx_1-1),k) = AIR_CONDUCTIVITY
                                CCSIGY(i,KIdx_1,k) = AIR_CONDUCTIVITY * LenRatio_CCSIGX(i,KIdx_1,K)+TAR_CONDUCTIVITY(2)*(1-LenRatio_CCSIGY(i,KIdx_1,k))
                            ENDIF
                        ENDIF
                    ENDIF
                    
                    IF(ii==mmy_per(Ry))THEN  !The last intersection point	
						IF (crosspoint_YY(ii,Ry)%Log_In) THEN !The intersection point is the penetration point, that is, the air penetrates into the stratum
                            DO j=1,Ny
                                Logic_1=(crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y > coordinates_y(j)).AND. &
                                      & (crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y < coordinates_y(j+1))
                                IF(Logic_1) THEN
                                    KIdx_1=j
                                    EXIT
                                ENDIF
                            ENDDO
                            CCSIGY(i,(KIdx_1+1):NY,k) = TAR_CONDUCTIVITY(2)
                        ELSE  !The intersection point is the exit point, that is, it penetrates into the air from the stratum
                            DO j=1,Ny
                                Logic_1=(crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y > coordinates_y(j)).AND. &
                                      & (crosspoint_YY(mmy_per(Ry),Ry)%Global_Coord%Coord_Y < coordinates_y(j+1))
                                IF(Logic_1) THEN
                                    KIdx_1=j
                                    EXIT
                                ENDIF
                            ENDDO
                            CCSIGY(i,(KIdx_1+1):NY,k) = AIR_CONDUCTIVITY
                        ENDIF
                    ENDIF
                ENDDO	
            ELSEIF(mmy_per(Ry)==0)THEN  !In the absence of an intersection point, directly compare the Z-direction positional relationship between the ray and the z-intersection point on the plane
                IF(coordinates_z(k)>crosspoint_ZZ(1,(i-1)*NXB+2)%Global_Coord%Coord_Z) THEN
                    CCSIGY(i,:,k) = TAR_CONDUCTIVITY(2)
                ELSE
                    CCSIGY(i,:,k) = AIR_CONDUCTIVITY
                ENDIF
            ENDIF
        ENDDO
    ENDDO
    !=========================================Terrain conformal in the z-direction==========================================
    !>The z-conductivity is being treated
    DO j=1, NYB       
        DO i=1, NXB    
            Rz=(j-1)*NXB+i
            IF(mmz_per(Rz)>0)THEN
                DO k=1,NZ
                    Logic_1=(crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_Z > coordinates_z(k)).AND. &
                          & (crosspoint_ZZ(mmz_per(Rz),Rz)%Global_Coord%Coord_Z < coordinates_z(k+1))
                    IF(Logic_1) THEN
                        KIdx_1=k
                        EXIT
                    ENDIF
                ENDDO
            ENDIF
           LenRatio_CCSIGZ(i,j,KIdx_1) = (crosspoint_ZZ(1,Rz)%Global_Coord%Coord_Z-coordinates_z(KIdx_1))/Cdelz(KIdx_1)
            CCSIGZ(i,j,1:(KIdx_1-1)) = AIR_CONDUCTIVITY
            CCSIGZ(i,j,KIdx_1) = LenRatio_CCSIGZ(i,j,KIdx_1)*AIR_CONDUCTIVITY+(1-LenRatio_CCSIGZ(i,j,KIdx_1))*TAR_CONDUCTIVITY(2)
            CCSIGZ(i,j,(KIdx_1+1):NZ) = TAR_CONDUCTIVITY(2)
        ENDDO    
    ENDDO
    END SUBROUTINE terrain_conformal