[cig-commits] r12672 - in seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta: setup src
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Mon Aug 18 10:04:31 PDT 2008
Author: dkomati1
Date: 2008-08-18 10:04:30 -0700 (Mon, 18 Aug 2008)
New Revision: 12672
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/setup/constants.h
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/read_compute_parameters.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/write_seismograms.F90
Log:
fixed instabilities (time steps too large) in the case of 3D models at low resolution;
added a flag called COMPUTE_STORE_SEISMOGRAMS in setup/constants.h to completely suppress seismograms (computing them and storing them to disk at the end of the run).
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/setup/constants.h
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/setup/constants.h 2008-08-18 02:29:55 UTC (rev 12671)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/setup/constants.h 2008-08-18 17:04:30 UTC (rev 12672)
@@ -38,9 +38,12 @@
! (much) faster detection of receivers at high resolution: use grid points only
logical, parameter :: FASTER_RECEIVERS_POINTS_ONLY = .true.
+! suppress calculation and storage of seismograms if needed
+ logical, parameter :: COMPUTE_STORE_SEISMOGRAMS = .true.
+
! decrease the number of MPI messages by 3 but increase the size
! of several MPI buffers by 3 but in order to do that
- logical, parameter :: FEWER_MESSAGES_LARGER_BUFFERS = .false.
+ logical, parameter :: FEWER_MESSAGES_LARGER_BUFFERS = .true.
!
! solver in single or double precision depending on the machine (4 or 8 bytes)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/read_compute_parameters.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/read_compute_parameters.F90 2008-08-18 02:29:55 UTC (rev 12671)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/read_compute_parameters.F90 2008-08-18 17:04:30 UTC (rev 12672)
@@ -731,26 +731,45 @@
endif
+! we cannot honor a two-layer crust with only one element in the crust
+ if(NER_CRUST == 1) ONE_CRUST = .true.
+
!----
!---- change some values in the case of regular PREM with two crustal layers or of 3D models
!----
! case of regular PREM with two crustal layers: change the time step for small meshes
! because of a different size of elements in the radial direction in the crust
- if (HONOR_1D_SPHERICAL_MOHO) then
- if (.not. ONE_CRUST) then
- ! case 1D + two crustal layers
- if (NER_CRUST<2) NER_CRUST=2
- if(NEX_MAX*multiplication_factor <= 160) then
- DT = 0.20d0
- else if(NEX_MAX*multiplication_factor <= 256) then
- DT = 0.20d0
- endif
+! case 1D + two crustal layers
+ if(HONOR_1D_SPHERICAL_MOHO .and. .not. ONE_CRUST .and. NER_CRUST < 2) NER_CRUST = 2
+
+ if (.not. ONE_CRUST) then
+ if(NEX_MAX*multiplication_factor <= 160) then
+ DT = 0.20d0
+ else if(NEX_MAX*multiplication_factor <= 256) then
+ DT = 0.20d0
endif
- else
- ! case 3D
- if (NER_CRUST<2) NER_CRUST=2
- if(NEX_MAX*multiplication_factor <= 160) then
+ endif
+
+!! DK DK I do not understand why we should impose two layers in the crust here
+!! DK DK therefore I got rid of it
+!! DK DK else
+!! DK DK ! case 3D
+!! DK DK if (NER_CRUST<2) NER_CRUST=2
+!! DK DK if(NEX_MAX*multiplication_factor <= 160) then
+!! DK DK DT = 0.15d0
+!! DK DK else if(NEX_MAX*multiplication_factor <= 256) then
+!! DK DK DT = 0.17d0
+!! DK DK else if(NEX_MAX*multiplication_factor <= 320) then
+!! DK DK DT = 0.155d0
+!! DK DK endif
+
+! because of distorted elements (for instance due to topography)
+! we need a smaller time steps at low resolution in the 3D case
+ if(CASE_3D) then
+ if(NEX_MAX*multiplication_factor <= 80) then
+ DT = 0.125d0
+ else if(NEX_MAX*multiplication_factor <= 160) then
DT = 0.15d0
else if(NEX_MAX*multiplication_factor <= 256) then
DT = 0.17d0
@@ -760,7 +779,7 @@
endif
! model 1066A has one very thin layer therefore we need to drastically reduce the time step
- if (REFERENCE_1D_MODEL == REFERENCE_MODEL_1066A) DT = DT*0.20d0
+ if (REFERENCE_1D_MODEL == REFERENCE_MODEL_1066A) DT = DT / 5.d0
!----
! if the chunk is smaller than full 90 degrees, call Brian Savage's set of auto_ner routines
@@ -777,15 +796,13 @@
endif ! of call to Brian Savage's set of auto_ner routines
- if (HONOR_1D_SPHERICAL_MOHO) then
- if (.not. ONE_CRUST) then
- ! case 1D + two crustal layers
- if (NER_CRUST<2) NER_CRUST=2
- endif
- else
- ! case 3D
- if (NER_CRUST<2) NER_CRUST=2
- endif
+! case 1D + two crustal layers
+ if(HONOR_1D_SPHERICAL_MOHO .and. .not. ONE_CRUST .and. NER_CRUST < 2) NER_CRUST = 2
+!! DK DK I do not understand why we should impose two layers in the crust here
+!! DK DK therefore I got rid of it
+!! DK DK else
+!! DK DK ! case 3D
+!! DK DK if (NER_CRUST<2) NER_CRUST=2
! determine suitable attenuation periods for the mesh
call auto_attenuation_periods(ANGULAR_WIDTH_XI_IN_DEGREES,NEX_MAX,MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.F90 2008-08-18 02:29:55 UTC (rev 12671)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.F90 2008-08-18 17:04:30 UTC (rev 12672)
@@ -1015,7 +1015,8 @@
call get_value_string(STATIONS, 'solver.STATIONS', rec_filename)
! locate receivers in the crust in the mesh
- call locate_receivers(myrank,DT,NSTEP,NSPEC_CRUST_MANTLE,NGLOB_CRUST_MANTLE,ibool_crust_mantle, &
+ if(COMPUTE_STORE_SEISMOGRAMS) &
+ call locate_receivers(myrank,DT,NSTEP,NSPEC_CRUST_MANTLE,NGLOB_CRUST_MANTLE,ibool_crust_mantle, &
xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle, &
xigll,yigll,zigll,trim(rec_filename), &
nrec,islice_selected_rec,ispec_selected_rec, &
@@ -1047,6 +1048,8 @@
!--- select local receivers
+ if(COMPUTE_STORE_SEISMOGRAMS) then
+
! count number of receivers located in this slice
nrec_local = 0
do irec = 1,nrec
@@ -1114,6 +1117,8 @@
endif
endif
+ endif ! of if(COMPUTE_STORE_SEISMOGRAMS)
+
!! DK DK end of section with sources and receivers excluded in the serial case
#else
t0 = 0.03d0
@@ -1215,9 +1220,15 @@
else
write(IMAIN,*) 'no general mantle anisotropy'
endif
+
write(IMAIN,*)
+ if(COMPUTE_STORE_SEISMOGRAMS) then
+ write(IMAIN,*) 'computing and storing seismograms'
+ else
+ write(IMAIN,*) 'not computing and storing seismograms'
+ endif
+
write(IMAIN,*)
- write(IMAIN,*)
endif
@@ -1545,7 +1556,7 @@
#ifdef USE_MPI
! allocate seismogram array
- if (nrec_local > 0) then
+ if(COMPUTE_STORE_SEISMOGRAMS .and. nrec_local > 0) then
allocate(uxdstore(nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS),stat=ier)
allocate(uydstore(nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS),stat=ier)
allocate(uzdstore(nrec_local,NTSTEP_BETWEEN_OUTPUT_SEISMOS),stat=ier)
@@ -2437,7 +2448,7 @@
! store the seismograms only if there is at least one receiver located in this slice
#ifdef USE_MPI
- if (nrec_local > 0) then
+ if(COMPUTE_STORE_SEISMOGRAMS .and. nrec_local > 0) then
do irec_local = 1,nrec_local
@@ -2486,7 +2497,7 @@
endif ! nrec_local
! write the current or final seismograms
- if(seismo_current == NTSTEP_BETWEEN_OUTPUT_SEISMOS .or. it == it_end) then
+ if(COMPUTE_STORE_SEISMOGRAMS .and. (seismo_current == NTSTEP_BETWEEN_OUTPUT_SEISMOS .or. it == it_end)) then
call write_seismograms(myrank,uxdstore,uydstore,uzdstore,number_receiver_global,station_name, &
network_name,stlat,stlon,stele,nrec,nrec_local,DT,t0,it_end, &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/write_seismograms.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/write_seismograms.F90 2008-08-18 02:29:55 UTC (rev 12671)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/write_seismograms.F90 2008-08-18 17:04:30 UTC (rev 12672)
@@ -104,6 +104,9 @@
logical SAVE_ALL_SEISMOS_IN_ONE_FILE
logical USE_BINARY_FOR_LARGE_FILE
+! suppress calculation and storage of seismograms if needed
+ if(.not. COMPUTE_STORE_SEISMOGRAMS) return
+
! check that the sum of the number of receivers in each slice is nrec
#ifdef USE_MPI
call MPI_REDUCE(nrec_local,nrec_tot_found,1,MPI_INTEGER,MPI_SUM,0,MPI_COMM_WORLD,ier)
More information about the cig-commits
mailing list