[cig-commits] r10444 - in seismo/2D/SPECFEM2D/branches/SEM_tangente_branch: . DATA
nlegoff at geodynamics.org
nlegoff at geodynamics.org
Fri Feb 15 01:48:23 PST 2008
Author: nlegoff
Date: 2008-02-15 01:48:21 -0800 (Fri, 15 Feb 2008)
New Revision: 10444
Modified:
seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file
seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file_Abel_Balanche_bathy_source_solid
seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file_M2_UPPA
seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file_Ronan_SV_30
seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file_canyon
seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file_no_canyon
seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file_unstruct
seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/locate_receivers.F90
seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/meshfem2D.F90
seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/specfem2D.F90
Log:
normal to the surface for the source and the receivers is now related to the mesh, and no longer to an external file containing a list of points delimiting the velocity model.
Modified: seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file
===================================================================
--- seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file 2007-03-26 18:32:04 UTC (rev 10443)
+++ seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file 2008-02-15 09:48:21 UTC (rev 10444)
@@ -10,7 +10,6 @@
materials_file = ./DATA/Mesh_canyon/canyon_materials_file # file containing the material number for each element
free_surface_file = ./DATA/Mesh_canyon/canyon_free_surface_file # file containing the free surface
absorbing_surface_file = ./DATA/Mesh_canyon/canyon_absorbing_surface_file # file containing the absorbing surface
-tangential_detection_curve_file = ./DATA/FAILLE_MIROIR/courbe_eros_nodes # file containing the curve delimiting the velocity model
# parameters concerning partitioning
nproc = 1 # number of processes
@@ -52,6 +51,7 @@
Mzz = 1. # Mzz component (for a moment tensor source only)
Mxz = 0. # Mxz component (for a moment tensor source only)
factor = 1.d10 # amplification factor
+force_normal_to_surface = .false. # angleforce normal to surface (source must be on a surface element)
# constants for attenuation
N_SLS = 2 # number of standard linear solids for attenuation
@@ -64,7 +64,7 @@
generate_STATIONS = .true. # creates a STATION file in ./DATA
nreceiverlines = 1 # number of receiver lines
anglerec = 0.d0 # angle to rotate components at receivers
-tangential_detection = .false. # detecting the normal at each receivers
+rotate_seismograms_RT = .false. # rotate the seismograms (normal to surface)
# first receiver line
nrec = 11 # number of receivers
Modified: seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file_Abel_Balanche_bathy_source_solid
===================================================================
--- seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file_Abel_Balanche_bathy_source_solid 2007-03-26 18:32:04 UTC (rev 10443)
+++ seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file_Abel_Balanche_bathy_source_solid 2008-02-15 09:48:21 UTC (rev 10444)
@@ -10,7 +10,6 @@
materials_file = ./DATA/Mesh_canyon/canyon_materials_file # file containing the material number for each element
free_surface_file = ./DATA/Mesh_canyon/canyon_free_surface_file # file containing the free surface
absorbing_surface_file = ./DATA/Mesh_canyon/canyon_absorbing_surface_file # file containing the absorbing surface
-tangential_detection_curve_file = ./DATA/FAILLE_MIROIR/courbe_eros_nodes # file containing the curve delimiting the velocity model
# parameters concerning partitioning
nproc = 1 # number of processes
@@ -51,6 +50,7 @@
Mzz = 1. # Mzz component (for a moment tensor source only)
Mxz = 0. # Mxz component (for a moment tensor source only)
factor = 1.d10 # amplification factor
+force_normal_to_surface = .false. # angleforce normal to surface (source must be on a surface element)
# constants for attenuation
N_SLS = 2 # number of standard linear solids for attenuation
@@ -63,7 +63,7 @@
generate_STATIONS = .true. # creates a STATION file in ./DATA
nreceiverlines = 2 # number of receiver lines
anglerec = 0.d0 # angle to rotate components at receivers
-tangential_detection = .false. # detecting the normal at each receivers
+rotate_seismograms_RT = .false. # rotate the seismograms (normal to surface)
# first receiver line (in the ocean)
nrec = 11 # number of receivers
Modified: seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file_M2_UPPA
===================================================================
--- seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file_M2_UPPA 2007-03-26 18:32:04 UTC (rev 10443)
+++ seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file_M2_UPPA 2008-02-15 09:48:21 UTC (rev 10444)
@@ -10,7 +10,6 @@
materials_file = ./DATA/Mesh_canyon/canyon_materials_file # file containing the material number for each element
free_surface_file = ./DATA/Mesh_canyon/canyon_free_surface_file # file containing the free surface
absorbing_surface_file = ./DATA/Mesh_canyon/canyon_absorbing_surface_file # file containing the absorbing surface
-tangential_detection_curve_file = ./DATA/FAILLE_MIROIR/courbe_eros_nodes # file containing the curve delimiting the velocity model
# parameters concerning partitioning
nproc = 1 # number of processes
@@ -52,6 +51,7 @@
Mzz = 1. # Mzz component (for a moment tensor source only)
Mxz = 0. # Mxz component (for a moment tensor source only)
factor = 1.d10 # amplification factor
+force_normal_to_surface = .false. # angleforce normal to surface (source must be on a surface element)
# constants for attenuation
N_SLS = 2 # number of standard linear solids for attenuation
@@ -64,7 +64,7 @@
generate_STATIONS = .true. # creates a STATION file in ./DATA
nreceiverlines = 1 # number of receiver lines
anglerec = 0.d0 # angle to rotate components at receivers
-tangential_detection = .false. # detecting the normal at each receivers
+rotate_seismograms_RT = .false. # rotate the seismograms (normal to surface)
# first receiver line
nrec = 11 # number of receivers
Modified: seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file_Ronan_SV_30
===================================================================
--- seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file_Ronan_SV_30 2007-03-26 18:32:04 UTC (rev 10443)
+++ seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file_Ronan_SV_30 2008-02-15 09:48:21 UTC (rev 10444)
@@ -10,7 +10,6 @@
materials_file = ./DATA/unstructured_fluide_solide_test/mat # file containing the material number for each element
free_surface_file = ./DATA/unstructured_fluide_solide_test/surface_free # file containing the free surface
absorbing_surface_file = ./DATA/unstructured_fluide_solide_test/surface_abs # file containing the absorbing surface
-tangential_detection_curve_file = ./DATA/FAILLE_MIROIR/courbe_eros_nodes # file containing the curve delimiting the velocity model
# parameters concerning partitioning
nproc = 1 # number of processes
@@ -51,6 +50,7 @@
Mzz = 1. # Mzz component (for a moment tensor source only)
Mxz = 0. # Mxz component (for a moment tensor source only)
factor = 1.d10 # amplification factor
+rotate_seismograms_RT = .false. # rotate the seismograms (normal to surface)
# constants for attenuation
N_SLS = 2 # number of standard linear solids for attenuation
@@ -63,7 +63,7 @@
generate_STATIONS = .true. # creates a STATION file in ./DATA
nreceiverlines = 1 # number of receiver lines
anglerec = 0.d0 # angle to rotate components at receivers
-tangential_detection = .false. # detecting the normal at each receivers
+rotate_seismograms_RT = .false. # rotate the seismograms (normal to surface)
# first receiver line
nrec = 11 # number of receivers
Modified: seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file_canyon
===================================================================
--- seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file_canyon 2007-03-26 18:32:04 UTC (rev 10443)
+++ seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file_canyon 2008-02-15 09:48:21 UTC (rev 10444)
@@ -10,7 +10,6 @@
materials_file = ./DATA/Mesh_canyon/canyon_materials_file # file containing the material number for each element
free_surface_file = ./DATA/Mesh_canyon/canyon_free_surface_file # file containing the free surface
absorbing_surface_file = ./DATA/Mesh_canyon/canyon_absorbing_surface_file # file containing the absorbing surface
-tangential_detection_curve_file = ./DATA/FAILLE_MIROIR/courbe_eros_nodes # file containing the curve delimiting the velocity model
# parameters concerning partitioning
nproc = 1 # number of processes
@@ -51,6 +50,7 @@
Mzz = 1. # Mzz component (for a moment tensor source only)
Mxz = 0. # Mxz component (for a moment tensor source only)
factor = 1.d10 # amplification factor
+force_normal_to_surface = .false. # angleforce normal to surface (source must be on a surface element)
# constants for attenuation
N_SLS = 2 # number of standard linear solids for attenuation
@@ -63,7 +63,7 @@
generate_STATIONS = .true. # creates a STATION file in ./DATA
nreceiverlines = 1 # number of receiver lines
anglerec = 0.d0 # angle to rotate components at receivers
-tangential_detection = .false. # detecting the normal at each receivers
+rotate_seismograms_RT = .false. # rotate the seismograms (normal to surface)
# first receiver line
nrec = 60 # number of receivers
Modified: seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file_no_canyon
===================================================================
--- seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file_no_canyon 2007-03-26 18:32:04 UTC (rev 10443)
+++ seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file_no_canyon 2008-02-15 09:48:21 UTC (rev 10444)
@@ -10,7 +10,6 @@
materials_file = ./DATA/Mesh_canyon/canyon_materials_file # file containing the material number for each element
free_surface_file = ./DATA/Mesh_canyon/canyon_free_surface_file # file containing the free surface
absorbing_surface_file = ./DATA/Mesh_canyon/canyon_absorbing_surface_file # file containing the absorbing surface
-tangential_detection_curve_file = ./DATA/FAILLE_MIROIR/courbe_eros_nodes # file containing the curve delimiting the velocity model
# parameters concerning partitioning
nproc = 1 # number of processes
@@ -51,6 +50,7 @@
Mzz = 1. # Mzz component (for a moment tensor source only)
Mxz = 0. # Mxz component (for a moment tensor source only)
factor = 1.d10 # amplification factor
+force_normal_to_surface = .false. # angleforce normal to surface (source must be on a surface element)
# constants for attenuation
N_SLS = 2 # number of standard linear solids for attenuation
@@ -63,7 +63,7 @@
generate_STATIONS = .false. # creates a STATION file in ./DATA
nreceiverlines = 1 # number of receiver lines
anglerec = 0.d0 # angle to rotate components at receivers
-tangential_detection = .false. # detecting the normal at each receivers
+rotate_seismograms_RT = .false. # rotate the seismograms (normal to surface)
# first receiver line
nrec = 60 # number of receivers
Modified: seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file_unstruct
===================================================================
--- seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file_unstruct 2007-03-26 18:32:04 UTC (rev 10443)
+++ seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/DATA/Par_file_unstruct 2008-02-15 09:48:21 UTC (rev 10444)
@@ -10,7 +10,6 @@
materials_file = ./DATA/unstructured_fluide_solide_test/mat # file containing the material number for each element
free_surface_file = ./DATA/unstructured_fluide_solide_test/surface_free # file containing the free surface
absorbing_surface_file = ./DATA/unstructured_fluide_solide_test/surface_abs # file containing the absorbing surface
-tangential_detection_curve_file = ./DATA/FAILLE_MIROIR/courbe_eros_nodes # file containing the curve delimiting the velocity model
# parameters concerning partitioning
nproc = 8 # number of processes
@@ -51,6 +50,7 @@
Mzz = 1. # Mzz component (for a moment tensor source only)
Mxz = 0. # Mxz component (for a moment tensor source only)
factor = 1.d10 # amplification factor
+force_normal_to_surface = .false. # angleforce normal to surface (source must be on a surface element)
# constants for attenuation
N_SLS = 2 # number of standard linear solids for attenuation
@@ -63,6 +63,7 @@
generate_STATIONS = .true. # creates a STATION file in ./DATA
nreceiverlines = 1 # number of receiver lines
anglerec = 0.d0 # angle to rotate components at receivers
+rotate_seismograms_RT = .false. # rotate the seismograms (normal to surface)
# first receiver line
nrec = 100 # number of receivers
Modified: seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/locate_receivers.F90
===================================================================
--- seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/locate_receivers.F90 2007-03-26 18:32:04 UTC (rev 10443)
+++ seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/locate_receivers.F90 2008-02-15 09:48:21 UTC (rev 10444)
@@ -47,7 +47,7 @@
subroutine locate_receivers(ibool,coord,nspec,npoin,xigll,zigll,nrec,nrecloc,recloc,which_proc_receiver,nproc,myrank, &
st_xval,st_zval,ispec_selected_rec, &
xi_receiver,gamma_receiver,station_name,network_name,x_source,z_source,coorg,knods,ngnod,npgeo, &
- x_final_receiver, z_final_receiver)
+ recloc_ispec)
implicit none
@@ -95,7 +95,7 @@
double precision, dimension(nrec) :: st_xval,st_zval
! tangential detection
- double precision, dimension(nrec) :: x_final_receiver, z_final_receiver
+ integer, dimension(nrec) :: recloc_ispec
double precision, dimension(nrec,nproc) :: gather_final_distance
@@ -212,9 +212,6 @@
! compute final distance between asked and found
final_distance(irec) = sqrt((st_xval(irec)-x)**2 + (st_zval(irec)-z)**2)
- x_final_receiver(irec) = x
- z_final_receiver(irec) = z
-
enddo
! close receiver file
@@ -257,6 +254,7 @@
if ( which_proc_receiver(irec) == myrank ) then
nrecloc = nrecloc + 1
recloc(nrecloc) = irec
+ recloc_ispec(nrecloc) = ispec_selected_rec(irec)
end if
end do
Modified: seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/meshfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/meshfem2D.F90 2007-03-26 18:32:04 UTC (rev 10443)
+++ seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/meshfem2D.F90 2008-02-15 09:48:21 UTC (rev 10444)
@@ -201,10 +201,7 @@
double precision :: f0_attenuation
! variables used for tangential detection
- logical :: tangential_detection
- character(len=256) :: tangential_detection_curve_file
- integer :: nnodes_tangential_curve
- double precision, dimension(:,:), allocatable :: nodes_tangential_curve
+ logical :: force_normal_to_surface,rotate_seismograms_RT
#if defined USE_METIS || defined USE_SCOTCH
integer :: edgecut
@@ -238,7 +235,6 @@
call read_value_string(IIN,IGNORE_JUNK,materials_file)
call read_value_string(IIN,IGNORE_JUNK,free_surface_file)
call read_value_string(IIN,IGNORE_JUNK,absorbing_surface_file)
- call read_value_string(IIN,IGNORE_JUNK,tangential_detection_curve_file)
! read info about partitioning
call read_value_integer(IIN,IGNORE_JUNK,nproc)
@@ -426,6 +422,7 @@
call read_value_double_precision(IIN,IGNORE_JUNK,Mzz)
call read_value_double_precision(IIN,IGNORE_JUNK,Mxz)
call read_value_double_precision(IIN,IGNORE_JUNK,factor)
+ call read_value_logical(IIN,IGNORE_JUNK,force_normal_to_surface)
! if Dirac source time function, use a very thin Gaussian instead
! if Heaviside source time function, use a very thin error function instead
@@ -466,7 +463,7 @@
call read_value_logical(IIN,IGNORE_JUNK,generate_STATIONS)
call read_value_integer(IIN,IGNORE_JUNK,nreceiverlines)
call read_value_double_precision(IIN,IGNORE_JUNK,anglerec)
- call read_value_logical(IIN,IGNORE_JUNK,tangential_detection)
+ call read_value_logical(IIN,IGNORE_JUNK,rotate_seismograms_RT)
if(nreceiverlines < 1) stop 'number of receiver lines must be greater than 1'
@@ -565,20 +562,6 @@
print *
enddo
-! tangential detection
- if (tangential_detection) then
- open(unit=IIN,file=tangential_detection_curve_file,status='old',action='read')
- read(IIN,*) nnodes_tangential_curve
- allocate(nodes_tangential_curve(2,nnodes_tangential_curve))
- do i = 1, nnodes_tangential_curve
- read(IIN,*) nodes_tangential_curve(1,i), nodes_tangential_curve(2,i)
- enddo
- close(IIN)
- else
- nnodes_tangential_curve = 0
- allocate(nodes_tangential_curve(2,1))
- endif
-
if ( read_external_mesh ) then
call read_mat(materials_file, nelmnts, num_material)
else
@@ -1131,8 +1114,8 @@
write(15,*) 'meshvect modelvect boundvect cutsnaps subsamp sizemax_arrows'
write(15,*) meshvect,modelvect,boundvect,cutsnaps,subsamp,sizemax_arrows
- write(15,*) 'anglerec'
- write(15,*) anglerec
+ write(15,*) 'anglerec force_normal_to_surface rotate_seismograms_RT'
+ write(15,*) anglerec,force_normal_to_surface,rotate_seismograms_RT
write(15,*) 'initialfield add_Bielak_conditions'
write(15,*) initialfield,add_Bielak_conditions
@@ -1181,8 +1164,8 @@
call write_fluidsolid_edges_database(15, nedges_coupled, nedges_coupled_loc, &
edges_coupled, glob2loc_elmnts, part, iproc, 1)
- write(15,*) 'nelemabs nelem_acoustic_surface num_fluid_solid_edges nnodes_tangential_curve'
- write(15,*) nelemabs_loc,nelem_acoustic_surface_loc,nedges_coupled_loc, nnodes_tangential_curve
+ write(15,*) 'nelemabs nelem_acoustic_surface num_fluid_solid_edges'
+ write(15,*) nelemabs_loc,nelem_acoustic_surface_loc,nedges_coupled_loc
write(15,*) 'Material sets (num 1 rho vp vs 0 0) or (num 2 rho c11 c13 c33 c44)'
@@ -1234,11 +1217,6 @@
call write_fluidsolid_edges_database(15, nedges_coupled, nedges_coupled_loc, &
edges_coupled, glob2loc_elmnts, part, iproc, 2)
- write(15,*) 'List of tangential detection curve nodes:'
- write(15,*) nnodes_tangential_curve
- do i = 1, nnodes_tangential_curve
- write(15,*) 1000*nodes_tangential_curve(1,i), 1000*nodes_tangential_curve(2,i)
- enddo
end do
Modified: seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/specfem2D.F90 2007-03-26 18:32:04 UTC (rev 10443)
+++ seismo/2D/SPECFEM2D/branches/SEM_tangente_branch/specfem2D.F90 2008-02-15 09:48:21 UTC (rev 10444)
@@ -361,15 +361,12 @@
! tangential detection
double precision, dimension(:), allocatable :: anglerec_irec
double precision, dimension(:), allocatable :: cosrot_irec, sinrot_irec
- double precision, dimension(:), allocatable :: x_final_receiver, z_final_receiver
- logical :: tangential_detection
- integer :: nnodes_tangential_curve, source_courbe_eros
- double precision, dimension(:,:), allocatable :: nodes_tangential_curve
- integer :: n1_tangential_detection_curve
- integer, dimension(4) :: n_tangential_detection_curve
- integer, dimension(:), allocatable :: rec_tangential_detection_curve
+ integer, dimension(:), allocatable :: recloc_ispec
+ integer, dimension(2,4) :: recloc_iedge_elect
+ logical :: force_normal_to_surface, rotate_seismograms_RT
+
double precision :: distmin, dist_current, angleforce_recv
- double precision, dimension(:), allocatable :: dist_tangential_detection_curve
+
double precision :: x_final_receiver_dummy, z_final_receiver_dummy
@@ -451,7 +448,7 @@
cutsnaps = cutsnaps / 100.d0
read(IIN,"(a80)") datlin
- read(IIN,*) anglerec
+ read(IIN,*) anglerec,force_normal_to_surface,rotate_seismograms_RT
read(IIN,"(a80)") datlin
read(IIN,*) initialfield,add_Bielak_conditions
@@ -549,7 +546,7 @@
read(IIN,"(a80)") datlin
read(IIN,*) numat,ngnod,nspec,pointsdisp,plot_lowerleft_corner_only
read(IIN,"(a80)") datlin
- read(IIN,*) nelemabs,nelem_acoustic_surface,num_fluid_solid_edges,nnodes_tangential_curve
+ read(IIN,*) nelemabs,nelem_acoustic_surface,num_fluid_solid_edges
!
!---- allocate arrays
@@ -751,21 +748,6 @@
end if
!
-!---- read tangential detection curve
-!
- read(IIN,"(a80)") datlin
- if (nnodes_tangential_curve > 0) then
- tangential_detection = .true.
- allocate(nodes_tangential_curve(2,nnodes_tangential_curve))
- do i = 1, nnodes_tangential_curve
- read(IIN,*) nodes_tangential_curve(1,i),nodes_tangential_curve(2,i)
- enddo
- else
- tangential_detection = .false.
- nnodes_tangential_curve = 0
- allocate(nodes_tangential_curve(2,1))
- endif
-!
!---- close input file
!
close(IIN)
@@ -1043,190 +1025,15 @@
call exit_MPI('incorrect source type')
endif
- allocate(x_final_receiver(nrec))
- allocate(z_final_receiver(nrec))
+ allocate(recloc_ispec(nrec))
! locate receivers in the mesh
call locate_receivers(ibool,coord,nspec,npoin,xigll,zigll,nrec,nrecloc,recloc,which_proc_receiver,nproc,myrank,&
st_xval,st_zval,ispec_selected_rec, &
xi_receiver,gamma_receiver,station_name,network_name,x_source,z_source,coorg,knods,ngnod,npgeo, &
- x_final_receiver, z_final_receiver)
+ recloc_ispec)
- allocate(anglerec_irec(nrecloc))
- anglerec_irec(:) = 0
- allocate(cosrot_irec(nrecloc))
- allocate(sinrot_irec(nrecloc))
- cosrot_irec(:) = cos(anglerec)
- sinrot_irec(:) = sin(anglerec)
- allocate(rec_tangential_detection_curve(nrecloc))
- if (tangential_detection) then
- irecloc = 0
- do irec = 1, nrec
- if (which_proc_receiver(irec) == myrank) then
- irecloc = irecloc + 1
- distmin = HUGEVAL
- do i = 1, nnodes_tangential_curve
- dist_current = sqrt((x_final_receiver(irec)-nodes_tangential_curve(1,i))**2 + &
- (z_final_receiver(irec)-nodes_tangential_curve(2,i))**2)
- if ( dist_current < distmin ) then
- n1_tangential_detection_curve = i
- distmin = dist_current
- endif
- enddo
-
- rec_tangential_detection_curve(irecloc) = n1_tangential_detection_curve
- call tri_quad(n_tangential_detection_curve, n1_tangential_detection_curve, nnodes_tangential_curve)
-
- call calcul_normale( anglerec_irec(irecloc), nodes_tangential_curve(1,n_tangential_detection_curve(1)), &
- nodes_tangential_curve(1,n_tangential_detection_curve(2)), &
- nodes_tangential_curve(1,n_tangential_detection_curve(3)), &
- nodes_tangential_curve(1,n_tangential_detection_curve(4)), &
- nodes_tangential_curve(2,n_tangential_detection_curve(1)), &
- nodes_tangential_curve(2,n_tangential_detection_curve(2)), &
- nodes_tangential_curve(2,n_tangential_detection_curve(3)), &
- nodes_tangential_curve(2,n_tangential_detection_curve(4)) )
- endif
-
- enddo
- cosrot_irec(:) = cos(anglerec_irec(:))
- sinrot_irec(:) = sin(anglerec_irec(:))
-
- distmin = HUGEVAL
- do i = 1, nnodes_tangential_curve
- dist_current = sqrt((coord(1,iglob_source)-nodes_tangential_curve(1,i))**2 + &
- (coord(2,iglob_source)-nodes_tangential_curve(2,i))**2)
- if ( dist_current < distmin ) then
- n1_tangential_detection_curve = i
- distmin = dist_current
-
- endif
- enddo
-
- call tri_quad(n_tangential_detection_curve, n1_tangential_detection_curve, nnodes_tangential_curve)
-
- call calcul_normale( angleforce, nodes_tangential_curve(1,n_tangential_detection_curve(1)), &
- nodes_tangential_curve(1,n_tangential_detection_curve(2)), &
- nodes_tangential_curve(1,n_tangential_detection_curve(3)), &
- nodes_tangential_curve(1,n_tangential_detection_curve(4)), &
- nodes_tangential_curve(2,n_tangential_detection_curve(1)), &
- nodes_tangential_curve(2,n_tangential_detection_curve(2)), &
- nodes_tangential_curve(2,n_tangential_detection_curve(3)), &
- nodes_tangential_curve(2,n_tangential_detection_curve(4)) )
-
- source_courbe_eros = n1_tangential_detection_curve
- if ( myrank == 0 .and. is_proc_source == 1 .and. nb_proc_source == 1 ) then
- source_courbe_eros = n1_tangential_detection_curve
- angleforce_recv = angleforce
-#ifdef USE_MPI
- else if ( myrank == 0 ) then
- do i = 1, nb_proc_source - is_proc_source
- call MPI_recv(source_courbe_eros,1,MPI_INTEGER,MPI_ANY_SOURCE,42,MPI_COMM_WORLD,request_mpi_status,ier)
- call MPI_recv(angleforce_recv,1,MPI_DOUBLE_PRECISION,MPI_ANY_SOURCE,43,MPI_COMM_WORLD,request_mpi_status,ier)
- end do
- else if ( is_proc_source == 1 ) then
- call MPI_send(n1_tangential_detection_curve,1,MPI_INTEGER,0,42,MPI_COMM_WORLD,ier)
- call MPI_send(angleforce,1,MPI_DOUBLE_PRECISION,0,43,MPI_COMM_WORLD,ier)
-#endif
- end if
-
-#ifdef USE_MPI
- call MPI_bcast(angleforce_recv,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
- angleforce = angleforce_recv
-#endif
- print *, 'POYOP', angleforce, angleforce_recv
- allocate(dist_tangential_detection_curve(nnodes_tangential_curve))
- dist_tangential_detection_curve(source_courbe_eros) = 0
- do i = source_courbe_eros+1, nnodes_tangential_curve
- dist_tangential_detection_curve(i) = dist_tangential_detection_curve(i-1) + &
- sqrt((nodes_tangential_curve(1,i)-nodes_tangential_curve(1,i-1))**2 + &
- (nodes_tangential_curve(2,i)-nodes_tangential_curve(2,i-1))**2)
- enddo
- dist_tangential_detection_curve(1) = dist_tangential_detection_curve(nnodes_tangential_curve) + &
- sqrt((nodes_tangential_curve(1,1)-nodes_tangential_curve(1,nnodes_tangential_curve))**2 + &
- (nodes_tangential_curve(2,1)-nodes_tangential_curve(2,nnodes_tangential_curve))**2)
- do i = 2, source_courbe_eros-1
- dist_tangential_detection_curve(i) = dist_tangential_detection_curve(i-1) + &
- sqrt((nodes_tangential_curve(1,i)-nodes_tangential_curve(1,i-1))**2 + &
- (nodes_tangential_curve(2,i)-nodes_tangential_curve(2,i-1))**2)
- enddo
- do i = source_courbe_eros-1, 1, -1
- dist_current = dist_tangential_detection_curve(i+1) + &
- sqrt((nodes_tangential_curve(1,i)-nodes_tangential_curve(1,i+1))**2 + &
- (nodes_tangential_curve(2,i)-nodes_tangential_curve(2,i+1))**2)
- if ( dist_current < dist_tangential_detection_curve(i) ) then
- dist_tangential_detection_curve(i) = dist_current
- endif
- enddo
- dist_current = dist_tangential_detection_curve(1) + &
- sqrt((nodes_tangential_curve(1,1)-nodes_tangential_curve(1,nnodes_tangential_curve))**2 + &
- (nodes_tangential_curve(2,1)-nodes_tangential_curve(2,nnodes_tangential_curve))**2)
- if ( dist_current < dist_tangential_detection_curve(nnodes_tangential_curve) ) then
- dist_tangential_detection_curve(nnodes_tangential_curve) = dist_current
- endif
- do i = nnodes_tangential_curve-1, source_courbe_eros+1, -1
- dist_current = dist_tangential_detection_curve(i+1) + &
- sqrt((nodes_tangential_curve(1,i)-nodes_tangential_curve(1,i+1))**2 + &
- (nodes_tangential_curve(2,i)-nodes_tangential_curve(2,i+1))**2)
- if ( dist_current < dist_tangential_detection_curve(i) ) then
- dist_tangential_detection_curve(i) = dist_current
- endif
- enddo
-
- if ( myrank == 0 ) then
- open(unit=11,file='OUTPUT_FILES/dist_rec_tangential_detection_curve', form='formatted', status='unknown')
- end if
- irecloc = 0
- do irec = 1,nrec
-
- if ( myrank == 0 ) then
- if ( which_proc_receiver(irec) == myrank ) then
- irecloc = irecloc + 1
- n1_tangential_detection_curve = rec_tangential_detection_curve(irecloc)
- x_final_receiver_dummy = x_final_receiver(irec)
- z_final_receiver_dummy = z_final_receiver(irec)
-#ifdef USE_MPI
- else
-
- call MPI_RECV(n1_tangential_detection_curve,1,MPI_INTEGER,&
- which_proc_receiver(irec),irec,MPI_COMM_WORLD,request_mpi_status,ier)
- call MPI_RECV(x_final_receiver_dummy,1,MPI_DOUBLE_PRECISION,&
- which_proc_receiver(irec),irec,MPI_COMM_WORLD,request_mpi_status,ier)
- call MPI_RECV(z_final_receiver_dummy,1,MPI_DOUBLE_PRECISION,&
- which_proc_receiver(irec),irec,MPI_COMM_WORLD,request_mpi_status,ier)
-
-#endif
- end if
-
-#ifdef USE_MPI
- else
- if ( which_proc_receiver(irec) == myrank ) then
- irecloc = irecloc + 1
- call MPI_SEND(rec_tangential_detection_curve(irecloc),1,MPI_INTEGER,0,irec,MPI_COMM_WORLD,ier)
- call MPI_SEND(x_final_receiver(irec),1,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ier)
- call MPI_SEND(z_final_receiver(irec),1,MPI_DOUBLE_PRECISION,0,irec,MPI_COMM_WORLD,ier)
- end if
-#endif
-
- end if
- if ( myrank == 0 ) then
- write(11,*) dist_tangential_detection_curve(n1_tangential_detection_curve)
- write(12,*) x_final_receiver_dummy
- write(13,*) z_final_receiver_dummy
- end if
-
- enddo
- if ( myrank == 0 ) then
- close(11)
- close(12)
- close(13)
- end if
-
- else
-
-
- endif
-
! allocate seismogram arrays
allocate(sisux(NTSTEP_BETWEEN_OUTPUT_SEISMO,nrecloc))
allocate(sisuz(NTSTEP_BETWEEN_OUTPUT_SEISMO,nrecloc))
@@ -1448,6 +1255,229 @@
#endif
+
+ allocate(anglerec_irec(nrecloc))
+ anglerec_irec(:) = anglerec
+ allocate(cosrot_irec(nrecloc))
+ allocate(sinrot_irec(nrecloc))
+ cosrot_irec(:) = cos(anglerec)
+ sinrot_irec(:) = sin(anglerec)
+
+ if (force_normal_to_surface) then
+ if (is_proc_source == 1) then
+ ispec = ispec_selected_source
+ recloc_iedge_elect(1,IBOTTOM) = ibool(3,1,ispec)
+ recloc_iedge_elect(1,IRIGHT) = ibool(5,3,ispec)
+ recloc_iedge_elect(1,ITOP) = ibool(3,5,ispec)
+ recloc_iedge_elect(1,ILEFT) = ibool(1,3,ispec)
+
+ recloc_iedge_elect(2,1:4) = 0
+
+ do ispec = 1, nspec
+ do j = 1, NGLLZ
+ do i = 1, NGLLX
+ iglob = ibool(i,j,ispec)
+ do inum = 1, 4
+ if (recloc_iedge_elect(1,inum) == iglob) then
+ recloc_iedge_elect(2,inum) = recloc_iedge_elect(2,inum)+1
+ endif
+ enddo
+#ifdef USE_MPI
+ do num_interface = 1, ninterface_elastic
+ do ipoin = 1, nibool_interfaces_elastic(num_interface)
+ do inum = 1, 4
+ if (recloc_iedge_elect(1,inum) == ibool_interfaces_elastic(ipoin,num_interface)) then
+ recloc_iedge_elect(2,inum) = recloc_iedge_elect(2,inum)+1
+ endif
+ enddo
+ enddo
+ enddo
+ do num_interface = 1, ninterface_acoustic
+ do ipoin = 1, nibool_interfaces_acoustic(num_interface)
+ do inum = 1, 4
+ if (recloc_iedge_elect(1,inum) == ibool_interfaces_acoustic(ipoin,num_interface)) then
+ recloc_iedge_elect(2,inum) = recloc_iedge_elect(2,inum)+1
+ endif
+ enddo
+ enddo
+ enddo
+#endif
+ enddo
+ enddo
+ enddo
+
+ nx = 0
+ nz = 1
+
+ if (recloc_iedge_elect(2,IBOTTOM) == 1) then
+ i = 3
+ j = 1
+
+ xxi = + gammaz(i,j,ispec_selected_source) * jacobian(i,j,ispec_selected_source)
+ zxi = - gammax(i,j,ispec_selected_source) * jacobian(i,j,ispec_selected_source)
+ jacobian1D = sqrt(xxi**2 + zxi**2)
+ nx = + zxi / jacobian1D
+ nz = - xxi / jacobian1D
+ endif
+
+ if (recloc_iedge_elect(2,IRIGHT) == 1) then
+ i = 5
+ j = 3
+
+ xgamma = - xiz(i,j,ispec_selected_source) * jacobian(i,j,ispec_selected_source)
+ zgamma = + xix(i,j,ispec_selected_source) * jacobian(i,j,ispec_selected_source)
+ jacobian1D = sqrt(xgamma**2 + zgamma**2)
+ nx = + zgamma / jacobian1D
+ nz = - xgamma / jacobian1D
+ endif
+
+ if (recloc_iedge_elect(2,ITOP) == 1) then
+ i = 3
+ j = 5
+
+ xxi = + gammaz(i,j,ispec_selected_source) * jacobian(i,j,ispec_selected_source)
+ zxi = - gammax(i,j,ispec_selected_source) * jacobian(i,j,ispec_selected_source)
+ jacobian1D = sqrt(xxi**2 + zxi**2)
+ nx = + zxi / jacobian1D
+ nz = - xxi / jacobian1D
+ endif
+
+ if (recloc_iedge_elect(2,ILEFT) == 1) then
+ i = 1
+ j = 3
+
+ xgamma = - xiz(i,j,ispec_selected_source) * jacobian(i,j,ispec_selected_source)
+ zgamma = + xix(i,j,ispec_selected_source) * jacobian(i,j,ispec_selected_source)
+ jacobian1D = sqrt(xgamma**2 + zgamma**2)
+ nx = + zgamma / jacobian1D
+ nz = - xgamma / jacobian1D
+ endif
+
+ angleforce = acos(nz)
+ angleforce = -sign(angleforce,nx)
+ angleforce = pi + angleforce
+
+ endif
+ endif
+
+ if (rotate_seismograms_RT) then
+ do irec = 1, nrecloc
+ ispec = recloc_ispec(irec)
+ recloc_iedge_elect(1,IBOTTOM) = ibool(3,1,ispec)
+ recloc_iedge_elect(1,IRIGHT) = ibool(5,3,ispec)
+ recloc_iedge_elect(1,ITOP) = ibool(3,5,ispec)
+ recloc_iedge_elect(1,ILEFT) = ibool(1,3,ispec)
+
+ recloc_iedge_elect(2,1:4) = 0
+
+ do ispec = 1, nspec
+ do j = 1, NGLLZ
+ do i = 1, NGLLX
+ iglob = ibool(i,j,ispec)
+ do inum = 1, 4
+ if (recloc_iedge_elect(1,inum) == iglob) then
+ recloc_iedge_elect(2,inum) = recloc_iedge_elect(2,inum)+1
+ endif
+ enddo
+#ifdef USE_MPI
+ do num_interface = 1, ninterface_elastic
+ do ipoin = 1, nibool_interfaces_elastic(num_interface)
+ do inum = 1, 4
+ if (recloc_iedge_elect(1,inum) == ibool_interfaces_elastic(ipoin,num_interface)) then
+ recloc_iedge_elect(2,inum) = recloc_iedge_elect(2,inum)+1
+ endif
+ enddo
+ enddo
+ enddo
+ do num_interface = 1, ninterface_acoustic
+ do ipoin = 1, nibool_interfaces_acoustic(num_interface)
+ do inum = 1, 4
+ if (recloc_iedge_elect(1,inum) == ibool_interfaces_acoustic(ipoin,num_interface)) then
+ recloc_iedge_elect(2,inum) = recloc_iedge_elect(2,inum)+1
+ endif
+ enddo
+ enddo
+ enddo
+#endif
+ enddo
+ enddo
+ enddo
+
+ nx = 0
+ nz = 1
+
+ if (recloc_iedge_elect(2,IBOTTOM) == 1) then
+ i = 3
+ j = 1
+
+ xxi = + gammaz(i,j,recloc_ispec(irec)) * jacobian(i,j,recloc_ispec(irec))
+ zxi = - gammax(i,j,recloc_ispec(irec)) * jacobian(i,j,recloc_ispec(irec))
+ jacobian1D = sqrt(xxi**2 + zxi**2)
+ nx = + zxi / jacobian1D
+ nz = - xxi / jacobian1D
+ endif
+
+ if (recloc_iedge_elect(2,IRIGHT) == 1) then
+ i = 5
+ j = 3
+
+ xgamma = - xiz(i,j,recloc_ispec(irec)) * jacobian(i,j,recloc_ispec(irec))
+ zgamma = + xix(i,j,recloc_ispec(irec)) * jacobian(i,j,recloc_ispec(irec))
+ jacobian1D = sqrt(xgamma**2 + zgamma**2)
+ nx = + zgamma / jacobian1D
+ nz = - xgamma / jacobian1D
+ endif
+
+ if (recloc_iedge_elect(2,ITOP) == 1) then
+ i = 3
+ j = 5
+
+ xxi = + gammaz(i,j,recloc_ispec(irec)) * jacobian(i,j,recloc_ispec(irec))
+ zxi = - gammax(i,j,recloc_ispec(irec)) * jacobian(i,j,recloc_ispec(irec))
+ jacobian1D = sqrt(xxi**2 + zxi**2)
+ nx = + zxi / jacobian1D
+ nz = - xxi / jacobian1D
+ endif
+
+ if (recloc_iedge_elect(2,ILEFT) == 1) then
+ i = 1
+ j = 3
+
+ xgamma = - xiz(i,j,recloc_ispec(irec)) * jacobian(i,j,recloc_ispec(irec))
+ zgamma = + xix(i,j,recloc_ispec(irec)) * jacobian(i,j,recloc_ispec(irec))
+ jacobian1D = sqrt(xgamma**2 + zgamma**2)
+ nx = + zgamma / jacobian1D
+ nz = - xgamma / jacobian1D
+ endif
+
+ anglerec_irec(irec) = acos(nz)
+ anglerec_irec(irec) = -sign(anglerec_irec(irec),nx)
+ anglerec_irec(irec) = pi+anglerec_irec(irec)
+ cosrot_irec(irec) = cos(anglerec_irec(irec))
+ sinrot_irec(irec) = sin(anglerec_irec(irec))
+
+ enddo
+ endif
+
+
+ if ( myrank == 0 .and. is_proc_source == 1 .and. nb_proc_source == 1 ) then
+ angleforce_recv = angleforce
+#ifdef USE_MPI
+ else if ( myrank == 0 ) then
+ do i = 1, nb_proc_source - is_proc_source
+ call MPI_recv(angleforce_recv,1,MPI_DOUBLE_PRECISION,MPI_ANY_SOURCE,43,MPI_COMM_WORLD,request_mpi_status,ier)
+ end do
+ else if ( is_proc_source == 1 ) then
+ call MPI_send(angleforce,1,MPI_DOUBLE_PRECISION,0,43,MPI_COMM_WORLD,ier)
+#endif
+ end if
+
+#ifdef USE_MPI
+ call MPI_bcast(angleforce_recv,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ angleforce = angleforce_recv
+#endif
+
+
! fill mass matrix with fictitious non-zero values to make sure it can be inverted globally
if(any_elastic) where(rmass_inverse_elastic <= 0._CUSTOM_REAL) rmass_inverse_elastic = 1._CUSTOM_REAL
if(any_acoustic) where(rmass_inverse_acoustic <= 0._CUSTOM_REAL) rmass_inverse_acoustic = 1._CUSTOM_REAL
@@ -3080,72 +3110,3 @@
end subroutine is_in_convex_quadrilateral
-subroutine tri_quad(n, n1, nnodes)
-
- implicit none
-
- integer :: n1, n2, n3, n4, nnodes
- integer, dimension(4) :: n
-
-
- n(2) = n1
-
- if ( n1 == 1 ) then
- n(1) = nnodes
- else
- n(1) = n1-1
- endif
-
- if ( n1 == nnodes ) then
- n(3) = 1
- else
- n(3) = n1+1
- endif
-
- if ( n(3) == nnodes ) then
- n(4) = 1
- else
- n(4) = n(3)+1
- endif
-
-
-end subroutine tri_quad
-
-
-subroutine calcul_normale( angle, n1_x, n2_x, n3_x, n4_x, n1_z, n2_z, n3_z, n4_z )
-
- implicit none
-
- include 'constants.h'
-
- double precision :: angle
- double precision :: n1_x, n2_x, n3_x, n4_x, n1_z, n2_z, n3_z, n4_z
-
- double precision :: theta1, theta2, theta3
- double precision :: costheta1, costheta2, costheta3
-
- if ( abs(n2_z - n1_z) < TINYVAL ) then
- costheta1 = 0
- else
- costheta1 = (n2_z - n1_z) / sqrt((n2_x - n1_x)**2 + (n2_z - n1_z)**2)
- endif
- if ( abs(n3_z - n2_z) < TINYVAL ) then
- costheta2 = 0
- else
- costheta2 = (n3_z - n2_z) / sqrt((n3_x - n2_x)**2 + (n3_z - n2_z)**2)
- endif
- if ( abs(n4_z - n3_z) < TINYVAL ) then
- costheta3 = 0
- else
- costheta3 = (n4_z - n3_z) / sqrt((n4_x - n3_x)**2 + (n4_z - n3_z)**2)
- endif
-
- theta1 = - sign(1.d0,n2_x - n1_x) * acos(costheta1)
- theta2 = - sign(1.d0,n3_x - n2_x) * acos(costheta2)
- theta3 = - sign(1.d0,n4_x - n3_x) * acos(costheta3)
-
-
- angle = ( theta1 + theta2 + theta3 ) / 3.d0 + PI/2.d0
- !angle = theta2
-
-end subroutine calcul_normale
More information about the cig-commits
mailing list