Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 10 additions & 5 deletions src/specfem3D/get_cmt.f90
Original file line number Diff line number Diff line change
Expand Up @@ -29,8 +29,8 @@ subroutine get_cmt(CMTSOLUTION,yr,jda,mo,da,ho,mi,sec, &
tshift_src,hdur,lat,long,depth,moment_tensor, &
DT,NSOURCES,min_tshift_src_original,user_source_time_function)

use constants, only: IIN,MAX_STRING_LEN,CUSTOM_REAL
use shared_parameters, only: USE_EXTERNAL_SOURCE_FILE,NSTEP_STF,NSOURCES_STF,NOISE_TOMOGRAPHY
use constants, only: IIN,MAX_STRING_LEN,CUSTOM_REAL,TINYVAL
use shared_parameters, only: USE_EXTERNAL_SOURCE_FILE,NSTEP_STF,NSOURCES_STF,NOISE_TOMOGRAPHY,USE_RICKER_TIME_FUNCTION

implicit none

Expand Down Expand Up @@ -364,9 +364,14 @@ subroutine get_cmt(CMTSOLUTION,yr,jda,mo,da,ho,mi,sec, &
endif

! checks half-duration
! null half-duration indicates a Heaviside
! replace with very short error function
if (hdur(isource) < 5.d0 * DT) hdur(isource) = 5.d0 * DT
if (USE_RICKER_TIME_FUNCTION) then
! hdur is dominant frequency f0; only guard against zero/negative
if (hdur(isource) < TINYVAL) hdur(isource) = TINYVAL
else
! null half-duration indicates a Heaviside
! replace with very short error function
if (hdur(isource) < 5.d0 * DT) hdur(isource) = 5.d0 * DT
endif

! reads USER EXTERNAL SOURCE if needed
if (USE_EXTERNAL_SOURCE_FILE) then
Expand Down
33 changes: 29 additions & 4 deletions src/specfem3D/get_elevation.f90
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,9 @@ subroutine get_elevation_and_z_coordinate_all(npoints,plon,plat,pbur,utm_x,utm_y

use constants, only: CUSTOM_REAL,HUGEVAL,ILONGLAT2UTM
use specfem_par, only: ibool,myrank,NSPEC_AB,NGLOB_AB,USE_SOURCES_RECEIVERS_Z, &
xstore,ystore,zstore,NPROC,num_free_surface_faces,free_surface_ispec,free_surface_ijk
xstore,ystore,zstore,NPROC,num_free_surface_faces,free_surface_ispec,free_surface_ijk, &
NSPEC2D_BOTTOM
use shared_parameters, only: BOTTOM_FREE_SURFACE
implicit none

integer, intent(in) :: npoints
Expand All @@ -52,6 +54,7 @@ subroutine get_elevation_and_z_coordinate_all(npoints,plon,plat,pbur,utm_x,utm_y
double precision, dimension(:,:), allocatable :: elevation_all,elevation_distmin_all
real(kind=CUSTOM_REAL) :: xloc,yloc,loc_ele,loc_distmin
integer :: iproc,ier
integer :: num_topo_faces, iface_offset

!debug
!print *,'get elevation: ',npoints
Expand Down Expand Up @@ -107,6 +110,16 @@ subroutine get_elevation_and_z_coordinate_all(npoints,plon,plat,pbur,utm_x,utm_y
elevation_all(:,:) = 0.d0
elevation_distmin_all(:,:) = HUGEVAL

! when the bottom is a free surface, its faces are stored first in free_surface_ispec
! (bottom faces are filled before top faces in get_absorbing_boundary.F90).
! skip those bottom faces so elevation is always taken from the top surface.
if (BOTTOM_FREE_SURFACE) then
iface_offset = NSPEC2D_BOTTOM
else
iface_offset = 0
endif
num_topo_faces = num_free_surface_faces - iface_offset

! gets elevation
do ipoin = 1, npoints
xloc = x_target(ipoin)
Expand All @@ -115,7 +128,7 @@ subroutine get_elevation_and_z_coordinate_all(npoints,plon,plat,pbur,utm_x,utm_y
! get approximate topography elevation at point x/y coordinates
call get_topo_elevation_free(xloc,yloc,loc_ele,loc_distmin, &
NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, &
num_free_surface_faces,free_surface_ispec,free_surface_ijk)
num_topo_faces,free_surface_ispec(iface_offset+1:),free_surface_ijk(:,:,iface_offset+1:))

elevation(ipoin) = loc_ele
elevation_distmin(ipoin) = loc_distmin
Expand Down Expand Up @@ -166,7 +179,9 @@ subroutine get_elevation_and_z_coordinate(lon,lat,utm_x,utm_y,z_target,elevation

use constants
use specfem_par, only: USE_SOURCES_RECEIVERS_Z,ibool,myrank,NSPEC_AB,NGLOB_AB, &
xstore,ystore,zstore,NPROC,num_free_surface_faces,free_surface_ispec,free_surface_ijk
xstore,ystore,zstore,NPROC,num_free_surface_faces,free_surface_ispec,free_surface_ijk, &
NSPEC2D_BOTTOM
use shared_parameters, only: BOTTOM_FREE_SURFACE

implicit none

Expand All @@ -179,6 +194,7 @@ subroutine get_elevation_and_z_coordinate(lon,lat,utm_x,utm_y,z_target,elevation
double precision,dimension(1) :: altitude_rec,distmin_ele
double precision,dimension(NPROC) :: distmin_ele_all,elevation_all
real(kind=CUSTOM_REAL) :: xloc,yloc,loc_ele,loc_distmin
integer :: num_topo_faces,iface_offset

! convert station location to UTM
call utm_geo(lon,lat,utm_x,utm_y,ILONGLAT2UTM)
Expand All @@ -194,10 +210,19 @@ subroutine get_elevation_and_z_coordinate(lon,lat,utm_x,utm_y,z_target,elevation
return
endif

! when the bottom is a free surface, its faces are stored first in free_surface_ispec;
! skip them so elevation is always taken from the top surface.
if (BOTTOM_FREE_SURFACE) then
iface_offset = NSPEC2D_BOTTOM
else
iface_offset = 0
endif
num_topo_faces = num_free_surface_faces - iface_offset

! get approximate topography elevation at point long/lat coordinates
call get_topo_elevation_free(xloc,yloc,loc_ele,loc_distmin, &
NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, &
num_free_surface_faces,free_surface_ispec,free_surface_ijk)
num_topo_faces,free_surface_ispec(iface_offset+1:),free_surface_ijk(:,:,iface_offset+1:))

altitude_rec(1) = loc_ele
distmin_ele(1) = loc_distmin
Expand Down
Loading