[cig-commits] [commit] devel: updates reading in CMTSOLUTION and station files to be less stringent about formatting and spacing (2daa4d3)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Sep 16 08:26:13 PDT 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branch : devel
Link : https://github.com/geodynamics/specfem3d_globe/compare/0aa2df069164153a157cf61edc5ea1bce1e70644...6034d3445ce173fc6a0a1cfd5c31dfd2558a2eaf
>---------------------------------------------------------------
commit 2daa4d3264b0dad63de52d77a99ccf76d0f3d0c9
Author: daniel peter <peterda at ethz.ch>
Date: Tue Sep 2 17:52:45 2014 +0200
updates reading in CMTSOLUTION and station files to be less stringent about formatting and spacing
>---------------------------------------------------------------
2daa4d3264b0dad63de52d77a99ccf76d0f3d0c9
src/auxiliaries/combine_AVS_DX.f90 | 4 +-
src/shared/calendar.f90 | 6 +-
src/specfem3D/get_cmt.f90 | 264 ++++++++++++++++++++++++++++----
src/specfem3D/initialize_simulation.f90 | 17 +-
src/specfem3D/locate_receivers.f90 | 72 +++++----
src/specfem3D/locate_sources.f90 | 7 +-
6 files changed, 295 insertions(+), 75 deletions(-)
diff --git a/src/auxiliaries/combine_AVS_DX.f90 b/src/auxiliaries/combine_AVS_DX.f90
index ac9b15e..feeb638 100644
--- a/src/auxiliaries/combine_AVS_DX.f90
+++ b/src/auxiliaries/combine_AVS_DX.f90
@@ -441,7 +441,7 @@
! get source information for frequency for number of points per lambda
print *,'reading source duration from the CMTSOLUTION file'
call get_cmt(yr,jda,ho,mi,sec,tshift_cmt,hdur,lat,long,depth,moment_tensor, &
- DT,1,min_tshift_cmt_original)
+ DT,1,min_tshift_cmt_original)
! set global element and point offsets to zero
iglobpointoffset = 0
@@ -939,7 +939,7 @@
! get source information
print *,'reading position of the source from the CMTSOLUTION file'
call get_cmt(yr,jda,ho,mi,sec,tshift_cmt,hdur,lat,long,depth,moment_tensor, &
- DT,1,min_tshift_cmt_original)
+ DT,1,min_tshift_cmt_original)
! convert geographic latitude lat (degrees) to geocentric colatitude theta (radians)
call lat_2_geocentric_colat_dble(lat(1),theta)
diff --git a/src/shared/calendar.f90 b/src/shared/calendar.f90
index c17cf38..46d2c7e 100644
--- a/src/shared/calendar.f90
+++ b/src/shared/calendar.f90
@@ -29,10 +29,10 @@
implicit none
- integer yr,mo,da
+ integer :: yr,mo,da
- integer mon(12)
- integer lpyr
+ integer :: mon(12)
+ integer :: lpyr
data mon /0,31,59,90,120,151,181,212,243,273,304,334/
julian_day = da + mon(mo)
diff --git a/src/specfem3D/get_cmt.f90 b/src/specfem3D/get_cmt.f90
index e6932b1..f895d87 100644
--- a/src/specfem3D/get_cmt.f90
+++ b/src/specfem3D/get_cmt.f90
@@ -26,9 +26,10 @@
!=====================================================================
subroutine get_cmt(yr,jda,ho,mi,sec,tshift_cmt,hdur,lat,long,depth,moment_tensor, &
- DT,NSOURCES,min_tshift_cmt_original)
+ DT,NSOURCES,min_tshift_cmt_original)
- use constants
+ use constants,only: IIN,IMAIN,USE_FORCE_POINT_SOURCE,EXTERNAL_SOURCE_TIME_FUNCTION, &
+ RHOAV,R_EARTH,PI,GRAV,TINYVAL
implicit none
@@ -42,13 +43,13 @@
double precision, dimension(NSOURCES), intent(out) :: tshift_cmt,hdur,lat,long,depth
double precision, dimension(6,NSOURCES), intent(out) :: moment_tensor
-!--- local variables below
-
- integer mo,da,julian_day,isource
- double precision scaleM
- double precision t_shift(NSOURCES)
- character(len=5) datasource
- character(len=256) string
+ ! local variables below
+ integer :: mo,da,julian_day,isource
+ integer :: i,itype,istart,iend,ier
+ double precision :: scaleM
+ double precision :: t_shift(NSOURCES)
+ !character(len=5) :: datasource
+ character(len=256) :: string
! initializes
lat(:) = 0.d0
@@ -58,76 +59,236 @@
tshift_cmt(:) = 0.d0
hdur(:) = 0.d0
moment_tensor(:,:) = 0.d0
- yr = 0
- jda = 0
- ho = 0
- mi = 0
- sec = 0.d0
!
!---- read hypocenter info
!
- open(unit = 1,file='DATA/CMTSOLUTION',status='old',action='read')
+ open(unit = IIN,file='DATA/CMTSOLUTION',status='old',action='read',iostat=ier)
+ if (ier /= 0) stop 'Error opening DATA/CMTSOLUTION file'
! read source number isource
do isource = 1,NSOURCES
- read(1,"(a256)") string
+ ! initializes
+ yr = 0
+ da = 0
+ ho = -1
+ mi = -1
+ sec = -1.d0
+
+ ! gets header line
+ read(IIN,"(a256)",iostat=ier) string
+ if (ier /= 0) then
+ write(IMAIN,*) 'Error reading header line in source ',isource
+ stop 'Error reading header line in station in CMTSOLUTION file'
+ endif
+
! skips empty lines
do while( len_trim(string) == 0 )
- read(1,"(a256)") string
+ read(IIN,"(a256)",iostat=ier) string
+ if (ier /= 0) then
+ write(IMAIN,*) 'Error reading header line in source ',isource
+ stop 'Error reading header line in station in CMTSOLUTION file'
+ endif
enddo
- ! read header with event information
- read(string,"(a4,i5,i3,i3,i3,i3,f6.2)") datasource,yr,mo,da,ho,mi,sec
- jda=julian_day(yr,mo,da)
+ ! debug
+ !print*,'line ----',string,'----'
+
+ ! read header line with event information
+ ! old line: read(string,"(a4,i5,i3,i3,i3,i3,f6.2)") datasource,yr,mo,da,ho,mi,sec
+
+ ! reads in chunks of the first line (which contain numbers)
+ ! to get rid of the first datasource qualifyer string like PDE,PDEQ,.. which can have variable length)
+ !
+ ! reads in year,month,day,hour,minutes,seconds
+ istart = 1
+ do itype = 1,6
+ ! determines where first number starts
+ do i = istart,len_trim(string)
+ if (is_numeric(string(i:i)) ) then
+ istart = i
+ exit
+ endif
+ enddo
+ if ( istart >= len_trim(string) ) stop 'Error determining datasource length in header line in CMTSOLUTION file'
+ if ( istart <= 1 ) stop 'Error determining datasource length in header line in CMTSOLUTION file'
+
+ ! determines end and length of number
+ iend = istart
+ do i = istart,len_trim(string)
+ if (itype /= 6 ) then
+ ! integer values
+ if (.not. is_numeric(string(i:i))) then
+ iend = i
+ exit
+ endif
+ else
+ ! seconds will have a digit number
+ ! digit numbers, e.g. 39.60, can contain '.'
+ if (.not. is_digit(string(i:i))) then
+ iend = i
+ exit
+ endif
+ endif
+ enddo
+ iend = iend-1
+ if ( iend >= len_trim(string) ) stop 'Error determining number length in header line in CMTSOLUTION file'
+ if ( iend < istart ) stop 'Error determining number with negative length in header line in CMTSOLUTION file'
+
+ ! debug
+ !print*,itype,'line ----',string(istart:iend),'----'
+
+ ! reads in event time information
+ select case( itype )
+ case( 1 )
+ ! year (as integer value)
+ read(string(istart:iend),*) yr
+ case( 2 )
+ ! month (as integer value)
+ read(string(istart:iend),*) mo
+ case( 3 )
+ ! day (as integer value)
+ read(string(istart:iend),*) da
+ case( 4 )
+ ! hour (as integer value)
+ read(string(istart:iend),*) ho
+ case( 5 )
+ ! minutes (as integer value)
+ read(string(istart:iend),*) mi
+ case( 6 )
+ ! seconds (as float value)
+ read(string(istart:iend),*) sec
+ end select
+
+ ! advances string
+ istart = iend + 1
+ enddo
+
+
+ ! checks time information
+ if (yr <= 0 .or. yr > 10000) then
+ write(IMAIN,*) 'Error reading year: ',yr,' in source ',isource,'is invalid'
+ stop 'Error reading year out of header line in CMTSOLUTION file'
+ endif
+ if (mo < 1 .or. mo > 12) then
+ write(IMAIN,*) 'Error reading month: ',mo,' in source ',isource,'is invalid'
+ stop 'Error reading month out of header line in CMTSOLUTION file'
+ endif
+ if (da < 1 .or. da > 31) then
+ write(IMAIN,*) 'Error reading day: ',da,' in source ',isource,'is invalid'
+ stop 'Error reading day of header line in CMTSOLUTION file'
+ endif
+ if (ho < 0 .or. ho > 24) then
+ write(IMAIN,*) 'Error reading hour: ',ho,' in source ',isource,'is invalid'
+ stop 'Error reading hour of header line in CMTSOLUTION file'
+ endif
+ if (mi < 0 .or. mi > 59) then
+ write(IMAIN,*) 'Error reading minute: ',mi,' in source ',isource,'is invalid'
+ stop 'Error reading minute of header line in CMTSOLUTION file'
+ endif
+ if (sec < 0.0 .or. sec >= 60.0) then
+ write(IMAIN,*) 'Error reading second: ',sec,' in source ',isource,'is invalid'
+ stop 'Error reading second of header line in CMTSOLUTION file'
+ endif
+
+ ! gets julian day number
+ jda = julian_day(yr,mo,da)
! ignore line with event name
- read(1,"(a)") string
+ read(IIN,"(a)",iostat=ier) string
+ if (ier /= 0) then
+ write(IMAIN,*) 'Error reading event name in source ',isource
+ stop 'Error reading event name in station in CMTSOLUTION file'
+ endif
! read time shift
- read(1,"(a)") string
+ read(IIN,"(a)",iostat=ier) string
+ if (ier /= 0) then
+ write(IMAIN,*) 'Error reading time shift in source ',isource
+ stop 'Error reading time shift in station in CMTSOLUTION file'
+ endif
!read(string(12:len_trim(string)),*) tshift_cmt(isource)
read(string(12:len_trim(string)),*) t_shift(isource)
! read half duration
- read(1,"(a)") string
+ read(IIN,"(a)",iostat=ier) string
+ if (ier /= 0) then
+ write(IMAIN,*) 'Error reading half duration in source ',isource
+ stop 'Error reading half duration in station in CMTSOLUTION file'
+ endif
read(string(15:len_trim(string)),*) hdur(isource)
! read latitude
- read(1,"(a)") string
+ read(IIN,"(a)",iostat=ier) string
+ if (ier /= 0) then
+ write(IMAIN,*) 'Error reading latitude in source ',isource
+ stop 'Error reading latitude in station in CMTSOLUTION file'
+ endif
read(string(10:len_trim(string)),*) lat(isource)
! read longitude
- read(1,"(a)") string
+ read(IIN,"(a)",iostat=ier) string
+ if (ier /= 0) then
+ write(IMAIN,*) 'Error reading longitude in source ',isource
+ stop 'Error reading longitude in station in CMTSOLUTION file'
+ endif
read(string(11:len_trim(string)),*) long(isource)
! read depth
- read(1,"(a)") string
+ read(IIN,"(a)",iostat=ier) string
+ if (ier /= 0) then
+ write(IMAIN,*) 'Error reading depth in source ',isource
+ stop 'Error reading depth in station in CMTSOLUTION file'
+ endif
read(string(7:len_trim(string)),*) depth(isource)
! read Mrr
- read(1,"(a)") string
+ read(IIN,"(a)",iostat=ier) string
+ if (ier /= 0) then
+ write(IMAIN,*) 'Error reading Mrr in source ',isource
+ stop 'Error reading Mrr in station in CMTSOLUTION file'
+ endif
read(string(5:len_trim(string)),*) moment_tensor(1,isource)
! read Mtt
- read(1,"(a)") string
+ read(IIN,"(a)",iostat=ier) string
+ if (ier /= 0) then
+ write(IMAIN,*) 'Error reading Mtt in source ',isource
+ stop 'Error reading Mtt in station in CMTSOLUTION file'
+ endif
read(string(5:len_trim(string)),*) moment_tensor(2,isource)
! read Mpp
- read(1,"(a)") string
+ read(IIN,"(a)",iostat=ier) string
+ if (ier /= 0) then
+ write(IMAIN,*) 'Error reading Mpp in source ',isource
+ stop 'Error reading Mpp in station in CMTSOLUTION file'
+ endif
read(string(5:len_trim(string)),*) moment_tensor(3,isource)
! read Mrt
- read(1,"(a)") string
+ read(IIN,"(a)",iostat=ier) string
+ if (ier /= 0) then
+ write(IMAIN,*) 'Error reading Mrt in source ',isource
+ stop 'Error reading Mrt in station in CMTSOLUTION file'
+ endif
read(string(5:len_trim(string)),*) moment_tensor(4,isource)
! read Mrp
- read(1,"(a)") string
+ read(IIN,"(a)",iostat=ier) string
+ if (ier /= 0) then
+ write(IMAIN,*) 'Error reading Mrp in source ',isource
+ stop 'Error reading Mrp in station in CMTSOLUTION file'
+ endif
read(string(5:len_trim(string)),*) moment_tensor(5,isource)
! read Mtp
- read(1,"(a)") string
+ read(IIN,"(a)",iostat=ier) string
+ if (ier /= 0) then
+ write(IMAIN,*) 'Error reading Mtp in source ',isource
+ stop 'Error reading Mtp in station in CMTSOLUTION file'
+ endif
read(string(5:len_trim(string)),*) moment_tensor(6,isource)
! checks half-duration
@@ -150,7 +311,7 @@
hdur(:) = 0.d0
end if
- close(1)
+ close(IIN)
! Sets tshift_cmt to zero to initiate the simulation!
if (NSOURCES == 1) then
@@ -172,5 +333,40 @@
scaleM = 1.d7 * RHOAV * (R_EARTH**5) * PI*GRAV*RHOAV
moment_tensor(:,:) = moment_tensor(:,:) / scaleM
- end subroutine get_cmt
+ contains
+
+ !--------------------------------------------------------------
+
+ logical function is_numeric(char)
+
+ ! returns .true. if input character is a number
+
+ implicit none
+ character(len=1), intent(in) :: char
+
+ is_numeric = .false.
+
+ if ( index('0123456789', char) /= 0) then
+ is_numeric = .true.
+ endif
+ end function
+
+ !--------------------------------------------------------------
+
+ logical function is_digit(char)
+
+ ! returns .true. if input character is a number or a '.'
+
+ implicit none
+ character(len=1), intent(in) :: char
+
+ is_digit = .false.
+
+ if ( index('0123456789.', char) /= 0) then
+ is_digit = .true.
+ endif
+
+ end function
+
+ end subroutine get_cmt
diff --git a/src/specfem3D/initialize_simulation.f90 b/src/specfem3D/initialize_simulation.f90
index 1d81eff..e3f4aff 100644
--- a/src/specfem3D/initialize_simulation.f90
+++ b/src/specfem3D/initialize_simulation.f90
@@ -34,7 +34,7 @@
! local parameters
integer :: sizeprocs
- integer :: ios
+ integer :: ier
character(len=150) :: dummystring
! sizeprocs returns number of processes started (should be equal to NPROCTOT).
@@ -78,8 +78,8 @@
! open main output file, only written to by process 0
if (myrank == 0 .and. IMAIN /= ISTANDARD_OUTPUT) then
- open(unit=IMAIN,file=trim(OUTPUT_FILES)//'/output_solver.txt',status='unknown',action='write',iostat=ios)
- if (ios /= 0 ) call exit_MPI(myrank,'Error opening file output_solver.txt for writing output info')
+ open(unit=IMAIN,file=trim(OUTPUT_FILES)//'/output_solver.txt',status='unknown',action='write',iostat=ier)
+ if (ier /= 0 ) call exit_MPI(myrank,'Error opening file output_solver.txt for writing output info')
endif
if (myrank == 0) then
@@ -228,11 +228,14 @@
! get total number of receivers
if (myrank == 0) then
- open(unit=IIN,file=STATIONS,iostat=ios,status='old',action='read')
+ open(unit=IIN,file=STATIONS,iostat=ier,status='old',action='read')
nrec = 0
- do while(ios == 0)
- read(IIN,"(a)",iostat=ios) dummystring
- if (ios == 0) nrec = nrec + 1
+ do while(ier == 0)
+ read(IIN,"(a)",iostat=ier) dummystring
+ if (ier == 0) then
+ ! excludes empty lines
+ if (len_trim(dummystring) > 0 ) nrec = nrec + 1
+ endif
enddo
close(IIN)
endif
diff --git a/src/specfem3D/locate_receivers.f90 b/src/specfem3D/locate_receivers.f90
index 44e6945..beb5dbb 100644
--- a/src/specfem3D/locate_receivers.f90
+++ b/src/specfem3D/locate_receivers.f90
@@ -139,6 +139,7 @@
logical :: located_target
character(len=2) :: bic
+ character(len=256) :: string
integer :: nrec_SUBSET_current_size,irec_in_this_subset,irec_already_done
double precision, allocatable, dimension(:) :: x_found_subset,y_found_subset,z_found_subset
@@ -186,27 +187,18 @@
! use 10 times the distance as a criterion for source detection
typical_size = 10. * typical_size
- if (myrank == 0) then
- write(IMAIN,*)
- write(IMAIN,*) '****************************'
- write(IMAIN,*) 'reading receiver information'
- write(IMAIN,*) '****************************'
- write(IMAIN,*)
- call flush_IMAIN()
- endif
-
! allocate memory for arrays using number of stations
allocate(epidist(nrec), &
- ix_initial_guess(nrec), &
- iy_initial_guess(nrec), &
- iz_initial_guess(nrec), &
- x_target(nrec), &
- y_target(nrec), &
- z_target(nrec), &
- x_found(nrec), &
- y_found(nrec), &
- z_found(nrec), &
- final_distance(nrec),stat=ier)
+ ix_initial_guess(nrec), &
+ iy_initial_guess(nrec), &
+ iz_initial_guess(nrec), &
+ x_target(nrec), &
+ y_target(nrec), &
+ z_target(nrec), &
+ x_found(nrec), &
+ y_found(nrec), &
+ z_found(nrec), &
+ final_distance(nrec),stat=ier)
if (ier /= 0) call exit_MPI(myrank,'Error allocating temporary receiver arrays')
! initializes search distances
@@ -214,12 +206,40 @@
! read that STATIONS file on the master
if (myrank == 0) then
+ ! user output
+ write(IMAIN,*) 'reading receiver information...'
+ write(IMAIN,*)
+ call flush_IMAIN()
+
+ ! opens station file STATIONS or STATIONS_ADJOINT
open(unit=IIN,file=trim(rec_filename),status='old',action='read',iostat=ier)
if (ier /= 0 ) call exit_MPI(myrank,'Error opening STATIONS file')
! loop on all the stations to read station information
do irec = 1,nrec
- read(IIN,*,iostat=ier) station_name(irec),network_name(irec),stlat(irec),stlon(irec),stele(irec),stbur(irec)
+
+ ! old line:
+ !read(IIN,*,iostat=ier) station_name(irec),network_name(irec),stlat(irec),stlon(irec),stele(irec),stbur(irec)
+
+ ! reads in line as string
+ read(IIN,"(a256)",iostat=ier) string
+ if (ier /= 0) then
+ write(IMAIN,*) 'Error reading in station ',irec
+ call exit_MPI(myrank,'Error reading in station in STATIONS file')
+ endif
+
+ ! skips empty lines
+ do while( len_trim(string) == 0 )
+ read(IIN,"(a256)",iostat=ier) string
+ if (ier /= 0) then
+ write(IMAIN,*) 'Error reading in station ',irec
+ call exit_MPI(myrank,'Error reading in station in STATIONS file')
+ endif
+ enddo
+
+ ! reads in station information
+ read(string(1:len_trim(string)),*,iostat=ier) station_name(irec),network_name(irec), &
+ stlat(irec),stlon(irec),stele(irec),stbur(irec)
if (ier /= 0) then
write(IMAIN,*) 'Error reading in station ',irec
call exit_MPI(myrank,'Error reading in station in STATIONS file')
@@ -227,8 +247,8 @@
! checks latitude
if (stlat(irec) < -90.d0 .or. stlat(irec) > 90.d0) then
- print*,'Error station ',trim(station_name(irec)),': latitude ',stlat(irec),' is invalid, please check STATIONS record'
- close(IIN)
+ write(IMAIN,*) 'Error station ',trim(station_name(irec)),': latitude ',stlat(irec), &
+ ' is invalid, please check STATIONS record'
call exit_MPI(myrank,'Error station latitude invalid')
endif
@@ -502,9 +522,8 @@
write(IMAIN,*) 'Stations sorted by epicentral distance:'
do i = 1,nrec
irec = irec_dist_ordered(i)
- write(IMAIN,*) 'Station #',irec,': ',station_name(irec)(1:len_trim(station_name(irec))), &
- '.',network_name(irec)(1:len_trim(network_name(irec))), &
- ' epicentral distance: ',sngl(epidist(irec)),' degrees'
+ write(IMAIN,'(a,i6,a,a24,a,f12.6,a)') ' Station #',irec,': ',trim(station_name(irec))//'.'//trim(network_name(irec)), &
+ ' epicentral distance: ',sngl(epidist(irec)),' degrees'
enddo
deallocate(irec_dist_ordered)
@@ -823,8 +842,7 @@
if (myrank == 0) then
! appends receiver locations to sr.vtk file
- open(IOUT_VTK,file='OUTPUT_FILES/sr_tmp.vtk', &
- position='append',status='old',iostat=ier)
+ open(IOUT_VTK,file='OUTPUT_FILES/sr_tmp.vtk',position='append',status='old',iostat=ier)
if (ier /= 0 ) call exit_MPI(myrank,'Error opening and appending receivers to file sr_tmp.vtk')
! chooses best receivers locations
diff --git a/src/specfem3D/locate_sources.f90 b/src/specfem3D/locate_sources.f90
index d34ec25..e2df265 100644
--- a/src/specfem3D/locate_sources.f90
+++ b/src/specfem3D/locate_sources.f90
@@ -143,8 +143,11 @@
final_distance_source(:) = HUGEVAL
! read all the sources
- if (myrank == 0) call get_cmt(yr,jda,ho,mi,sec,tshift_cmt,hdur,lat,long,depth,moment_tensor, &
- DT,NSOURCES,min_tshift_cmt_original)
+ if (myrank == 0) then
+ ! only master process reads in CMTSOLUTION file
+ call get_cmt(yr,jda,ho,mi,sec,tshift_cmt,hdur,lat,long,depth,moment_tensor, &
+ DT,NSOURCES,min_tshift_cmt_original)
+ endif
! broadcast the information read on the master to the nodes
call bcast_all_dp(tshift_cmt,NSOURCES)
More information about the CIG-COMMITS
mailing list