From f3c26a36cdaf3db827faf44aaee857594439c5f0 Mon Sep 17 00:00:00 2001 From: Daniel Peter Date: Wed, 3 Jul 2019 00:29:22 +0200 Subject: [PATCH] adds myrank to constants module; comments out unused routines; adds initial openmp statements --- Makefile.in | 5 +- src/auxiliaries/combine_vol_data.F90 | 4 +- .../combine_vol_data_vtk_binary.F90 | 2 +- .../create_movie_shakemap_AVS_DX_GMT.f90 | 4 +- ...t_and_combine_vol_data_on_regular_grid.f90 | 2 +- src/decompose_mesh/program_decompose_mesh.f90 | 2 +- .../program_decompose_mesh_mpi.f90 | 2 +- .../create_regions_mesh.f90 | 46 +- .../fault_generate_databases.f90 | 27 +- .../generate_databases_par.F90 | 21 +- .../get_absorbing_boundary.F90 | 60 +- .../get_coupling_surfaces.f90 | 30 +- src/generate_databases/get_model.F90 | 6 +- src/generate_databases/get_perm_color.f90 | 48 +- src/generate_databases/memory_eval.f90 | 12 +- src/generate_databases/model_coupled.f90 | 11 +- .../model_external_values.F90 | 4 +- .../model_salton_trough.f90 | 6 +- src/generate_databases/read_parameters.f90 | 2 +- src/generate_databases/setup_color_perm.f90 | 3 +- src/generate_databases/setup_mesh.f90 | 10 +- src/inverse_problem_for_model/rules.mk | 1 + .../specfem_interface_mod.F90 | 4 +- src/meshfem3D/check_mesh_quality.f90 | 4 +- src/meshfem3D/create_CPML_regions.f90 | 4 +- src/meshfem3D/create_interfaces_mesh.f90 | 4 +- src/meshfem3D/create_meshfem_mesh.f90 | 36 +- src/meshfem3D/define_subregions.f90 | 3 +- src/meshfem3D/define_subregions_heuristic.f90 | 19 +- src/meshfem3D/determine_cavity.f90 | 4 +- src/meshfem3D/earth_chunk.f90 | 6 +- src/meshfem3D/get_MPI_cutplanes_eta.f90 | 300 +++---- src/meshfem3D/get_MPI_cutplanes_xi.f90 | 291 +++---- src/meshfem3D/get_flags_boundaries.f90 | 1 + src/meshfem3D/meshfem3D.F90 | 2 +- src/meshfem3D/meshfem3D_adios_stubs.f90 | 22 +- src/meshfem3D/meshfem3D_par.f90 | 2 +- src/meshfem3D/save_databases.F90 | 3 +- src/meshfem3D/save_databases_adios.F90 | 7 +- src/meshfem3D/store_boundaries.f90 | 610 +++++++------- src/shared/check_mesh_resolution.f90 | 20 +- src/shared/detect_surface.f90 | 22 +- src/shared/exit_mpi.f90 | 6 +- src/shared/get_attenuation_model.f90 | 16 +- src/shared/get_jacobian_boundaries.f90 | 66 +- src/shared/get_shape2D.f90 | 20 +- src/shared/get_shape3D.f90 | 26 +- src/shared/read_parameter_file.F90 | 8 +- src/shared/rules.mk | 1 + src/shared/serial.f90 | 2 +- src/shared/shared_par.F90 | 3 + src/specfem3D/assemble_MPI_vector.f90 | 25 +- .../compute_add_sources_viscoelastic.F90 | 24 +- src/specfem3D/compute_forces_viscoelastic.f90 | 5 +- ...te_forces_viscoelastic_calling_routine.F90 | 52 +- src/specfem3D/compute_kernels.f90 | 760 +++++++++--------- src/specfem3D/compute_stacey_viscoelastic.F90 | 243 +++--- src/specfem3D/detect_mesh_surfaces.f90 | 20 +- src/specfem3D/initialize_simulation.F90 | 15 +- src/specfem3D/noise_tomography.f90 | 72 +- src/specfem3D/prepare_noise.f90 | 7 +- src/specfem3D/rules.mk | 1 + src/specfem3D/save_adjoint_kernels.f90 | 4 +- src/specfem3D/setup_GLL_points.f90 | 6 +- src/specfem3D/setup_sources_receivers.f90 | 4 +- src/specfem3D/specfem3D_par.F90 | 11 +- src/specfem3D/update_displacement_scheme.f90 | 76 +- src/specfem3D/vtk_window.F90 | 4 +- src/specfem3D/write_movie_output.F90 | 12 +- src/tomography/add_model_iso.f90 | 2 +- src/tomography/model_update.f90 | 12 +- .../clip_sem.f90 | 2 +- .../combine_sem.f90 | 2 +- .../smooth_sem.F90 | 8 +- src/tomography/rules.mk | 1 + src/tomography/sum_kernels.f90 | 2 +- src/tomography/sum_preconditioned_kernels.f90 | 2 +- tests/decompose_mesh/test_read.f90 | 2 +- 78 files changed, 1614 insertions(+), 1580 deletions(-) diff --git a/Makefile.in b/Makefile.in index 9db5a0d9f..cd7dd0ba5 100644 --- a/Makefile.in +++ b/Makefile.in @@ -225,10 +225,7 @@ GENCODE_70 = -gencode=arch=compute_70,code=\"sm_70,compute_70\" @COND_OMP_TRUE@OPENMP = yes @COND_OMP_FALSE@OPENMP = no -@COND_OMP_TRUE@OMP_FLAGS = @OMP_FCFLAGS@ $(FC_DEFINE)USE_OPENMP -@COND_OMP_FALSE@OMP_FLAGS = - -FCFLAGS += $(OMP_FLAGS) +@COND_OMP_TRUE@FCFLAGS += $(FC_DEFINE)USE_OPENMP @OMP_FCFLAGS@ @COND_OMP_TRUE@OMP_LIBS = $(OMP_LIB) @COND_OMP_FALSE@OMP_LIBS = diff --git a/src/auxiliaries/combine_vol_data.F90 b/src/auxiliaries/combine_vol_data.F90 index aa9b96705..d51e7b13c 100644 --- a/src/auxiliaries/combine_vol_data.F90 +++ b/src/auxiliaries/combine_vol_data.F90 @@ -72,7 +72,7 @@ program combine_vol_data character(len=MAX_STRING_LEN*2) :: mesh_file,local_data_file character(len=MAX_STRING_LEN) :: data_array_name logical :: HIGH_RESOLUTION_MESH,BROADCAST_AFTER_READ - integer :: ires,myrank + integer :: ires ! Variables to read ADIOS files integer :: sizeprocs @@ -96,7 +96,7 @@ program combine_vol_data ! needs local_path for mesh files myrank = 0 BROADCAST_AFTER_READ = .false. - call read_parameter_file(myrank,BROADCAST_AFTER_READ) + call read_parameter_file(BROADCAST_AFTER_READ) ! reads in arguments do i = 1, command_argument_count() diff --git a/src/auxiliaries/combine_vol_data_vtk_binary.F90 b/src/auxiliaries/combine_vol_data_vtk_binary.F90 index cb0ef13cd..2b7ab981e 100644 --- a/src/auxiliaries/combine_vol_data_vtk_binary.F90 +++ b/src/auxiliaries/combine_vol_data_vtk_binary.F90 @@ -89,7 +89,7 @@ program combine_vol_data ! needs local_path for mesh files myrank = 0 BROADCAST_AFTER_READ = .false. - call read_parameter_file(myrank,BROADCAST_AFTER_READ) + call read_parameter_file(BROADCAST_AFTER_READ) ! reads in arguments do i = 1, command_argument_count() diff --git a/src/auxiliaries/create_movie_shakemap_AVS_DX_GMT.f90 b/src/auxiliaries/create_movie_shakemap_AVS_DX_GMT.f90 index 1c3ea880f..4bc7937fb 100644 --- a/src/auxiliaries/create_movie_shakemap_AVS_DX_GMT.f90 +++ b/src/auxiliaries/create_movie_shakemap_AVS_DX_GMT.f90 @@ -105,7 +105,7 @@ program create_movie_shakemap ! order of points representing the 2D square element integer,dimension(NGNOD2D_FOUR_CORNERS_AVS_DX),parameter :: iorder = (/1,3,2,4/) - integer :: NSPEC_SURFACE_EXT_MESH,myrank + integer :: NSPEC_SURFACE_EXT_MESH logical :: BROADCAST_AFTER_READ ! ************** PROGRAM STARTS HERE ************** @@ -123,7 +123,7 @@ program create_movie_shakemap BROADCAST_AFTER_READ = .false. ! read the parameter file - call read_parameter_file(myrank,BROADCAST_AFTER_READ) + call read_parameter_file(BROADCAST_AFTER_READ) ! only one global array for movie data, but stored for all surfaces defined ! in file 'surface_from_mesher.h' diff --git a/src/auxiliaries/project_and_combine_vol_data_on_regular_grid.f90 b/src/auxiliaries/project_and_combine_vol_data_on_regular_grid.f90 index 92d4c40bf..b8dd44fd3 100644 --- a/src/auxiliaries/project_and_combine_vol_data_on_regular_grid.f90 +++ b/src/auxiliaries/project_and_combine_vol_data_on_regular_grid.f90 @@ -70,7 +70,7 @@ program project_and_combine_vol_data_on_regular_grid ! needs local_path for mesh files BROADCAST_AFTER_READ = .true. - call read_parameter_file(myrank,BROADCAST_AFTER_READ) + call read_parameter_file(BROADCAST_AFTER_READ) ! parse command line arguments if (command_argument_count() /= NARGS) then diff --git a/src/decompose_mesh/program_decompose_mesh.f90 b/src/decompose_mesh/program_decompose_mesh.f90 index c943e9e62..8295c76ad 100644 --- a/src/decompose_mesh/program_decompose_mesh.f90 +++ b/src/decompose_mesh/program_decompose_mesh.f90 @@ -66,7 +66,7 @@ program xdecompose_mesh ! needs local_path for mesh files myrank = 0 BROADCAST_AFTER_READ = .false. - call read_parameter_file(myrank,BROADCAST_AFTER_READ) + call read_parameter_file(BROADCAST_AFTER_READ) ! checks adios parameters if (ADIOS_FOR_DATABASES) then diff --git a/src/decompose_mesh/program_decompose_mesh_mpi.f90 b/src/decompose_mesh/program_decompose_mesh_mpi.f90 index fa7e4b2f5..61804d23c 100644 --- a/src/decompose_mesh/program_decompose_mesh_mpi.f90 +++ b/src/decompose_mesh/program_decompose_mesh_mpi.f90 @@ -58,7 +58,7 @@ program xdecompose_mesh_mpi call world_rank(myrank) ! 1/ read Par_file - call read_parameter_file(myrank, BROADCAST_AFTER_READ) + call read_parameter_file(BROADCAST_AFTER_READ) ! safety check, won't work yet with local-time stepping if (LTS_MODE) stop 'LTS_MODE not supported yet with parallel mesh decomposer!' diff --git a/src/generate_databases/create_regions_mesh.f90 b/src/generate_databases/create_regions_mesh.f90 index 8b9da2f5f..4c67ecf02 100644 --- a/src/generate_databases/create_regions_mesh.f90 +++ b/src/generate_databases/create_regions_mesh.f90 @@ -76,7 +76,7 @@ subroutine create_regions_mesh() nodes_coords_ext_mesh,nnodes_ext_mesh,elmnts_ext_mesh,nelmnts_ext_mesh,ANY_FAULT_IN_THIS_PROC) ! if faults exist this reads nodes_coords_open - call fault_read_input(prname,myrank) + call fault_read_input(prname) ! fills location and weights for Gauss-Lobatto-Legendre points, shape and derivations, ! returns jacobianstore,xixstore,...gammazstore @@ -129,8 +129,8 @@ subroutine create_regions_mesh() endif ! at this point (xyz)store_dummy are still open if (.not. PARALLEL_FAULT) then - call fault_setup (ibool,nnodes_ext_mesh,nodes_coords_ext_mesh, & - xstore,ystore,zstore,nspec,nglob,myrank) + call fault_setup(ibool,nnodes_ext_mesh,nodes_coords_ext_mesh, & + xstore,ystore,zstore,nspec,nglob) endif ! this closes (xyz)store_dummy endif @@ -155,8 +155,8 @@ subroutine create_regions_mesh() if (PARALLEL_FAULT .and. ANY_FAULT) then call synchronize_all() !at this point (xyz)store_dummy are still open - call fault_setup (ibool,nnodes_ext_mesh,nodes_coords_ext_mesh, & - xstore,ystore,zstore,nspec,nglob,myrank) + call fault_setup(ibool,nnodes_ext_mesh,nodes_coords_ext_mesh, & + xstore,ystore,zstore,nspec,nglob) ! this closes (xyz)store_dummy endif @@ -330,13 +330,13 @@ subroutine create_regions_mesh() if (POROELASTIC_SIMULATION) then !chris: check for poro: At the moment cpI & cpII are for eta=0 - call check_mesh_resolution_poro(myrank,nspec,nglob_dummy,ibool, & + call check_mesh_resolution_poro(nspec,nglob_dummy,ibool, & xstore_dummy,ystore_dummy,zstore_dummy, & -1.0d0, model_speed_max,min_resolved_period, & phistore,tortstore,rhoarraystore,rho_vpI,rho_vpII,rho_vsI, & LOCAL_PATH,SAVE_MESH_FILES ) else - call check_mesh_resolution(myrank,nspec,nglob_dummy, & + call check_mesh_resolution(nspec,nglob_dummy, & ibool,xstore_dummy,ystore_dummy,zstore_dummy, & kappastore,mustore,rho_vp,rho_vs, & -1.0d0,model_speed_max,min_resolved_period, & @@ -351,7 +351,7 @@ subroutine create_regions_mesh() write(IMAIN,*) ' ...saving attenuation databases' call flush_IMAIN() endif - call get_attenuation_model(myrank,nspec,USE_OLSEN_ATTENUATION,OLSEN_ATTENUATION_RATIO, & + call get_attenuation_model(nspec,USE_OLSEN_ATTENUATION,OLSEN_ATTENUATION_RATIO, & mustore,rho_vs,kappastore,rho_vp,qkappa_attenuation_store,qmu_attenuation_store, & ispec_is_elastic,min_resolved_period,prname,ATTENUATION_f0_REFERENCE) endif @@ -821,13 +821,13 @@ subroutine crm_ext_setup_jacobian(myrank, & call zwgljd(zigll,wzgll,NGLLZ,GAUSSALPHA,GAUSSBETA) ! get the 3-D shape functions - call get_shape3D(myrank,shape3D,dershape3D,xigll,yigll,zigll,NGNOD) + call get_shape3D(shape3D,dershape3D,xigll,yigll,zigll,NGNOD) ! get the 2-D shape functions - call get_shape2D(myrank,shape2D_x,dershape2D_x,yigll,zigll,NGLLY,NGLLZ,NGNOD,NGNOD2D) - call get_shape2D(myrank,shape2D_y,dershape2D_y,xigll,zigll,NGLLX,NGLLZ,NGNOD,NGNOD2D) - call get_shape2D(myrank,shape2D_bottom,dershape2D_bottom,xigll,yigll,NGLLX,NGLLY,NGNOD,NGNOD2D) - call get_shape2D(myrank,shape2D_top,dershape2D_top,xigll,yigll,NGLLX,NGLLY,NGNOD,NGNOD2D) + call get_shape2D(shape2D_x,dershape2D_x,yigll,zigll,NGLLY,NGLLZ,NGNOD,NGNOD2D) + call get_shape2D(shape2D_y,dershape2D_y,xigll,zigll,NGLLX,NGLLZ,NGNOD,NGNOD2D) + call get_shape2D(shape2D_bottom,dershape2D_bottom,xigll,yigll,NGLLX,NGLLY,NGNOD,NGNOD2D) + call get_shape2D(shape2D_top,dershape2D_top,xigll,yigll,NGLLX,NGLLY,NGNOD,NGNOD2D) ! 2D weights do j=1,NGLLY @@ -1111,11 +1111,11 @@ subroutine crm_setup_moho( myrank,nspec, & call get_element_face_gll_indices(iface,ijk_face,NGLLX,NGLLZ) ! weighted jacobian and normal - call get_jacobian_boundary_face(myrank,nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & - dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & - wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & - ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ,NGNOD2D) + call get_jacobian_boundary_face(nspec, & + xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & + dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & + wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & + ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ,NGNOD2D) ! normal convention: points away from element ! switch normal direction if necessary @@ -1204,11 +1204,11 @@ subroutine crm_setup_moho( myrank,nspec, & ! re-computes face infos ! weighted jacobian and normal - call get_jacobian_boundary_face(myrank,nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & - dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & - wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & - ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ,NGNOD2D) + call get_jacobian_boundary_face(nspec, & + xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & + dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & + wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & + ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ,NGNOD2D) ! normal convention: points away from element ! switch normal direction if necessary diff --git a/src/generate_databases/fault_generate_databases.f90 b/src/generate_databases/fault_generate_databases.f90 index e42aaeaef..360524b2a 100644 --- a/src/generate_databases/fault_generate_databases.f90 +++ b/src/generate_databases/fault_generate_databases.f90 @@ -90,12 +90,12 @@ module fault_generate_databases contains !================================================================================================================= - subroutine fault_read_input(prname,myrank) + subroutine fault_read_input(prname) - use constants, only: MAX_STRING_LEN, IN_DATA_FILES + use constants, only: MAX_STRING_LEN, IN_DATA_FILES,myrank + implicit none character(len=MAX_STRING_LEN), intent(in) :: prname - integer, intent(in) :: myrank integer :: nb,i,iflt,ier,nspec,dummy_node integer, parameter :: IIN_PAR = 100 @@ -188,12 +188,13 @@ end subroutine fault_read_input !================================================================================================================== subroutine fault_setup(ibool,nnodes_ext_mesh,nodes_coords_ext_mesh, & - xstore,ystore,zstore,nspec,nglob,myrank) + xstore,ystore,zstore,nspec,nglob) + + implicit none integer, intent(in) :: nspec double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec), intent(inout) :: xstore,ystore,zstore integer, dimension(NGLLX,NGLLY,NGLLZ,nspec), intent(in) :: ibool - integer, intent(in) :: myrank integer, intent(in) :: nglob integer, intent(in) :: nnodes_ext_mesh double precision, dimension(NDIM,nnodes_ext_mesh),intent(in) :: nodes_coords_ext_mesh @@ -226,7 +227,7 @@ subroutine fault_setup(ibool,nnodes_ext_mesh,nodes_coords_ext_mesh, & call save_fault_xyzcoord_ibulk(fault_db(iflt)) - call setup_normal_jacobian(fault_db(iflt),ibool,nspec,nglob,myrank) + call setup_normal_jacobian(fault_db(iflt),ibool,nspec,nglob) enddo @@ -497,7 +498,7 @@ end subroutine save_fault_xyzcoord_ibulk !================================================================================= - subroutine setup_normal_jacobian(fdb,ibool,nspec,nglob,myrank) + subroutine setup_normal_jacobian(fdb,ibool,nspec,nglob) use create_regions_mesh_ext_par, only: xstore_dummy,ystore_dummy,zstore_dummy, & dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & @@ -509,8 +510,6 @@ subroutine setup_normal_jacobian(fdb,ibool,nspec,nglob,myrank) type(fault_db_type), intent(inout) :: fdb integer, intent(in) :: nspec,nglob, ibool(NGLLX,NGLLY,NGLLZ,nspec) - integer, intent(in) :: myrank - ! (assumes NGLLX=NGLLY=NGLLZ) real(kind=CUSTOM_REAL),dimension(NGNOD2D) :: xcoord,ycoord,zcoord real(kind=CUSTOM_REAL) :: jacobian2Dw_face(NGLLX,NGLLY) @@ -545,11 +544,11 @@ subroutine setup_normal_jacobian(fdb,ibool,nspec,nglob,myrank) enddo ! gets face GLL 2Djacobian, weighted from element face - call get_jacobian_boundary_face(myrank,nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob, & - dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & - wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & - ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D) + call get_jacobian_boundary_face(nspec, & + xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob, & + dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & + wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & + ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D) ! normal convention: points away from domain1, reference element. do j=1,NGLLY diff --git a/src/generate_databases/generate_databases_par.F90 b/src/generate_databases/generate_databases_par.F90 index e173f2056..f93385778 100644 --- a/src/generate_databases/generate_databases_par.F90 +++ b/src/generate_databases/generate_databases_par.F90 @@ -27,24 +27,7 @@ module generate_databases_par - use constants, only: NGLLX,NGLLY,NGLLZ,NGLLSQUARE,NDIM,NDIM2D,NGNOD2D_FOUR_CORNERS,N_SLS, & - CUSTOM_REAL,SIZE_REAL,SIZE_DOUBLE, & - IMAIN,IIN,IOUT,ISTANDARD_OUTPUT, & - ZERO,ONE,TWO,FOUR_THIRDS,PI,HUGEVAL,GAUSSALPHA,GAUSSBETA, & - SMALLVAL_TOL,TINYVAL,HUGEVAL,R_EARTH, & - MAX_STRING_LEN,ATTENUATION_COMP_MAXIMUM, & - MINIMUM_THICKNESS_3D_OCEANS,RHO_APPROXIMATE_OCEAN_LOAD, & - CPML_X_ONLY,CPML_Y_ONLY,CPML_Z_ONLY, & - CPML_XY_ONLY,CPML_XZ_ONLY,CPML_YZ_ONLY,CPML_XYZ, & - NPOWER,CPML_Rcoef, & - IMODEL_DEFAULT,IMODEL_GLL,IMODEL_1D_PREM,IMODEL_1D_CASCADIA,IMODEL_1D_SOCAL, & - IMODEL_SALTON_TROUGH,IMODEL_TOMO,IMODEL_USER_EXTERNAL,IMODEL_IPATI,IMODEL_IPATI_WATER, & - IMODEL_1D_PREM_PB,IMODEL_SEP,IMODEL_COUPLED, & - IDOMAIN_ACOUSTIC,IDOMAIN_ELASTIC,IDOMAIN_POROELASTIC, & - OUTPUT_FILES, & - NX_TOPO_FILE,NY_TOPO_FILE, & - USE_MESH_COLORING_GPU,MAX_NUMBER_OF_COLORS, & - ADIOS_TRANSPORT_METHOD + use constants use shared_parameters @@ -60,7 +43,7 @@ module generate_databases_par double precision, dimension(:,:,:,:), allocatable :: xstore,ystore,zstore ! proc numbers for MPI - integer :: myrank,sizeprocs,ier + integer :: sizeprocs,ier integer :: NX_TOPO,NY_TOPO integer, dimension(:,:), allocatable :: itopo_bathy diff --git a/src/generate_databases/get_absorbing_boundary.F90 b/src/generate_databases/get_absorbing_boundary.F90 index 5059c9c00..ba9440a6d 100644 --- a/src/generate_databases/get_absorbing_boundary.F90 +++ b/src/generate_databases/get_absorbing_boundary.F90 @@ -159,11 +159,11 @@ subroutine get_absorbing_boundary(myrank,nspec,ibool, & call get_element_face_gll_indices(iface,ijk_face,NGLLX,NGLLZ) ! weighted jacobian and normal - call get_jacobian_boundary_face(myrank,nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & - dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & - wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & - ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ,NGNOD2D) + call get_jacobian_boundary_face(nspec, & + xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & + dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & + wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & + ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ,NGNOD2D) ! normal convention: points away from element ! switch normal direction if necessary @@ -233,11 +233,11 @@ subroutine get_absorbing_boundary(myrank,nspec,ibool, & call get_element_face_gll_indices(iface,ijk_face,NGLLX,NGLLZ) ! weighted jacobian and normal - call get_jacobian_boundary_face(myrank,nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & - dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & - wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & - ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ,NGNOD2D) + call get_jacobian_boundary_face(nspec, & + xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & + dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & + wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & + ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLZ,NGNOD2D) ! normal convention: points away from element ! switch normal direction if necessary @@ -301,11 +301,11 @@ subroutine get_absorbing_boundary(myrank,nspec,ibool, & call get_element_face_gll_indices(iface,ijk_face,NGLLY,NGLLZ) ! weighted jacobian and normal - call get_jacobian_boundary_face(myrank,nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & - dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & - wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & - ispec,iface,jacobian2Dw_face,normal_face,NGLLY,NGLLZ,NGNOD2D) + call get_jacobian_boundary_face(nspec, & + xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & + dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & + wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & + ispec,iface,jacobian2Dw_face,normal_face,NGLLY,NGLLZ,NGNOD2D) ! normal convention: points away from element ! switch normal direction if necessary @@ -369,11 +369,11 @@ subroutine get_absorbing_boundary(myrank,nspec,ibool, & call get_element_face_gll_indices(iface,ijk_face,NGLLY,NGLLZ) ! weighted jacobian and normal - call get_jacobian_boundary_face(myrank,nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & - dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & - wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & - ispec,iface,jacobian2Dw_face,normal_face,NGLLY,NGLLZ,NGNOD2D) + call get_jacobian_boundary_face(nspec, & + xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & + dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & + wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & + ispec,iface,jacobian2Dw_face,normal_face,NGLLY,NGLLZ,NGNOD2D) ! normal convention: points away from element ! switch normal direction if necessary @@ -437,11 +437,11 @@ subroutine get_absorbing_boundary(myrank,nspec,ibool, & call get_element_face_gll_indices(iface,ijk_face,NGLLX,NGLLY) ! weighted jacobian and normal - call get_jacobian_boundary_face(myrank,nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & - dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & - wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & - ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D) + call get_jacobian_boundary_face(nspec, & + xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & + dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & + wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & + ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D) ! normal convention: points away from element ! switch normal direction if necessary @@ -524,11 +524,11 @@ subroutine get_absorbing_boundary(myrank,nspec,ibool, & call get_element_face_gll_indices(iface,ijk_face,NGLLX,NGLLY) ! weighted jacobian and normal - call get_jacobian_boundary_face(myrank,nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & - dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & - wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & - ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D) + call get_jacobian_boundary_face(nspec, & + xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & + dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & + wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & + ispec,iface,jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D) ! normal convention: points away from element ! switch normal direction if necessary diff --git a/src/generate_databases/get_coupling_surfaces.f90 b/src/generate_databases/get_coupling_surfaces.f90 index 64be9be54..e55938095 100644 --- a/src/generate_databases/get_coupling_surfaces.f90 +++ b/src/generate_databases/get_coupling_surfaces.f90 @@ -299,11 +299,11 @@ subroutine get_coupling_surfaces_ac_el(myrank,nspec,ibool,elastic_flag) (elastic_flag(iglob_midpoint) >= 1)) then ! gets face GLL 2Djacobian, weighted from element face - call get_jacobian_boundary_face(myrank,nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & - dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & - wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & - ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D) + call get_jacobian_boundary_face(nspec, & + xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & + dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & + wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & + ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D) ! normal convention: points away from acoustic, reference element ! switch normal direction if necessary @@ -487,11 +487,11 @@ subroutine get_coupling_surfaces_ac_poro(myrank,nspec,ibool,acoustic_flag) call get_element_face_gll_indices(iface_ref,ijk_face,NGLLX,NGLLY) ! gets face GLL 2Djacobian, weighted from element face - call get_jacobian_boundary_face(myrank,nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & - dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & - wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & - ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D) + call get_jacobian_boundary_face(nspec, & + xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & + dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & + wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & + ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D) ! normal convention: points away from acoustic, reference element ! switch normal direction if necessary @@ -692,11 +692,11 @@ subroutine get_coupling_surfaces_el_poro(myrank,nspec,ibool,elastic_flag) call get_element_face_gll_indices(iface_ref_el,ijk_face_el,NGLLX,NGLLY) ! gets face GLL 2Djacobian, weighted from element face - call get_jacobian_boundary_face(myrank,nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & - dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & - wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & - ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D) + call get_jacobian_boundary_face(nspec, & + xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob_dummy, & + dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & + wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & + ispec,iface_ref,jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D) ! normal convention: points away from poroelastic, reference element do j=1,NGLLY diff --git a/src/generate_databases/get_model.F90 b/src/generate_databases/get_model.F90 index a8245171f..1027077f0 100644 --- a/src/generate_databases/get_model.F90 +++ b/src/generate_databases/get_model.F90 @@ -79,11 +79,11 @@ subroutine get_model(myrank) ! prepares external model values if needed select case (IMODEL) case (IMODEL_USER_EXTERNAL) - call model_external_broadcast(myrank) + call model_external_broadcast() case (IMODEL_SALTON_TROUGH) - call model_salton_trough_broadcast(myrank) + call model_salton_trough_broadcast() case (IMODEL_COUPLED) - call model_coupled_broadcast(myrank) + call model_coupled_broadcast() end select ! DP: not sure if this here should check if (COUPLE_WITH_INJECTION_TECHNIQUE .or. MESH_A_CHUNK_OF_THE_EARTH) .. diff --git a/src/generate_databases/get_perm_color.f90 b/src/generate_databases/get_perm_color.f90 index 5e932c78c..7bd1a7d28 100644 --- a/src/generate_databases/get_perm_color.f90 +++ b/src/generate_databases/get_perm_color.f90 @@ -41,8 +41,7 @@ subroutine get_perm_color_faster(is_on_a_slice_edge,ispec_is_d, & nspec,nglob, & nb_colors_outer_elements,nb_colors_inner_elements, & nspec_outer,nspec_inner,nspec_domain, & - first_elem_number_in_this_color, & - myrank) + first_elem_number_in_this_color) use constants @@ -60,14 +59,13 @@ subroutine get_perm_color_faster(is_on_a_slice_edge,ispec_is_d, & integer, intent(out) :: nb_colors_outer_elements,nb_colors_inner_elements integer, intent(out) :: nspec_outer,nspec_inner,nspec_domain - integer, intent(in) :: myrank ! local variables integer :: nb_colors ! coloring algorithm w/ Droux call get_color_faster(ibool, is_on_a_slice_edge, ispec_is_d, & - myrank, nspec, nglob, & + nspec, nglob, & color, nb_colors_outer_elements, nb_colors_inner_elements, & nspec_outer,nspec_inner,nspec_domain) @@ -100,20 +98,20 @@ end subroutine get_perm_color_faster ! subroutine get_color_faster(ibool, is_on_a_slice_edge, ispec_is_d, & - myrank, nspec, nglob, & - color, nb_colors_outer_elements, nb_colors_inner_elements, & - nspec_outer,nspec_inner,nspec_domain) + nspec, nglob, & + color, nb_colors_outer_elements, nb_colors_inner_elements, & + nspec_outer,nspec_inner,nspec_domain) use constants implicit none - integer nspec,nglob + integer :: nspec,nglob logical, dimension(nspec) :: is_on_a_slice_edge,ispec_is_d integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool integer, dimension(nspec) :: color - integer :: nb_colors_outer_elements,nb_colors_inner_elements,myrank + integer :: nb_colors_outer_elements,nb_colors_inner_elements integer :: nspec_outer,nspec_inner,nspec_domain @@ -188,8 +186,8 @@ subroutine get_color_faster(ibool, is_on_a_slice_edge, ispec_is_d, & if (DISPLAY_MIN_POSSIBLE_COLORS .or. try_Droux_coloring) then ! gets maximum values of valence for inner and outer element points call count_mesh_valence(ibool,is_on_a_slice_edge,ispec_is_d, & - myrank, nspec, nglob, & - maxval_count_ibool_outer,maxval_count_ibool_inner) + nspec, nglob, & + maxval_count_ibool_outer,maxval_count_ibool_inner) endif ! allocates mask @@ -349,7 +347,7 @@ subroutine get_color_faster(ibool, is_on_a_slice_edge, ispec_is_d, & ! tries to find a balanced coloring call balance_colors_Droux(ibool,is_on_a_slice_edge,ispec_is_d, & - myrank, nspec, nglob, & + nspec, nglob, & color,nb_colors_outer_elements,nb_colors_inner_elements, & nspec_outer,nspec_inner,maxval_count_ibool_inner, & mask_ibool,fail_safe) @@ -365,9 +363,9 @@ subroutine get_color_faster(ibool, is_on_a_slice_edge, ispec_is_d, & ! balances colors using a simple algorithm (if Droux was not used) if (BALANCE_COLORS_SIMPLE_ALGO) then call balance_colors_simple(ibool,is_on_a_slice_edge,ispec_is_d, & - myrank, nspec, nglob, & - color,nb_colors_outer_elements,nb_colors_inner_elements, & - nspec_outer,nspec_inner,mask_ibool) + nspec, nglob, & + color,nb_colors_outer_elements,nb_colors_inner_elements, & + nspec_outer,nspec_inner,mask_ibool) endif ! checks that all the color sets are independent @@ -425,7 +423,7 @@ end subroutine get_color_faster ! subroutine count_mesh_valence(ibool,is_on_a_slice_edge,ispec_is_d, & - myrank, nspec, nglob, & + nspec, nglob, & maxval_count_ibool_outer,maxval_count_ibool_inner) use constants @@ -433,12 +431,9 @@ subroutine count_mesh_valence(ibool,is_on_a_slice_edge,ispec_is_d, & implicit none integer :: nspec,nglob - integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool logical, dimension(nspec) :: is_on_a_slice_edge,ispec_is_d - - integer :: myrank integer :: maxval_count_ibool_outer,maxval_count_ibool_inner ! local parameters @@ -531,7 +526,7 @@ end subroutine count_mesh_valence ! subroutine balance_colors_Droux(ibool,is_on_a_slice_edge,ispec_is_d, & - myrank, nspec, nglob, & + nspec, nglob, & color, nb_colors_outer_elements, nb_colors_inner_elements, & nspec_outer,nspec_inner,maxval_count_ibool_inner, & mask_ibool,fail_safe) @@ -541,15 +536,12 @@ subroutine balance_colors_Droux(ibool,is_on_a_slice_edge,ispec_is_d, & implicit none integer :: nspec,nglob - integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool logical, dimension(nspec) :: is_on_a_slice_edge,ispec_is_d logical, dimension(nglob) :: mask_ibool integer, dimension(nspec) :: color - - integer :: myrank integer :: nb_colors_outer_elements,nb_colors_inner_elements integer :: nspec_outer,nspec_inner @@ -721,24 +713,21 @@ end subroutine balance_colors_Droux ! subroutine balance_colors_simple(ibool,is_on_a_slice_edge,ispec_is_d, & - myrank, nspec, nglob, & - color, nb_colors_outer_elements, nb_colors_inner_elements, & - nspec_outer,nspec_inner,mask_ibool) + nspec, nglob, & + color, nb_colors_outer_elements, nb_colors_inner_elements, & + nspec_outer,nspec_inner,mask_ibool) use constants implicit none integer :: nspec,nglob - integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool logical, dimension(nspec) :: is_on_a_slice_edge,ispec_is_d logical, dimension(nglob) :: mask_ibool integer, dimension(nspec) :: color - - integer :: myrank integer :: nb_colors_outer_elements,nb_colors_inner_elements integer :: nspec_outer,nspec_inner @@ -1042,7 +1031,6 @@ subroutine get_final_perm(color,perm,first_elem_number_in_this_color, & integer, intent(inout) :: first_elem_number_in_this_color(nb_colors) logical,dimension(nspec),intent(in) :: ispec_is_d - integer,intent(in) :: nb_colors_outer_elements,nspec_domain ! local parameters diff --git a/src/generate_databases/memory_eval.f90 b/src/generate_databases/memory_eval.f90 index 9dbde056e..203211119 100644 --- a/src/generate_databases/memory_eval.f90 +++ b/src/generate_databases/memory_eval.f90 @@ -224,18 +224,18 @@ end subroutine memory_eval ! compute the approximate amount of memory needed to run the mesher - subroutine memory_eval_mesher(myrank,nspec,npointot,nnodes_ext_mesh, & - nelmnts_ext_mesh,nmat_ext_mesh,num_interfaces_ext_mesh, & - max_interface_size_ext_mesh,nspec2D_xmin,nspec2D_xmax, & - nspec2D_ymin,nspec2D_ymax,nspec2D_bottom,nspec2D_top, & - memory_size_request) + subroutine memory_eval_mesher(nspec,npointot,nnodes_ext_mesh, & + nelmnts_ext_mesh,nmat_ext_mesh,num_interfaces_ext_mesh, & + max_interface_size_ext_mesh,nspec2D_xmin,nspec2D_xmax, & + nspec2D_ymin,nspec2D_ymax,nspec2D_bottom,nspec2D_top, & + memory_size_request) use constants use generate_databases_par, only: NGNOD,NGNOD2D implicit none - integer,intent(in) :: myrank,nspec,npointot,nnodes_ext_mesh,nelmnts_ext_mesh, & + integer,intent(in) :: nspec,npointot,nnodes_ext_mesh,nelmnts_ext_mesh, & nmat_ext_mesh,num_interfaces_ext_mesh, & max_interface_size_ext_mesh,nspec2D_xmin,nspec2D_xmax, & nspec2D_ymin,nspec2D_ymax,nspec2D_bottom,nspec2D_top diff --git a/src/generate_databases/model_coupled.f90 b/src/generate_databases/model_coupled.f90 index 6417a7967..c71a46943 100644 --- a/src/generate_databases/model_coupled.f90 +++ b/src/generate_databases/model_coupled.f90 @@ -56,7 +56,7 @@ end module model_coupled_par !------------------------------------------------------------------------------------------------- ! - subroutine model_coupled_broadcast(myrank) + subroutine model_coupled_broadcast() ! standard routine to setup model @@ -68,15 +68,13 @@ subroutine model_coupled_broadcast(myrank) implicit none - integer,intent(in) :: myrank - ! safety check if (.not. (COUPLE_WITH_INJECTION_TECHNIQUE .or. MESH_A_CHUNK_OF_THE_EARTH)) then print *,'Error: model coupling requires coupling with injection technique or mesh a chunk of the earth' stop 'Error model coupling' endif - call read_model_for_coupling_or_chunk(myrank) + call read_model_for_coupling_or_chunk() end subroutine model_coupled_broadcast @@ -84,14 +82,13 @@ end subroutine model_coupled_broadcast !------------------------------------------------------------------------------------------------- ! - subroutine read_model_for_coupling_or_chunk(myrank) + subroutine read_model_for_coupling_or_chunk() - use constants, only: IMAIN,IN_DATA_FILES + use constants, only: IMAIN,IN_DATA_FILES,myrank use model_coupled_par !! VM VM custom subroutine for coupling with DSM implicit none - integer, intent(in) :: myrank ! local parameters character(len=256) :: filename diff --git a/src/generate_databases/model_external_values.F90 b/src/generate_databases/model_external_values.F90 index 7f32a9b0a..d17baf9ba 100644 --- a/src/generate_databases/model_external_values.F90 +++ b/src/generate_databases/model_external_values.F90 @@ -59,7 +59,7 @@ end module external_model !------------------------------------------------------------------------------------------------- ! - subroutine model_external_broadcast(myrank) + subroutine model_external_broadcast() ! standard routine to setup model @@ -71,8 +71,6 @@ subroutine model_external_broadcast(myrank) implicit none - integer :: myrank - ! local parameters integer :: idummy diff --git a/src/generate_databases/model_salton_trough.f90 b/src/generate_databases/model_salton_trough.f90 index f0b9c032c..1459d9b70 100644 --- a/src/generate_databases/model_salton_trough.f90 +++ b/src/generate_databases/model_salton_trough.f90 @@ -53,18 +53,16 @@ end module salton_trough_par !------------------------------------------------------------------------------------------------- ! - subroutine model_salton_trough_broadcast(myrank) + subroutine model_salton_trough_broadcast() use constants use salton_trough_par - implicit none - integer :: myrank + implicit none ! local parameters integer :: ier - allocate(vp_array(GOCAD_ST_NU,GOCAD_ST_NV,GOCAD_ST_NW),stat=ier) if (ier /= 0) call exit_MPI_without_rank('error allocating array 867') if (ier /= 0) call exit_mpi(myrank,'error allocating vp_array for salton') diff --git a/src/generate_databases/read_parameters.f90 b/src/generate_databases/read_parameters.f90 index dd1743b09..7dbf43e95 100644 --- a/src/generate_databases/read_parameters.f90 +++ b/src/generate_databases/read_parameters.f90 @@ -37,7 +37,7 @@ subroutine read_parameters() ! reads Par_file BROADCAST_AFTER_READ = .true. - call read_parameter_file(myrank,BROADCAST_AFTER_READ) + call read_parameter_file(BROADCAST_AFTER_READ) ! check that the code is running with the requested nb of processes if (sizeprocs /= NPROC) then diff --git a/src/generate_databases/setup_color_perm.f90 b/src/generate_databases/setup_color_perm.f90 index a5e47f75c..bed6dba20 100644 --- a/src/generate_databases/setup_color_perm.f90 +++ b/src/generate_databases/setup_color_perm.f90 @@ -201,8 +201,7 @@ subroutine setup_color(myrank,nspec,nglob,ibool,perm,ispec_is_d,idomain, & nspec,nglob, & nb_colors_outer_elements,nb_colors_inner_elements, & nspec_outer,nspec_inner,nspec_domain, & - first_elem_number_in_this_color, & - myrank) + first_elem_number_in_this_color) ! for the last color, the next color is fictitious and its first (fictitious) element number is nspec + 1 first_elem_number_in_this_color(nb_colors_outer_elements + nb_colors_inner_elements + 1) & diff --git a/src/generate_databases/setup_mesh.f90 b/src/generate_databases/setup_mesh.f90 index a55b103a1..50cd12233 100644 --- a/src/generate_databases/setup_mesh.f90 +++ b/src/generate_databases/setup_mesh.f90 @@ -52,11 +52,11 @@ subroutine setup_mesh if (ier /= 0) call exit_MPI_without_rank('error allocating array 608') if (ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays') - call memory_eval_mesher(myrank,NSPEC_AB,npointot,nnodes_ext_mesh, & - nelmnts_ext_mesh,nmat_ext_mesh,num_interfaces_ext_mesh, & - max_interface_size_ext_mesh,nspec2D_xmin,nspec2D_xmax, & - nspec2D_ymin,nspec2D_ymax,nspec2D_bottom,nspec2D_top, & - max_memory_size_request) + call memory_eval_mesher(NSPEC_AB,npointot,nnodes_ext_mesh, & + nelmnts_ext_mesh,nmat_ext_mesh,num_interfaces_ext_mesh, & + max_interface_size_ext_mesh,nspec2D_xmin,nspec2D_xmax, & + nspec2D_ymin,nspec2D_ymax,nspec2D_bottom,nspec2D_top, & + max_memory_size_request) max_memory_size = max_memory_size_request diff --git a/src/inverse_problem_for_model/rules.mk b/src/inverse_problem_for_model/rules.mk index 63ff9d0de..1f2d55cba 100644 --- a/src/inverse_problem_for_model/rules.mk +++ b/src/inverse_problem_for_model/rules.mk @@ -214,6 +214,7 @@ inverse_problem_for_model_SHARED_OBJECTS = \ $O/gll_library.shared.o \ $O/heap_sort.shared.o \ $O/hex_nodes.shared.o \ + $O/init_openmp.shared.o \ $O/lagrange_poly.shared.o \ $O/netlib_specfun_erf.shared.o \ $O/param_reader.cc.o \ diff --git a/src/inverse_problem_for_model/specfem_interface/specfem_interface_mod.F90 b/src/inverse_problem_for_model/specfem_interface/specfem_interface_mod.F90 index 481c94fc8..e7f638682 100644 --- a/src/inverse_problem_for_model/specfem_interface/specfem_interface_mod.F90 +++ b/src/inverse_problem_for_model/specfem_interface/specfem_interface_mod.F90 @@ -747,7 +747,7 @@ subroutine InitSpecfemForOneRun(acqui_simu, ievent, inversion_param, iter_invers !! info on mesh and parameters --------------- if (SIMULATION_TYPE == 1) then if (ELASTIC_SIMULATION) then - call check_mesh_resolution(myrank,NSPEC_AB,NGLOB_AB, & + call check_mesh_resolution(NSPEC_AB,NGLOB_AB, & ibool,xstore,ystore,zstore, & kappastore,mustore,rho_vp,rho_vs, & DT_dble,model_speed_max,min_resolved_period, & @@ -761,7 +761,7 @@ subroutine InitSpecfemForOneRun(acqui_simu, ievent, inversion_param, iter_invers if (ier /= 0) stop 'error allocating array rho_vs' rho_vp = sqrt( kappastore / rhostore ) * rhostore rho_vs = 0.0_CUSTOM_REAL - call check_mesh_resolution(myrank,NSPEC_AB,NGLOB_AB, & + call check_mesh_resolution(NSPEC_AB,NGLOB_AB, & ibool,xstore,ystore,zstore, & kappastore,mustore,rho_vp,rho_vs, & DT_dble,model_speed_max,min_resolved_period, & diff --git a/src/meshfem3D/check_mesh_quality.f90 b/src/meshfem3D/check_mesh_quality.f90 index 118fcdbc6..6c03b5037 100644 --- a/src/meshfem3D/check_mesh_quality.f90 +++ b/src/meshfem3D/check_mesh_quality.f90 @@ -35,7 +35,7 @@ !! DK DK this routine could be improved by computing the mean in addition to min and max of ratios !! DK DK - subroutine check_mesh_quality(myrank,VP_MAX,NGLOB,NSPEC,x,y,z,ibool, & + subroutine check_mesh_quality(VP_MAX,NGLOB,NSPEC,x,y,z,ibool, & CREATE_VTK_FILES,prname) use constants @@ -43,8 +43,6 @@ subroutine check_mesh_quality(myrank,VP_MAX,NGLOB,NSPEC,x,y,z,ibool, & implicit none - integer :: myrank - double precision :: VP_MAX ! maximum vp in volume block id 3 integer :: NGLOB ! number of nodes diff --git a/src/meshfem3D/create_CPML_regions.f90 b/src/meshfem3D/create_CPML_regions.f90 index 8beb9309b..096f89d2c 100644 --- a/src/meshfem3D/create_CPML_regions.f90 +++ b/src/meshfem3D/create_CPML_regions.f90 @@ -27,7 +27,7 @@ subroutine create_CPML_regions(nspec,nglob,nodes_coords) - use meshfem3D_par, only: myrank,ibool,prname, & + use meshfem3D_par, only: ibool,prname, & nspec_CPML,is_CPML,CPML_to_spec,CPML_regions, & THICKNESS_OF_X_PML,THICKNESS_OF_Y_PML,THICKNESS_OF_Z_PML, & SUPPRESS_UTM_PROJECTION,CREATE_VTK_FILES @@ -35,7 +35,7 @@ subroutine create_CPML_regions(nspec,nglob,nodes_coords) ! create the different regions of the mesh use constants, only: IMAIN,CUSTOM_REAL,SMALL_PERCENTAGE_TOLERANCE, & CPML_X_ONLY,CPML_Y_ONLY,CPML_Z_ONLY,CPML_XY_ONLY,CPML_XZ_ONLY,CPML_YZ_ONLY,CPML_XYZ, & - PI,IOVTK,MAX_STRING_LEN,NDIM + PI,IOVTK,MAX_STRING_LEN,NDIM,myrank ! CPML use shared_parameters, only: PML_CONDITIONS,PML_INSTEAD_OF_FREE_SURFACE diff --git a/src/meshfem3D/create_interfaces_mesh.f90 b/src/meshfem3D/create_interfaces_mesh.f90 index 0c88278d4..25cd5b7ed 100644 --- a/src/meshfem3D/create_interfaces_mesh.f90 +++ b/src/meshfem3D/create_interfaces_mesh.f90 @@ -456,11 +456,11 @@ end subroutine create_interfaces_mesh subroutine get_interfaces_mesh_count() - use meshfem3D_par, only: myrank,INTERFACES_FILE, & + use meshfem3D_par, only: INTERFACES_FILE, & number_of_interfaces,max_npx_interface,max_npy_interface, & number_of_layers,ner_layer - use constants, only: IMAIN,IIN,MF_IN_DATA_FILES,DONT_IGNORE_JUNK,MAX_STRING_LEN + use constants, only: IMAIN,IIN,MF_IN_DATA_FILES,DONT_IGNORE_JUNK,MAX_STRING_LEN,myrank implicit none diff --git a/src/meshfem3D/create_meshfem_mesh.f90 b/src/meshfem3D/create_meshfem_mesh.f90 index ceada7e9a..92a53c888 100644 --- a/src/meshfem3D/create_meshfem_mesh.f90 +++ b/src/meshfem3D/create_meshfem_mesh.f90 @@ -62,7 +62,7 @@ subroutine create_meshfem_mesh() ! create the different regions of the mesh - use constants, only: IMAIN + use constants, only: IMAIN,myrank use meshfem3D_par, only: & ibool, & @@ -71,7 +71,7 @@ subroutine create_meshfem_mesh() NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, & NPROC_XI,NPROC_ETA, & NMATERIALS,material_properties, & - myrank, sizeprocs, prname, & + sizeprocs, prname, & LOCAL_PATH, & CREATE_ABAQUS_FILES,CREATE_DX_FILES,CREATE_VTK_FILES, & ADIOS_ENABLED, ADIOS_FOR_DATABASES, & @@ -131,22 +131,22 @@ subroutine create_meshfem_mesh() prname,nodes_coords,ibool,ispec_material_id) ! stores boundary informations - call store_boundaries(myrank,iboun,nspec, & - ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, & - nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, & - NSPEC2D_BOTTOM,NSPEC2D_TOP, & - NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX) + call store_boundaries(iboun,nspec, & + ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, & + nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, & + NSPEC2D_BOTTOM,NSPEC2D_TOP, & + NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX) ! checks mesh resolution VP_MAX = maxval(material_properties(:,2)) - call check_mesh_quality(myrank,VP_MAX,nglob,nspec, & - nodes_coords(:,1),nodes_coords(:,2),nodes_coords(:,3),ibool, & - CREATE_VTK_FILES,prname) + call check_mesh_quality(VP_MAX,nglob,nspec, & + nodes_coords(:,1),nodes_coords(:,2),nodes_coords(:,3),ibool, & + CREATE_VTK_FILES,prname) ! saves mesh as databases file if (ADIOS_FOR_DATABASES) then - call save_databases_adios(LOCAL_PATH, myrank, sizeprocs, & + call save_databases_adios(LOCAL_PATH,sizeprocs, & nspec,nglob,iproc_xi_current,iproc_eta_current, & NPROC_XI,NPROC_ETA,addressing,iMPIcut_xi,iMPIcut_eta, & ibool,nodes_coords,ispec_material_id, & @@ -189,11 +189,11 @@ end subroutine create_meshfem_mesh subroutine cmm_allocate_arrays() - use constants, only: IMAIN + use constants, only: IMAIN,myrank use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M,NGLLCUBE_M use meshfem3D_par, only: NSPEC_AB,nspec,ibool, & - xstore,ystore,zstore,npointot,myrank, & + xstore,ystore,zstore,npointot, & NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP use create_meshfem_par @@ -273,13 +273,13 @@ end subroutine cmm_allocate_arrays subroutine cmm_create_mesh_elements() - use constants, only: CUSTOM_REAL,IMAIN,NGNOD_EIGHT_CORNERS,HUGEVAL + use constants, only: CUSTOM_REAL,IMAIN,NGNOD_EIGHT_CORNERS,HUGEVAL,myrank use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M, & NGLOB_DOUBLING_SUPERBRICK,NSPEC_DOUBLING_SUPERBRICK, & IFLAG_BASEMENT_TOPO,IFLAG_ONE_LAYER_TOPOGRAPHY - use meshfem3D_par, only: myrank,xstore,ystore,zstore, & + use meshfem3D_par, only: xstore,ystore,zstore, & xgrid,ygrid,zgrid, & UTM_X_MIN,UTM_X_MAX,UTM_Y_MIN,UTM_Y_MAX,Z_DEPTH_BLOCK, & NEX_XI,NEX_ETA,NPROC_XI,NPROC_ETA, & @@ -408,7 +408,7 @@ subroutine cmm_create_mesh_elements() endif ! define shape of elements - call define_mesh_regions(myrank,USE_REGULAR_MESH,isubregion,NER,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, & + call define_mesh_regions(USE_REGULAR_MESH,isubregion,NER,NEX_PER_PROC_XI,NEX_PER_PROC_ETA, & iproc_xi_current,iproc_eta_current, & NDOUBLINGS,ner_doublings, & iaddx,iaddy,iaddz,ix1,ix2,dix,iy1,iy2,diy,ir1,ir2,dir,iax,iay,iar) @@ -651,10 +651,10 @@ end subroutine cmm_create_mesh_elements subroutine cmm_create_addressing(nglob) - use constants, only: NDIM,IOVTK,IMAIN + use constants, only: NDIM,IOVTK,IMAIN,myrank use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M,NGLLCUBE_M - use meshfem3D_par, only: myrank,prname,ibool,xstore,ystore,zstore, & + use meshfem3D_par, only: prname,ibool,xstore,ystore,zstore, & UTM_X_MIN,UTM_X_MAX,NGLOB_AB,nspec,npointot use create_meshfem_par diff --git a/src/meshfem3D/define_subregions.f90 b/src/meshfem3D/define_subregions.f90 index 62e41a71e..13c08adf0 100644 --- a/src/meshfem3D/define_subregions.f90 +++ b/src/meshfem3D/define_subregions.f90 @@ -88,7 +88,7 @@ end subroutine define_model_regions !------------------------------------------------------------------------------------ ! - subroutine define_mesh_regions(myrank,USE_REGULAR_MESH,isubregion,NER, & + subroutine define_mesh_regions(USE_REGULAR_MESH,isubregion,NER, & NEX_PER_PROC_XI,NEX_PER_PROC_ETA,iproc_xi,iproc_eta, & NDOUBLINGS,ner_doublings, & iaddx,iaddy,iaddz,ix1,ix2,dix,iy1,iy2,diy,ir1,ir2,dir,iax,iay,iar) @@ -97,7 +97,6 @@ subroutine define_mesh_regions(myrank,USE_REGULAR_MESH,isubregion,NER, & implicit none - integer,intent(in) :: myrank logical,intent(in) :: USE_REGULAR_MESH integer,intent(in) :: isubregion diff --git a/src/meshfem3D/define_subregions_heuristic.f90 b/src/meshfem3D/define_subregions_heuristic.f90 index 5452366b9..56a4c953c 100644 --- a/src/meshfem3D/define_subregions_heuristic.f90 +++ b/src/meshfem3D/define_subregions_heuristic.f90 @@ -28,7 +28,7 @@ ! note: ! this routine is not used anymore, but kept for reference... - subroutine define_subregions_heuristic(myrank,isubregion,iaddx,iaddy,iaddz, & + subroutine define_subregions_heuristic(isubregion,iaddx,iaddy,iaddz, & ix1,ix2,dix,iy1,iy2,diy,ir1,ir2,dir,iax,iay,iar, & itype_element,npx,npy, & NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM) @@ -41,18 +41,17 @@ subroutine define_subregions_heuristic(myrank,isubregion,iaddx,iaddy,iaddz, & implicit none - integer myrank - integer ix1,ix2,dix,iy1,iy2,diy,ir1,ir2,dir - integer iax,iay,iar - integer isubregion,itype_element - integer npx,npy + integer :: ix1,ix2,dix,iy1,iy2,diy,ir1,ir2,dir + integer :: iax,iay,iar + integer :: isubregion,itype_element + integer :: npx,npy - integer NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM + integer :: NER_BOTTOM_MOHO,NER_MOHO_16,NER_16_BASEMENT,NER_BASEMENT_SEDIM ! topology of the elements - integer iaddx(NGNOD_EIGHT_CORNERS) - integer iaddy(NGNOD_EIGHT_CORNERS) - integer iaddz(NGNOD_EIGHT_CORNERS) + integer :: iaddx(NGNOD_EIGHT_CORNERS) + integer :: iaddy(NGNOD_EIGHT_CORNERS) + integer :: iaddz(NGNOD_EIGHT_CORNERS) ! ************** diff --git a/src/meshfem3D/determine_cavity.f90 b/src/meshfem3D/determine_cavity.f90 index 04e17060d..76816b254 100644 --- a/src/meshfem3D/determine_cavity.f90 +++ b/src/meshfem3D/determine_cavity.f90 @@ -28,12 +28,12 @@ subroutine cmm_determine_cavity(nglob) - use constants, only: MF_IN_DATA_FILES,MAX_STRING_LEN,IMAIN,HUGEVAL,TINYVAL,NDIM + use constants, only: MF_IN_DATA_FILES,MAX_STRING_LEN,IMAIN,HUGEVAL,TINYVAL,NDIM,myrank use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M use create_meshfem_par, only: nodes_coords,ispec_material_id,iboun,iMPIcut_xi,iMPIcut_eta - use meshfem3D_par, only: ibool,xstore,ystore,zstore,nspec,NPROC_XI,NPROC_ETA,myrank,CAVITY_FILE + use meshfem3D_par, only: ibool,xstore,ystore,zstore,nspec,NPROC_XI,NPROC_ETA,CAVITY_FILE implicit none diff --git a/src/meshfem3D/earth_chunk.f90 b/src/meshfem3D/earth_chunk.f90 index 58f2feae9..cf30455f9 100644 --- a/src/meshfem3D/earth_chunk.f90 +++ b/src/meshfem3D/earth_chunk.f90 @@ -55,7 +55,6 @@ subroutine earth_chunk_HEX8_Mesher(NGNOD) !--- Parameters ! - integer, parameter :: myrank = 0 integer, parameter :: nlayer = 12 !! (number of layer in the model iasp91, or ak135, or prem (one more layer than the model) double precision, parameter :: GAUSSALPHA = 0.d0, GAUSSBETA = 0.d0 @@ -267,7 +266,7 @@ subroutine earth_chunk_HEX8_Mesher(NGNOD) ! !--- get the 3-D shape functions ! - call get_shape3D(myrank,shape3D,dershape3D,xigll,yigll,zigll,NGNOD) + call get_shape3D(shape3D,dershape3D,xigll,yigll,zigll,NGNOD) ! !--- rotation matrix to switch to the geographical coordinates @@ -780,7 +779,6 @@ subroutine earth_chunk_HEX27_Mesher(NGNOD) !--- Parameters ! - integer, parameter :: myrank = 0 integer, parameter :: nlayer = 12 !! (number of layer in the model iasp91, or ak135, or prem (one more layer than the model) double precision, parameter :: GAUSSALPHA = 0.d0, GAUSSBETA = 0.d0 @@ -1071,7 +1069,7 @@ subroutine earth_chunk_HEX27_Mesher(NGNOD) ! !--- get the 3-D shape functions ! - call get_shape3D(myrank,shape3D,dershape3D,xigll,yigll,zigll,NGNOD) + call get_shape3D(shape3D,dershape3D,xigll,yigll,zigll,NGNOD) ! !--- rotation matrix to switch to the geographical coordinates diff --git a/src/meshfem3D/get_MPI_cutplanes_eta.f90 b/src/meshfem3D/get_MPI_cutplanes_eta.f90 index ac8105d16..05b7a7c3f 100644 --- a/src/meshfem3D/get_MPI_cutplanes_eta.f90 +++ b/src/meshfem3D/get_MPI_cutplanes_eta.f90 @@ -28,158 +28,158 @@ ! not used any more, remains here for reference... ! note: MPI setup is now done in xgenerate_databases - subroutine get_MPI_cutplanes_eta(myrank,prname,nspec,iMPIcut_eta,ibool, & - xstore,ystore,zstore,mask_ibool,npointot, & - NSPEC2D_A_XI,NSPEC2D_B_XI) - -! this routine detects cut planes along eta -! In principle the left cut plane of the first slice -! and the right cut plane of the last slice are not used -! in the solver except if we want to have periodic conditions - - use constants - use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M - - implicit none - - integer nspec,myrank - integer NSPEC2D_A_XI,NSPEC2D_B_XI - - logical iMPIcut_eta(2,nspec) - - integer ibool(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) - - double precision xstore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) - double precision ystore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) - double precision zstore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) - -! logical mask used to create arrays iboolleft_eta and iboolright_eta - integer npointot - logical mask_ibool(npointot) - -! global element numbering - integer ispec - -! MPI cut-plane element numbering - integer ispecc1,ispecc2,npoin2D_eta,ix,iy,iz - integer nspec2Dtheor1,nspec2Dtheor2 - -! processor identification - character(len=MAX_STRING_LEN) :: prname - -! theoretical number of surface elements in the buffers -! cut planes along eta=constant correspond to XI faces - nspec2Dtheor1 = NSPEC2D_A_XI - nspec2Dtheor2 = NSPEC2D_B_XI - -! write the MPI buffers for the left and right edges of the slice -! and the position of the points to check that the buffers are fine - +! subroutine get_MPI_cutplanes_eta(prname,nspec,iMPIcut_eta,ibool, & +! xstore,ystore,zstore,mask_ibool,npointot, & +! NSPEC2D_A_XI,NSPEC2D_B_XI) ! -! determine if the element falls on the left MPI cut plane +!! this routine detects cut planes along eta +!! In principle the left cut plane of the first slice +!! and the right cut plane of the last slice are not used +!! in the solver except if we want to have periodic conditions ! - -! global point number and coordinates left MPI cut-plane - open(unit=10,file=prname(1:len_trim(prname))//'iboolleft_eta.txt',status='unknown') - -! erase the logical mask used to mark points already found - mask_ibool(:) = .false. - -! nb of global points shared with the other slice - npoin2D_eta = 0 - -! nb of elements in this cut-plane - ispecc1=0 - - do ispec=1,nspec - if (iMPIcut_eta(1,ispec)) then - - ispecc1=ispecc1+1 - - ! loop on all the points in that 2-D element, including edges - iy = 1 - do ix=1,NGLLX_M - do iz=1,NGLLZ_M - - ! select point, if not already selected - if (.not. mask_ibool(ibool(ix,iy,iz,ispec))) then - mask_ibool(ibool(ix,iy,iz,ispec)) = .true. - npoin2D_eta = npoin2D_eta + 1 - - write(10,*) ibool(ix,iy,iz,ispec),xstore(ix,iy,iz,ispec), & - ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec) - endif - - enddo - enddo - - endif - enddo - -! put flag to indicate end of the list of points - write(10,*) '0 0 0. 0. 0.' - -! write total number of points - write(10,*) npoin2D_eta - - close(10) - -! compare number of surface elements detected to analytical value - if (ispecc1 /= nspec2Dtheor1 .and. ispecc1 /= nspec2Dtheor2) & - call exit_MPI(myrank,'error MPI cut-planes detection in eta=left') - +! use constants +! use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M ! -! determine if the element falls on the right MPI cut plane +! implicit none ! - -! global point number and coordinates right MPI cut-plane - open(unit=10,file=prname(1:len_trim(prname))//'iboolright_eta.txt',status='unknown') - -! erase the logical mask used to mark points already found - mask_ibool(:) = .false. - -! nb of global points shared with the other slice - npoin2D_eta = 0 - -! nb of elements in this cut-plane - ispecc2=0 - - do ispec=1,nspec - if (iMPIcut_eta(2,ispec)) then - - ispecc2=ispecc2+1 - - ! loop on all the points in that 2-D element, including edges - iy = NGLLY_M - do ix=1,NGLLX_M - do iz=1,NGLLZ_M - - ! select point, if not already selected - if (.not. mask_ibool(ibool(ix,iy,iz,ispec))) then - mask_ibool(ibool(ix,iy,iz,ispec)) = .true. - npoin2D_eta = npoin2D_eta + 1 - - write(10,*) ibool(ix,iy,iz,ispec),xstore(ix,iy,iz,ispec), & - ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec) - endif - - enddo - enddo - - endif - enddo - -! put flag to indicate end of the list of points - write(10,*) '0 0 0. 0. 0.' - -! write total number of points - write(10,*) npoin2D_eta - - close(10) - -! compare number of surface elements detected to analytical value - if (ispecc2 /= nspec2Dtheor1 .and. ispecc2 /= nspec2Dtheor2) & - call exit_MPI(myrank,'error MPI cut-planes detection in eta=right') - - end subroutine get_MPI_cutplanes_eta +! integer :: nspec +! integer :: NSPEC2D_A_XI,NSPEC2D_B_XI +! +! logical :: iMPIcut_eta(2,nspec) +! +! integer :: ibool(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) +! +! double precision :: xstore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) +! double precision :: ystore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) +! double precision :: zstore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) +! +!! logical mask used to create arrays iboolleft_eta and iboolright_eta +! integer :: npointot +! logical :: mask_ibool(npointot) +! +!! global element numbering +! integer :: ispec +! +!! MPI cut-plane element numbering +! integer :: ispecc1,ispecc2,npoin2D_eta,ix,iy,iz +! integer :: nspec2Dtheor1,nspec2Dtheor2 +! +!! processor identification +! character(len=MAX_STRING_LEN) :: prname +! +!! theoretical number of surface elements in the buffers +!! cut planes along eta=constant correspond to XI faces +! nspec2Dtheor1 = NSPEC2D_A_XI +! nspec2Dtheor2 = NSPEC2D_B_XI +! +!! write the MPI buffers for the left and right edges of the slice +!! and the position of the points to check that the buffers are fine +! +!! +!! determine if the element falls on the left MPI cut plane +!! +! +!! global point number and coordinates left MPI cut-plane +! open(unit=10,file=prname(1:len_trim(prname))//'iboolleft_eta.txt',status='unknown') +! +!! erase the logical mask used to mark points already found +! mask_ibool(:) = .false. +! +!! nb of global points shared with the other slice +! npoin2D_eta = 0 +! +!! nb of elements in this cut-plane +! ispecc1=0 +! +! do ispec=1,nspec +! if (iMPIcut_eta(1,ispec)) then +! +! ispecc1=ispecc1+1 +! +! ! loop on all the points in that 2-D element, including edges +! iy = 1 +! do ix=1,NGLLX_M +! do iz=1,NGLLZ_M +! +! ! select point, if not already selected +! if (.not. mask_ibool(ibool(ix,iy,iz,ispec))) then +! mask_ibool(ibool(ix,iy,iz,ispec)) = .true. +! npoin2D_eta = npoin2D_eta + 1 +! +! write(10,*) ibool(ix,iy,iz,ispec),xstore(ix,iy,iz,ispec), & +! ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec) +! endif +! +! enddo +! enddo +! +! endif +! enddo +! +!! put flag to indicate end of the list of points +! write(10,*) '0 0 0. 0. 0.' +! +!! write total number of points +! write(10,*) npoin2D_eta +! +! close(10) +! +!! compare number of surface elements detected to analytical value +! if (ispecc1 /= nspec2Dtheor1 .and. ispecc1 /= nspec2Dtheor2) & +! call exit_MPI(myrank,'error MPI cut-planes detection in eta=left') +! +!! +!! determine if the element falls on the right MPI cut plane +!! +! +!! global point number and coordinates right MPI cut-plane +! open(unit=10,file=prname(1:len_trim(prname))//'iboolright_eta.txt',status='unknown') +! +!! erase the logical mask used to mark points already found +! mask_ibool(:) = .false. +! +!! nb of global points shared with the other slice +! npoin2D_eta = 0 +! +!! nb of elements in this cut-plane +! ispecc2=0 +! +! do ispec=1,nspec +! if (iMPIcut_eta(2,ispec)) then +! +! ispecc2=ispecc2+1 +! +! ! loop on all the points in that 2-D element, including edges +! iy = NGLLY_M +! do ix=1,NGLLX_M +! do iz=1,NGLLZ_M +! +! ! select point, if not already selected +! if (.not. mask_ibool(ibool(ix,iy,iz,ispec))) then +! mask_ibool(ibool(ix,iy,iz,ispec)) = .true. +! npoin2D_eta = npoin2D_eta + 1 +! +! write(10,*) ibool(ix,iy,iz,ispec),xstore(ix,iy,iz,ispec), & +! ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec) +! endif +! +! enddo +! enddo +! +! endif +! enddo +! +!! put flag to indicate end of the list of points +! write(10,*) '0 0 0. 0. 0.' +! +!! write total number of points +! write(10,*) npoin2D_eta +! +! close(10) +! +!! compare number of surface elements detected to analytical value +! if (ispecc2 /= nspec2Dtheor1 .and. ispecc2 /= nspec2Dtheor2) & +! call exit_MPI(myrank,'error MPI cut-planes detection in eta=right') +! +! end subroutine get_MPI_cutplanes_eta diff --git a/src/meshfem3D/get_MPI_cutplanes_xi.f90 b/src/meshfem3D/get_MPI_cutplanes_xi.f90 index 8593f67df..537986fa7 100644 --- a/src/meshfem3D/get_MPI_cutplanes_xi.f90 +++ b/src/meshfem3D/get_MPI_cutplanes_xi.f90 @@ -28,153 +28,154 @@ ! not used any more, remains here for reference... ! note: MPI setup is now done in xgenerate_databases - subroutine get_MPI_cutplanes_xi(myrank,prname,nspec,iMPIcut_xi,ibool, & - xstore,ystore,zstore,mask_ibool,npointot, & - NSPEC2D_A_ETA,NSPEC2D_B_ETA) - -! this routine detects cut planes along xi -! In principle the left cut plane of the first slice -! and the right cut plane of the last slice are not used -! in the solver except if we want to have periodic conditions - - use constants - use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M - - implicit none - - integer nspec,myrank - integer NSPEC2D_A_ETA,NSPEC2D_B_ETA - - logical iMPIcut_xi(2,nspec) - - integer ibool(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) - - double precision xstore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) - double precision ystore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) - double precision zstore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) - -! logical mask used to create arrays iboolleft_xi and iboolright_xi - integer npointot - logical mask_ibool(npointot) - -! global element numbering - integer ispec - -! MPI cut-plane element numbering - integer ispecc1,ispecc2,npoin2D_xi,ix,iy,iz - integer nspec2Dtheor1,nspec2Dtheor2 - -! processor identification - character(len=MAX_STRING_LEN) :: prname - -! theoretical number of surface elements in the buffers -! cut planes along xi=constant correspond to ETA faces - nspec2Dtheor1 = NSPEC2D_A_ETA - nspec2Dtheor2 = NSPEC2D_B_ETA - -! write the MPI buffers for the left and right edges of the slice -! and the position of the points to check that the buffers are fine - +! subroutine get_MPI_cutplanes_xi(prname,nspec,iMPIcut_xi,ibool, & +! xstore,ystore,zstore,mask_ibool,npointot, & +! NSPEC2D_A_ETA,NSPEC2D_B_ETA) ! -! determine if the element falls on the left MPI cut plane +!! this routine detects cut planes along xi +!! In principle the left cut plane of the first slice +!! and the right cut plane of the last slice are not used +!! in the solver except if we want to have periodic conditions ! - -! global point number and coordinates left MPI cut-plane - open(unit=10,file=prname(1:len_trim(prname))//'iboolleft_xi.txt',status='unknown') - -! erase the logical mask used to mark points already found - mask_ibool(:) = .false. - -! nb of global points shared with the other slice - npoin2D_xi = 0 - -! nb of elements in this cut-plane - ispecc1=0 - - do ispec=1,nspec - if (iMPIcut_xi(1,ispec)) then - - ispecc1=ispecc1+1 - - ! loop on all the points in that 2-D element, including edges - ix = 1 - do iy=1,NGLLY_M - do iz=1,NGLLZ_M - - ! select point, if not already selected - if (.not. mask_ibool(ibool(ix,iy,iz,ispec))) then - mask_ibool(ibool(ix,iy,iz,ispec)) = .true. - npoin2D_xi = npoin2D_xi + 1 - - write(10,*) ibool(ix,iy,iz,ispec),xstore(ix,iy,iz,ispec), & - ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec) - endif - enddo - enddo - endif - enddo - -! put flag to indicate end of the list of points - write(10,*) '0 0 0. 0. 0.' - -! write total number of points - write(10,*) npoin2D_xi - - close(10) - -! compare number of surface elements detected to analytical value - if (ispecc1 /= nspec2Dtheor1 .and. ispecc1 /= nspec2Dtheor2) & - call exit_MPI(myrank,'error MPI cut-planes detection in xi=left') - +! use constants +! use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M ! -! determine if the element falls on the right MPI cut plane +! implicit none ! +! integer :: nspec +! integer :: NSPEC2D_A_ETA,NSPEC2D_B_ETA +! +! logical :: iMPIcut_xi(2,nspec) +! +! integer :: ibool(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) +! +! double precision :: xstore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) +! double precision :: ystore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) +! double precision :: zstore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) +! +!! logical mask used to create arrays iboolleft_xi and iboolright_xi +! integer :: npointot +! logical :: mask_ibool(npointot) +! +!! global element numbering +! integer :: ispec +! +!! MPI cut-plane element numbering +! integer :: ispecc1,ispecc2,npoin2D_xi,ix,iy,iz +! integer :: nspec2Dtheor1,nspec2Dtheor2 +! +!! processor identification +! character(len=MAX_STRING_LEN) :: prname +! +!! theoretical number of surface elements in the buffers +!! cut planes along xi=constant correspond to ETA faces +! nspec2Dtheor1 = NSPEC2D_A_ETA +! nspec2Dtheor2 = NSPEC2D_B_ETA +! +!! write the MPI buffers for the left and right edges of the slice +!! and the position of the points to check that the buffers are fine +! +!! +!! determine if the element falls on the left MPI cut plane +!! +! +!! global point number and coordinates left MPI cut-plane +! open(unit=10,file=prname(1:len_trim(prname))//'iboolleft_xi.txt',status='unknown') +! +!! erase the logical mask used to mark points already found +! mask_ibool(:) = .false. +! +!! nb of global points shared with the other slice +! npoin2D_xi = 0 +! +!! nb of elements in this cut-plane +! ispecc1=0 +! +! do ispec=1,nspec +! if (iMPIcut_xi(1,ispec)) then +! +! ispecc1=ispecc1+1 +! +! ! loop on all the points in that 2-D element, including edges +! ix = 1 +! do iy=1,NGLLY_M +! do iz=1,NGLLZ_M +! +! ! select point, if not already selected +! if (.not. mask_ibool(ibool(ix,iy,iz,ispec))) then +! mask_ibool(ibool(ix,iy,iz,ispec)) = .true. +! npoin2D_xi = npoin2D_xi + 1 +! +! write(10,*) ibool(ix,iy,iz,ispec),xstore(ix,iy,iz,ispec), & +! ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec) +! endif +! enddo +! enddo +! endif +! enddo +! +!! put flag to indicate end of the list of points +! write(10,*) '0 0 0. 0. 0.' +! +!! write total number of points +! write(10,*) npoin2D_xi +! +! close(10) +! +!! compare number of surface elements detected to analytical value +! if (ispecc1 /= nspec2Dtheor1 .and. ispecc1 /= nspec2Dtheor2) & +! call exit_MPI(myrank,'error MPI cut-planes detection in xi=left') +! +!! +!! determine if the element falls on the right MPI cut plane +!! +! +!! global point number and coordinates right MPI cut-plane +! open(unit=10,file=prname(1:len_trim(prname))//'iboolright_xi.txt',status='unknown') +! +!! erase the logical mask used to mark points already found +! mask_ibool(:) = .false. +! +!! nb of global points shared with the other slice +! npoin2D_xi = 0 +! +!! nb of elements in this cut-plane +! ispecc2=0 +! +! do ispec=1,nspec +! if (iMPIcut_xi(2,ispec)) then +! +! ispecc2=ispecc2+1 +! +! ! loop on all the points in that 2-D element, including edges +! ix = NGLLX_M +! do iy=1,NGLLY_M +! do iz=1,NGLLZ_M +! +! ! select point, if not already selected +! if (.not. mask_ibool(ibool(ix,iy,iz,ispec))) then +! mask_ibool(ibool(ix,iy,iz,ispec)) = .true. +! npoin2D_xi = npoin2D_xi + 1 +! +! write(10,*) ibool(ix,iy,iz,ispec),xstore(ix,iy,iz,ispec), & +! ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec) +! endif +! enddo +! enddo +! endif +! enddo +! +!! put flag to indicate end of the list of points +! write(10,*) '0 0 0. 0. 0.' +! +!! write total number of points +! write(10,*) npoin2D_xi +! +! close(10) +! +!! compare number of surface elements detected to analytical value +! if (ispecc2 /= nspec2Dtheor1 .and. ispecc2 /= nspec2Dtheor2) & +! call exit_MPI(myrank,'error MPI cut-planes detection in xi=right') +! +! end subroutine get_MPI_cutplanes_xi -! global point number and coordinates right MPI cut-plane - open(unit=10,file=prname(1:len_trim(prname))//'iboolright_xi.txt',status='unknown') - -! erase the logical mask used to mark points already found - mask_ibool(:) = .false. - -! nb of global points shared with the other slice - npoin2D_xi = 0 - -! nb of elements in this cut-plane - ispecc2=0 - - do ispec=1,nspec - if (iMPIcut_xi(2,ispec)) then - - ispecc2=ispecc2+1 - - ! loop on all the points in that 2-D element, including edges - ix = NGLLX_M - do iy=1,NGLLY_M - do iz=1,NGLLZ_M - - ! select point, if not already selected - if (.not. mask_ibool(ibool(ix,iy,iz,ispec))) then - mask_ibool(ibool(ix,iy,iz,ispec)) = .true. - npoin2D_xi = npoin2D_xi + 1 - - write(10,*) ibool(ix,iy,iz,ispec),xstore(ix,iy,iz,ispec), & - ystore(ix,iy,iz,ispec),zstore(ix,iy,iz,ispec) - endif - enddo - enddo - endif - enddo - -! put flag to indicate end of the list of points - write(10,*) '0 0 0. 0. 0.' - -! write total number of points - write(10,*) npoin2D_xi - - close(10) - -! compare number of surface elements detected to analytical value - if (ispecc2 /= nspec2Dtheor1 .and. ispecc2 /= nspec2Dtheor2) & - call exit_MPI(myrank,'error MPI cut-planes detection in xi=right') - - end subroutine get_MPI_cutplanes_xi diff --git a/src/meshfem3D/get_flags_boundaries.f90 b/src/meshfem3D/get_flags_boundaries.f90 index a5a357cdb..fbd68bedc 100644 --- a/src/meshfem3D/get_flags_boundaries.f90 +++ b/src/meshfem3D/get_flags_boundaries.f90 @@ -25,6 +25,7 @@ ! !===================================================================== + subroutine get_flags_boundaries(nspec,iproc_xi,iproc_eta,ispec,idoubling, & xstore,ystore,zstore,iboun,iMPIcut_xi,iMPIcut_eta, & NPROC_XI,NPROC_ETA, & diff --git a/src/meshfem3D/meshfem3D.F90 b/src/meshfem3D/meshfem3D.F90 index d448f80b3..fbda0f718 100644 --- a/src/meshfem3D/meshfem3D.F90 +++ b/src/meshfem3D/meshfem3D.F90 @@ -284,7 +284,7 @@ program xmeshfem3D ! read the parameter file (DATA/Par_file) BROADCAST_AFTER_READ = .true. - call read_parameter_file(myrank,BROADCAST_AFTER_READ) + call read_parameter_file(BROADCAST_AFTER_READ) ! make sure everybody is synchronized call synchronize_all() diff --git a/src/meshfem3D/meshfem3D_adios_stubs.f90 b/src/meshfem3D/meshfem3D_adios_stubs.f90 index a6c38cba8..48feecd7a 100644 --- a/src/meshfem3D/meshfem3D_adios_stubs.f90 +++ b/src/meshfem3D/meshfem3D_adios_stubs.f90 @@ -34,15 +34,16 @@ !! \author MPBL !============================================================================== -subroutine save_databases_adios(LOCAL_PATH, myrank, sizeprocs, & - nspec,nglob,iproc_xi,iproc_eta, & - NPROC_XI,NPROC_ETA,addressing,iMPIcut_xi,iMPIcut_eta, & - ibool,nodes_coords,ispec_material_id, & - nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, & - NSPEC2D_BOTTOM,NSPEC2D_TOP, NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX, & - ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, & - NMATERIALS,material_properties, & - nspec_CPML,CPML_to_spec,CPML_regions,is_CPML) + + subroutine save_databases_adios(LOCAL_PATH,sizeprocs, & + nspec,nglob,iproc_xi,iproc_eta, & + NPROC_XI,NPROC_ETA,addressing,iMPIcut_xi,iMPIcut_eta, & + ibool,nodes_coords,ispec_material_id, & + nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, & + NSPEC2D_BOTTOM,NSPEC2D_TOP, NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX, & + ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, & + NMATERIALS,material_properties, & + nspec_CPML,CPML_to_spec,CPML_regions,is_CPML) use constants use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M,NUMBER_OF_MATERIAL_PROPERTIES @@ -52,7 +53,7 @@ subroutine save_databases_adios(LOCAL_PATH, myrank, sizeprocs, & implicit none character(len=MAX_STRING_LEN) :: LOCAL_PATH - integer :: myrank, sizeprocs + integer :: sizeprocs integer :: nspec integer :: nglob integer :: iproc_xi,iproc_eta @@ -81,7 +82,6 @@ subroutine save_databases_adios(LOCAL_PATH, myrank, sizeprocs, & double precision :: unused_dp unused_i4 = len_trim(LOCAL_PATH) - unused_i4 = myrank unused_i4 = sizeprocs unused_i4 = iproc_xi unused_i4 = iproc_eta diff --git a/src/meshfem3D/meshfem3D_par.f90 b/src/meshfem3D/meshfem3D_par.f90 index 2b5a5e36c..abfca5ce1 100644 --- a/src/meshfem3D/meshfem3D_par.f90 +++ b/src/meshfem3D/meshfem3D_par.f90 @@ -63,7 +63,7 @@ module meshfem3D_par double precision, dimension(:,:,:,:), allocatable :: xstore,ystore,zstore ! proc numbers for MPI - integer :: myrank,sizeprocs + integer :: sizeprocs ! mesh point steps for interfaces integer :: npx_element_steps,npy_element_steps diff --git a/src/meshfem3D/save_databases.F90 b/src/meshfem3D/save_databases.F90 index cceb4f38c..34b79e67f 100644 --- a/src/meshfem3D/save_databases.F90 +++ b/src/meshfem3D/save_databases.F90 @@ -609,7 +609,6 @@ subroutine save_output_mesh_files_for_coupled_model(nspec, & double precision :: z_bottom ! for axisem coupling case ( only serial case for mesher use scotch after) - integer, parameter :: myrank = 0 integer, parameter :: nlayer = 12 !! (number of layer in the model iasp91, or ak135, or prem (one more layer than the model) double precision, parameter :: GAUSSALPHA = 0.d0, GAUSSBETA = 0.d0 double precision :: rotation_matrix(3,3) @@ -684,7 +683,7 @@ subroutine save_output_mesh_files_for_coupled_model(nspec, & ! allocate(shape3D(NGNOD,NGLLX,NGLLY,NGLLZ),dershape3D(NDIM,NGNOD,NGLLX,NGLLY,NGLLZ),stat=ier) if (ier /= 0) call exit_MPI_without_rank('error allocating array 1351') - call get_shape3D(myrank,shape3D,dershape3D,xigll,yigll,zigll,NGNOD) + call get_shape3D(shape3D,dershape3D,xigll,yigll,zigll,NGNOD) ! !! reading parameters for coupling diff --git a/src/meshfem3D/save_databases_adios.F90 b/src/meshfem3D/save_databases_adios.F90 index a663eb38b..9b504cd84 100644 --- a/src/meshfem3D/save_databases_adios.F90 +++ b/src/meshfem3D/save_databases_adios.F90 @@ -34,7 +34,7 @@ !============================================================================== -subroutine save_databases_adios(LOCAL_PATH, myrank, sizeprocs, & +subroutine save_databases_adios(LOCAL_PATH,sizeprocs, & nspec,nglob,iproc_xi,iproc_eta, & NPROC_XI,NPROC_ETA,addressing,iMPIcut_xi,iMPIcut_eta, & ibool,nodes_coords,ispec_material_id, & @@ -45,7 +45,8 @@ subroutine save_databases_adios(LOCAL_PATH, myrank, sizeprocs, & nspec_CPML,CPML_to_spec,CPML_regions,is_CPML) use constants, only: MAX_STRING_LEN,IDOMAIN_ACOUSTIC,IDOMAIN_ELASTIC,ADIOS_TRANSPORT_METHOD, & - NGLLX,NGLLY,NGLLZ,NDIM + NGLLX,NGLLY,NGLLZ,NDIM,myrank + use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M,NUMBER_OF_MATERIAL_PROPERTIES use adios_helpers_mod, only: define_adios_global_array1d,define_adios_scalar, & @@ -57,7 +58,7 @@ subroutine save_databases_adios(LOCAL_PATH, myrank, sizeprocs, & implicit none ! MPI variables - integer :: myrank, sizeprocs + integer :: sizeprocs ! number of spectral elements in each block integer nspec diff --git a/src/meshfem3D/store_boundaries.f90 b/src/meshfem3D/store_boundaries.f90 index 27501bd96..d2d75211d 100644 --- a/src/meshfem3D/store_boundaries.f90 +++ b/src/meshfem3D/store_boundaries.f90 @@ -25,7 +25,7 @@ ! !===================================================================== - subroutine store_boundaries(myrank,iboun,nspec, & + subroutine store_boundaries(iboun,nspec, & ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, & nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, & NSPEC2D_BOTTOM,NSPEC2D_TOP, & @@ -35,21 +35,21 @@ subroutine store_boundaries(myrank,iboun,nspec, & implicit none - integer nspec,myrank - integer NSPEC2D_BOTTOM,NSPEC2D_TOP,NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX + integer :: nspec + integer :: NSPEC2D_BOTTOM,NSPEC2D_TOP,NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX - integer nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax - integer ibelm_xmin(NSPEC2DMAX_XMIN_XMAX),ibelm_xmax(NSPEC2DMAX_XMIN_XMAX) - integer ibelm_ymin(NSPEC2DMAX_YMIN_YMAX),ibelm_ymax(NSPEC2DMAX_YMIN_YMAX) - integer ibelm_bottom(NSPEC2D_BOTTOM),ibelm_top(NSPEC2D_TOP) + integer :: nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax + integer :: ibelm_xmin(NSPEC2DMAX_XMIN_XMAX),ibelm_xmax(NSPEC2DMAX_XMIN_XMAX) + integer :: ibelm_ymin(NSPEC2DMAX_YMIN_YMAX),ibelm_ymax(NSPEC2DMAX_YMIN_YMAX) + integer :: ibelm_bottom(NSPEC2D_BOTTOM),ibelm_top(NSPEC2D_TOP) - logical iboun(6,nspec) + logical :: iboun(6,nspec) ! global element numbering - integer ispec + integer :: ispec ! counters to keep track of number of elements on each of the boundaries - integer ispecb1,ispecb2,ispecb3,ispecb4,ispecb5,ispecb6 + integer :: ispecb1,ispecb2,ispecb3,ispecb4,ispecb5,ispecb6 ! initializes ispecb1 = 0 @@ -123,300 +123,306 @@ end subroutine store_boundaries !------------------------------------------------------------------------------------------------- ! - subroutine get_jacobian_boundaries(myrank,iboun,nspec,xstore,ystore,zstore, & - dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & - ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, & - nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, & - jacobian2D_xmin,jacobian2D_xmax, & - jacobian2D_ymin,jacobian2D_ymax, & - jacobian2D_bottom,jacobian2D_top, & - normal_xmin,normal_xmax, & - normal_ymin,normal_ymax, & - normal_bottom,normal_top, & - NSPEC2D_BOTTOM,NSPEC2D_TOP, & - NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX) - - use constants - use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M - - implicit none - - integer nspec,myrank - integer NSPEC2D_BOTTOM,NSPEC2D_TOP,NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX - - integer nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax - integer ibelm_xmin(NSPEC2DMAX_XMIN_XMAX),ibelm_xmax(NSPEC2DMAX_XMIN_XMAX) - integer ibelm_ymin(NSPEC2DMAX_YMIN_YMAX),ibelm_ymax(NSPEC2DMAX_YMIN_YMAX) - integer ibelm_bottom(NSPEC2D_BOTTOM),ibelm_top(NSPEC2D_TOP) - - logical iboun(6,nspec) - - double precision xstore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) - double precision ystore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) - double precision zstore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) - - real(kind=CUSTOM_REAL) jacobian2D_xmin(NGLLY_M,NGLLZ_M,NSPEC2DMAX_XMIN_XMAX) - real(kind=CUSTOM_REAL) jacobian2D_xmax(NGLLY_M,NGLLZ_M,NSPEC2DMAX_XMIN_XMAX) - real(kind=CUSTOM_REAL) jacobian2D_ymin(NGLLX_M,NGLLZ_M,NSPEC2DMAX_YMIN_YMAX) - real(kind=CUSTOM_REAL) jacobian2D_ymax(NGLLX_M,NGLLZ_M,NSPEC2DMAX_YMIN_YMAX) - real(kind=CUSTOM_REAL) jacobian2D_bottom(NGLLX_M,NGLLY_M,NSPEC2D_BOTTOM) - real(kind=CUSTOM_REAL) jacobian2D_top(NGLLX_M,NGLLY_M,NSPEC2D_TOP) - - real(kind=CUSTOM_REAL) normal_xmin(NDIM,NGLLY_M,NGLLZ_M,NSPEC2DMAX_XMIN_XMAX) - real(kind=CUSTOM_REAL) normal_xmax(NDIM,NGLLY_M,NGLLZ_M,NSPEC2DMAX_XMIN_XMAX) - real(kind=CUSTOM_REAL) normal_ymin(NDIM,NGLLX_M,NGLLZ_M,NSPEC2DMAX_YMIN_YMAX) - real(kind=CUSTOM_REAL) normal_ymax(NDIM,NGLLX_M,NGLLZ_M,NSPEC2DMAX_YMIN_YMAX) - real(kind=CUSTOM_REAL) normal_bottom(NDIM,NGLLX_M,NGLLY_M,NSPEC2D_BOTTOM) - real(kind=CUSTOM_REAL) normal_top(NDIM,NGLLX_M,NGLLY_M,NSPEC2D_TOP) - - double precision dershape2D_x(NDIM2D,NGNOD2D_FOUR_CORNERS,NGLLY_M,NGLLZ_M) - double precision dershape2D_y(NDIM2D,NGNOD2D_FOUR_CORNERS,NGLLX_M,NGLLZ_M) - double precision dershape2D_bottom(NDIM2D,NGNOD2D_FOUR_CORNERS,NGLLX_M,NGLLY_M) - double precision dershape2D_top(NDIM2D,NGNOD2D_FOUR_CORNERS,NGLLX_M,NGLLY_M) - -! global element numbering - integer ispec - -! counters to keep track of number of elements on each of the boundaries - integer ispecb1,ispecb2,ispecb3,ispecb4,ispecb5,ispecb6 - - double precision xelm(NGNOD2D_FOUR_CORNERS),yelm(NGNOD2D_FOUR_CORNERS),zelm(NGNOD2D_FOUR_CORNERS) - - ispecb1 = 0 - ispecb2 = 0 - ispecb3 = 0 - ispecb4 = 0 - ispecb5 = 0 - ispecb6 = 0 - - do ispec=1,nspec - -! determine if the element falls on a boundary - -! on boundary: xmin - - if (iboun(1,ispec)) then - - ispecb1 = ispecb1+1 - ibelm_xmin(ispecb1) = ispec - -! specify the 4 nodes for the 2-D boundary element - xelm(1)=xstore(1,1,1,ispec) - yelm(1)=ystore(1,1,1,ispec) - zelm(1)=zstore(1,1,1,ispec) - xelm(2)=xstore(1,NGLLY_M,1,ispec) - yelm(2)=ystore(1,NGLLY_M,1,ispec) - zelm(2)=zstore(1,NGLLY_M,1,ispec) - xelm(3)=xstore(1,NGLLY_M,NGLLZ_M,ispec) - yelm(3)=ystore(1,NGLLY_M,NGLLZ_M,ispec) - zelm(3)=zstore(1,NGLLY_M,NGLLZ_M,ispec) - xelm(4)=xstore(1,1,NGLLZ_M,ispec) - yelm(4)=ystore(1,1,NGLLZ_M,ispec) - zelm(4)=zstore(1,1,NGLLZ_M,ispec) - - call compute_jacobian_2D(myrank,ispecb1,xelm,yelm,zelm,dershape2D_x, & - jacobian2D_xmin,normal_xmin,NGLLY_M,NGLLZ_M,NSPEC2DMAX_XMIN_XMAX) - - endif - -! on boundary: xmax - - if (iboun(2,ispec)) then - - ispecb2 = ispecb2+1 - ibelm_xmax(ispecb2) = ispec - -! specify the 4 nodes for the 2-D boundary element - xelm(1)=xstore(NGLLX_M,1,1,ispec) - yelm(1)=ystore(NGLLX_M,1,1,ispec) - zelm(1)=zstore(NGLLX_M,1,1,ispec) - xelm(2)=xstore(NGLLX_M,NGLLY_M,1,ispec) - yelm(2)=ystore(NGLLX_M,NGLLY_M,1,ispec) - zelm(2)=zstore(NGLLX_M,NGLLY_M,1,ispec) - xelm(3)=xstore(NGLLX_M,NGLLY_M,NGLLZ_M,ispec) - yelm(3)=ystore(NGLLX_M,NGLLY_M,NGLLZ_M,ispec) - zelm(3)=zstore(NGLLX_M,NGLLY_M,NGLLZ_M,ispec) - xelm(4)=xstore(NGLLX_M,1,NGLLZ_M,ispec) - yelm(4)=ystore(NGLLX_M,1,NGLLZ_M,ispec) - zelm(4)=zstore(NGLLX_M,1,NGLLZ_M,ispec) - - call compute_jacobian_2D(myrank,ispecb2,xelm,yelm,zelm,dershape2D_x, & - jacobian2D_xmax,normal_xmax,NGLLY_M,NGLLZ_M,NSPEC2DMAX_XMIN_XMAX) - - endif - -! on boundary: ymin - - if (iboun(3,ispec)) then - - ispecb3 = ispecb3+1 - ibelm_ymin(ispecb3) = ispec - -! specify the 4 nodes for the 2-D boundary element - xelm(1)=xstore(1,1,1,ispec) - yelm(1)=ystore(1,1,1,ispec) - zelm(1)=zstore(1,1,1,ispec) - xelm(2)=xstore(NGLLX_M,1,1,ispec) - yelm(2)=ystore(NGLLX_M,1,1,ispec) - zelm(2)=zstore(NGLLX_M,1,1,ispec) - xelm(3)=xstore(NGLLX_M,1,NGLLZ_M,ispec) - yelm(3)=ystore(NGLLX_M,1,NGLLZ_M,ispec) - zelm(3)=zstore(NGLLX_M,1,NGLLZ_M,ispec) - xelm(4)=xstore(1,1,NGLLZ_M,ispec) - yelm(4)=ystore(1,1,NGLLZ_M,ispec) - zelm(4)=zstore(1,1,NGLLZ_M,ispec) - - call compute_jacobian_2D(myrank,ispecb3,xelm,yelm,zelm,dershape2D_y, & - jacobian2D_ymin,normal_ymin,NGLLX_M,NGLLZ_M,NSPEC2DMAX_YMIN_YMAX) - - endif - -! on boundary: ymax - - if (iboun(4,ispec)) then - - ispecb4 = ispecb4+1 - ibelm_ymax(ispecb4) = ispec - -! specify the 4 nodes for the 2-D boundary element - xelm(1)=xstore(1,NGLLY_M,1,ispec) - yelm(1)=ystore(1,NGLLY_M,1,ispec) - zelm(1)=zstore(1,NGLLY_M,1,ispec) - xelm(2)=xstore(NGLLX_M,NGLLY_M,1,ispec) - yelm(2)=ystore(NGLLX_M,NGLLY_M,1,ispec) - zelm(2)=zstore(NGLLX_M,NGLLY_M,1,ispec) - xelm(3)=xstore(NGLLX_M,NGLLY_M,NGLLZ_M,ispec) - yelm(3)=ystore(NGLLX_M,NGLLY_M,NGLLZ_M,ispec) - zelm(3)=zstore(NGLLX_M,NGLLY_M,NGLLZ_M,ispec) - xelm(4)=xstore(1,NGLLY_M,NGLLZ_M,ispec) - yelm(4)=ystore(1,NGLLY_M,NGLLZ_M,ispec) - zelm(4)=zstore(1,NGLLY_M,NGLLZ_M,ispec) - - call compute_jacobian_2D(myrank,ispecb4,xelm,yelm,zelm,dershape2D_y, & - jacobian2D_ymax,normal_ymax,NGLLX_M,NGLLZ_M,NSPEC2DMAX_YMIN_YMAX) - - endif - -! on boundary: bottom - - if (iboun(5,ispec)) then - - ispecb5 = ispecb5+1 - ibelm_bottom(ispecb5) = ispec - - xelm(1)=xstore(1,1,1,ispec) - yelm(1)=ystore(1,1,1,ispec) - zelm(1)=zstore(1,1,1,ispec) - xelm(2)=xstore(NGLLX_M,1,1,ispec) - yelm(2)=ystore(NGLLX_M,1,1,ispec) - zelm(2)=zstore(NGLLX_M,1,1,ispec) - xelm(3)=xstore(NGLLX_M,NGLLY_M,1,ispec) - yelm(3)=ystore(NGLLX_M,NGLLY_M,1,ispec) - zelm(3)=zstore(NGLLX_M,NGLLY_M,1,ispec) - xelm(4)=xstore(1,NGLLY_M,1,ispec) - yelm(4)=ystore(1,NGLLY_M,1,ispec) - zelm(4)=zstore(1,NGLLY_M,1,ispec) - - call compute_jacobian_2D(myrank,ispecb5,xelm,yelm,zelm,dershape2D_bottom, & - jacobian2D_bottom,normal_bottom,NGLLX_M,NGLLY_M,NSPEC2D_BOTTOM) - - endif - -! on boundary: top - - if (iboun(6,ispec)) then - - ispecb6 = ispecb6+1 - ibelm_top(ispecb6) = ispec - - xelm(1)=xstore(1,1,NGLLZ_M,ispec) - yelm(1)=ystore(1,1,NGLLZ_M,ispec) - zelm(1)=zstore(1,1,NGLLZ_M,ispec) - xelm(2)=xstore(NGLLX_M,1,NGLLZ_M,ispec) - yelm(2)=ystore(NGLLX_M,1,NGLLZ_M,ispec) - zelm(2)=zstore(NGLLX_M,1,NGLLZ_M,ispec) - xelm(3)=xstore(NGLLX_M,NGLLY_M,NGLLZ_M,ispec) - yelm(3)=ystore(NGLLX_M,NGLLY_M,NGLLZ_M,ispec) - zelm(3)=zstore(NGLLX_M,NGLLY_M,NGLLZ_M,ispec) - xelm(4)=xstore(1,NGLLY_M,NGLLZ_M,ispec) - yelm(4)=ystore(1,NGLLY_M,NGLLZ_M,ispec) - zelm(4)=zstore(1,NGLLY_M,NGLLZ_M,ispec) - - call compute_jacobian_2D(myrank,ispecb6,xelm,yelm,zelm,dershape2D_top, & - jacobian2D_top,normal_top,NGLLX_M,NGLLY_M,NSPEC2D_TOP) - - endif - - enddo - -! check theoretical value of elements at the bottom - if (ispecb5 /= NSPEC2D_BOTTOM) call exit_MPI(myrank,'ispecb5 should equal NSPEC2D_BOTTOM') - -! check theoretical value of elements at the top - if (ispecb6 /= NSPEC2D_TOP) call exit_MPI(myrank,'ispecb6 should equal NSPEC2D_TOP') - - nspec2D_xmin = ispecb1 - nspec2D_xmax = ispecb2 - nspec2D_ymin = ispecb3 - nspec2D_ymax = ispecb4 - - end subroutine get_jacobian_boundaries +! unused routine, left for reference... + +! subroutine get_jacobian_boundaries(iboun,nspec,xstore,ystore,zstore, & +! dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & +! ibelm_xmin,ibelm_xmax,ibelm_ymin,ibelm_ymax,ibelm_bottom,ibelm_top, & +! nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax, & +! jacobian2D_xmin,jacobian2D_xmax, & +! jacobian2D_ymin,jacobian2D_ymax, & +! jacobian2D_bottom,jacobian2D_top, & +! normal_xmin,normal_xmax, & +! normal_ymin,normal_ymax, & +! normal_bottom,normal_top, & +! NSPEC2D_BOTTOM,NSPEC2D_TOP, & +! NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX) +! +! use constants +! use constants_meshfem3D, only: NGLLX_M,NGLLY_M,NGLLZ_M +! +! implicit none +! +! integer :: nspec +! integer :: NSPEC2D_BOTTOM,NSPEC2D_TOP,NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX +! +! integer :: nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax +! integer :: ibelm_xmin(NSPEC2DMAX_XMIN_XMAX),ibelm_xmax(NSPEC2DMAX_XMIN_XMAX) +! integer :: ibelm_ymin(NSPEC2DMAX_YMIN_YMAX),ibelm_ymax(NSPEC2DMAX_YMIN_YMAX) +! integer :: ibelm_bottom(NSPEC2D_BOTTOM),ibelm_top(NSPEC2D_TOP) +! +! logical :: iboun(6,nspec) +! +! double precision :: xstore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) +! double precision :: ystore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) +! double precision :: zstore(NGLLX_M,NGLLY_M,NGLLZ_M,nspec) +! +! real(kind=CUSTOM_REAL) :: jacobian2D_xmin(NGLLY_M,NGLLZ_M,NSPEC2DMAX_XMIN_XMAX) +! real(kind=CUSTOM_REAL) :: jacobian2D_xmax(NGLLY_M,NGLLZ_M,NSPEC2DMAX_XMIN_XMAX) +! real(kind=CUSTOM_REAL) :: jacobian2D_ymin(NGLLX_M,NGLLZ_M,NSPEC2DMAX_YMIN_YMAX) +! real(kind=CUSTOM_REAL) :: jacobian2D_ymax(NGLLX_M,NGLLZ_M,NSPEC2DMAX_YMIN_YMAX) +! real(kind=CUSTOM_REAL) :: jacobian2D_bottom(NGLLX_M,NGLLY_M,NSPEC2D_BOTTOM) +! real(kind=CUSTOM_REAL) :: jacobian2D_top(NGLLX_M,NGLLY_M,NSPEC2D_TOP) +! +! real(kind=CUSTOM_REAL) :: normal_xmin(NDIM,NGLLY_M,NGLLZ_M,NSPEC2DMAX_XMIN_XMAX) +! real(kind=CUSTOM_REAL) :: normal_xmax(NDIM,NGLLY_M,NGLLZ_M,NSPEC2DMAX_XMIN_XMAX) +! real(kind=CUSTOM_REAL) :: normal_ymin(NDIM,NGLLX_M,NGLLZ_M,NSPEC2DMAX_YMIN_YMAX) +! real(kind=CUSTOM_REAL) :: normal_ymax(NDIM,NGLLX_M,NGLLZ_M,NSPEC2DMAX_YMIN_YMAX) +! real(kind=CUSTOM_REAL) :: normal_bottom(NDIM,NGLLX_M,NGLLY_M,NSPEC2D_BOTTOM) +! real(kind=CUSTOM_REAL) :: normal_top(NDIM,NGLLX_M,NGLLY_M,NSPEC2D_TOP) +! +! double precision :: dershape2D_x(NDIM2D,NGNOD2D_FOUR_CORNERS,NGLLY_M,NGLLZ_M) +! double precision :: dershape2D_y(NDIM2D,NGNOD2D_FOUR_CORNERS,NGLLX_M,NGLLZ_M) +! double precision :: dershape2D_bottom(NDIM2D,NGNOD2D_FOUR_CORNERS,NGLLX_M,NGLLY_M) +! double precision :: dershape2D_top(NDIM2D,NGNOD2D_FOUR_CORNERS,NGLLX_M,NGLLY_M) +! +!! global element numbering +! integer :: ispec +! +!! counters to keep track of number of elements on each of the boundaries +! integer :: ispecb1,ispecb2,ispecb3,ispecb4,ispecb5,ispecb6 +! +! double precision :: xelm(NGNOD2D_FOUR_CORNERS),yelm(NGNOD2D_FOUR_CORNERS),zelm(NGNOD2D_FOUR_CORNERS) +! +! ispecb1 = 0 +! ispecb2 = 0 +! ispecb3 = 0 +! ispecb4 = 0 +! ispecb5 = 0 +! ispecb6 = 0 +! +! do ispec=1,nspec +! +!! determine if the element falls on a boundary +! +!! on boundary: xmin +! +! if (iboun(1,ispec)) then +! +! ispecb1 = ispecb1+1 +! ibelm_xmin(ispecb1) = ispec +! +!! specify the 4 nodes for the 2-D boundary element +! xelm(1)=xstore(1,1,1,ispec) +! yelm(1)=ystore(1,1,1,ispec) +! zelm(1)=zstore(1,1,1,ispec) +! xelm(2)=xstore(1,NGLLY_M,1,ispec) +! yelm(2)=ystore(1,NGLLY_M,1,ispec) +! zelm(2)=zstore(1,NGLLY_M,1,ispec) +! xelm(3)=xstore(1,NGLLY_M,NGLLZ_M,ispec) +! yelm(3)=ystore(1,NGLLY_M,NGLLZ_M,ispec) +! zelm(3)=zstore(1,NGLLY_M,NGLLZ_M,ispec) +! xelm(4)=xstore(1,1,NGLLZ_M,ispec) +! yelm(4)=ystore(1,1,NGLLZ_M,ispec) +! zelm(4)=zstore(1,1,NGLLZ_M,ispec) +! +! call compute_jacobian_2D(myrank,ispecb1,xelm,yelm,zelm,dershape2D_x, & +! jacobian2D_xmin,normal_xmin,NGLLY_M,NGLLZ_M,NSPEC2DMAX_XMIN_XMAX) +! +! endif +! +!! on boundary: xmax +! +! if (iboun(2,ispec)) then +! +! ispecb2 = ispecb2+1 +! ibelm_xmax(ispecb2) = ispec +! +!! specify the 4 nodes for the 2-D boundary element +! xelm(1)=xstore(NGLLX_M,1,1,ispec) +! yelm(1)=ystore(NGLLX_M,1,1,ispec) +! zelm(1)=zstore(NGLLX_M,1,1,ispec) +! xelm(2)=xstore(NGLLX_M,NGLLY_M,1,ispec) +! yelm(2)=ystore(NGLLX_M,NGLLY_M,1,ispec) +! zelm(2)=zstore(NGLLX_M,NGLLY_M,1,ispec) +! xelm(3)=xstore(NGLLX_M,NGLLY_M,NGLLZ_M,ispec) +! yelm(3)=ystore(NGLLX_M,NGLLY_M,NGLLZ_M,ispec) +! zelm(3)=zstore(NGLLX_M,NGLLY_M,NGLLZ_M,ispec) +! xelm(4)=xstore(NGLLX_M,1,NGLLZ_M,ispec) +! yelm(4)=ystore(NGLLX_M,1,NGLLZ_M,ispec) +! zelm(4)=zstore(NGLLX_M,1,NGLLZ_M,ispec) +! +! call compute_jacobian_2D(myrank,ispecb2,xelm,yelm,zelm,dershape2D_x, & +! jacobian2D_xmax,normal_xmax,NGLLY_M,NGLLZ_M,NSPEC2DMAX_XMIN_XMAX) +! +! endif +! +!! on boundary: ymin +! +! if (iboun(3,ispec)) then +! +! ispecb3 = ispecb3+1 +! ibelm_ymin(ispecb3) = ispec +! +!! specify the 4 nodes for the 2-D boundary element +! xelm(1)=xstore(1,1,1,ispec) +! yelm(1)=ystore(1,1,1,ispec) +! zelm(1)=zstore(1,1,1,ispec) +! xelm(2)=xstore(NGLLX_M,1,1,ispec) +! yelm(2)=ystore(NGLLX_M,1,1,ispec) +! zelm(2)=zstore(NGLLX_M,1,1,ispec) +! xelm(3)=xstore(NGLLX_M,1,NGLLZ_M,ispec) +! yelm(3)=ystore(NGLLX_M,1,NGLLZ_M,ispec) +! zelm(3)=zstore(NGLLX_M,1,NGLLZ_M,ispec) +! xelm(4)=xstore(1,1,NGLLZ_M,ispec) +! yelm(4)=ystore(1,1,NGLLZ_M,ispec) +! zelm(4)=zstore(1,1,NGLLZ_M,ispec) +! +! call compute_jacobian_2D(myrank,ispecb3,xelm,yelm,zelm,dershape2D_y, & +! jacobian2D_ymin,normal_ymin,NGLLX_M,NGLLZ_M,NSPEC2DMAX_YMIN_YMAX) +! +! endif +! +!! on boundary: ymax +! +! if (iboun(4,ispec)) then +! +! ispecb4 = ispecb4+1 +! ibelm_ymax(ispecb4) = ispec +! +!! specify the 4 nodes for the 2-D boundary element +! xelm(1)=xstore(1,NGLLY_M,1,ispec) +! yelm(1)=ystore(1,NGLLY_M,1,ispec) +! zelm(1)=zstore(1,NGLLY_M,1,ispec) +! xelm(2)=xstore(NGLLX_M,NGLLY_M,1,ispec) +! yelm(2)=ystore(NGLLX_M,NGLLY_M,1,ispec) +! zelm(2)=zstore(NGLLX_M,NGLLY_M,1,ispec) +! xelm(3)=xstore(NGLLX_M,NGLLY_M,NGLLZ_M,ispec) +! yelm(3)=ystore(NGLLX_M,NGLLY_M,NGLLZ_M,ispec) +! zelm(3)=zstore(NGLLX_M,NGLLY_M,NGLLZ_M,ispec) +! xelm(4)=xstore(1,NGLLY_M,NGLLZ_M,ispec) +! yelm(4)=ystore(1,NGLLY_M,NGLLZ_M,ispec) +! zelm(4)=zstore(1,NGLLY_M,NGLLZ_M,ispec) +! +! call compute_jacobian_2D(myrank,ispecb4,xelm,yelm,zelm,dershape2D_y, & +! jacobian2D_ymax,normal_ymax,NGLLX_M,NGLLZ_M,NSPEC2DMAX_YMIN_YMAX) +! +! endif +! +!! on boundary: bottom +! +! if (iboun(5,ispec)) then +! +! ispecb5 = ispecb5+1 +! ibelm_bottom(ispecb5) = ispec +! +! xelm(1)=xstore(1,1,1,ispec) +! yelm(1)=ystore(1,1,1,ispec) +! zelm(1)=zstore(1,1,1,ispec) +! xelm(2)=xstore(NGLLX_M,1,1,ispec) +! yelm(2)=ystore(NGLLX_M,1,1,ispec) +! zelm(2)=zstore(NGLLX_M,1,1,ispec) +! xelm(3)=xstore(NGLLX_M,NGLLY_M,1,ispec) +! yelm(3)=ystore(NGLLX_M,NGLLY_M,1,ispec) +! zelm(3)=zstore(NGLLX_M,NGLLY_M,1,ispec) +! xelm(4)=xstore(1,NGLLY_M,1,ispec) +! yelm(4)=ystore(1,NGLLY_M,1,ispec) +! zelm(4)=zstore(1,NGLLY_M,1,ispec) +! +! call compute_jacobian_2D(myrank,ispecb5,xelm,yelm,zelm,dershape2D_bottom, & +! jacobian2D_bottom,normal_bottom,NGLLX_M,NGLLY_M,NSPEC2D_BOTTOM) +! +! endif +! +!! on boundary: top +! +! if (iboun(6,ispec)) then +! +! ispecb6 = ispecb6+1 +! ibelm_top(ispecb6) = ispec +! +! xelm(1)=xstore(1,1,NGLLZ_M,ispec) +! yelm(1)=ystore(1,1,NGLLZ_M,ispec) +! zelm(1)=zstore(1,1,NGLLZ_M,ispec) +! xelm(2)=xstore(NGLLX_M,1,NGLLZ_M,ispec) +! yelm(2)=ystore(NGLLX_M,1,NGLLZ_M,ispec) +! zelm(2)=zstore(NGLLX_M,1,NGLLZ_M,ispec) +! xelm(3)=xstore(NGLLX_M,NGLLY_M,NGLLZ_M,ispec) +! yelm(3)=ystore(NGLLX_M,NGLLY_M,NGLLZ_M,ispec) +! zelm(3)=zstore(NGLLX_M,NGLLY_M,NGLLZ_M,ispec) +! xelm(4)=xstore(1,NGLLY_M,NGLLZ_M,ispec) +! yelm(4)=ystore(1,NGLLY_M,NGLLZ_M,ispec) +! zelm(4)=zstore(1,NGLLY_M,NGLLZ_M,ispec) +! +! call compute_jacobian_2D(myrank,ispecb6,xelm,yelm,zelm,dershape2D_top, & +! jacobian2D_top,normal_top,NGLLX_M,NGLLY_M,NSPEC2D_TOP) +! +! endif +! +! enddo +! +!! check theoretical value of elements at the bottom +! if (ispecb5 /= NSPEC2D_BOTTOM) call exit_MPI(myrank,'ispecb5 should equal NSPEC2D_BOTTOM') +! +!! check theoretical value of elements at the top +! if (ispecb6 /= NSPEC2D_TOP) call exit_MPI(myrank,'ispecb6 should equal NSPEC2D_TOP') +! +! nspec2D_xmin = ispecb1 +! nspec2D_xmax = ispecb2 +! nspec2D_ymin = ispecb3 +! nspec2D_ymax = ispecb4 +! +! end subroutine get_jacobian_boundaries +! ! ------------------------------------------------------- +! - subroutine compute_jacobian_2D(myrank,ispecb,xelm,yelm,zelm,dershape2D,jacobian2D,normal,NGLLA,NGLLB,NSPEC2DMAX_AB) - - use constants - - implicit none - -! generic routine that accepts any polynomial degree in each direction - - integer ispecb,NGLLA,NGLLB,NSPEC2DMAX_AB,myrank - - double precision xelm(NGNOD2D_FOUR_CORNERS),yelm(NGNOD2D_FOUR_CORNERS),zelm(NGNOD2D_FOUR_CORNERS) - double precision dershape2D(NDIM2D,NGNOD2D_FOUR_CORNERS,NGLLA,NGLLB) - - real(kind=CUSTOM_REAL) jacobian2D(NGLLA,NGLLB,NSPEC2DMAX_AB) - real(kind=CUSTOM_REAL) normal(3,NGLLA,NGLLB,NSPEC2DMAX_AB) - - integer i,j,ia - double precision xxi,xeta,yxi,yeta,zxi,zeta - double precision unx,uny,unz,jacobian - - do j = 1,NGLLB - do i=1,NGLLA - xxi = ZERO - xeta = ZERO - yxi = ZERO - yeta = ZERO - zxi = ZERO - zeta = ZERO - do ia = 1,NGNOD2D_FOUR_CORNERS - xxi = xxi+dershape2D(1,ia,i,j)*xelm(ia) - xeta = xeta+dershape2D(2,ia,i,j)*xelm(ia) - yxi = yxi+dershape2D(1,ia,i,j)*yelm(ia) - yeta = yeta+dershape2D(2,ia,i,j)*yelm(ia) - zxi = zxi+dershape2D(1,ia,i,j)*zelm(ia) - zeta = zeta+dershape2D(2,ia,i,j)*zelm(ia) - enddo - - ! calculate the unnormalized normal to the boundary - unx = yxi*zeta-yeta*zxi - uny = zxi*xeta-zeta*xxi - unz = xxi*yeta-xeta*yxi - jacobian = dsqrt(unx**2+uny**2+unz**2) - if (jacobian <= ZERO) call exit_MPI(myrank,'2D Jacobian undefined') - - ! normalize normal vector and store surface jacobian - ! distinguish if single or double precision for reals - jacobian2D(i,j,ispecb) = real(jacobian,kind=CUSTOM_REAL) - normal(1,i,j,ispecb) = real(unx/jacobian,kind=CUSTOM_REAL) - normal(2,i,j,ispecb) = real(uny/jacobian,kind=CUSTOM_REAL) - normal(3,i,j,ispecb) = real(unz/jacobian,kind=CUSTOM_REAL) - - enddo - enddo +! unused routine, left for reference... - end subroutine compute_jacobian_2D +! subroutine compute_jacobian_2D(ispecb,xelm,yelm,zelm,dershape2D,jacobian2D,normal,NGLLA,NGLLB,NSPEC2DMAX_AB) +! +! use constants +! +! implicit none +! +!! generic routine that accepts any polynomial degree in each direction +! +! integer :: ispecb,NGLLA,NGLLB,NSPEC2DMAX_AB +! +! double precision :: xelm(NGNOD2D_FOUR_CORNERS),yelm(NGNOD2D_FOUR_CORNERS),zelm(NGNOD2D_FOUR_CORNERS) +! double precision :: dershape2D(NDIM2D,NGNOD2D_FOUR_CORNERS,NGLLA,NGLLB) +! +! real(kind=CUSTOM_REAL) :: jacobian2D(NGLLA,NGLLB,NSPEC2DMAX_AB) +! real(kind=CUSTOM_REAL) :: normal(3,NGLLA,NGLLB,NSPEC2DMAX_AB) +! +! integer :: i,j,ia +! double precision :: xxi,xeta,yxi,yeta,zxi,zeta +! double precision :: unx,uny,unz,jacobian +! +! do j = 1,NGLLB +! do i=1,NGLLA +! xxi = ZERO +! xeta = ZERO +! yxi = ZERO +! yeta = ZERO +! zxi = ZERO +! zeta = ZERO +! do ia = 1,NGNOD2D_FOUR_CORNERS +! xxi = xxi+dershape2D(1,ia,i,j)*xelm(ia) +! xeta = xeta+dershape2D(2,ia,i,j)*xelm(ia) +! yxi = yxi+dershape2D(1,ia,i,j)*yelm(ia) +! yeta = yeta+dershape2D(2,ia,i,j)*yelm(ia) +! zxi = zxi+dershape2D(1,ia,i,j)*zelm(ia) +! zeta = zeta+dershape2D(2,ia,i,j)*zelm(ia) +! enddo +! +! ! calculate the unnormalized normal to the boundary +! unx = yxi*zeta-yeta*zxi +! uny = zxi*xeta-zeta*xxi +! unz = xxi*yeta-xeta*yxi +! jacobian = dsqrt(unx**2+uny**2+unz**2) +! if (jacobian <= ZERO) call exit_MPI(myrank,'2D Jacobian undefined') +! +! ! normalize normal vector and store surface jacobian +! ! distinguish if single or double precision for reals +! jacobian2D(i,j,ispecb) = real(jacobian,kind=CUSTOM_REAL) +! normal(1,i,j,ispecb) = real(unx/jacobian,kind=CUSTOM_REAL) +! normal(2,i,j,ispecb) = real(uny/jacobian,kind=CUSTOM_REAL) +! normal(3,i,j,ispecb) = real(unz/jacobian,kind=CUSTOM_REAL) +! +! enddo +! enddo +! +! end subroutine compute_jacobian_2D diff --git a/src/shared/check_mesh_resolution.f90 b/src/shared/check_mesh_resolution.f90 index a3e398ef0..c85853e29 100644 --- a/src/shared/check_mesh_resolution.f90 +++ b/src/shared/check_mesh_resolution.f90 @@ -25,10 +25,10 @@ ! !===================================================================== - subroutine check_mesh_resolution(myrank,NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, & - kappastore,mustore,rho_vp,rho_vs, & - DT, model_speed_max,min_resolved_period, & - LOCAL_PATH,SAVE_MESH_FILES) + subroutine check_mesh_resolution(NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, & + kappastore,mustore,rho_vp,rho_vs, & + DT, model_speed_max,min_resolved_period, & + LOCAL_PATH,SAVE_MESH_FILES) ! check the mesh, stability and resolved period for acoustic / elastic domains ! @@ -62,7 +62,6 @@ subroutine check_mesh_resolution(myrank,NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zs logical:: DT_PRESENT - integer :: myrank integer :: NSPEC_AB_global_min,NSPEC_AB_global_max,NSPEC_AB_global_sum integer :: NGLOB_AB_global_min,NGLOB_AB_global_max,NGLOB_AB_global_sum integer :: ispec @@ -466,11 +465,11 @@ end subroutine check_mesh_resolution !------------------------------------------------------------------------------------------------- ! - subroutine check_mesh_resolution_poro(myrank,NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, & - DT, model_speed_max,min_resolved_period, & - phistore,tortstore,rhoarraystore, & - rho_vpI,rho_vpII,rho_vsI, & - LOCAL_PATH,SAVE_MESH_FILES) + subroutine check_mesh_resolution_poro(NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, & + DT, model_speed_max,min_resolved_period, & + phistore,tortstore,rhoarraystore, & + rho_vpI,rho_vpII,rho_vsI, & + LOCAL_PATH,SAVE_MESH_FILES) ! check the mesh, stability and resolved period for poroelastic domains ! @@ -502,7 +501,6 @@ subroutine check_mesh_resolution_poro(myrank,NSPEC_AB,NGLOB_AB,ibool,xstore,ysto logical:: DT_PRESENT - integer :: myrank integer :: NSPEC_AB_global_min,NSPEC_AB_global_max,NSPEC_AB_global_sum integer :: NGLOB_AB_global_min,NGLOB_AB_global_max,NGLOB_AB_global_sum integer :: ispec diff --git a/src/shared/detect_surface.f90 b/src/shared/detect_surface.f90 index cccdfc4d2..aa2023279 100644 --- a/src/shared/detect_surface.f90 +++ b/src/shared/detect_surface.f90 @@ -214,16 +214,16 @@ end subroutine ds_set_surface_flags ! subroutine detect_surface_cross_section(NPROC,nglob,nspec,ibool, & - ispec_is_surface_cross_section, & - iglob_is_surface_cross_section, & - nfaces_surface, & - num_interfaces_ext_mesh, & - max_nibool_interfaces_ext_mesh, & - nibool_interfaces_ext_mesh, & - my_neighbors_ext_mesh, & - ibool_interfaces_ext_mesh, & - x_section,y_section,z_section, & - xstore,ystore,zstore,myrank) + ispec_is_surface_cross_section, & + iglob_is_surface_cross_section, & + nfaces_surface, & + num_interfaces_ext_mesh, & + max_nibool_interfaces_ext_mesh, & + nibool_interfaces_ext_mesh, & + my_neighbors_ext_mesh, & + ibool_interfaces_ext_mesh, & + x_section,y_section,z_section, & + xstore,ystore,zstore) ! instead of surface of model, this returns cross-section surfaces through model ! at specified x,y,z - coordinates @@ -238,7 +238,7 @@ subroutine detect_surface_cross_section(NPROC,nglob,nspec,ibool, & implicit none ! global indexing - integer :: NPROC,nglob,nspec,myrank + integer :: NPROC,nglob,nspec integer, dimension(NGLLX,NGLLY,NGLLZ,nspec):: ibool ! surface diff --git a/src/shared/exit_mpi.f90 b/src/shared/exit_mpi.f90 index 87a686c69..182b6ce39 100644 --- a/src/shared/exit_mpi.f90 +++ b/src/shared/exit_mpi.f90 @@ -29,15 +29,15 @@ subroutine exit_MPI(myrank,error_msg) - use constants + use constants, only: MAX_STRING_LEN,IMAIN,ISTANDARD_OUTPUT,OUTPUT_FILES implicit none ! identifier for error message file integer, parameter :: IERROR = 30 - integer myrank - character(len=*) :: error_msg + integer, intent(in) :: myrank + character(len=*),intent(in) :: error_msg character(len=MAX_STRING_LEN) :: outputname diff --git a/src/shared/get_attenuation_model.f90 b/src/shared/get_attenuation_model.f90 index a18feed87..845c32b39 100644 --- a/src/shared/get_attenuation_model.f90 +++ b/src/shared/get_attenuation_model.f90 @@ -124,7 +124,7 @@ end subroutine get_attenuation_model_olsen !------------------------------------------------------------------------------------------------- ! - subroutine get_attenuation_model(myrank,nspec,USE_OLSEN_ATTENUATION,OLSEN_ATTENUATION_RATIO, & + subroutine get_attenuation_model(nspec,USE_OLSEN_ATTENUATION,OLSEN_ATTENUATION_RATIO, & mustore,rho_vs,kappastore,rho_vp,qkappa_attenuation_store,qmu_attenuation_store, & ispec_is_elastic,min_resolved_period,prname,ATTENUATION_f0_REFERENCE) @@ -137,7 +137,7 @@ subroutine get_attenuation_model(myrank,nspec,USE_OLSEN_ATTENUATION,OLSEN_ATTENU implicit none double precision,intent(in) :: OLSEN_ATTENUATION_RATIO,ATTENUATION_f0_REFERENCE - integer,intent(in) :: myrank,nspec + integer,intent(in) :: nspec logical,intent(in) :: USE_OLSEN_ATTENUATION real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec),intent(in) :: mustore @@ -347,7 +347,7 @@ subroutine get_attenuation_model(myrank,nspec,USE_OLSEN_ATTENUATION,OLSEN_ATTENU ! gets beta, on_minus_sum_beta and factor_scale ! based on calculation of strain relaxation times tau_eps - call get_attenuation_factors(myrank,Q_mu,MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD, & + call get_attenuation_factors(Q_mu,MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD, & f_c_source,tau_sigma_dble, & beta_dble,one_minus_sum_beta_dble,factor_scale_dble, & Q_kappa,beta_dble_kappa,one_minus_sum_beta_dble_kappa,factor_scale_dble_kappa, & @@ -492,7 +492,7 @@ end subroutine get_attenuation_constants !------------------------------------------------------------------------------------------------- ! - subroutine get_attenuation_factors(myrank,Q_mu,MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD, & + subroutine get_attenuation_factors(Q_mu,MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD, & f_c_source,tau_sigma, & beta,one_minus_sum_beta,factor_scale, & Q_kappa,beta_kappa,one_minus_sum_beta_kappa,factor_scale_kappa,ATTENUATION_f0_REFERENCE) @@ -510,7 +510,6 @@ subroutine get_attenuation_factors(myrank,Q_mu,MIN_ATTENUATION_PERIOD,MAX_ATTENU implicit none - integer:: myrank double precision :: MIN_ATTENUATION_PERIOD,MAX_ATTENUATION_PERIOD,ATTENUATION_f0_REFERENCE double precision :: f_c_source,Q_mu,Q_kappa double precision, dimension(N_SLS) :: tau_sigma @@ -528,7 +527,7 @@ subroutine get_attenuation_factors(myrank,Q_mu,MIN_ATTENUATION_PERIOD,MAX_ATTENU call get_attenuation_property_values(tau_sigma,tau_eps_kappa,beta_kappa,one_minus_sum_beta_kappa) ! determines the "scale factor" - call get_attenuation_scale_factor(myrank,f_c_source,tau_eps_kappa,tau_sigma,Q_kappa,factor_scale_kappa,ATTENUATION_f0_REFERENCE) + call get_attenuation_scale_factor(f_c_source,tau_eps_kappa,tau_sigma,Q_kappa,factor_scale_kappa,ATTENUATION_f0_REFERENCE) ! uncomment this to print the constants to use in the 3D viscoelastic analytical code for validation purposes ! if (myrank == 0) & ! print *,'for Q_Kappa,tau_eps_kappa,tau_sigma,factor_scale_kappa = ',Q_Kappa,tau_eps_kappa(:),tau_sigma(:),factor_scale_kappa @@ -540,7 +539,7 @@ subroutine get_attenuation_factors(myrank,Q_mu,MIN_ATTENUATION_PERIOD,MAX_ATTENU call get_attenuation_property_values(tau_sigma,tau_eps,beta,one_minus_sum_beta) ! determines the "scale factor" - call get_attenuation_scale_factor(myrank,f_c_source,tau_eps,tau_sigma,Q_mu,factor_scale,ATTENUATION_f0_REFERENCE) + call get_attenuation_scale_factor(f_c_source,tau_eps,tau_sigma,Q_mu,factor_scale,ATTENUATION_f0_REFERENCE) ! uncomment this to print the constants to use in the 3D viscoelastic analytical code for validation purposes ! if (myrank == 0) & ! print *,'for Q_mu,tau_eps,tau_sigma,factor_scale = ',Q_mu,tau_eps(:),tau_sigma(:),factor_scale @@ -587,7 +586,7 @@ end subroutine get_attenuation_property_values !------------------------------------------------------------------------------------------------- ! - subroutine get_attenuation_scale_factor(myrank,f_c_source,tau_eps,tau_sigma,Q_val,scale_factor,ATTENUATION_f0_REFERENCE) + subroutine get_attenuation_scale_factor(f_c_source,tau_eps,tau_sigma,Q_val,scale_factor,ATTENUATION_f0_REFERENCE) ! returns: physical dispersion scaling factor scale_factor @@ -595,7 +594,6 @@ subroutine get_attenuation_scale_factor(myrank,f_c_source,tau_eps,tau_sigma,Q_va implicit none - integer :: myrank double precision, intent(in) :: ATTENUATION_f0_REFERENCE double precision :: scale_factor, Q_val, f_c_source ! strain and stress relaxation times diff --git a/src/shared/get_jacobian_boundaries.f90 b/src/shared/get_jacobian_boundaries.f90 index 7639ef397..9f95d824f 100644 --- a/src/shared/get_jacobian_boundaries.f90 +++ b/src/shared/get_jacobian_boundaries.f90 @@ -25,11 +25,11 @@ ! !===================================================================== - subroutine get_jacobian_boundary_face(myrank,nspec, & - xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob, & - dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & - wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & - ispec,iface,jacobian2Dw_face,normal_face,NGLLA,NGLLB,NGNOD2D) + subroutine get_jacobian_boundary_face(nspec, & + xstore_dummy,ystore_dummy,zstore_dummy,ibool,nglob, & + dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top, & + wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, & + ispec,iface,jacobian2Dw_face,normal_face,NGLLA,NGLLB,NGNOD2D) ! returns jacobian2Dw_face and normal_face (pointing outwards of element) @@ -37,7 +37,7 @@ subroutine get_jacobian_boundary_face(myrank,nspec, & implicit none - integer :: nspec,myrank,nglob + integer :: nspec,nglob ! arrays with the mesh integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: ibool @@ -100,9 +100,9 @@ subroutine get_jacobian_boundary_face(myrank,nspec, & zelm(9)=zstore_dummy( ibool(1,(NGLLY+1)/2,(NGLLZ+1)/2,ispec) ) endif - call compute_jacobian_2D_face(myrank,xelm,yelm,zelm, & - dershape2D_x,wgllwgll_yz, & - jacobian2Dw_face,normal_face,NGLLY,NGLLZ,NGNOD2D) + call compute_jacobian_2D_face(xelm,yelm,zelm, & + dershape2D_x,wgllwgll_yz, & + jacobian2Dw_face,normal_face,NGLLY,NGLLZ,NGNOD2D) ! on boundary: xmax case (2) @@ -137,9 +137,9 @@ subroutine get_jacobian_boundary_face(myrank,nspec, & zelm(9)=zstore_dummy( ibool(NGLLX,(NGLLY+1)/2,(NGLLZ+1)/2,ispec) ) endif - call compute_jacobian_2D_face(myrank,xelm,yelm,zelm, & - dershape2D_x,wgllwgll_yz, & - jacobian2Dw_face,normal_face,NGLLY,NGLLZ,NGNOD2D) + call compute_jacobian_2D_face(xelm,yelm,zelm, & + dershape2D_x,wgllwgll_yz, & + jacobian2Dw_face,normal_face,NGLLY,NGLLZ,NGNOD2D) ! on boundary: ymin case (3) @@ -174,9 +174,9 @@ subroutine get_jacobian_boundary_face(myrank,nspec, & zelm(9)=zstore_dummy( ibool((NGLLX+1)/2,1,(NGLLZ+1)/2,ispec) ) endif - call compute_jacobian_2D_face(myrank,xelm,yelm,zelm, & - dershape2D_y,wgllwgll_xz, & - jacobian2Dw_face,normal_face,NGLLX,NGLLZ,NGNOD2D) + call compute_jacobian_2D_face(xelm,yelm,zelm, & + dershape2D_y,wgllwgll_xz, & + jacobian2Dw_face,normal_face,NGLLX,NGLLZ,NGNOD2D) ! on boundary: ymax case (4) @@ -211,9 +211,9 @@ subroutine get_jacobian_boundary_face(myrank,nspec, & zelm(9)=zstore_dummy( ibool((NGLLX+1)/2,NGLLY,(NGLLZ+1)/2,ispec) ) endif - call compute_jacobian_2D_face(myrank,xelm,yelm,zelm, & - dershape2D_y, wgllwgll_xz, & - jacobian2Dw_face,normal_face,NGLLX,NGLLZ,NGNOD2D) + call compute_jacobian_2D_face(xelm,yelm,zelm, & + dershape2D_y, wgllwgll_xz, & + jacobian2Dw_face,normal_face,NGLLX,NGLLZ,NGNOD2D) ! on boundary: bottom @@ -249,9 +249,9 @@ subroutine get_jacobian_boundary_face(myrank,nspec, & zelm(9)=zstore_dummy( ibool((NGLLX+1)/2,(NGLLY+1)/2,1,ispec) ) endif - call compute_jacobian_2D_face(myrank,xelm,yelm,zelm, & - dershape2D_bottom,wgllwgll_xy, & - jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D) + call compute_jacobian_2D_face(xelm,yelm,zelm, & + dershape2D_bottom,wgllwgll_xy, & + jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D) ! on boundary: top case (6) @@ -286,9 +286,9 @@ subroutine get_jacobian_boundary_face(myrank,nspec, & zelm(9)=zstore_dummy( ibool((NGLLX+1)/2,(NGLLY+1)/2,NGLLZ,ispec) ) endif - call compute_jacobian_2D_face(myrank,xelm,yelm,zelm, & - dershape2D_top, wgllwgll_xy, & - jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D) + call compute_jacobian_2D_face(xelm,yelm,zelm, & + dershape2D_top, wgllwgll_xy, & + jacobian2Dw_face,normal_face,NGLLX,NGLLY,NGNOD2D) case default stop 'error 2D jacobian' @@ -300,9 +300,9 @@ end subroutine get_jacobian_boundary_face ! ------------------------------------------------------------------------------------------------ ! - subroutine compute_jacobian_2D_face(myrank,xelm,yelm,zelm, & - dershape2D,wgllwgll, & - jacobian2Dw_face,normal_face,NGLLA,NGLLB,NGNOD2D) + subroutine compute_jacobian_2D_face(xelm,yelm,zelm, & + dershape2D,wgllwgll, & + jacobian2Dw_face,normal_face,NGLLA,NGLLB,NGNOD2D) use constants @@ -311,18 +311,18 @@ subroutine compute_jacobian_2D_face(myrank,xelm,yelm,zelm, & ! generic routine that accepts any polynomial degree in each direction ! returns 2D jacobian and normal for this face only - integer NGLLA,NGLLB,NGNOD2D,myrank + integer :: NGLLA,NGLLB,NGNOD2D - double precision xelm(NGNOD2D),yelm(NGNOD2D),zelm(NGNOD2D) - double precision dershape2D(NDIM2D,NGNOD2D,NGLLA,NGLLB) + double precision :: xelm(NGNOD2D),yelm(NGNOD2D),zelm(NGNOD2D) + double precision :: dershape2D(NDIM2D,NGNOD2D,NGLLA,NGLLB) double precision, dimension(NGLLA,NGLLB) :: wgllwgll real(kind=CUSTOM_REAL), dimension(NGLLA,NGLLB) :: jacobian2Dw_face real(kind=CUSTOM_REAL), dimension(NDIM,NGLLA,NGLLB) :: normal_face - integer i,j,ia - double precision xxi,xeta,yxi,yeta,zxi,zeta - double precision unx,uny,unz,jacobian + integer :: i,j,ia + double precision :: xxi,xeta,yxi,yeta,zxi,zeta + double precision :: unx,uny,unz,jacobian do j = 1,NGLLB do i = 1,NGLLA diff --git a/src/shared/get_shape2D.f90 b/src/shared/get_shape2D.f90 index 6aeac4f27..4f60f5f37 100644 --- a/src/shared/get_shape2D.f90 +++ b/src/shared/get_shape2D.f90 @@ -25,7 +25,7 @@ ! !===================================================================== - subroutine get_shape2D(myrank,shape2D,dershape2D,xigll,yigll,NGLLA,NGLLB,NGNOD,NGNOD2D) + subroutine get_shape2D(shape2D,dershape2D,xigll,yigll,NGLLA,NGLLB,NGNOD,NGNOD2D) use constants @@ -33,23 +33,23 @@ subroutine get_shape2D(myrank,shape2D,dershape2D,xigll,yigll,NGLLA,NGLLB,NGNOD,N ! generic routine that accepts any polynomial degree in each direction - integer NGLLA,NGLLB,NGNOD,NGNOD2D,myrank + integer :: NGLLA,NGLLB,NGNOD,NGNOD2D - double precision xigll(NGLLA) - double precision yigll(NGLLB) + double precision :: xigll(NGLLA) + double precision :: yigll(NGLLB) ! 2D shape functions and their derivatives - double precision shape2D(NGNOD2D,NGLLA,NGLLB) - double precision dershape2D(NDIM2D,NGNOD2D,NGLLA,NGLLB) + double precision :: shape2D(NGNOD2D,NGLLA,NGLLB) + double precision :: dershape2D(NDIM2D,NGNOD2D,NGLLA,NGLLB) - integer i,j,ia + integer :: i,j,ia ! location of the nodes of the 2D quadrilateral elements - double precision xi,eta - double precision xi_map,eta_map + double precision :: xi,eta + double precision :: xi_map,eta_map ! for checking the 2D shape functions - double precision sumshape,sumdershapexi,sumdershapeeta + double precision :: sumshape,sumdershapexi,sumdershapeeta ! check that the parameter file is correct if (NGNOD /= 8 .and. NGNOD /= 27) & diff --git a/src/shared/get_shape3D.f90 b/src/shared/get_shape3D.f90 index 37fc27677..0686addf7 100644 --- a/src/shared/get_shape3D.f90 +++ b/src/shared/get_shape3D.f90 @@ -27,31 +27,31 @@ ! 3D shape functions for 8-node or 27-node element - subroutine get_shape3D(myrank,shape3D,dershape3D,xigll,yigll,zigll,NGNOD) + subroutine get_shape3D(shape3D,dershape3D,xigll,yigll,zigll,NGNOD) use constants implicit none - integer NGNOD,myrank + integer :: NGNOD ! Gauss-Lobatto-Legendre points of integration - double precision xigll(NGLLX) - double precision yigll(NGLLY) - double precision zigll(NGLLZ) + double precision :: xigll(NGLLX) + double precision :: yigll(NGLLY) + double precision :: zigll(NGLLZ) ! 3D shape functions and their derivatives - double precision shape3D(NGNOD,NGLLX,NGLLY,NGLLZ) - double precision dershape3D(NDIM,NGNOD,NGLLX,NGLLY,NGLLZ) + double precision :: shape3D(NGNOD,NGLLX,NGLLY,NGLLZ) + double precision :: dershape3D(NDIM,NGNOD,NGLLX,NGLLY,NGLLZ) - integer i,j,k,ia + integer :: i,j,k,ia ! location of the nodes of the 3D hexahedra elements - double precision xi,eta,gamma - double precision ra1,ra2,rb1,rb2,rc1,rc2 + double precision :: xi,eta,gamma + double precision :: ra1,ra2,rb1,rb2,rc1,rc2 ! for checking the 3D shape functions - double precision sumshape,sumdershapexi,sumdershapeeta,sumdershapegamma + double precision :: sumshape,sumdershapexi,sumdershapeeta,sumdershapegamma double precision, parameter :: ONE_EIGHTH = 0.125d0 @@ -168,13 +168,13 @@ end subroutine get_shape3D ! 3D shape functions for given, single xi/eta/gamma location - subroutine eval_shape3D_single(myrank,shape3D,xi,eta,gamma,NGNOD) + subroutine eval_shape3D_single(shape3D,xi,eta,gamma,NGNOD) use constants implicit none - integer :: myrank,NGNOD + integer :: NGNOD ! 3D shape functions double precision :: shape3D(NGNOD) diff --git a/src/shared/read_parameter_file.F90 b/src/shared/read_parameter_file.F90 index a18a7b9da..7ceef0c2f 100644 --- a/src/shared/read_parameter_file.F90 +++ b/src/shared/read_parameter_file.F90 @@ -25,14 +25,16 @@ ! !===================================================================== - subroutine read_parameter_file(myrank,BROADCAST_AFTER_READ) + subroutine read_parameter_file(BROADCAST_AFTER_READ) - use constants + use constants, only: & + myrank, & + INJECTION_TECHNIQUE_IS_AXISEM,INJECTION_TECHNIQUE_IS_DSM,INJECTION_TECHNIQUE_IS_FK + use shared_parameters implicit none - integer, intent(in) :: myrank logical, intent(in) :: BROADCAST_AFTER_READ ! local variables diff --git a/src/shared/rules.mk b/src/shared/rules.mk index 2929e0cde..1e7f4230d 100644 --- a/src/shared/rules.mk +++ b/src/shared/rules.mk @@ -54,6 +54,7 @@ shared_OBJECTS = \ $O/gll_library.shared.o \ $O/hex_nodes.shared.o \ $O/heap_sort.shared.o \ + $O/init_openmp.shared.o \ $O/lagrange_poly.shared.o \ $O/merge_sort.shared.o \ $O/netlib_specfun_erf.shared.o \ diff --git a/src/shared/serial.f90 b/src/shared/serial.f90 index 1a938e0bf..db625d947 100644 --- a/src/shared/serial.f90 +++ b/src/shared/serial.f90 @@ -489,7 +489,7 @@ subroutine init_mpi() ! we need to make sure that NUMBER_OF_SIMULTANEOUS_RUNS is read, thus read the parameter file myrank = 0 BROADCAST_AFTER_READ = .false. - call read_parameter_file(myrank,BROADCAST_AFTER_READ) + call read_parameter_file(BROADCAST_AFTER_READ) if (NUMBER_OF_SIMULTANEOUS_RUNS <= 0) stop 'NUMBER_OF_SIMULTANEOUS_RUNS <= 0 makes no sense' diff --git a/src/shared/shared_par.F90 b/src/shared/shared_par.F90 index da1859077..12f4897b2 100644 --- a/src/shared/shared_par.F90 +++ b/src/shared/shared_par.F90 @@ -31,6 +31,9 @@ module constants include "constants.h" + ! proc number for MPI process + integer :: myrank + ! a negative initial value is a convention that indicates that groups (i.e. sub-communicators, one per run) are off by default integer :: mygroup = -1 diff --git a/src/specfem3D/assemble_MPI_vector.f90 b/src/specfem3D/assemble_MPI_vector.f90 index e01532b42..b41ad580e 100644 --- a/src/specfem3D/assemble_MPI_vector.f90 +++ b/src/specfem3D/assemble_MPI_vector.f90 @@ -133,9 +133,10 @@ end subroutine assemble_MPI_vector_blocking !------------------------------------------------------------------------------------------------- ! subroutine synchronize_MPI_vector_blocking_ord(NPROC,NGLOB_AB,array_val, & - num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, & - nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, & - my_neighbors_ext_mesh,myrank) + num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, & + nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, & + my_neighbors_ext_mesh) + ! kbai added this subroutine to synchronize a vector field ! to ensure that its values at nodes on MPI interfaces stay equal on all processors that share the node. ! Synchronize by setting the value to that of the processor with highest rank @@ -153,7 +154,7 @@ subroutine synchronize_MPI_vector_blocking_ord(NPROC,NGLOB_AB,array_val, & ! array to assemble real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: array_val - integer :: num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh,myrank + integer :: num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh integer, dimension(num_interfaces_ext_mesh) :: nibool_interfaces_ext_mesh,my_neighbors_ext_mesh integer, dimension(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh) :: ibool_interfaces_ext_mesh ! local parameters @@ -165,9 +166,7 @@ subroutine synchronize_MPI_vector_blocking_ord(NPROC,NGLOB_AB,array_val, & integer, dimension(:), allocatable :: request_send_vector integer, dimension(:), allocatable :: request_recv_vector - integer ipoin,iinterface,ier,iglob - - + integer :: ipoin,iinterface,ier,iglob ! setting the value to that of the processor with highest rank @@ -310,11 +309,11 @@ end subroutine assemble_MPI_vector_async_send ! subroutine assemble_MPI_vector_async_w_ord(NPROC,NGLOB_AB,array_val, & - buffer_recv_vector_ext_mesh,num_interfaces_ext_mesh, & - max_nibool_interfaces_ext_mesh, & - nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, & - request_send_vector_ext_mesh,request_recv_vector_ext_mesh, & - my_neighbors_ext_mesh,myrank) + buffer_recv_vector_ext_mesh,num_interfaces_ext_mesh, & + max_nibool_interfaces_ext_mesh, & + nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, & + request_send_vector_ext_mesh,request_recv_vector_ext_mesh, & + my_neighbors_ext_mesh) ! waits for data to receive and assembles @@ -338,7 +337,7 @@ subroutine assemble_MPI_vector_async_w_ord(NPROC,NGLOB_AB,array_val, & ! array to assemble real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: array_val - integer :: num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh,myrank + integer :: num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh real(kind=CUSTOM_REAL), dimension(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh) :: & buffer_recv_vector_ext_mesh diff --git a/src/specfem3D/compute_add_sources_viscoelastic.F90 b/src/specfem3D/compute_add_sources_viscoelastic.F90 index cbf03897d..b8b6c917a 100644 --- a/src/specfem3D/compute_add_sources_viscoelastic.F90 +++ b/src/specfem3D/compute_add_sources_viscoelastic.F90 @@ -66,6 +66,11 @@ subroutine compute_add_sources_viscoelastic() ! forward simulations if (SIMULATION_TYPE == 1 .and. NOISE_TOMOGRAPHY == 0 .and. nsources_local > 0) then +! openmp solver +!$OMP PARALLEL if (NSOURCES > 100) & +!$OMP DEFAULT(SHARED) & +!$OMP PRIVATE(isource,time_source_dble,iglob,stf_used,stf,ispec,i,j,k) +!$OMP DO do isource = 1,NSOURCES ! add the source (only if this proc carries the source) @@ -97,7 +102,12 @@ subroutine compute_add_sources_viscoelastic() do j=1,NGLLY do i=1,NGLLX iglob = ibool(i,j,k,ispec) - accel(:,iglob) = accel(:,iglob) + sourcearrays(isource,:,i,j,k)*stf_used +!$OMP ATOMIC + accel(1,iglob) = accel(1,iglob) + sourcearrays(isource,1,i,j,k)*stf_used +!$OMP ATOMIC + accel(2,iglob) = accel(2,iglob) + sourcearrays(isource,2,i,j,k)*stf_used +!$OMP ATOMIC + accel(3,iglob) = accel(3,iglob) + sourcearrays(isource,3,i,j,k)*stf_used enddo enddo enddo @@ -105,6 +115,9 @@ subroutine compute_add_sources_viscoelastic() endif ! ispec_is_elastic endif ! myrank enddo ! NSOURCES +!$OMP ENDDO +!$OMP END PARALLEL + endif ! forward ! NOTE: adjoint sources and backward wavefield timing: @@ -216,11 +229,10 @@ subroutine compute_add_sources_viscoelastic() ! hence, instead of a moment tensor 'sourcearrays', a 'noise_sourcearray' for a point force is needed. ! furthermore, the CMTSOLUTION needs to be zero, i.e., no earthquakes. ! now this must be manually set in DATA/CMTSOLUTION, by USERS. - call add_source_master_rec_noise(myrank,nrec, & - NSTEP,accel,noise_sourcearray, & - ibool,islice_selected_rec,ispec_selected_rec, & - it,irec_master_noise, & - NSPEC_AB,NGLOB_AB) + call add_source_master_rec_noise(nrec,NSTEP,accel,noise_sourcearray, & + ibool,islice_selected_rec,ispec_selected_rec, & + it,irec_master_noise, & + NSPEC_AB,NGLOB_AB) else if (NOISE_TOMOGRAPHY == 2) then ! second step of noise tomography, i.e., read the surface movie saved at every timestep ! use the movie to drive the ensemble forward wavefield diff --git a/src/specfem3D/compute_forces_viscoelastic.f90 b/src/specfem3D/compute_forces_viscoelastic.f90 index c9e7b298c..0f9f635fa 100644 --- a/src/specfem3D/compute_forces_viscoelastic.f90 +++ b/src/specfem3D/compute_forces_viscoelastic.f90 @@ -188,9 +188,7 @@ subroutine compute_forces_viscoelastic(iphase, & !ZN ZN Do we need a stop statement somewhere in case of !ZN ZN SIMULATION_TYPE == 3 and allocated(Kelvin_Voigt_eta) = .true. ?? if (allocated(Kelvin_Voigt_eta)) then - eta = Kelvin_Voigt_eta(ispec) - if (is_CPML(ispec) .and. eta /= 0._CUSTOM_REAL) stop 'you cannot put a fault inside a PML layer' do k=1,NGLLZ @@ -203,7 +201,9 @@ subroutine compute_forces_viscoelastic(iphase, & enddo enddo enddo + else + do k=1,NGLLZ do j=1,NGLLY do i=1,NGLLX @@ -214,6 +214,7 @@ subroutine compute_forces_viscoelastic(iphase, & enddo enddo enddo + endif if (is_CPML(ispec)) then diff --git a/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 b/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 index 0e6df6688..0cab7ec18 100644 --- a/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 +++ b/src/specfem3D/compute_forces_viscoelastic_calling_routine.F90 @@ -47,13 +47,13 @@ subroutine compute_forces_viscoelastic_calling() ! Do this only for dynamic rupture simulations if (SIMULATION_TYPE_DYN) then call synchronize_MPI_vector_blocking_ord(NPROC,NGLOB_AB,displ, & - num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, & - nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, & - my_neighbors_ext_mesh,myrank) + num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, & + nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, & + my_neighbors_ext_mesh) call synchronize_MPI_vector_blocking_ord(NPROC,NGLOB_AB,veloc, & - num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, & - nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, & - my_neighbors_ext_mesh,myrank) + num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, & + nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, & + my_neighbors_ext_mesh) endif ! distinguishes two runs: for elements in contact with MPI interfaces, and elements within the partitions @@ -147,25 +147,27 @@ subroutine compute_forces_viscoelastic_calling() else ! waits for send/receive requests to be completed and assembles values call assemble_MPI_vector_async_w_ord(NPROC,NGLOB_AB,accel, & - buffer_recv_vector_ext_mesh,num_interfaces_ext_mesh, & - max_nibool_interfaces_ext_mesh, & - nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, & - request_send_vector_ext_mesh,request_recv_vector_ext_mesh, & - my_neighbors_ext_mesh,myrank) + buffer_recv_vector_ext_mesh,num_interfaces_ext_mesh, & + max_nibool_interfaces_ext_mesh, & + nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, & + request_send_vector_ext_mesh,request_recv_vector_ext_mesh, & + my_neighbors_ext_mesh) endif enddo -! Percy, Fault boundary term B*tau is added to the assembled forces -! which at this point are stored in the array 'accel' + ! Percy, Fault boundary term B*tau is added to the assembled forces + ! which at this point are stored in the array 'accel' if (SIMULATION_TYPE_DYN) call bc_dynflt_set3d_all(accel,veloc,displ) if (SIMULATION_TYPE_KIN) call bc_kinflt_set_all(accel,veloc,displ) - ! multiplies with inverse of mass matrix (note: rmass has been inverted already) - accel(1,:) = accel(1,:)*rmassx(:) - accel(2,:) = accel(2,:)*rmassy(:) - accel(3,:) = accel(3,:)*rmassz(:) + ! multiplies with inverse of mass matrix (note: rmass has been inverted already) + do iglob = 1,NGLOB_AB + accel(1,iglob) = accel(1,iglob)*rmassx(iglob) + accel(2,iglob) = accel(2,iglob)*rmassy(iglob) + accel(3,iglob) = accel(3,iglob)*rmassz(iglob) + enddo -! updates acceleration with ocean load term + ! updates acceleration with ocean load term if (APPROXIMATE_OCEAN_LOAD) then call compute_coupling_ocean(NSPEC_AB,NGLOB_AB, & ibool,rmassx,rmassy,rmassz, & @@ -174,7 +176,7 @@ subroutine compute_forces_viscoelastic_calling() num_free_surface_faces) endif -! impose Dirichlet conditions on the outer edges of the C-PML layers + ! impose Dirichlet conditions on the outer edges of the C-PML layers if (PML_CONDITIONS) then do iface = 1,num_abs_boundary_faces ispec = abs_boundary_ispec(iface) @@ -188,7 +190,7 @@ subroutine compute_forces_viscoelastic_calling() j = abs_boundary_ijk(2,igll,iface) k = abs_boundary_ijk(3,igll,iface) - iglob=ibool(i,j,k,ispec) + iglob = ibool(i,j,k,ispec) accel(:,iglob) = 0._CUSTOM_REAL veloc(:,iglob) = 0._CUSTOM_REAL @@ -331,11 +333,11 @@ subroutine compute_forces_viscoelastic_backward_calling() ! waits for send/receive requests to be completed and assembles values ! adjoint simulations call assemble_MPI_vector_async_w_ord(NPROC,NGLOB_ADJOINT,b_accel, & - b_buffer_recv_vector_ext_mesh,num_interfaces_ext_mesh, & - max_nibool_interfaces_ext_mesh, & - nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, & - b_request_send_vector_ext_mesh,b_request_recv_vector_ext_mesh, & - my_neighbors_ext_mesh,myrank) + b_buffer_recv_vector_ext_mesh,num_interfaces_ext_mesh, & + max_nibool_interfaces_ext_mesh, & + nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh, & + b_request_send_vector_ext_mesh,b_request_recv_vector_ext_mesh, & + my_neighbors_ext_mesh) endif enddo diff --git a/src/specfem3D/compute_kernels.f90 b/src/specfem3D/compute_kernels.f90 index a81040308..c4a68bfda 100644 --- a/src/specfem3D/compute_kernels.f90 +++ b/src/specfem3D/compute_kernels.f90 @@ -361,6 +361,7 @@ subroutine compute_kernels_po() endif end subroutine compute_kernels_po + ! !------------------------------------------------------------------------------------------------- ! @@ -472,14 +473,15 @@ subroutine compute_kernels_Hessian() end subroutine compute_kernels_Hessian +! !------------------------------------------------------------------------------------------------- ! + ! subroutine to compute the kernels for the 21 elastic coefficients ! Last modified 19/04/2007 -!------------------------------------------------------------------- subroutine compute_strain_product(prod,eps_trace_over_3,epsdev, & - b_eps_trace_over_3,b_epsdev) + b_eps_trace_over_3,b_epsdev) ! Purpose : compute the 21 strain products at a grid point ! (ispec,i,j,k fixed) and at a time t to compute then the kernels cij_kl (Voigt notation) @@ -532,85 +534,91 @@ subroutine compute_strain_product(prod,eps_trace_over_3,epsdev, & end subroutine compute_strain_product +! +!------------------------------------------------------------------------------------------------- +! + subroutine compute_anisotropic_kernels_for_velocity_data(vsem_fwd,asem_fwd,dsem_adj,vsem_adj) - use specfem_par, only: ngllx, nglly, ngllz, ibool, hprime_xx, hprime_yy, hprime_zz, nspec_ab, & - xix, xiy, xiz, etax, etay, etaz, gammax, gammay, gammaz, CUSTOM_REAL, & - nglob_ab, ndim, deltat, irregular_element_number, xix_regular - use specfem_par_elastic, only: ispec_is_elastic, rho_kl, cijkl_kl + use specfem_par, only: ngllx, nglly, ngllz, ibool, hprime_xx, hprime_yy, hprime_zz, nspec_ab, & + xix, xiy, xiz, etax, etay, etaz, gammax, gammay, gammaz, CUSTOM_REAL, & + nglob_ab, ndim, deltat, irregular_element_number, xix_regular + use specfem_par_elastic, only: ispec_is_elastic, rho_kl, cijkl_kl + implicit none - real(kind=CUSTOM_REAL), dimension(ndim,nglob_ab), intent(in) :: vsem_fwd, asem_fwd - real(kind=CUSTOM_REAL), dimension(ndim,nglob_ab), intent(in) :: dsem_adj, vsem_adj + real(kind=CUSTOM_REAL), dimension(ndim,nglob_ab), intent(in) :: vsem_fwd, asem_fwd + real(kind=CUSTOM_REAL), dimension(ndim,nglob_ab), intent(in) :: dsem_adj, vsem_adj - integer :: i, j, k, ielem, iglob, ielem_irreg ! ,l + ! local parameters + integer :: i, j, k, ielem, iglob, ielem_irreg ! ,l - real(kind=CUSTOM_REAL) :: dxil_dxl, dxil_dyl, dxil_dzl - real(kind=CUSTOM_REAL) :: detal_dxl, detal_dyl, detal_dzl - real(kind=CUSTOM_REAL) :: dgaml_dxl, dgaml_dyl, dgaml_dzl + real(kind=CUSTOM_REAL) :: dxil_dxl, dxil_dyl, dxil_dzl + real(kind=CUSTOM_REAL) :: detal_dxl, detal_dyl, detal_dzl + real(kind=CUSTOM_REAL) :: dgaml_dxl, dgaml_dyl, dgaml_dzl - real(kind=CUSTOM_REAL) :: dux_dxil, dux_detal, dux_dgaml - real(kind=CUSTOM_REAL) :: duy_dxil, duy_detal, duy_dgaml - real(kind=CUSTOM_REAL) :: duz_dxil, duz_detal, duz_dgaml + real(kind=CUSTOM_REAL) :: dux_dxil, dux_detal, dux_dgaml + real(kind=CUSTOM_REAL) :: duy_dxil, duy_detal, duy_dgaml + real(kind=CUSTOM_REAL) :: duz_dxil, duz_detal, duz_dgaml - real(kind=CUSTOM_REAL) :: dvx_dxil, dvx_detal, dvx_dgaml - real(kind=CUSTOM_REAL) :: dvy_dxil, dvy_detal, dvy_dgaml - real(kind=CUSTOM_REAL) :: dvz_dxil, dvz_detal, dvz_dgaml + real(kind=CUSTOM_REAL) :: dvx_dxil, dvx_detal, dvx_dgaml + real(kind=CUSTOM_REAL) :: dvy_dxil, dvy_detal, dvy_dgaml + real(kind=CUSTOM_REAL) :: dvz_dxil, dvz_detal, dvz_dgaml - real(kind=CUSTOM_REAL) :: dux_dxl, dux_dyl, dux_dzl - real(kind=CUSTOM_REAL) :: duy_dxl, duy_dyl, duy_dzl - real(kind=CUSTOM_REAL) :: duz_dxl, duz_dyl, duz_dzl - real(kind=CUSTOM_REAL) :: dux_dyl_plus_duy_dxl - real(kind=CUSTOM_REAL) :: duz_dxl_plus_dux_dzl, duz_dyl_plus_duy_dzl + real(kind=CUSTOM_REAL) :: dux_dxl, dux_dyl, dux_dzl + real(kind=CUSTOM_REAL) :: duy_dxl, duy_dyl, duy_dzl + real(kind=CUSTOM_REAL) :: duz_dxl, duz_dyl, duz_dzl + real(kind=CUSTOM_REAL) :: dux_dyl_plus_duy_dxl + real(kind=CUSTOM_REAL) :: duz_dxl_plus_dux_dzl, duz_dyl_plus_duy_dzl - real(kind=CUSTOM_REAL) :: dvx_dxl, dvx_dyl, dvx_dzl - real(kind=CUSTOM_REAL) :: dvy_dxl, dvy_dyl, dvy_dzl - real(kind=CUSTOM_REAL) :: dvz_dxl, dvz_dyl, dvz_dzl - real(kind=CUSTOM_REAL) :: dvx_dyl_plus_dvy_dxl - real(kind=CUSTOM_REAL) :: dvz_dxl_plus_dvx_dzl, dvz_dyl_plus_dvy_dzl + real(kind=CUSTOM_REAL) :: dvx_dxl, dvx_dyl, dvx_dzl + real(kind=CUSTOM_REAL) :: dvy_dxl, dvy_dyl, dvy_dzl + real(kind=CUSTOM_REAL) :: dvz_dxl, dvz_dyl, dvz_dzl + real(kind=CUSTOM_REAL) :: dvx_dyl_plus_dvy_dxl + real(kind=CUSTOM_REAL) :: dvz_dxl_plus_dvx_dzl, dvz_dyl_plus_dvy_dzl - real(kind=CUSTOM_REAL) :: fac + real(kind=CUSTOM_REAL) :: fac - real(kind=CUSTOM_REAL), dimension(3,ngllx,nglly,ngllz) :: vadj_gll, afwd_gll - real(kind=CUSTOM_REAL), dimension(3,ngllx,nglly,ngllz) :: vsem_fwd_gll, dsem_adj_gll + real(kind=CUSTOM_REAL), dimension(3,ngllx,nglly,ngllz) :: vadj_gll, afwd_gll + real(kind=CUSTOM_REAL), dimension(3,ngllx,nglly,ngllz) :: vsem_fwd_gll, dsem_adj_gll - !*** Loop over GLL points - do ielem=1,nspec_ab + !*** Loop over GLL points + do ielem = 1,nspec_ab - if (ispec_is_elastic(ielem)) then + if (ispec_is_elastic(ielem)) then - !*** Real first loop for optim - do k=1,ngllz - do j=1,nglly - do i=1,ngllx + !*** Real first loop for optim + do k = 1,ngllz + do j = 1,nglly + do i = 1,ngllx - iglob = ibool(i,j,k,ielem) ! find global index + iglob = ibool(i,j,k,ielem) ! find global index - vadj_gll(1,i,j,k) = vsem_adj(1,iglob) - vadj_gll(2,i,j,k) = vsem_adj(2,iglob) - vadj_gll(3,i,j,k) = vsem_adj(3,iglob) + vadj_gll(1,i,j,k) = vsem_adj(1,iglob) + vadj_gll(2,i,j,k) = vsem_adj(2,iglob) + vadj_gll(3,i,j,k) = vsem_adj(3,iglob) - afwd_gll(1,i,j,k) = asem_fwd(1,iglob) - afwd_gll(2,i,j,k) = asem_fwd(2,iglob) - afwd_gll(3,i,j,k) = asem_fwd(3,iglob) + afwd_gll(1,i,j,k) = asem_fwd(1,iglob) + afwd_gll(2,i,j,k) = asem_fwd(2,iglob) + afwd_gll(3,i,j,k) = asem_fwd(3,iglob) - dsem_adj_gll(1,i,j,k) = dsem_adj(1,iglob) - dsem_adj_gll(2,i,j,k) = dsem_adj(2,iglob) - dsem_adj_gll(3,i,j,k) = dsem_adj(3,iglob) + dsem_adj_gll(1,i,j,k) = dsem_adj(1,iglob) + dsem_adj_gll(2,i,j,k) = dsem_adj(2,iglob) + dsem_adj_gll(3,i,j,k) = dsem_adj(3,iglob) - vsem_fwd_gll(1,i,j,k) = vsem_fwd(1,iglob) - vsem_fwd_gll(2,i,j,k) = vsem_fwd(2,iglob) - vsem_fwd_gll(3,i,j,k) = vsem_fwd(3,iglob) + vsem_fwd_gll(1,i,j,k) = vsem_fwd(1,iglob) + vsem_fwd_gll(2,i,j,k) = vsem_fwd(2,iglob) + vsem_fwd_gll(3,i,j,k) = vsem_fwd(3,iglob) - enddo - enddo - enddo + enddo + enddo + enddo - ielem_irreg = irregular_element_number(ielem) + ielem_irreg = irregular_element_number(ielem) - do k=1,ngllz - do j=1,nglly - do i=1,ngllx + do k = 1,ngllz + do j = 1,nglly + do i = 1,ngllx !! DK DK Oct 2018: we could (and should) use the Deville matrix products instead here !! DK DK Oct 2018: we could (and should) use the Deville matrix products instead here @@ -623,325 +631,325 @@ subroutine compute_anisotropic_kernels_for_velocity_data(vsem_fwd,asem_fwd,dsem_ !! DK DK Oct 2018: we could (and should) use the Deville matrix products instead here !! DK DK Oct 2018: we could (and should) use the Deville matrix products instead here - !================================================================ - ! Compute strain related terms (adjoit strain and time der of normal strain) - !*** Init derivatives to 0 - !* Normal - dvx_dxil = 0._CUSTOM_REAL - dvx_detal = 0._CUSTOM_REAL - dvx_dgaml = 0._CUSTOM_REAL - dvy_dxil = 0._CUSTOM_REAL - dvy_detal = 0._CUSTOM_REAL - dvy_dgaml = 0._CUSTOM_REAL - dvz_dxil = 0._CUSTOM_REAL - dvz_detal = 0._CUSTOM_REAL - dvz_dgaml = 0._CUSTOM_REAL - - !* Adjoint - dux_dxil = 0._CUSTOM_REAL - dux_detal = 0._CUSTOM_REAL - dux_dgaml = 0._CUSTOM_REAL - duy_dxil = 0._CUSTOM_REAL - duy_detal = 0._CUSTOM_REAL - duy_dgaml = 0._CUSTOM_REAL - duz_dxil = 0._CUSTOM_REAL - duz_detal = 0._CUSTOM_REAL - duz_dgaml = 0._CUSTOM_REAL - - !*** Field derivatives wrt xi - fac = hprime_xx(1,i) ! derivative of local lagrange polynomials - dvx_dxil = dvx_dxil + vsem_fwd_gll(1,1,j,k) * fac - dvy_dxil = dvy_dxil + vsem_fwd_gll(2,1,j,k) * fac - dvz_dxil = dvz_dxil + vsem_fwd_gll(3,1,j,k) * fac - dux_dxil = dux_dxil + dsem_adj_gll(1,1,j,k) * fac - duy_dxil = duy_dxil + dsem_adj_gll(2,1,j,k) * fac - duz_dxil = duz_dxil + dsem_adj_gll(3,1,j,k) * fac - - fac = hprime_xx(2,i) ! derivative of local lagrange polynomials - dvx_dxil = dvx_dxil + vsem_fwd_gll(1,2,j,k) * fac - dvy_dxil = dvy_dxil + vsem_fwd_gll(2,2,j,k) * fac - dvz_dxil = dvz_dxil + vsem_fwd_gll(3,2,j,k) * fac - dux_dxil = dux_dxil + dsem_adj_gll(1,2,j,k) * fac - duy_dxil = duy_dxil + dsem_adj_gll(2,2,j,k) * fac - duz_dxil = duz_dxil + dsem_adj_gll(3,2,j,k) * fac - - fac = hprime_xx(3,i) ! derivative of local lagrange polynomials - dvx_dxil = dvx_dxil + vsem_fwd_gll(1,3,j,k) * fac - dvy_dxil = dvy_dxil + vsem_fwd_gll(2,3,j,k) * fac - dvz_dxil = dvz_dxil + vsem_fwd_gll(3,3,j,k) * fac - dux_dxil = dux_dxil + dsem_adj_gll(1,3,j,k) * fac - duy_dxil = duy_dxil + dsem_adj_gll(2,3,j,k) * fac - duz_dxil = duz_dxil + dsem_adj_gll(3,3,j,k) * fac - - fac = hprime_xx(4,i) ! derivative of local lagrange polynomials - dvx_dxil = dvx_dxil + vsem_fwd_gll(1,4,j,k) * fac - dvy_dxil = dvy_dxil + vsem_fwd_gll(2,4,j,k) * fac - dvz_dxil = dvz_dxil + vsem_fwd_gll(3,4,j,k) * fac - dux_dxil = dux_dxil + dsem_adj_gll(1,4,j,k) * fac - duy_dxil = duy_dxil + dsem_adj_gll(2,4,j,k) * fac - duz_dxil = duz_dxil + dsem_adj_gll(3,4,j,k) * fac - - fac = hprime_xx(5,i) ! derivative of local lagrange polynomials - dvx_dxil = dvx_dxil + vsem_fwd_gll(1,5,j,k) * fac - dvy_dxil = dvy_dxil + vsem_fwd_gll(2,5,j,k) * fac - dvz_dxil = dvz_dxil + vsem_fwd_gll(3,5,j,k) * fac - dux_dxil = dux_dxil + dsem_adj_gll(1,5,j,k) * fac - duy_dxil = duy_dxil + dsem_adj_gll(2,5,j,k) * fac - duz_dxil = duz_dxil + dsem_adj_gll(3,5,j,k) * fac - - !*** Field derivatives wrt eta - fac = hprime_yy(1,j) ! derivative of local lageange polynomials - dvx_detal = dvx_detal + vsem_fwd_gll(1,i,1,k) * fac - dvy_detal = dvy_detal + vsem_fwd_gll(2,i,1,k) * fac - dvz_detal = dvz_detal + vsem_fwd_gll(3,i,1,k) * fac - dux_detal = dux_detal + dsem_adj_gll(1,i,1,k) * fac - duy_detal = duy_detal + dsem_adj_gll(2,i,1,k) * fac - duz_detal = duz_detal + dsem_adj_gll(3,i,1,k) * fac - - fac = hprime_yy(2,j) ! derivative of local lageange polynomials - dvx_detal = dvx_detal + vsem_fwd_gll(1,i,2,k) * fac - dvy_detal = dvy_detal + vsem_fwd_gll(2,i,2,k) * fac - dvz_detal = dvz_detal + vsem_fwd_gll(3,i,2,k) * fac - dux_detal = dux_detal + dsem_adj_gll(1,i,2,k) * fac - duy_detal = duy_detal + dsem_adj_gll(2,i,2,k) * fac - duz_detal = duz_detal + dsem_adj_gll(3,i,2,k) * fac - - fac = hprime_yy(3,j) ! derivative of local lageange polynomials - dvx_detal = dvx_detal + vsem_fwd_gll(1,i,3,k) * fac - dvy_detal = dvy_detal + vsem_fwd_gll(2,i,3,k) * fac - dvz_detal = dvz_detal + vsem_fwd_gll(3,i,3,k) * fac - dux_detal = dux_detal + dsem_adj_gll(1,i,3,k) * fac - duy_detal = duy_detal + dsem_adj_gll(2,i,3,k) * fac - duz_detal = duz_detal + dsem_adj_gll(3,i,3,k) * fac - - fac = hprime_yy(4,j) ! derivative of local lageange polynomials - dvx_detal = dvx_detal + vsem_fwd_gll(1,i,4,k) * fac - dvy_detal = dvy_detal + vsem_fwd_gll(2,i,4,k) * fac - dvz_detal = dvz_detal + vsem_fwd_gll(3,i,4,k) * fac - dux_detal = dux_detal + dsem_adj_gll(1,i,4,k) * fac - duy_detal = duy_detal + dsem_adj_gll(2,i,4,k) * fac - duz_detal = duz_detal + dsem_adj_gll(3,i,4,k) * fac - - fac = hprime_yy(5,j) ! derivative of local lageange polynomials - dvx_detal = dvx_detal + vsem_fwd_gll(1,i,5,k) * fac - dvy_detal = dvy_detal + vsem_fwd_gll(2,i,5,k) * fac - dvz_detal = dvz_detal + vsem_fwd_gll(3,i,5,k) * fac - dux_detal = dux_detal + dsem_adj_gll(1,i,5,k) * fac - duy_detal = duy_detal + dsem_adj_gll(2,i,5,k) * fac - duz_detal = duz_detal + dsem_adj_gll(3,i,5,k) * fac - - !*** Field derivatives wrt gamma - fac = hprime_zz(1,k) ! derivative of local lagange polynomials - dvx_dgaml = dvx_dgaml + vsem_fwd_gll(1,i,j,1) * fac - dvy_dgaml = dvy_dgaml + vsem_fwd_gll(2,i,j,1) * fac - dvz_dgaml = dvz_dgaml + vsem_fwd_gll(3,i,j,1) * fac - dux_dgaml = dux_dgaml + dsem_adj_gll(1,i,j,1) * fac - duy_dgaml = duy_dgaml + dsem_adj_gll(2,i,j,1) * fac - duz_dgaml = duz_dgaml + dsem_adj_gll(3,i,j,1) * fac - - fac = hprime_zz(2,k) ! derivative of local lagange polynomials - dvx_dgaml = dvx_dgaml + vsem_fwd_gll(1,i,j,2) * fac - dvy_dgaml = dvy_dgaml + vsem_fwd_gll(2,i,j,2) * fac - dvz_dgaml = dvz_dgaml + vsem_fwd_gll(3,i,j,2) * fac - dux_dgaml = dux_dgaml + dsem_adj_gll(1,i,j,2) * fac - duy_dgaml = duy_dgaml + dsem_adj_gll(2,i,j,2) * fac - duz_dgaml = duz_dgaml + dsem_adj_gll(3,i,j,2) * fac - - fac = hprime_zz(3,k) ! derivative of local lagange polynomials - dvx_dgaml = dvx_dgaml + vsem_fwd_gll(1,i,j,3) * fac - dvy_dgaml = dvy_dgaml + vsem_fwd_gll(2,i,j,3) * fac - dvz_dgaml = dvz_dgaml + vsem_fwd_gll(3,i,j,3) * fac - dux_dgaml = dux_dgaml + dsem_adj_gll(1,i,j,3) * fac - duy_dgaml = duy_dgaml + dsem_adj_gll(2,i,j,3) * fac - duz_dgaml = duz_dgaml + dsem_adj_gll(3,i,j,3) * fac - - fac = hprime_zz(4,k) ! derivative of local lagange polynomials - dvx_dgaml = dvx_dgaml + vsem_fwd_gll(1,i,j,4) * fac - dvy_dgaml = dvy_dgaml + vsem_fwd_gll(2,i,j,4) * fac - dvz_dgaml = dvz_dgaml + vsem_fwd_gll(3,i,j,4) * fac - dux_dgaml = dux_dgaml + dsem_adj_gll(1,i,j,4) * fac - duy_dgaml = duy_dgaml + dsem_adj_gll(2,i,j,4) * fac - duz_dgaml = duz_dgaml + dsem_adj_gll(3,i,j,4) * fac - - fac = hprime_zz(5,k) ! derivative of local lagange polynomials - dvx_dgaml = dvx_dgaml + vsem_fwd_gll(1,i,j,5) * fac - dvy_dgaml = dvy_dgaml + vsem_fwd_gll(2,i,j,5) * fac - dvz_dgaml = dvz_dgaml + vsem_fwd_gll(3,i,j,5) * fac - dux_dgaml = dux_dgaml + dsem_adj_gll(1,i,j,5) * fac - duy_dgaml = duy_dgaml + dsem_adj_gll(2,i,j,5) * fac - duz_dgaml = duz_dgaml + dsem_adj_gll(3,i,j,5) * fac - - if (ielem_irreg /= 0) then ! irregular element - - !*** Get local derivatives of ref square coord wrt Cartesian ones (jacobian) - dxil_dxl = xix(i,j,k,ielem_irreg) - dxil_dyl = xiy(i,j,k,ielem_irreg) - dxil_dzl = xiz(i,j,k,ielem_irreg) - detal_dxl = etax(i,j,k,ielem_irreg) - detal_dyl = etay(i,j,k,ielem_irreg) - detal_dzl = etaz(i,j,k,ielem_irreg) - dgaml_dxl = gammax(i,j,k,ielem_irreg) - dgaml_dyl = gammay(i,j,k,ielem_irreg) - dgaml_dzl = gammaz(i,j,k,ielem_irreg) - - !*** Strain - !* Normal state normal strain - dvx_dxl = dvx_dxil * dxil_dxl + dvx_detal * detal_dxl + dvx_dgaml * dgaml_dxl - dvx_dyl = dvx_dxil * dxil_dyl + dvx_detal * detal_dyl + dvx_dgaml * dgaml_dyl - dvx_dzl = dvx_dxil * dxil_dzl + dvx_detal * detal_dzl + dvx_dgaml * dgaml_dzl - dvy_dxl = dvy_dxil * dxil_dxl + dvy_detal * detal_dxl + dvy_dgaml * dgaml_dxl - dvy_dyl = dvy_dxil * dxil_dyl + dvy_detal * detal_dyl + dvy_dgaml * dgaml_dyl - dvy_dzl = dvy_dxil * dxil_dzl + dvy_detal * detal_dzl + dvy_dgaml * dgaml_dzl - dvz_dxl = dvz_dxil * dxil_dxl + dvz_detal * detal_dxl + dvz_dgaml * dgaml_dxl - dvz_dyl = dvz_dxil * dxil_dyl + dvz_detal * detal_dyl + dvz_dgaml * dgaml_dyl - dvz_dzl = dvz_dxil * dxil_dzl + dvz_detal * detal_dzl + dvz_dgaml * dgaml_dzl - - !* Adjoint state normal strain - dux_dxl = dux_dxil * dxil_dxl + dux_detal * detal_dxl + dux_dgaml * dgaml_dxl - dux_dyl = dux_dxil * dxil_dyl + dux_detal * detal_dyl + dux_dgaml * dgaml_dyl - dux_dzl = dux_dxil * dxil_dzl + dux_detal * detal_dzl + dux_dgaml * dgaml_dzl - duy_dxl = duy_dxil * dxil_dxl + duy_detal * detal_dxl + duy_dgaml * dgaml_dxl - duy_dyl = duy_dxil * dxil_dyl + duy_detal * detal_dyl + duy_dgaml * dgaml_dyl - duy_dzl = duy_dxil * dxil_dzl + duy_detal * detal_dzl + duy_dgaml * dgaml_dzl - duz_dxl = duz_dxil * dxil_dxl + duz_detal * detal_dxl + duz_dgaml * dgaml_dxl - duz_dyl = duz_dxil * dxil_dyl + duz_detal * detal_dyl + duz_dgaml * dgaml_dyl - duz_dzl = duz_dxil * dxil_dzl + duz_detal * detal_dzl + duz_dgaml * dgaml_dzl - - else ! regular element - - !*** Strain - !* Normal state normal strain - dvx_dxl = dvx_dxil * xix_regular - dvx_dyl = dvx_detal * xix_regular - dvx_dzl = dvx_dgaml * xix_regular - dvy_dxl = dvy_dxil * xix_regular - dvy_dyl = dvy_detal * xix_regular - dvy_dzl = dvy_dgaml * xix_regular - dvz_dxl = dvz_dxil * xix_regular - dvz_dyl = dvz_detal * xix_regular - dvz_dzl = dvz_dgaml * xix_regular - - !* Adjoint state normal strain - dux_dxl = dux_dxil * xix_regular - dux_dyl = dux_detal * xix_regular - dux_dzl = dux_dgaml * xix_regular - duy_dxl = duy_dxil * xix_regular - duy_dyl = duy_detal * xix_regular - duy_dzl = duy_dgaml * xix_regular - duz_dxl = duz_dxil * xix_regular - duz_dyl = duz_detal * xix_regular - duz_dzl = duz_dgaml * xix_regular - - endif ! element regularity - - !* Normal state non normal strain - dvx_dyl_plus_dvy_dxl = dvx_dyl + dvy_dxl - dvz_dxl_plus_dvx_dzl = dvz_dxl + dvx_dzl - dvz_dyl_plus_dvy_dzl = dvz_dyl + dvy_dzl - - !* Adjoint state non normal strain - dux_dyl_plus_duy_dxl = dux_dyl + duy_dxl - duz_dxl_plus_dux_dzl = duz_dxl + dux_dzl - duz_dyl_plus_duy_dzl = duz_dyl + duy_dzl - - - !=========================================================================== - ! Gradient - rho_kl(i,j,k,ielem) = rho_kl(i,j,k,ielem) & - + (afwd_gll(1,i,j,k) * vadj_gll(1,i,j,k) & - + afwd_gll(2,i,j,k) * vadj_gll(2,i,j,k) & - + afwd_gll(3,i,j,k) * vadj_gll(3,i,j,k))*deltat - - !*** Gradient wrt c_ij - !* c11 - cijkl_kl(1,i,j,k,ielem) = cijkl_kl(1,i,j,k,ielem) + dvx_dxl*dux_dxl*deltat - - !* c12 - cijkl_kl(2,i,j,k,ielem) = cijkl_kl(2,i,j,k,ielem) & - +(dvx_dxl * duy_dyl + dvy_dyl * dux_dxl)*deltat - - !* c13 - cijkl_kl(3,i,j,k,ielem) = cijkl_kl(3,i,j,k,ielem) & - +(dvx_dxl * duz_dzl + dvz_dzl * dux_dxl)*deltat - - !* c14 - cijkl_kl(4,i,j,k,ielem) = cijkl_kl(4,i,j,k,ielem) & - +(dvx_dxl * duz_dyl_plus_duy_dzl + dvz_dyl_plus_dvy_dzl * dux_dxl)*deltat - - !* c15 - cijkl_kl(5,i,j,k,ielem) = cijkl_kl(5,i,j,k,ielem) & - +(dvx_dxl * duz_dxl_plus_dux_dzl + dvz_dxl_plus_dvx_dzl * dux_dxl)*deltat - - !* c16 - cijkl_kl(6,i,j,k,ielem) = cijkl_kl(6,i,j,k,ielem) & - +(dvx_dxl * dux_dyl_plus_duy_dxl + dvx_dyl_plus_dvy_dxl * dux_dxl)*deltat - - !* c22 - cijkl_kl(7,i,j,k,ielem) = cijkl_kl(7,i,j,k,ielem) + dvy_dyl*duy_dyl*deltat - - !* c23 - cijkl_kl(8,i,j,k,ielem) = cijkl_kl(8,i,j,k,ielem) & - +(dvy_dyl * duz_dzl + dvz_dzl * duy_dyl)*deltat - - !* c24 - cijkl_kl(9,i,j,k,ielem) = cijkl_kl(9,i,j,k,ielem) & - +(dvy_dyl * duz_dyl_plus_duy_dzl + dvz_dyl_plus_dvy_dzl * duy_dyl)*deltat - - !* c25 - cijkl_kl(10,i,j,k,ielem) = cijkl_kl(10,i,j,k,ielem) & - +(dvy_dyl * duz_dxl_plus_dux_dzl + dvz_dxl_plus_dvx_dzl * duy_dyl)*deltat - - !* c26 - cijkl_kl(11,i,j,k,ielem) = cijkl_kl(11,i,j,k,ielem) & - +(dvy_dyl * dux_dyl_plus_duy_dxl + dvx_dyl_plus_dvy_dxl * duy_dyl)*deltat - - !* c33 - cijkl_kl(12,i,j,k,ielem) = cijkl_kl(12,i,j,k,ielem) + dvz_dzl*duz_dzl*deltat - - !* c34 - cijkl_kl(13,i,j,k,ielem) = cijkl_kl(13,i,j,k,ielem) & - +(dvz_dzl * duz_dyl_plus_duy_dzl + dvz_dyl_plus_dvy_dzl * duz_dzl)*deltat - - !* c35 - cijkl_kl(14,i,j,k,ielem) = cijkl_kl(14,i,j,k,ielem) & - +(dvz_dzl * duz_dxl_plus_dux_dzl + dvz_dxl_plus_dvx_dzl * duz_dzl)*deltat - - !* c36 - cijkl_kl(15,i,j,k,ielem) = cijkl_kl(15,i,j,k,ielem) & - +(dvz_dzl * dux_dyl_plus_duy_dxl + dvx_dyl_plus_dvy_dxl * duz_dzl)*deltat - - !* c44 - cijkl_kl(16,i,j,k,ielem) = cijkl_kl(16,i,j,k,ielem) & - + dvz_dyl_plus_dvy_dzl * duz_dyl_plus_duy_dzl*deltat - - !* c45 - cijkl_kl(17,i,j,k,ielem) = cijkl_kl(17,i,j,k,ielem) & - +(dvz_dyl_plus_dvy_dzl * duz_dxl_plus_dux_dzl & - + dvz_dxl_plus_dvx_dzl * duz_dyl_plus_duy_dzl)*deltat - - !* c46 - cijkl_kl(18,i,j,k,ielem) = cijkl_kl(18,i,j,k,ielem) & - +(dvz_dyl_plus_dvy_dzl * dux_dyl_plus_duy_dxl & - + dvx_dyl_plus_dvy_dxl * duz_dyl_plus_duy_dzl)*deltat - - !* c55 - cijkl_kl(19,i,j,k,ielem) = cijkl_kl(19,i,j,k,ielem) & - + dvz_dxl_plus_dvx_dzl * duz_dxl_plus_dux_dzl*deltat - - !* c56 - cijkl_kl(20,i,j,k,ielem) = cijkl_kl(20,i,j,k,ielem) & - +(dvz_dxl_plus_dvx_dzl * dux_dyl_plus_duy_dxl & - + dvx_dyl_plus_dvy_dxl * duz_dxl_plus_dux_dzl)*deltat - - !* c66 - cijkl_kl(21,i,j,k,ielem) = cijkl_kl(21,i,j,k,ielem) & - + dvx_dyl_plus_dvy_dxl * dux_dyl_plus_duy_dxl*deltat - enddo - enddo - enddo + !================================================================ + ! Compute strain related terms (adjoit strain and time der of normal strain) + !*** Init derivatives to 0 + !* Normal + dvx_dxil = 0._CUSTOM_REAL + dvx_detal = 0._CUSTOM_REAL + dvx_dgaml = 0._CUSTOM_REAL + dvy_dxil = 0._CUSTOM_REAL + dvy_detal = 0._CUSTOM_REAL + dvy_dgaml = 0._CUSTOM_REAL + dvz_dxil = 0._CUSTOM_REAL + dvz_detal = 0._CUSTOM_REAL + dvz_dgaml = 0._CUSTOM_REAL + + !* Adjoint + dux_dxil = 0._CUSTOM_REAL + dux_detal = 0._CUSTOM_REAL + dux_dgaml = 0._CUSTOM_REAL + duy_dxil = 0._CUSTOM_REAL + duy_detal = 0._CUSTOM_REAL + duy_dgaml = 0._CUSTOM_REAL + duz_dxil = 0._CUSTOM_REAL + duz_detal = 0._CUSTOM_REAL + duz_dgaml = 0._CUSTOM_REAL + + !*** Field derivatives wrt xi + fac = hprime_xx(1,i) ! derivative of local lagrange polynomials + dvx_dxil = dvx_dxil + vsem_fwd_gll(1,1,j,k) * fac + dvy_dxil = dvy_dxil + vsem_fwd_gll(2,1,j,k) * fac + dvz_dxil = dvz_dxil + vsem_fwd_gll(3,1,j,k) * fac + dux_dxil = dux_dxil + dsem_adj_gll(1,1,j,k) * fac + duy_dxil = duy_dxil + dsem_adj_gll(2,1,j,k) * fac + duz_dxil = duz_dxil + dsem_adj_gll(3,1,j,k) * fac + + fac = hprime_xx(2,i) ! derivative of local lagrange polynomials + dvx_dxil = dvx_dxil + vsem_fwd_gll(1,2,j,k) * fac + dvy_dxil = dvy_dxil + vsem_fwd_gll(2,2,j,k) * fac + dvz_dxil = dvz_dxil + vsem_fwd_gll(3,2,j,k) * fac + dux_dxil = dux_dxil + dsem_adj_gll(1,2,j,k) * fac + duy_dxil = duy_dxil + dsem_adj_gll(2,2,j,k) * fac + duz_dxil = duz_dxil + dsem_adj_gll(3,2,j,k) * fac + + fac = hprime_xx(3,i) ! derivative of local lagrange polynomials + dvx_dxil = dvx_dxil + vsem_fwd_gll(1,3,j,k) * fac + dvy_dxil = dvy_dxil + vsem_fwd_gll(2,3,j,k) * fac + dvz_dxil = dvz_dxil + vsem_fwd_gll(3,3,j,k) * fac + dux_dxil = dux_dxil + dsem_adj_gll(1,3,j,k) * fac + duy_dxil = duy_dxil + dsem_adj_gll(2,3,j,k) * fac + duz_dxil = duz_dxil + dsem_adj_gll(3,3,j,k) * fac + + fac = hprime_xx(4,i) ! derivative of local lagrange polynomials + dvx_dxil = dvx_dxil + vsem_fwd_gll(1,4,j,k) * fac + dvy_dxil = dvy_dxil + vsem_fwd_gll(2,4,j,k) * fac + dvz_dxil = dvz_dxil + vsem_fwd_gll(3,4,j,k) * fac + dux_dxil = dux_dxil + dsem_adj_gll(1,4,j,k) * fac + duy_dxil = duy_dxil + dsem_adj_gll(2,4,j,k) * fac + duz_dxil = duz_dxil + dsem_adj_gll(3,4,j,k) * fac + + fac = hprime_xx(5,i) ! derivative of local lagrange polynomials + dvx_dxil = dvx_dxil + vsem_fwd_gll(1,5,j,k) * fac + dvy_dxil = dvy_dxil + vsem_fwd_gll(2,5,j,k) * fac + dvz_dxil = dvz_dxil + vsem_fwd_gll(3,5,j,k) * fac + dux_dxil = dux_dxil + dsem_adj_gll(1,5,j,k) * fac + duy_dxil = duy_dxil + dsem_adj_gll(2,5,j,k) * fac + duz_dxil = duz_dxil + dsem_adj_gll(3,5,j,k) * fac + + !*** Field derivatives wrt eta + fac = hprime_yy(1,j) ! derivative of local lageange polynomials + dvx_detal = dvx_detal + vsem_fwd_gll(1,i,1,k) * fac + dvy_detal = dvy_detal + vsem_fwd_gll(2,i,1,k) * fac + dvz_detal = dvz_detal + vsem_fwd_gll(3,i,1,k) * fac + dux_detal = dux_detal + dsem_adj_gll(1,i,1,k) * fac + duy_detal = duy_detal + dsem_adj_gll(2,i,1,k) * fac + duz_detal = duz_detal + dsem_adj_gll(3,i,1,k) * fac + + fac = hprime_yy(2,j) ! derivative of local lageange polynomials + dvx_detal = dvx_detal + vsem_fwd_gll(1,i,2,k) * fac + dvy_detal = dvy_detal + vsem_fwd_gll(2,i,2,k) * fac + dvz_detal = dvz_detal + vsem_fwd_gll(3,i,2,k) * fac + dux_detal = dux_detal + dsem_adj_gll(1,i,2,k) * fac + duy_detal = duy_detal + dsem_adj_gll(2,i,2,k) * fac + duz_detal = duz_detal + dsem_adj_gll(3,i,2,k) * fac + + fac = hprime_yy(3,j) ! derivative of local lageange polynomials + dvx_detal = dvx_detal + vsem_fwd_gll(1,i,3,k) * fac + dvy_detal = dvy_detal + vsem_fwd_gll(2,i,3,k) * fac + dvz_detal = dvz_detal + vsem_fwd_gll(3,i,3,k) * fac + dux_detal = dux_detal + dsem_adj_gll(1,i,3,k) * fac + duy_detal = duy_detal + dsem_adj_gll(2,i,3,k) * fac + duz_detal = duz_detal + dsem_adj_gll(3,i,3,k) * fac + + fac = hprime_yy(4,j) ! derivative of local lageange polynomials + dvx_detal = dvx_detal + vsem_fwd_gll(1,i,4,k) * fac + dvy_detal = dvy_detal + vsem_fwd_gll(2,i,4,k) * fac + dvz_detal = dvz_detal + vsem_fwd_gll(3,i,4,k) * fac + dux_detal = dux_detal + dsem_adj_gll(1,i,4,k) * fac + duy_detal = duy_detal + dsem_adj_gll(2,i,4,k) * fac + duz_detal = duz_detal + dsem_adj_gll(3,i,4,k) * fac + + fac = hprime_yy(5,j) ! derivative of local lageange polynomials + dvx_detal = dvx_detal + vsem_fwd_gll(1,i,5,k) * fac + dvy_detal = dvy_detal + vsem_fwd_gll(2,i,5,k) * fac + dvz_detal = dvz_detal + vsem_fwd_gll(3,i,5,k) * fac + dux_detal = dux_detal + dsem_adj_gll(1,i,5,k) * fac + duy_detal = duy_detal + dsem_adj_gll(2,i,5,k) * fac + duz_detal = duz_detal + dsem_adj_gll(3,i,5,k) * fac + + !*** Field derivatives wrt gamma + fac = hprime_zz(1,k) ! derivative of local lagange polynomials + dvx_dgaml = dvx_dgaml + vsem_fwd_gll(1,i,j,1) * fac + dvy_dgaml = dvy_dgaml + vsem_fwd_gll(2,i,j,1) * fac + dvz_dgaml = dvz_dgaml + vsem_fwd_gll(3,i,j,1) * fac + dux_dgaml = dux_dgaml + dsem_adj_gll(1,i,j,1) * fac + duy_dgaml = duy_dgaml + dsem_adj_gll(2,i,j,1) * fac + duz_dgaml = duz_dgaml + dsem_adj_gll(3,i,j,1) * fac + + fac = hprime_zz(2,k) ! derivative of local lagange polynomials + dvx_dgaml = dvx_dgaml + vsem_fwd_gll(1,i,j,2) * fac + dvy_dgaml = dvy_dgaml + vsem_fwd_gll(2,i,j,2) * fac + dvz_dgaml = dvz_dgaml + vsem_fwd_gll(3,i,j,2) * fac + dux_dgaml = dux_dgaml + dsem_adj_gll(1,i,j,2) * fac + duy_dgaml = duy_dgaml + dsem_adj_gll(2,i,j,2) * fac + duz_dgaml = duz_dgaml + dsem_adj_gll(3,i,j,2) * fac + + fac = hprime_zz(3,k) ! derivative of local lagange polynomials + dvx_dgaml = dvx_dgaml + vsem_fwd_gll(1,i,j,3) * fac + dvy_dgaml = dvy_dgaml + vsem_fwd_gll(2,i,j,3) * fac + dvz_dgaml = dvz_dgaml + vsem_fwd_gll(3,i,j,3) * fac + dux_dgaml = dux_dgaml + dsem_adj_gll(1,i,j,3) * fac + duy_dgaml = duy_dgaml + dsem_adj_gll(2,i,j,3) * fac + duz_dgaml = duz_dgaml + dsem_adj_gll(3,i,j,3) * fac + + fac = hprime_zz(4,k) ! derivative of local lagange polynomials + dvx_dgaml = dvx_dgaml + vsem_fwd_gll(1,i,j,4) * fac + dvy_dgaml = dvy_dgaml + vsem_fwd_gll(2,i,j,4) * fac + dvz_dgaml = dvz_dgaml + vsem_fwd_gll(3,i,j,4) * fac + dux_dgaml = dux_dgaml + dsem_adj_gll(1,i,j,4) * fac + duy_dgaml = duy_dgaml + dsem_adj_gll(2,i,j,4) * fac + duz_dgaml = duz_dgaml + dsem_adj_gll(3,i,j,4) * fac + + fac = hprime_zz(5,k) ! derivative of local lagange polynomials + dvx_dgaml = dvx_dgaml + vsem_fwd_gll(1,i,j,5) * fac + dvy_dgaml = dvy_dgaml + vsem_fwd_gll(2,i,j,5) * fac + dvz_dgaml = dvz_dgaml + vsem_fwd_gll(3,i,j,5) * fac + dux_dgaml = dux_dgaml + dsem_adj_gll(1,i,j,5) * fac + duy_dgaml = duy_dgaml + dsem_adj_gll(2,i,j,5) * fac + duz_dgaml = duz_dgaml + dsem_adj_gll(3,i,j,5) * fac + + if (ielem_irreg /= 0) then ! irregular element + + !*** Get local derivatives of ref square coord wrt Cartesian ones (jacobian) + dxil_dxl = xix(i,j,k,ielem_irreg) + dxil_dyl = xiy(i,j,k,ielem_irreg) + dxil_dzl = xiz(i,j,k,ielem_irreg) + detal_dxl = etax(i,j,k,ielem_irreg) + detal_dyl = etay(i,j,k,ielem_irreg) + detal_dzl = etaz(i,j,k,ielem_irreg) + dgaml_dxl = gammax(i,j,k,ielem_irreg) + dgaml_dyl = gammay(i,j,k,ielem_irreg) + dgaml_dzl = gammaz(i,j,k,ielem_irreg) + + !*** Strain + !* Normal state normal strain + dvx_dxl = dvx_dxil * dxil_dxl + dvx_detal * detal_dxl + dvx_dgaml * dgaml_dxl + dvx_dyl = dvx_dxil * dxil_dyl + dvx_detal * detal_dyl + dvx_dgaml * dgaml_dyl + dvx_dzl = dvx_dxil * dxil_dzl + dvx_detal * detal_dzl + dvx_dgaml * dgaml_dzl + dvy_dxl = dvy_dxil * dxil_dxl + dvy_detal * detal_dxl + dvy_dgaml * dgaml_dxl + dvy_dyl = dvy_dxil * dxil_dyl + dvy_detal * detal_dyl + dvy_dgaml * dgaml_dyl + dvy_dzl = dvy_dxil * dxil_dzl + dvy_detal * detal_dzl + dvy_dgaml * dgaml_dzl + dvz_dxl = dvz_dxil * dxil_dxl + dvz_detal * detal_dxl + dvz_dgaml * dgaml_dxl + dvz_dyl = dvz_dxil * dxil_dyl + dvz_detal * detal_dyl + dvz_dgaml * dgaml_dyl + dvz_dzl = dvz_dxil * dxil_dzl + dvz_detal * detal_dzl + dvz_dgaml * dgaml_dzl + + !* Adjoint state normal strain + dux_dxl = dux_dxil * dxil_dxl + dux_detal * detal_dxl + dux_dgaml * dgaml_dxl + dux_dyl = dux_dxil * dxil_dyl + dux_detal * detal_dyl + dux_dgaml * dgaml_dyl + dux_dzl = dux_dxil * dxil_dzl + dux_detal * detal_dzl + dux_dgaml * dgaml_dzl + duy_dxl = duy_dxil * dxil_dxl + duy_detal * detal_dxl + duy_dgaml * dgaml_dxl + duy_dyl = duy_dxil * dxil_dyl + duy_detal * detal_dyl + duy_dgaml * dgaml_dyl + duy_dzl = duy_dxil * dxil_dzl + duy_detal * detal_dzl + duy_dgaml * dgaml_dzl + duz_dxl = duz_dxil * dxil_dxl + duz_detal * detal_dxl + duz_dgaml * dgaml_dxl + duz_dyl = duz_dxil * dxil_dyl + duz_detal * detal_dyl + duz_dgaml * dgaml_dyl + duz_dzl = duz_dxil * dxil_dzl + duz_detal * detal_dzl + duz_dgaml * dgaml_dzl + + else ! regular element + + !*** Strain + !* Normal state normal strain + dvx_dxl = dvx_dxil * xix_regular + dvx_dyl = dvx_detal * xix_regular + dvx_dzl = dvx_dgaml * xix_regular + dvy_dxl = dvy_dxil * xix_regular + dvy_dyl = dvy_detal * xix_regular + dvy_dzl = dvy_dgaml * xix_regular + dvz_dxl = dvz_dxil * xix_regular + dvz_dyl = dvz_detal * xix_regular + dvz_dzl = dvz_dgaml * xix_regular + + !* Adjoint state normal strain + dux_dxl = dux_dxil * xix_regular + dux_dyl = dux_detal * xix_regular + dux_dzl = dux_dgaml * xix_regular + duy_dxl = duy_dxil * xix_regular + duy_dyl = duy_detal * xix_regular + duy_dzl = duy_dgaml * xix_regular + duz_dxl = duz_dxil * xix_regular + duz_dyl = duz_detal * xix_regular + duz_dzl = duz_dgaml * xix_regular + + endif ! element regularity + + !* Normal state non normal strain + dvx_dyl_plus_dvy_dxl = dvx_dyl + dvy_dxl + dvz_dxl_plus_dvx_dzl = dvz_dxl + dvx_dzl + dvz_dyl_plus_dvy_dzl = dvz_dyl + dvy_dzl + + !* Adjoint state non normal strain + dux_dyl_plus_duy_dxl = dux_dyl + duy_dxl + duz_dxl_plus_dux_dzl = duz_dxl + dux_dzl + duz_dyl_plus_duy_dzl = duz_dyl + duy_dzl + + + !=========================================================================== + ! Gradient + rho_kl(i,j,k,ielem) = rho_kl(i,j,k,ielem) & + + (afwd_gll(1,i,j,k) * vadj_gll(1,i,j,k) & + + afwd_gll(2,i,j,k) * vadj_gll(2,i,j,k) & + + afwd_gll(3,i,j,k) * vadj_gll(3,i,j,k))*deltat + + !*** Gradient wrt c_ij + !* c11 + cijkl_kl(1,i,j,k,ielem) = cijkl_kl(1,i,j,k,ielem) + dvx_dxl*dux_dxl*deltat + + !* c12 + cijkl_kl(2,i,j,k,ielem) = cijkl_kl(2,i,j,k,ielem) & + +(dvx_dxl * duy_dyl + dvy_dyl * dux_dxl)*deltat + + !* c13 + cijkl_kl(3,i,j,k,ielem) = cijkl_kl(3,i,j,k,ielem) & + +(dvx_dxl * duz_dzl + dvz_dzl * dux_dxl)*deltat + + !* c14 + cijkl_kl(4,i,j,k,ielem) = cijkl_kl(4,i,j,k,ielem) & + +(dvx_dxl * duz_dyl_plus_duy_dzl + dvz_dyl_plus_dvy_dzl * dux_dxl)*deltat + + !* c15 + cijkl_kl(5,i,j,k,ielem) = cijkl_kl(5,i,j,k,ielem) & + +(dvx_dxl * duz_dxl_plus_dux_dzl + dvz_dxl_plus_dvx_dzl * dux_dxl)*deltat + + !* c16 + cijkl_kl(6,i,j,k,ielem) = cijkl_kl(6,i,j,k,ielem) & + +(dvx_dxl * dux_dyl_plus_duy_dxl + dvx_dyl_plus_dvy_dxl * dux_dxl)*deltat + + !* c22 + cijkl_kl(7,i,j,k,ielem) = cijkl_kl(7,i,j,k,ielem) + dvy_dyl*duy_dyl*deltat + + !* c23 + cijkl_kl(8,i,j,k,ielem) = cijkl_kl(8,i,j,k,ielem) & + +(dvy_dyl * duz_dzl + dvz_dzl * duy_dyl)*deltat + + !* c24 + cijkl_kl(9,i,j,k,ielem) = cijkl_kl(9,i,j,k,ielem) & + +(dvy_dyl * duz_dyl_plus_duy_dzl + dvz_dyl_plus_dvy_dzl * duy_dyl)*deltat + + !* c25 + cijkl_kl(10,i,j,k,ielem) = cijkl_kl(10,i,j,k,ielem) & + +(dvy_dyl * duz_dxl_plus_dux_dzl + dvz_dxl_plus_dvx_dzl * duy_dyl)*deltat + + !* c26 + cijkl_kl(11,i,j,k,ielem) = cijkl_kl(11,i,j,k,ielem) & + +(dvy_dyl * dux_dyl_plus_duy_dxl + dvx_dyl_plus_dvy_dxl * duy_dyl)*deltat + + !* c33 + cijkl_kl(12,i,j,k,ielem) = cijkl_kl(12,i,j,k,ielem) + dvz_dzl*duz_dzl*deltat + + !* c34 + cijkl_kl(13,i,j,k,ielem) = cijkl_kl(13,i,j,k,ielem) & + +(dvz_dzl * duz_dyl_plus_duy_dzl + dvz_dyl_plus_dvy_dzl * duz_dzl)*deltat + + !* c35 + cijkl_kl(14,i,j,k,ielem) = cijkl_kl(14,i,j,k,ielem) & + +(dvz_dzl * duz_dxl_plus_dux_dzl + dvz_dxl_plus_dvx_dzl * duz_dzl)*deltat + + !* c36 + cijkl_kl(15,i,j,k,ielem) = cijkl_kl(15,i,j,k,ielem) & + +(dvz_dzl * dux_dyl_plus_duy_dxl + dvx_dyl_plus_dvy_dxl * duz_dzl)*deltat + + !* c44 + cijkl_kl(16,i,j,k,ielem) = cijkl_kl(16,i,j,k,ielem) & + + dvz_dyl_plus_dvy_dzl * duz_dyl_plus_duy_dzl*deltat + + !* c45 + cijkl_kl(17,i,j,k,ielem) = cijkl_kl(17,i,j,k,ielem) & + +(dvz_dyl_plus_dvy_dzl * duz_dxl_plus_dux_dzl & + + dvz_dxl_plus_dvx_dzl * duz_dyl_plus_duy_dzl)*deltat + + !* c46 + cijkl_kl(18,i,j,k,ielem) = cijkl_kl(18,i,j,k,ielem) & + +(dvz_dyl_plus_dvy_dzl * dux_dyl_plus_duy_dxl & + + dvx_dyl_plus_dvy_dxl * duz_dyl_plus_duy_dzl)*deltat + + !* c55 + cijkl_kl(19,i,j,k,ielem) = cijkl_kl(19,i,j,k,ielem) & + + dvz_dxl_plus_dvx_dzl * duz_dxl_plus_dux_dzl*deltat + + !* c56 + cijkl_kl(20,i,j,k,ielem) = cijkl_kl(20,i,j,k,ielem) & + +(dvz_dxl_plus_dvx_dzl * dux_dyl_plus_duy_dxl & + + dvx_dyl_plus_dvy_dxl * duz_dxl_plus_dux_dzl)*deltat + + !* c66 + cijkl_kl(21,i,j,k,ielem) = cijkl_kl(21,i,j,k,ielem) & + + dvx_dyl_plus_dvy_dxl * dux_dyl_plus_duy_dxl*deltat + enddo + enddo + enddo - endif ! ispec_is_elastic + endif ! ispec_is_elastic - enddo + enddo end subroutine compute_anisotropic_kernels_for_velocity_data diff --git a/src/specfem3D/compute_stacey_viscoelastic.F90 b/src/specfem3D/compute_stacey_viscoelastic.F90 index 54dd07e6d..5ea5ed54c 100644 --- a/src/specfem3D/compute_stacey_viscoelastic.F90 +++ b/src/specfem3D/compute_stacey_viscoelastic.F90 @@ -117,87 +117,94 @@ subroutine compute_stacey_viscoelastic(NSPEC_AB,NGLOB_AB,accel, & ! ********************************************************************************* !! CD modif. : begin (implemented by VM) !! For coupling with DSM - integer :: kaxisem, ip + integer,dimension(NGLLSQUARE,num_abs_boundary_faces) :: ipt_table + ! only add these contributions in first pass + if (iphase /= 1) return + ! injecting boundary if (COUPLE_WITH_INJECTION_TECHNIQUE) then - ipt = 0 + ! note: only applies boundary condition in first phase if (INJECTION_TECHNIQUE_TYPE == INJECTION_TECHNIQUE_IS_DSM) then - if (old_DSM_coupling_from_Vadim) then - if (iphase == 1) then - if (mod(it_dsm,Ntime_step_dsm+1) == 0 .or. it == 1) then - call read_dsm_file(Veloc_dsm_boundary,Tract_dsm_boundary,num_abs_boundary_faces,it_dsm) - endif + if (mod(it_dsm,Ntime_step_dsm+1) == 0 .or. it == 1) then + call read_dsm_file(Veloc_dsm_boundary,Tract_dsm_boundary,num_abs_boundary_faces,it_dsm) endif else !! MODIFS DE NOBU 2D endif - else if (INJECTION_TECHNIQUE_TYPE == INJECTION_TECHNIQUE_IS_AXISEM) then - - if (iphase == 1) then - call read_axisem_file(Veloc_axisem,Tract_axisem,num_abs_boundary_faces*NGLLSQUARE) - - !! CD CD add this - if (RECIPROCITY_AND_KH_INTEGRAL) Tract_axisem_time(:,:,it) = Tract_axisem(:,:) - endif - - else if ( INJECTION_TECHNIQUE_TYPE == INJECTION_TECHNIQUE_IS_FK) then - if (iphase == 1) then - - !! find indices - ! example: np_resamp = 2 and it = 1,2,3,4,5,6, .. - ! --> ii = 1,1,2,2,3,3,.. - ii = floor( real(it + NP_RESAMP - 1) / real( NP_RESAMP)) - ! example: --> kk = 1,2,1,2,1,2,,.. - kk = it - (ii-1) * NP_RESAMP - - w = dble(kk-1) / dble(NP_RESAMP) - - ! Cubic spline values - cs(4) = w*w*w/6.d0 - cs(1) = 1.d0/6.d0+w*(w-1.d0)/2.d0-cs(4) - cs(3) = w+cs(1)-2.d0*cs(4) - cs(2) = 1.d0-cs(1)-cs(3)-cs(4) - - cs_single(:) = sngl(cs(:)) - - iim1 = ii-1 - iip1 = ii+1 - iip2 = ii+2 - - do ip=1, npt - - vx_FK(ip) = cs_single(1)* VX_t(ip,iim1) + cs_single(2)* VX_t(ip,ii) + cs_single(3)* VX_t(ip,iip1) + & - cs_single(4)* VX_t(ip,iip2) - vy_FK(ip) = cs_single(1)* VY_t(ip,iim1) + cs_single(2)* VY_t(ip,ii) + cs_single(3)* VY_t(ip,iip1) + & - cs_single(4)* VY_t(ip,iip2) - vz_FK(ip) = cs_single(1)* VZ_t(ip,iim1) + cs_single(2)* VZ_t(ip,ii) + cs_single(3)* VZ_t(ip,iip1) + & - cs_single(4)* VZ_t(ip,iip2) - tx_FK(ip) = cs_single(1)* TX_t(ip,iim1) + cs_single(2)* TX_t(ip,ii) + cs_single(3)* TX_t(ip,iip1) + & - cs_single(4)* TX_t(ip,iip2) - ty_FK(ip) = cs_single(1)* TY_t(ip,iim1) + cs_single(2)* TY_t(ip,ii) + cs_single(3)* TY_t(ip,iip1) + & - cs_single(4)* TY_t(ip,iip2) - tz_FK(ip) = cs_single(1)* TZ_t(ip,iim1) + cs_single(2)* TZ_t(ip,ii) + cs_single(3)* TZ_t(ip,iip1) + & - cs_single(4)* TZ_t(ip,iip2) - - enddo - it_fk=it_fk+1 - - endif - endif - -endif - - ! only add these contributions in first pass - if (iphase /= 1) return + else if (INJECTION_TECHNIQUE_TYPE == INJECTION_TECHNIQUE_IS_AXISEM) then + call read_axisem_file(Veloc_axisem,Tract_axisem,num_abs_boundary_faces*NGLLSQUARE) + + !! CD CD add this + if (RECIPROCITY_AND_KH_INTEGRAL) Tract_axisem_time(:,:,it) = Tract_axisem(:,:) + + else if ( INJECTION_TECHNIQUE_TYPE == INJECTION_TECHNIQUE_IS_FK) then + !! find indices + ! example: np_resamp = 2 and it = 1,2,3,4,5,6, .. + ! --> ii = 1,1,2,2,3,3,.. + ii = floor( real(it + NP_RESAMP - 1) / real( NP_RESAMP)) + ! example: --> kk = 1,2,1,2,1,2,,.. + kk = it - (ii-1) * NP_RESAMP + + w = dble(kk-1) / dble(NP_RESAMP) + + ! Cubic spline values + cs(4) = w*w*w/6.d0 + cs(1) = 1.d0/6.d0+w*(w-1.d0)/2.d0-cs(4) + cs(3) = w+cs(1)-2.d0*cs(4) + cs(2) = 1.d0-cs(1)-cs(3)-cs(4) + + cs_single(:) = sngl(cs(:)) + + iim1 = ii-1 + iip1 = ii+1 + iip2 = ii+2 + + do ip=1, npt + vx_FK(ip) = cs_single(1)* VX_t(ip,iim1) + cs_single(2)* VX_t(ip,ii) + cs_single(3)* VX_t(ip,iip1) + & + cs_single(4)* VX_t(ip,iip2) + vy_FK(ip) = cs_single(1)* VY_t(ip,iim1) + cs_single(2)* VY_t(ip,ii) + cs_single(3)* VY_t(ip,iip1) + & + cs_single(4)* VY_t(ip,iip2) + vz_FK(ip) = cs_single(1)* VZ_t(ip,iim1) + cs_single(2)* VZ_t(ip,ii) + cs_single(3)* VZ_t(ip,iip1) + & + cs_single(4)* VZ_t(ip,iip2) + tx_FK(ip) = cs_single(1)* TX_t(ip,iim1) + cs_single(2)* TX_t(ip,ii) + cs_single(3)* TX_t(ip,iip1) + & + cs_single(4)* TX_t(ip,iip2) + ty_FK(ip) = cs_single(1)* TY_t(ip,iim1) + cs_single(2)* TY_t(ip,ii) + cs_single(3)* TY_t(ip,iip1) + & + cs_single(4)* TY_t(ip,iip2) + tz_FK(ip) = cs_single(1)* TZ_t(ip,iim1) + cs_single(2)* TZ_t(ip,ii) + cs_single(3)* TZ_t(ip,iip1) + & + cs_single(4)* TZ_t(ip,iip2) + enddo + it_fk=it_fk+1 + + ! prepares ipt table + ! (needed also for openmp threads accessing correct ipt index) + ipt = 0 + do iface = 1,num_abs_boundary_faces + ispec = abs_boundary_ispec(iface) + if (ispec_is_elastic(ispec)) then + do igll = 1,NGLLSQUARE + ! counter + ipt = ipt + 1 + ipt_table(igll,iface) = ipt + enddo + endif + enddo + endif + endif ! COUPLE_WITH_INJECTION_TECHNIQUE ! checks if anything to do if (num_abs_boundary_faces == 0) return ! absorbs absorbing-boundary surface using Stacey condition (Clayton and Enquist) +! openmp solver +!$OMP PARALLEL if (num_abs_boundary_faces > 100) & +!$OMP DEFAULT(SHARED) & +!$OMP PRIVATE(iface,ispec,igll,i,j,k,iglob,vx,vy,vz,vn,nx,ny,nz,tx,ty,tz,jacobianw, & +!$OMP kaxisem,ixglob,ipt) +!$OMP DO do iface = 1,num_abs_boundary_faces ispec = abs_boundary_ispec(iface) @@ -218,29 +225,24 @@ subroutine compute_stacey_viscoelastic(NSPEC_AB,NGLOB_AB,accel, & vy = veloc(2,iglob) vz = veloc(3,iglob) - !! CD CD !! For coupling with EXTERNAL CODE - if (COUPLE_WITH_INJECTION_TECHNIQUE) then - - if (INJECTION_TECHNIQUE_TYPE == INJECTION_TECHNIQUE_IS_DSM) then !! To verify for NOBU version - vx = vx - Veloc_dsm_boundary(1,it_dsm,igll,iface) - vy = vy - Veloc_dsm_boundary(2,it_dsm,igll,iface) - vz = vz - Veloc_dsm_boundary(3,it_dsm,igll,iface) - - else if (INJECTION_TECHNIQUE_TYPE == INJECTION_TECHNIQUE_IS_AXISEM) then !! VM VM add AxiSEM - kaxisem = igll + NGLLSQUARE*(iface - 1) - vx = vx - Veloc_axisem(1,kaxisem) - vy = vy - Veloc_axisem(2,kaxisem) - vz = vz - Veloc_axisem(3,kaxisem) - endif - - endif - -! ********************************************************************************* -! added by Ping Tong (TP / Tong Ping) for the FK3D calculation - if (COUPLE_WITH_INJECTION_TECHNIQUE .and. INJECTION_TECHNIQUE_TYPE == INJECTION_TECHNIQUE_IS_FK) then - ipt = ipt + 1 - -!! DEBUGVM pour eviter de stocker pour profiler la vitesse de FK + !! CD CD !! For coupling with EXTERNAL CODE + if (COUPLE_WITH_INJECTION_TECHNIQUE) then + if (INJECTION_TECHNIQUE_TYPE == INJECTION_TECHNIQUE_IS_DSM) then !! To verify for NOBU version + vx = vx - Veloc_dsm_boundary(1,it_dsm,igll,iface) + vy = vy - Veloc_dsm_boundary(2,it_dsm,igll,iface) + vz = vz - Veloc_dsm_boundary(3,it_dsm,igll,iface) + + else if (INJECTION_TECHNIQUE_TYPE == INJECTION_TECHNIQUE_IS_AXISEM) then !! VM VM add AxiSEM + kaxisem = igll + NGLLSQUARE*(iface - 1) + vx = vx - Veloc_axisem(1,kaxisem) + vy = vy - Veloc_axisem(2,kaxisem) + vz = vz - Veloc_axisem(3,kaxisem) + else if (INJECTION_TECHNIQUE_TYPE == INJECTION_TECHNIQUE_IS_FK) then + ! added by Ping Tong (TP / Tong Ping) for the FK3D calculation + ! point index using table lookup + !ipt = ipt + 1 + ipt = ipt_table(igll,iface) + !! DEBUGVM pour eviter de stocker pour profiler la vitesse de FK !vx_FK = vxbd(it_fk,ipt) !vy_FK = vybd(it_fk,ipt) !vz_FK = vzbd(it_fk,ipt) @@ -254,7 +256,7 @@ subroutine compute_stacey_viscoelastic(NSPEC_AB,NGLOB_AB,accel, & vy = vy - vy_FK(ipt) vz = vz - vz_FK(ipt) endif -! ********************************************************************************* + endif ! gets associated normal nx = abs_boundary_normal(1,igll,iface) @@ -269,37 +271,33 @@ subroutine compute_stacey_viscoelastic(NSPEC_AB,NGLOB_AB,accel, & ty = rho_vp(i,j,k,ispec)*vn*ny + rho_vs(i,j,k,ispec)*(vy-vn*ny) tz = rho_vp(i,j,k,ispec)*vn*nz + rho_vs(i,j,k,ispec)*(vz-vn*nz) - !! CD CD !! For coupling with DSM - if (COUPLE_WITH_INJECTION_TECHNIQUE) then - - if (INJECTION_TECHNIQUE_TYPE == INJECTION_TECHNIQUE_IS_DSM) then !! To verify for NOBU version - tx = tx - Tract_dsm_boundary(1,it_dsm,igll,iface) - ty = ty - Tract_dsm_boundary(2,it_dsm,igll,iface) - tz = tz - Tract_dsm_boundary(3,it_dsm,igll,iface) - - else if (INJECTION_TECHNIQUE_TYPE == INJECTION_TECHNIQUE_IS_AXISEM) then - tx = tx - Tract_axisem(1,kaxisem) - ty = ty - Tract_axisem(2,kaxisem) - tz = tz - Tract_axisem(3,kaxisem) - endif - - endif - -! ********************************************************************************* -! added by Ping Tong (TP / Tong Ping) for the FK3D calculation - if (COUPLE_WITH_INJECTION_TECHNIQUE .and. INJECTION_TECHNIQUE_TYPE == INJECTION_TECHNIQUE_IS_FK) then + !! CD CD !! For coupling with DSM + if (COUPLE_WITH_INJECTION_TECHNIQUE) then + if (INJECTION_TECHNIQUE_TYPE == INJECTION_TECHNIQUE_IS_DSM) then !! To verify for NOBU version + tx = tx - Tract_dsm_boundary(1,it_dsm,igll,iface) + ty = ty - Tract_dsm_boundary(2,it_dsm,igll,iface) + tz = tz - Tract_dsm_boundary(3,it_dsm,igll,iface) + else if (INJECTION_TECHNIQUE_TYPE == INJECTION_TECHNIQUE_IS_AXISEM) then + tx = tx - Tract_axisem(1,kaxisem) + ty = ty - Tract_axisem(2,kaxisem) + tz = tz - Tract_axisem(3,kaxisem) + else if (INJECTION_TECHNIQUE_TYPE == INJECTION_TECHNIQUE_IS_FK) then + ! added by Ping Tong (TP / Tong Ping) for the FK3D calculation tx = tx - tx_FK(ipt) ty = ty - ty_FK(ipt) tz = tz - tz_FK(ipt) endif -! ********************************************************************************* + endif ! gets associated, weighted jacobian jacobianw = abs_boundary_jacobian2Dw(igll,iface) ! adds stacey term (weak form) +!$OMP ATOMIC accel(1,iglob) = accel(1,iglob) - tx*jacobianw +!$OMP ATOMIC accel(2,iglob) = accel(2,iglob) - ty*jacobianw +!$OMP ATOMIC accel(3,iglob) = accel(3,iglob) - tz*jacobianw ! adjoint simulations @@ -309,15 +307,30 @@ subroutine compute_stacey_viscoelastic(NSPEC_AB,NGLOB_AB,accel, & b_absorb_field(3,igll,iface) = tz*jacobianw endif !adjoint - !! CD CD added this - if (SAVE_RUN_BOUN_FOR_KH_INTEGRAL) then - write(237) b_absorb_field(1,igll,iface), b_absorb_field(2,igll,iface), b_absorb_field(3,igll,iface) - write(238) displ(1,iglob), displ(2,iglob), displ(3,iglob) - endif - enddo endif ! ispec_is_elastic enddo +!$OMP ENDDO +!$OMP END PARALLEL + + !! CD CD added this + if (SAVE_RUN_BOUN_FOR_KH_INTEGRAL) then + do iface = 1,num_abs_boundary_faces + ispec = abs_boundary_ispec(iface) + if (ispec_is_elastic(ispec)) then + ! reference GLL points on boundary face + do igll = 1,NGLLSQUARE + ! gets local indices for GLL point + i = abs_boundary_ijk(1,igll,iface) + j = abs_boundary_ijk(2,igll,iface) + k = abs_boundary_ijk(3,igll,iface) + iglob = ibool(i,j,k,ispec) + write(237) b_absorb_field(1,igll,iface), b_absorb_field(2,igll,iface), b_absorb_field(3,igll,iface) + write(238) displ(1,iglob), displ(2,iglob), displ(3,iglob) + enddo + endif + enddo + endif ! adjoint simulations: stores absorbed wavefield part if (SIMULATION_TYPE == 1 .and. SAVE_FORWARD) then diff --git a/src/specfem3D/detect_mesh_surfaces.f90 b/src/specfem3D/detect_mesh_surfaces.f90 index 42faed7c3..015c37803 100644 --- a/src/specfem3D/detect_mesh_surfaces.f90 +++ b/src/specfem3D/detect_mesh_surfaces.f90 @@ -46,16 +46,16 @@ subroutine detect_mesh_surfaces() if (MOVIE_SURFACE .or. CREATE_SHAKEMAP) then if (MOVIE_TYPE == 2 .and. PLOT_CROSS_SECTIONS) then call detect_surface_cross_section(NPROC,NGLOB_AB,NSPEC_AB,ibool, & - ispec_is_surface_external_mesh, & - iglob_is_surface_external_mesh, & - nfaces_surface, & - num_interfaces_ext_mesh, & - max_nibool_interfaces_ext_mesh, & - nibool_interfaces_ext_mesh, & - my_neighbors_ext_mesh, & - ibool_interfaces_ext_mesh, & - CROSS_SECTION_X,CROSS_SECTION_Y,CROSS_SECTION_Z, & - xstore,ystore,zstore,myrank) + ispec_is_surface_external_mesh, & + iglob_is_surface_external_mesh, & + nfaces_surface, & + num_interfaces_ext_mesh, & + max_nibool_interfaces_ext_mesh, & + nibool_interfaces_ext_mesh, & + my_neighbors_ext_mesh, & + ibool_interfaces_ext_mesh, & + CROSS_SECTION_X,CROSS_SECTION_Y,CROSS_SECTION_Z, & + xstore,ystore,zstore) endif endif diff --git a/src/specfem3D/initialize_simulation.F90 b/src/specfem3D/initialize_simulation.F90 index 9d9954194..5fdc7bad3 100644 --- a/src/specfem3D/initialize_simulation.F90 +++ b/src/specfem3D/initialize_simulation.F90 @@ -40,8 +40,8 @@ subroutine initialize_simulation() include 'version.fh' - integer ier - logical BROADCAST_AFTER_READ + integer :: ier + logical :: BROADCAST_AFTER_READ ! myrank is the rank of each process, between 0 and NPROC-1. ! as usual in MPI, process 0 is in charge of coordinating everything @@ -54,7 +54,7 @@ subroutine initialize_simulation() ! read the parameter file BROADCAST_AFTER_READ = .true. - call read_parameter_file(myrank,BROADCAST_AFTER_READ) + call read_parameter_file(BROADCAST_AFTER_READ) ! checks flags call initialize_simulation_check() @@ -89,12 +89,13 @@ subroutine initialize_simulation() else write(IMAIN,*) 'using double precision for the calculations' endif + if (FORCE_VECTORIZATION_VAL) write(IMAIN,*) 'using force vectorization' write(IMAIN,*) write(IMAIN,*) 'smallest and largest possible floating-point numbers are: ', & tiny(1._CUSTOM_REAL),huge(1._CUSTOM_REAL) write(IMAIN,*) - write(IMAIN,'(a)',advance='no') ' velocity model: ' + write(IMAIN,'(a)',advance='no') ' velocity model: ' select case (IMODEL) case (IMODEL_DEFAULT) write(IMAIN,'(a)',advance='yes') ' default ' @@ -257,6 +258,12 @@ subroutine initialize_simulation() ! initializes GPU cards if (GPU_MODE) call initialize_GPU() + ! output info for possible OpenMP + call init_openmp() + + ! synchronizes processes + call synchronize_all() + end subroutine initialize_simulation ! diff --git a/src/specfem3D/noise_tomography.f90 b/src/specfem3D/noise_tomography.f90 index 2085693c0..c77eadb46 100644 --- a/src/specfem3D/noise_tomography.f90 +++ b/src/specfem3D/noise_tomography.f90 @@ -234,7 +234,7 @@ end module user_noise_distribution ! ============================================================================================================= ! read parameters - subroutine read_parameters_noise(myrank,nrec,NSTEP,nmovie_points, & + subroutine read_parameters_noise(nrec,NSTEP,nmovie_points, & islice_selected_rec,xi_receiver,eta_receiver,gamma_receiver,nu, & noise_sourcearray,xigll,yigll,zigll, & ibool, & @@ -250,7 +250,7 @@ subroutine read_parameters_noise(myrank,nrec,NSTEP,nmovie_points, & implicit none ! input parameters - integer :: myrank, nrec, NSTEP, nmovie_points + integer :: nrec, NSTEP, nmovie_points integer :: NSPEC_AB_VAL,NGLOB_AB_VAL integer, dimension(nrec) :: islice_selected_rec @@ -305,9 +305,11 @@ subroutine read_parameters_noise(myrank,nrec,NSTEP,nmovie_points, & ! compute source arrays for "ensemble forward source", which is source of "ensemble forward wavefield" if (myrank == islice_selected_rec(irec_master_noise) .or. myrank == 0) then ! myrank == 0 is used for output only - call compute_arrays_source_noise(myrank, & - xi_receiver(irec_master_noise),eta_receiver(irec_master_noise),gamma_receiver(irec_master_noise), & - nu(:,:,irec_master_noise),noise_sourcearray, xigll,yigll,zigll,NSTEP) + call compute_arrays_source_noise(xi_receiver(irec_master_noise), & + eta_receiver(irec_master_noise), & + gamma_receiver(irec_master_noise), & + nu(:,:,irec_master_noise),noise_sourcearray, & + xigll,yigll,zigll,NSTEP) endif ! noise distribution and noise direction @@ -354,15 +356,15 @@ end subroutine read_parameters_noise ! ============================================================================================================= ! check for consistency of the parameters - subroutine check_parameters_noise(myrank,NOISE_TOMOGRAPHY,SIMULATION_TYPE,SAVE_FORWARD, & + subroutine check_parameters_noise(NOISE_TOMOGRAPHY,SIMULATION_TYPE,SAVE_FORWARD, & LOCAL_PATH,NSPEC_TOP,NSTEP) - use constants, only: CUSTOM_REAL,NDIM,NGLLSQUARE,MAX_STRING_LEN,IOUT_NOISE,OUTPUT_FILES + use constants, only: CUSTOM_REAL,NDIM,NGLLSQUARE,MAX_STRING_LEN,IOUT_NOISE,OUTPUT_FILES,myrank implicit none ! input parameters - integer :: myrank,NOISE_TOMOGRAPHY,SIMULATION_TYPE,NSPEC_TOP,NSTEP + integer :: NOISE_TOMOGRAPHY,SIMULATION_TYPE,NSPEC_TOP,NSTEP character(len=MAX_STRING_LEN) :: LOCAL_PATH logical :: SAVE_FORWARD ! local parameters @@ -456,16 +458,17 @@ end subroutine check_parameters_noise ! read and construct the "source" (source time function based upon noise spectrum) ! for "ensemble forward source" - subroutine compute_arrays_source_noise(myrank, & - xi_noise,eta_noise,gamma_noise,nu_single,noise_sourcearray, & + subroutine compute_arrays_source_noise(xi_noise,eta_noise,gamma_noise, & + nu_single,noise_sourcearray, & xigll,yigll,zigll,NSTEP) - use constants, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,NDIM,MAX_STRING_LEN,IIN_NOISE,IOUT_NOISE,OUTPUT_FILES + use constants, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,NDIM,MAX_STRING_LEN, & + IIN_NOISE,IOUT_NOISE,OUTPUT_FILES,myrank implicit none ! input parameters - integer :: myrank, NSTEP + integer :: NSTEP double precision, dimension(NGLLX) :: xigll double precision, dimension(NGLLY) :: yigll double precision, dimension(NGLLZ) :: zigll @@ -557,11 +560,10 @@ end subroutine compute_arrays_source_noise ! step 1: calculate the "ensemble forward source" ! add noise spectrum to the location of master receiver - subroutine add_source_master_rec_noise(myrank,nrec, & - NSTEP,accel,noise_sourcearray, & - ibool,islice_selected_rec,ispec_selected_rec, & - it,irec_master_noise, & - NSPEC_AB_VAL,NGLOB_AB_VAL) + subroutine add_source_master_rec_noise(nrec,NSTEP,accel,noise_sourcearray, & + ibool,islice_selected_rec,ispec_selected_rec, & + it,irec_master_noise, & + NSPEC_AB_VAL,NGLOB_AB_VAL) use constants, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,NDIM @@ -609,12 +611,11 @@ end subroutine add_source_master_rec_noise ! step 1: calculate the "ensemble forward source" ! save surface movie (displacement) at every time steps, for step 2 & 3. - subroutine noise_save_surface_movie(displ, & - ibool, & - noise_surface_movie,it, & - NSPEC_AB_VAL,NGLOB_AB_VAL, & - num_free_surface_faces,free_surface_ispec,free_surface_ijk, & - Mesh_pointer,GPU_MODE) + subroutine noise_save_surface_movie(displ,ibool, & + noise_surface_movie,it, & + NSPEC_AB_VAL,NGLOB_AB_VAL, & + num_free_surface_faces,free_surface_ispec,free_surface_ijk, & + Mesh_pointer,GPU_MODE) use constants, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,NDIM,NGLLSQUARE @@ -678,10 +679,10 @@ end subroutine noise_save_surface_movie ! in step 2, call noise_read_add_surface_movie(..., NSTEP-it+1 ,...) ! in step 3, call noise_read_add_surface_movie(..., it ,...) subroutine noise_read_add_surface_movie(nmovie_points,accel, & - normal_x_noise,normal_y_noise,normal_z_noise,mask_noise, & - ibool,noise_surface_movie,it,NSPEC_AB_VAL,NGLOB_AB_VAL, & - num_free_surface_faces,free_surface_ispec,free_surface_ijk, & - free_surface_jacobian2Dw) + normal_x_noise,normal_y_noise,normal_z_noise,mask_noise, & + ibool,noise_surface_movie,it,NSPEC_AB_VAL,NGLOB_AB_VAL, & + num_free_surface_faces,free_surface_ispec,free_surface_ijk, & + free_surface_jacobian2Dw) use constants @@ -788,12 +789,12 @@ end subroutine noise_read_add_surface_movie_GPU ! step 3: constructing noise source strength kernel subroutine compute_kernels_strength_noise(nmovie_points,ibool, & - sigma_kl,displ,deltat,it, & - normal_x_noise,normal_y_noise,normal_z_noise, & - noise_surface_movie, & - NSPEC_AB_VAL,NGLOB_AB_VAL, & - num_free_surface_faces,free_surface_ispec,free_surface_ijk, & - GPU_MODE,Mesh_pointer) + sigma_kl,displ,deltat,it, & + normal_x_noise,normal_y_noise,normal_z_noise, & + noise_surface_movie, & + NSPEC_AB_VAL,NGLOB_AB_VAL, & + num_free_surface_faces,free_surface_ispec,free_surface_ijk, & + GPU_MODE,Mesh_pointer) use constants @@ -877,14 +878,13 @@ end subroutine compute_kernels_strength_noise ! ============================================================================================================= ! step 3: save noise source strength kernel - subroutine save_kernels_strength_noise(myrank,LOCAL_PATH,sigma_kl,NSPEC_AB_VAL) + subroutine save_kernels_strength_noise(LOCAL_PATH,sigma_kl,NSPEC_AB_VAL) - use constants, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,MAX_STRING_LEN,IOUT_NOISE + use constants, only: CUSTOM_REAL,NGLLX,NGLLY,NGLLZ,MAX_STRING_LEN,IOUT_NOISE,myrank implicit none ! input parameters - integer myrank integer :: NSPEC_AB_VAL real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB_VAL) :: sigma_kl character(len=MAX_STRING_LEN) :: LOCAL_PATH diff --git a/src/specfem3D/prepare_noise.f90 b/src/specfem3D/prepare_noise.f90 index d2e35d06e..7d1937fdf 100644 --- a/src/specfem3D/prepare_noise.f90 +++ b/src/specfem3D/prepare_noise.f90 @@ -86,7 +86,7 @@ subroutine prepare_noise() noise_surface_movie(:,:,:) = 0._CUSTOM_REAL ! sets up noise source for master receiver station - call read_parameters_noise(myrank,nrec,NSTEP,NGLLSQUARE*num_free_surface_faces, & + call read_parameters_noise(nrec,NSTEP,NGLLSQUARE*num_free_surface_faces, & islice_selected_rec,xi_receiver,eta_receiver,gamma_receiver,nu, & noise_sourcearray,xigll,yigll,zigll, & ibool, & @@ -97,9 +97,8 @@ subroutine prepare_noise() ispec_is_acoustic) ! checks flags for noise simulation - call check_parameters_noise(myrank,NOISE_TOMOGRAPHY,SIMULATION_TYPE,SAVE_FORWARD, & - LOCAL_PATH, & - num_free_surface_faces,NSTEP) + call check_parameters_noise(NOISE_TOMOGRAPHY,SIMULATION_TYPE,SAVE_FORWARD, & + LOCAL_PATH,num_free_surface_faces,NSTEP) endif ! synchonizes diff --git a/src/specfem3D/rules.mk b/src/specfem3D/rules.mk index 6ea14e58c..49b504e68 100644 --- a/src/specfem3D/rules.mk +++ b/src/specfem3D/rules.mk @@ -143,6 +143,7 @@ specfem3D_SHARED_OBJECTS = \ $O/gll_library.shared.o \ $O/heap_sort.shared.o \ $O/hex_nodes.shared.o \ + $O/init_openmp.shared.o \ $O/lagrange_poly.shared.o \ $O/netlib_specfun_erf.shared.o \ $O/param_reader.cc.o \ diff --git a/src/specfem3D/save_adjoint_kernels.f90 b/src/specfem3D/save_adjoint_kernels.f90 index 0c239e9b7..83e2745d6 100644 --- a/src/specfem3D/save_adjoint_kernels.f90 +++ b/src/specfem3D/save_adjoint_kernels.f90 @@ -41,7 +41,7 @@ subroutine save_adjoint_kernels() use constants, only: CUSTOM_REAL, NGLLX, NGLLY, NGLLZ - use specfem_par, only: LOCAL_PATH, myrank, sigma_kl, NSPEC_AB, ADIOS_FOR_KERNELS, NOISE_TOMOGRAPHY, NSPEC_ADJOINT, & + use specfem_par, only: LOCAL_PATH, sigma_kl, NSPEC_AB, ADIOS_FOR_KERNELS, NOISE_TOMOGRAPHY, NSPEC_ADJOINT, & APPROXIMATE_HESS_KL, ANISOTROPIC_KL, SAVE_TRANSVERSE_KL use specfem_par_acoustic, only: ACOUSTIC_SIMULATION @@ -157,7 +157,7 @@ end subroutine save_kernels_elastic ! for noise simulations --- noise strength kernel if (NOISE_TOMOGRAPHY == 3) then - call save_kernels_strength_noise(myrank,LOCAL_PATH,sigma_kl,NSPEC_AB) + call save_kernels_strength_noise(LOCAL_PATH,sigma_kl,NSPEC_AB) endif ! for preconditioner diff --git a/src/specfem3D/setup_GLL_points.f90 b/src/specfem3D/setup_GLL_points.f90 index 818411256..ed09a786f 100644 --- a/src/specfem3D/setup_GLL_points.f90 +++ b/src/specfem3D/setup_GLL_points.f90 @@ -55,7 +55,7 @@ subroutine setup_GLL_points() ! checks Courant criteria on mesh if (ELASTIC_SIMULATION) then - call check_mesh_resolution(myrank,NSPEC_AB,NGLOB_AB, & + call check_mesh_resolution(NSPEC_AB,NGLOB_AB, & ibool,xstore,ystore,zstore, & kappastore,mustore,rho_vp,rho_vs, & DT,model_speed_max,min_resolved_period, & @@ -68,7 +68,7 @@ subroutine setup_GLL_points() if (ier /= 0) call exit_MPI_without_rank('error allocating array 2422') rho_vp = 0.0_CUSTOM_REAL rho_vs = 0.0_CUSTOM_REAL - call check_mesh_resolution_poro(myrank,NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, & + call check_mesh_resolution_poro(NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, & DT,model_speed_max,min_resolved_period, & phistore,tortstore,rhoarraystore,rho_vpI,rho_vpII,rho_vsI, & LOCAL_PATH,SAVE_MESH_FILES) @@ -82,7 +82,7 @@ subroutine setup_GLL_points() if (ier /= 0) stop 'error allocating array rho_vs' rho_vp = sqrt( kappastore / rhostore ) * rhostore rho_vs = 0.0_CUSTOM_REAL - call check_mesh_resolution(myrank,NSPEC_AB,NGLOB_AB, & + call check_mesh_resolution(NSPEC_AB,NGLOB_AB, & ibool,xstore,ystore,zstore, & kappastore,mustore,rho_vp,rho_vs, & DT,model_speed_max,min_resolved_period, & diff --git a/src/specfem3D/setup_sources_receivers.f90 b/src/specfem3D/setup_sources_receivers.f90 index ccf31c4ac..69ebe830e 100644 --- a/src/specfem3D/setup_sources_receivers.f90 +++ b/src/specfem3D/setup_sources_receivers.f90 @@ -1499,7 +1499,7 @@ subroutine setup_sources_receivers_VTKfile() xil = xi_source(isource) etal = eta_source(isource) gammal = gamma_source(isource) - call eval_shape3D_single(myrank,shape3D,xil,etal,gammal,NGNOD) + call eval_shape3D_single(shape3D,xil,etal,gammal,NGNOD) ! interpolates source locations xmesh = 0.0 @@ -1544,7 +1544,7 @@ subroutine setup_sources_receivers_VTKfile() xil = xi_receiver(irec) etal = eta_receiver(irec) gammal = gamma_receiver(irec) - call eval_shape3D_single(myrank,shape3D,xil,etal,gammal,NGNOD) + call eval_shape3D_single(shape3D,xil,etal,gammal,NGNOD) ! interpolates receiver locations xmesh = 0.0 diff --git a/src/specfem3D/specfem3D_par.F90 b/src/specfem3D/specfem3D_par.F90 index db6074585..70fab482a 100644 --- a/src/specfem3D/specfem3D_par.F90 +++ b/src/specfem3D/specfem3D_par.F90 @@ -105,7 +105,7 @@ module specfem_par real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_ocean_load ! time scheme - real(kind=CUSTOM_REAL) deltat,deltatover2,deltatsqover2 + real(kind=CUSTOM_REAL) :: deltat,deltatover2,deltatsqover2 ! LDDRK time scheme integer :: NSTAGE_TIME_SCHEME,istage @@ -186,7 +186,7 @@ module specfem_par real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: wgllwgll_xy_3D,wgllwgll_xz_3D,wgllwgll_yz_3D ! proc numbers for MPI - integer :: myrank, sizeprocs + integer :: sizeprocs ! timer MPI double precision :: time_start @@ -302,6 +302,13 @@ module specfem_par double precision, dimension(NTOTAL_OBSERVATION) :: x_observation,y_observation,z_observation, & g_x,g_y,g_z,G_xx,G_yy,G_zz,G_xy,G_xz,G_yz,temporary_array_for_sum + ! force vectorization +#ifdef FORCE_VECTORIZATION + logical, parameter :: FORCE_VECTORIZATION_VAL = .true. +#else + logical, parameter :: FORCE_VECTORIZATION_VAL = .false. +#endif + #ifdef VTK_VIS ! VTK window mode, default is off logical :: VTK_MODE = .false. diff --git a/src/specfem3D/update_displacement_scheme.f90 b/src/specfem3D/update_displacement_scheme.f90 index 1d1f13265..1825a7eea 100644 --- a/src/specfem3D/update_displacement_scheme.f90 +++ b/src/specfem3D/update_displacement_scheme.f90 @@ -200,39 +200,59 @@ subroutine update_displacement_elastic() ! updates elastic displacement and velocity if (PML_CONDITIONS .and. NSPEC_CPML > 0) then - do ispec_cpml=1,NSPEC_CPML - ispec = CPML_to_spec(ispec_cpml) - do k = 1, NGLLZ - do j = 1, NGLLY - do i = 1, NGLLX - iglob = ibool(i,j,k,ispec) - PML_displ_old(:,i,j,k,ispec_cpml)=displ(:,iglob) & - + deltatover2 * (1._CUSTOM_REAL - 2._CUSTOM_REAL*theta) * veloc(:,iglob) & - + deltatsqover2 * (1._CUSTOM_REAL - theta) * accel(:,iglob) - - enddo - enddo - enddo + do ispec_cpml=1,NSPEC_CPML + ispec = CPML_to_spec(ispec_cpml) + do k = 1, NGLLZ + do j = 1, NGLLY + do i = 1, NGLLX + iglob = ibool(i,j,k,ispec) + PML_displ_old(:,i,j,k,ispec_cpml) = displ(:,iglob) & + + deltatover2 * (1._CUSTOM_REAL - 2._CUSTOM_REAL*theta) * veloc(:,iglob) & + + deltatsqover2 * (1._CUSTOM_REAL - theta) * accel(:,iglob) + + enddo + enddo enddo + enddo endif - displ(:,:) = displ(:,:) + deltat * veloc(:,:) + deltatsqover2 * accel(:,:) - veloc(:,:) = veloc(:,:) + deltatover2 * accel(:,:) + if (SIMULATION_TYPE /= 1) accel_adj_coupling(:,:) = accel(:,:) - accel(:,:) = 0._CUSTOM_REAL + + if (FORCE_VECTORIZATION_VAL) then +! openmp solver +!$OMP PARALLEL DEFAULT(NONE) & +!$OMP SHARED( NGLOB_AB, displ, veloc, accel, & +!$OMP deltat, deltatsqover2, deltatover2 ) & +!$OMP PRIVATE(i) +!$OMP DO + do i = 1,NDIM * NGLOB_AB + displ(i,1) = displ(i,1) + deltat * veloc(i,1) + deltatsqover2 * accel(i,1) + veloc(i,1) = veloc(i,1) + deltatover2 * accel(i,1) + accel(i,1) = 0._CUSTOM_REAL + enddo +!$OMP ENDDO +!$OMP END PARALLEL + else + do i = 1,NGLOB_AB + displ(:,i) = displ(:,i) + deltat * veloc(:,i) + deltatsqover2 * accel(:,i) + veloc(:,i) = veloc(:,i) + deltatover2 * accel(:,i) + accel(:,i) = 0._CUSTOM_REAL + enddo + endif + if (PML_CONDITIONS .and. NSPEC_CPML > 0) then - do ispec_cpml=1,NSPEC_CPML - ispec = CPML_to_spec(ispec_cpml) - do k = 1, NGLLZ - do j = 1, NGLLY - do i = 1, NGLLX - iglob = ibool(i,j,k,ispec) - PML_displ_new(:,i,j,k,ispec_cpml) = displ(:,iglob)& - + deltatover2 * (1._CUSTOM_REAL - 2.0_CUSTOM_REAL*theta) * veloc(:,iglob) - - enddo - enddo - enddo + do ispec_cpml=1,NSPEC_CPML + ispec = CPML_to_spec(ispec_cpml) + do k = 1, NGLLZ + do j = 1, NGLLY + do i = 1, NGLLX + iglob = ibool(i,j,k,ispec) + PML_displ_new(:,i,j,k,ispec_cpml) = displ(:,iglob) & + + deltatover2 * (1._CUSTOM_REAL - 2.0_CUSTOM_REAL*theta) * veloc(:,iglob) + enddo + enddo enddo + enddo endif ! adjoint simulations diff --git a/src/specfem3D/vtk_window.F90 b/src/specfem3D/vtk_window.F90 index 80ece9ac4..4086a2956 100644 --- a/src/specfem3D/vtk_window.F90 +++ b/src/specfem3D/vtk_window.F90 @@ -183,7 +183,7 @@ subroutine vtk_window_prepare_source() xil = xi_source(isource) etal = eta_source(isource) gammal = gamma_source(isource) - call eval_shape3D_single(myrank,shape3D,xil,etal,gammal,NGNOD) + call eval_shape3D_single(shape3D,xil,etal,gammal,NGNOD) ! interpolates source locations xmesh = 0.0 @@ -281,7 +281,7 @@ subroutine vtk_window_prepare_receivers() xil = xi_receiver(irec) etal = eta_receiver(irec) gammal = gamma_receiver(irec) - call eval_shape3D_single(myrank,shape3D,xil,etal,gammal,NGNOD) + call eval_shape3D_single(shape3D,xil,etal,gammal,NGNOD) ! interpolates receiver locations xmesh = 0.0 diff --git a/src/specfem3D/write_movie_output.F90 b/src/specfem3D/write_movie_output.F90 index 15c762ce9..acf8bc9cc 100644 --- a/src/specfem3D/write_movie_output.F90 +++ b/src/specfem3D/write_movie_output.F90 @@ -697,14 +697,14 @@ subroutine wmo_movie_div_curl(veloc, & valence(:) = 0 ! loops over elements - do ispec=1,NSPEC_AB + do ispec = 1,NSPEC_AB if (.not. ispec_is(ispec)) cycle ispec_irreg = irregular_element_number(ispec) ! calculates divergence and curl of velocity field - do k=1,NGLLZ - do j=1,NGLLY - do i=1,NGLLX + do k = 1,NGLLZ + do j = 1,NGLLY + do i = 1,NGLLX tempx1l = 0._CUSTOM_REAL tempx2l = 0._CUSTOM_REAL tempx3l = 0._CUSTOM_REAL @@ -741,7 +741,7 @@ subroutine wmo_movie_div_curl(veloc, & !! DK DK Oct 2018: however this curl and div calculation routine for movies is almost never called !! DK DK Oct 2018: however this curl and div calculation routine for movies is almost never called - do l=1,NGLLX + do l = 1,NGLLX hp1 = hprime_xx(i,l) iglob = ibool(l,j,k,ispec) tempx1l = tempx1l + veloc(1,iglob)*hp1 @@ -830,7 +830,7 @@ subroutine wmo_movie_div_curl(veloc, & enddo enddo !NSPEC_AB - do i=1,NGLOB_AB + do i = 1,NGLOB_AB ! checks if point has a contribution ! note: might not be the case for points in acoustic elements if (valence(i) /= 0) then diff --git a/src/tomography/add_model_iso.f90 b/src/tomography/add_model_iso.f90 index c35fe8edf..87189b509 100644 --- a/src/tomography/add_model_iso.f90 +++ b/src/tomography/add_model_iso.f90 @@ -221,7 +221,7 @@ subroutine initialize() ! reads the parameter file BROADCAST_AFTER_READ = .true. - call read_parameter_file(myrank,BROADCAST_AFTER_READ) + call read_parameter_file(BROADCAST_AFTER_READ) if (ADIOS_ENABLED) stop 'Flag ADIOS_ENABLED not supported yet for xadd_model, please rerun program...' diff --git a/src/tomography/model_update.f90 b/src/tomography/model_update.f90 index 4b6c80566..1eb176a47 100644 --- a/src/tomography/model_update.f90 +++ b/src/tomography/model_update.f90 @@ -582,7 +582,7 @@ subroutine initialize() ! reads the parameter file BROADCAST_AFTER_READ = .true. - call read_parameter_file(myrank,BROADCAST_AFTER_READ) + call read_parameter_file(BROADCAST_AFTER_READ) if (ADIOS_ENABLED) stop 'Flag ADIOS_ENABLED set to .true. not supported yet for xmodel_update, please rerun program...' @@ -755,7 +755,7 @@ subroutine get_external_mesh() open(unit=IMAIN,file=trim(OUTPUT_MODEL_DIR)//'/output_mesh_resolution_initial.txt',status='unknown') if (ELASTIC_SIMULATION) then - call check_mesh_resolution(myrank,NSPEC_AB,NGLOB_AB, & + call check_mesh_resolution(NSPEC_AB,NGLOB_AB, & ibool,xstore,ystore,zstore, & kappastore,mustore,rho_vp,rho_vs, & DT,model_speed_max,min_resolved_period, & @@ -768,7 +768,7 @@ subroutine get_external_mesh() if (ier /= 0) call exit_MPI_without_rank('error allocating array 948') rho_vp = 0.0_CUSTOM_REAL rho_vs = 0.0_CUSTOM_REAL - call check_mesh_resolution_poro(myrank,NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, & + call check_mesh_resolution_poro(NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, & DT,model_speed_max,min_resolved_period, & phistore,tortstore,rhoarraystore,rho_vpI,rho_vpII,rho_vsI, & LOCAL_PATH,SAVE_MESH_FILES) @@ -782,7 +782,7 @@ subroutine get_external_mesh() if (ier /= 0) stop 'Error allocating array rho_vs' rho_vp = sqrt( kappastore / rhostore ) * rhostore rho_vs = 0.0_CUSTOM_REAL - call check_mesh_resolution(myrank,NSPEC_AB,NGLOB_AB, & + call check_mesh_resolution(NSPEC_AB,NGLOB_AB, & ibool,xstore,ystore,zstore, & kappastore,mustore,rho_vp,rho_vs, & DT,model_speed_max,min_resolved_period, & @@ -953,7 +953,7 @@ subroutine save_new_databases() ! calculate min_resolved_period (needed for attenuation model) if (ELASTIC_SIMULATION) then - call check_mesh_resolution(myrank,NSPEC_AB,NGLOB_AB, & + call check_mesh_resolution(NSPEC_AB,NGLOB_AB, & ibool,xstore,ystore,zstore, & kappastore_new,mustore_new,rho_vp_new,rho_vs_new, & -1.0d0,model_speed_max,min_resolved_period, & @@ -1070,7 +1070,7 @@ subroutine save_new_databases() call synchronize_all() ! calculates and stores attenuation arrays - call get_attenuation_model(myrank,NSPEC_AB,USE_OLSEN_ATTENUATION,OLSEN_ATTENUATION_RATIO, & + call get_attenuation_model(NSPEC_AB,USE_OLSEN_ATTENUATION,OLSEN_ATTENUATION_RATIO, & mustore_new,rho_vs_new,kappastore_new,rho_vp_new, & qkappa_attenuation_store,qmu_attenuation_store, & ispec_is_elastic,min_resolved_period,prname_new,ATTENUATION_f0_REFERENCE) diff --git a/src/tomography/postprocess_sensitivity_kernels/clip_sem.f90 b/src/tomography/postprocess_sensitivity_kernels/clip_sem.f90 index 3bae0b394..58d9fc964 100644 --- a/src/tomography/postprocess_sensitivity_kernels/clip_sem.f90 +++ b/src/tomography/postprocess_sensitivity_kernels/clip_sem.f90 @@ -117,7 +117,7 @@ program clip_sem ! read simulation parameters BROADCAST_AFTER_READ = .true. - call read_parameter_file(myrank,BROADCAST_AFTER_READ) + call read_parameter_file(BROADCAST_AFTER_READ) ! checks number of MPI processes if (sizeprocs /= NPROC) then diff --git a/src/tomography/postprocess_sensitivity_kernels/combine_sem.f90 b/src/tomography/postprocess_sensitivity_kernels/combine_sem.f90 index 3c937ba45..fc57de9e2 100644 --- a/src/tomography/postprocess_sensitivity_kernels/combine_sem.f90 +++ b/src/tomography/postprocess_sensitivity_kernels/combine_sem.f90 @@ -124,7 +124,7 @@ program combine_sem ! read simulation parameters BROADCAST_AFTER_READ = .true. - call read_parameter_file(myrank,BROADCAST_AFTER_READ) + call read_parameter_file(BROADCAST_AFTER_READ) ! checks number of MPI processes if (sizeprocs /= NPROC) then diff --git a/src/tomography/postprocess_sensitivity_kernels/smooth_sem.F90 b/src/tomography/postprocess_sensitivity_kernels/smooth_sem.F90 index 3a62e64e6..69f25fb33 100644 --- a/src/tomography/postprocess_sensitivity_kernels/smooth_sem.F90 +++ b/src/tomography/postprocess_sensitivity_kernels/smooth_sem.F90 @@ -230,7 +230,7 @@ program smooth_sem ! reads the parameter file BROADCAST_AFTER_READ = .true. - call read_parameter_file(myrank,BROADCAST_AFTER_READ) + call read_parameter_file(BROADCAST_AFTER_READ) if (ADIOS_ENABLED) stop 'Flag ADIOS_ENABLED not supported yet for smoothing, please rerun program...' @@ -353,7 +353,7 @@ program smooth_sem endif if (ELASTIC_SIMULATION) then - call check_mesh_resolution(myrank,NSPEC_AB,NGLOB_AB, & + call check_mesh_resolution(NSPEC_AB,NGLOB_AB, & ibool,xstore,ystore,zstore, & kappastore,mustore,rho_vp,rho_vs, & DT,model_speed_max,min_resolved_period, & @@ -366,7 +366,7 @@ program smooth_sem if (ier /= 0) call exit_MPI_without_rank('error allocating array 1010') rho_vp = 0.0_CUSTOM_REAL rho_vs = 0.0_CUSTOM_REAL - call check_mesh_resolution_poro(myrank,NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, & + call check_mesh_resolution_poro(NSPEC_AB,NGLOB_AB,ibool,xstore,ystore,zstore, & DT,model_speed_max,min_resolved_period, & phistore,tortstore,rhoarraystore,rho_vpI,rho_vpII,rho_vsI, & LOCAL_PATH,SAVE_MESH_FILES) @@ -380,7 +380,7 @@ program smooth_sem if (ier /= 0) stop 'Error allocating array rho_vs' rho_vp = sqrt( kappastore / rhostore ) * rhostore rho_vs = 0.0_CUSTOM_REAL - call check_mesh_resolution(myrank,NSPEC_AB,NGLOB_AB, & + call check_mesh_resolution(NSPEC_AB,NGLOB_AB, & ibool,xstore,ystore,zstore, & kappastore,mustore,rho_vp,rho_vs, & DT,model_speed_max,min_resolved_period, & diff --git a/src/tomography/rules.mk b/src/tomography/rules.mk index c660d1a65..5237f7e75 100644 --- a/src/tomography/rules.mk +++ b/src/tomography/rules.mk @@ -173,6 +173,7 @@ xmodel_update_SHARED_OBJECTS = \ $O/exit_mpi.shared.o \ $O/get_attenuation_model.shared.o \ $O/gll_library.shared.o \ + $O/init_openmp.shared.o \ $O/param_reader.cc.o \ $O/read_parameter_file.shared.o \ $O/read_value_parameters.shared.o \ diff --git a/src/tomography/sum_kernels.f90 b/src/tomography/sum_kernels.f90 index c52a95696..b586eeab7 100644 --- a/src/tomography/sum_kernels.f90 +++ b/src/tomography/sum_kernels.f90 @@ -101,7 +101,7 @@ program sum_kernels ! needs local_path for mesh files BROADCAST_AFTER_READ = .true. - call read_parameter_file(myrank,BROADCAST_AFTER_READ) + call read_parameter_file(BROADCAST_AFTER_READ) ! checks if number of MPI process as specified if (sizeprocs /= NPROC) then diff --git a/src/tomography/sum_preconditioned_kernels.f90 b/src/tomography/sum_preconditioned_kernels.f90 index bb9cf7ee2..db1d17635 100644 --- a/src/tomography/sum_preconditioned_kernels.f90 +++ b/src/tomography/sum_preconditioned_kernels.f90 @@ -100,7 +100,7 @@ program sum_preconditioned_kernels ! needs local_path for mesh files BROADCAST_AFTER_READ = .true. - call read_parameter_file(myrank,BROADCAST_AFTER_READ) + call read_parameter_file(BROADCAST_AFTER_READ) ! checks if number of MPI process as specified if (sizeprocs /= NPROC) then diff --git a/tests/decompose_mesh/test_read.f90 b/tests/decompose_mesh/test_read.f90 index 77c507160..42d3b3247 100644 --- a/tests/decompose_mesh/test_read.f90 +++ b/tests/decompose_mesh/test_read.f90 @@ -22,7 +22,7 @@ program test_read ! reads ../DATA/Par_file myrank = 0 BROADCAST_AFTER_READ = .false. - call read_parameter_file(myrank,BROADCAST_AFTER_READ) + call read_parameter_file(BROADCAST_AFTER_READ) ! punctual check of values for given default Par_file in SPECFEM3D/DATA/ directory print *,'NPROC = ',NPROC