!Copyright (c) 2013 by tdem.org under guide of Xiu Li(lixiu@chd.edu.cn)
!written by Huaifeng Sun(sunhuaifeng@gmail.com) and Xushan Lu(luxushan@gmail.com)
!Code distribution @ tdem.org or sunhuaifeng.com

SUBROUTINE GETDATA
  USE CONSTANTPARAMETERS
  !this line is added by Huaifeng Sun to get the dir 2016-10-30
  IMPLICIT NONE
  LOGICAL ALIVE
  INTEGER III,JJJ
  !this following lines 10-21 are added by Huaifeng Sun to get the dir 2016-10-30
  CHARACTER(255) InputFileName
  
  InputFileName='input.dat'
  INQUIRE(FILE=InputFileName, EXIST=ALIVE)
  IF(.NOT. ALIVE) THEN
     WRITE(10005,*) "input.dat DOES NOT EXIST."
     STOP
  ELSE
     OPEN(234,FILE=InputFileName,STATUS='OLD')
     READ(234,*)CAL_TYPE   !This is the calculation type, possible values are shown below.
     ! IF(CAL_TYPE=='TUNNEL' .OR. CAL_TYPE=='tunnel')THEN
        ! WRITE(10005,*)'The tunnel model calculation switch is set correctly!'
     ! ELSEIF(CAL_TYPE=='SEMI' .OR. CAL_TYPE=='semi')THEN
        ! WRITE(10005,*)'The SEMI-AIRBORNE compute switch is set correctly!'
     IF(CAL_TYPE==1)THEN
        WRITE(10005,*)'The ground model calculation switch is set correctly!'
     ELSEIF(CAL_TYPE==2)THEN
        WRITE(10005,*)'The SEMI-AIRBORNE compute switch is set correctly!'
	 ! ELSEIF(CAL_TYPE==3)THEN
        ! WRITE(10005,*)'The tunnel model calculation switch is set correctly!'
     ELSE
        WRITE(10005,*)'The model calculation switch is not set correctly. Please determine whether to use the ground model or the tunnel model!'
        STOP
     ENDIF
     READ(234,*)SourceLength   
     !The length of source, unit of which is meter, and you are supposed to set SourceLengh/GridSize as an odd number for the consideration of there will exist a central point within the source loop.
     READ(234,*)NX,NY,NZ        !The value of Nx, Ny and Nz varies from model to model.
     READ(234,*)Logic_PML       !Boundary condition switch: 1=CPML absorbing boundary, 0=original Dirichlet (zero field) boundary on the non-uniform grid
     READ(234,*)PML_X,PML_Y,PML_Z !PML thickness in x, y and z directions, only valid when Logic_PML=1
     IF(Logic_PML==1)THEN
        IF(PML_X<1 .OR. PML_Y<1 .OR. PML_Z<1)THEN
           WRITE(10005,*)'Error: PML thickness must be at least 1 when the CPML boundary is enabled!'
           WRITE(*,*)'Error: PML thickness must be at least 1 when the CPML boundary is enabled!'
           STOP
        ENDIF
        PML_X1=PML_X; PML_X2=PML_X
        PML_Y1=PML_Y; PML_Y2=PML_Y
        PML_Z1=PML_Z; PML_Z2=PML_Z
     ELSE
        PML_X1=0;  PML_X2=0
        PML_Y1=0;  PML_Y2=0
        PML_Z1=0;  PML_Z2=0
     ENDIF
     READ(234,*)UniGridNumX1,UniGridNumX2
     READ(234,*)UniGridNumY1,UniGridNumY2
     READ(234,*)UniGridNumZ1,UniGridNumZ2
	 READ(234,*)GridSize        !Most commonly used value is 10m
     READ(234,*)BACKGROUND_CONDUCTIVITY !Most commonly used value is 1e-2
     READ(234,*)TEMP_II                 !It depends on your model, and it should be set to 0 if you are doing homogeneous model calculation.
     ALLOCATE(TAR_X1(TEMP_II))
     ALLOCATE(TAR_X2(TEMP_II))
     ALLOCATE(TAR_Y1(TEMP_II))
     ALLOCATE(TAR_Y2(TEMP_II))
     ALLOCATE(TAR_Z1(TEMP_II))
     ALLOCATE(TAR_Z2(TEMP_II))
     ALLOCATE(TAR_CONDUCTIVITY(TEMP_II))
     DO III=1,TEMP_II
        READ(234,*)TAR_X1(III),TAR_X2(III)
        READ(234,*)TAR_Y1(III),TAR_Y2(III)
        READ(234,*)TAR_Z1(III),TAR_Z2(III)
        READ(234,*)TAR_CONDUCTIVITY(III)
     ENDDO
     READ(234,*)NSTOP           !The maximum iteration number.
     READ(234,*)MAX_OFF_TIME    !The maximum computation time, unit of which is ms
     READ(234,*)RAISETIME,RAISESTEP !Most commonly used value is: Raisetime=1e-6, Raisestep=1e-9
     READ(234,*)WAVE !,WAVESTEP
     READ(234,*)RAMP,RAMPSTEP   !Most commonly used value is: Ramp=1e-6, Rampstep=1e-9
     READ(234,*)TIMESTEP        !Most commonly used value is 1e-7
     READ(234,*)AMP             !It denotes the value of amplitude of transmitting source.
	 READ(234,*)tao_abnormal    !The electrical conductivity of the abnormal body
     !read(234,*)NumRecHeights    !It is determined by your recording configuration
     !allocate(FlightHeight(NumRecHeights),GridNumHeight(NumRecHeights),Nzs_Air(NumRecHeights))
     !READ(234,*)(FlightHeight(iii),iii=1,NumRecHeights)
     READ(234,'(a12)')SOURCE_TYPE !Currently the only possible value of Source_type is 'TIXING_UPCOS'
     READ(234,*)Point_Num                 !It depends on your problem, Number of measured points
	 IF(Point_Num>0) THEN
		ALLOCATE(Points_Observer(Point_Num))
        DO JJJ=1,Point_Num
			READ(234,*)Points_Observer(JJJ)%Idx_Num
			READ(234,*)Points_Observer(JJJ)%Local_Coord_To_Source%Coord_X,Points_Observer(JJJ)%Local_Coord_To_Source%Coord_Y,&
					   Points_Observer(JJJ)%Local_Coord_To_Source%Coord_Z                      !What is read here is the coordinates of the observation point with respect to the center of the source, and the positive and negative values are related to the positive direction of the axis
		ENDDO
	ENDIF
	 !read(234,'(a2)')RecFlag      !Possible values are 'HE' and 'Hz'
     !READ(234,*)NumRecLines
     !read(234,*)RecPointMin,RecPointMax
     !NumRecPoints=RecPointMax-RecPointMin+1
     ! IF(NumRecLines .EQ. 0)THEN
        ! WRITE(10005,*)'No additional receiving points were set, and the program continued to run！'
     ! ELSEIF(NumRecLines .GT. 0)THEN
        ! ALLOCATE(RecLine(NumRecLines),RecPoint(NumRecPoints))
     ! ELSE
        ! WRITE(10005,*)'Extra reception point settings are incorrect. Please refer to the input data file format description. Program terminated abnormally！'
        ! STOP
     ! ENDIF 
     CLOSE(234)
  ENDIF
  !>Detect the anomalous body surface mesh file. Two input formats are supported:
  !! 1. Complex_anomalous.dat - the original text format used by this program;
  !! 2. Complex_anomalous.stl - the ASCII STL format (e.g. exported from GiD).
  !! Only one of the two files should exist in the working folder, and the reading
  !! mode is chosen here. The variables filled later (Vert, Triangular_face_element,
  !! n_point, n_face) keep the same names in both formats.
  INQUIRE(FILE='Complex_anomalous.dat', EXIST=Logic_AnomalousDat)
  INQUIRE(FILE='Complex_anomalous.stl', EXIST=Logic_AnomalousStl)
  IF(Logic_AnomalousDat .AND. Logic_AnomalousStl)THEN
     WRITE(10005,*)'Both Complex_anomalous.dat and Complex_anomalous.stl exist! The .dat format takes precedence, the .stl file is ignored.'
     WRITE(*,*)'Both Complex_anomalous.dat and Complex_anomalous.stl exist! The .dat format takes precedence, the .stl file is ignored.'
  ELSEIF(Logic_AnomalousStl)THEN
     WRITE(10005,*)'Complex_anomalous.stl found, the anomalous body will be read in STL format.'
     WRITE(*,*)'Complex_anomalous.stl found, the anomalous body will be read in STL format.'
  ELSEIF(Logic_AnomalousDat)THEN
     WRITE(10005,*)'Complex_anomalous.dat found, the anomalous body will be read in the original format.'
     WRITE(*,*)'Complex_anomalous.dat found, the anomalous body will be read in the original format.'
  ELSE
     WRITE(10005,*)'Warning: no anomalous body mesh file (Complex_anomalous.dat / Complex_anomalous.stl) is found! The model is treated as homogeneous.'
     WRITE(*,*)'Warning: no anomalous body mesh file (Complex_anomalous.dat / Complex_anomalous.stl) is found! The model is treated as homogeneous.'
  ENDIF
  !>Detect the terrain surface mesh file. Two input formats are supported:
  !! 1. Complex_Terrain.dat - the original text format used by this program;
  !! 2. Complex_Terrain.stl - the ASCII STL format (e.g. exported from GiD).
  !! Only one of the two files should exist in the working folder, and the reading
  !! mode is chosen here. The variables filled later (Node_Label, CoordinatesX/Y/Z,
  !! Element_Label, Element_Node1/2/3, n_point, n_face) keep the same names in both formats.
  INQUIRE(FILE='Complex_Terrain.dat', EXIST=Logic_TerrainDat)
  INQUIRE(FILE='Complex_Terrain.stl', EXIST=Logic_TerrainStl)
  IF(Logic_TerrainDat .AND. Logic_TerrainStl)THEN
     WRITE(10005,*)'Both Complex_Terrain.dat and Complex_Terrain.stl exist! The .dat format takes precedence, the .stl file is ignored.'
     WRITE(*,*)'Both Complex_Terrain.dat and Complex_Terrain.stl exist! The .dat format takes precedence, the .stl file is ignored.'
  ELSEIF(Logic_TerrainStl)THEN
     WRITE(10005,*)'Complex_Terrain.stl found, the terrain will be read in STL format.'
     WRITE(*,*)'Complex_Terrain.stl found, the terrain will be read in STL format.'
  ELSEIF(Logic_TerrainDat)THEN
     WRITE(10005,*)'Complex_Terrain.dat found, the terrain will be read in the original format.'
     WRITE(*,*)'Complex_Terrain.dat found, the terrain will be read in the original format.'
  ELSE
     WRITE(10005,*)'Warning: no terrain mesh file (Complex_Terrain.dat / Complex_Terrain.stl) is found! The model is treated as without terrain.'
     WRITE(*,*)'Warning: no terrain mesh file (Complex_Terrain.dat / Complex_Terrain.stl) is found! The model is treated as without terrain.'
  ENDIF
  ! do iii=1,NumRecHeights
     ! GridNumHeight(iii)=FlightHeight(iii)/GridSize
  ! end do
  ! do iii=1,NumRecPoints,1
     ! RecPoint(iii)=iii+RecPointMin-1
  ! end do
  !Calculate other constants in CONSTANTPARAMETERS.
  NXB=NX+1
  NYB=NY+1
  NZB=NZ+1
  
  SourceGridNum=NINT(SourceLength/GridSize)
    IF(ABS(MOD(SourceGridNum,2))==1) THEN
        Logi_Sourcelenth=.TRUE.        !The number of grids occupied by the source is odd
		print*,'The number of grids in the core area is odd'
    ELSE
        Logi_Sourcelenth=.FALSE.        !The number of grids occupied by the source is even
		print*,'The number of grids in the core area is even'
    ENDIF
    IF(Logi_Sourcelenth) THEN
        NXS=(NX+1)/2
        NYS=(NY+1)/2
        NZS=NZ/2
    ELSE
        NXS=NX/2
        NYS=NY/2
        NZS=NZ/2
    ENDIF
  ! do iii=1,NumRecHeights
     ! NZS_AIR(iii)=NZS-GridNumHeight(iii)
  ! end do
  ! do iii=1,NumRecLines,1
     ! RecLine(iii)=nxs-(NumRecLines-1)/2+iii-1
  ! end do
  !Convert current into current density
  AMP=AMP/(GridSize*GridSize)
  SourceGridNum=int(SourceLength/GridSize)
  ALLOCATE(SOURCE(NSTOP))
ENDSUBROUTINE GETDATA