diff --git a/src/inverse_problem_for_model/elastic_tensor_tools_mod.f90 b/src/inverse_problem_for_model/elastic_tensor_tools_mod.f90 index 1e8161b57..cc938fa72 100644 --- a/src/inverse_problem_for_model/elastic_tensor_tools_mod.f90 +++ b/src/inverse_problem_for_model/elastic_tensor_tools_mod.f90 @@ -71,12 +71,15 @@ subroutine define_rotation_matrix(a0,b0,c0,rotmat,rotmat_t) rotmat_t = transpose(rotmat) end subroutine define_rotation_matrix - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Define bond rotation matrices to rotate a voigt tensor ! (see e.g. Auld 1973, for bond matrix definition secion D pages 73-76) ! HERE THE STRESS ONE => can be used to rotate the stress vector or stiffness tensor + function define_bond_stress_matrix(rotmat) result(bond) real(kind=dp), dimension(3,3), intent(in) :: rotmat @@ -131,12 +134,15 @@ function define_bond_stress_matrix(rotmat) result(bond) bond(6,6) = rotmat(1,1)*rotmat(2,2) + rotmat(1,2)*rotmat(2,1) end function define_bond_stress_matrix - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Define bond rotation matrices to rotate a voigt tensor ! (see e.g. Auld 1973, for bond matrix definition secion D pages 73-76) ! HERE THE STRAIN ONE => can be used to rotate the strain vector or compliance tensor + function define_bond_strain_matrix(rotmat) result(bond) real(kind=dp), dimension(3,3), intent(in) :: rotmat @@ -191,10 +197,13 @@ function define_bond_strain_matrix(rotmat) result(bond) bond(6,6) = rotmat(1,1)*rotmat(2,2) + rotmat(1,2)*rotmat(2,1) end function define_bond_strain_matrix - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Rotation of first order tensor (vector actually) + subroutine rotate_vector(rotmat,vi,vi_r) real(kind=dp), dimension(3), intent(in) :: vi @@ -212,10 +221,13 @@ subroutine rotate_vector(rotmat,vi,vi_r) enddo end subroutine rotate_vector - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Rotation of second order tensor (not efficient but corresponds to definition) + function rotate_second_order_tensor(rotmat,cij) result(cij_r) real(kind=dp), dimension(3,3), intent(in) :: cij @@ -237,10 +249,13 @@ function rotate_second_order_tensor(rotmat,cij) result(cij_r) enddo end function rotate_second_order_tensor - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Rotation of second order tensor (very not efficient but corresponds to definition) + subroutine rotate_fourth_order_tensor(rotmat,cijkl,cijkl_r) real(kind=dp), dimension(3,3,3,3), intent(in) :: cijkl @@ -272,11 +287,14 @@ subroutine rotate_fourth_order_tensor(rotmat,cijkl,cijkl_r) enddo end subroutine rotate_fourth_order_tensor - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Rotation of fourth order tensor in voigt matrix with bond matrix ! (see e.g. Auld 1973, for bond matrix definition) + function rotate_tensor_with_bond_matrix(bond,tensor) result(tensor_r) real(kind=dp), dimension(6,6), intent(in) :: tensor @@ -310,11 +328,14 @@ function rotate_tensor_with_bond_matrix(bond,tensor) result(tensor_r) enddo end function rotate_tensor_with_bond_matrix - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Get dilatational stiffness tensor according to Browaeys and Chevrot (2004) ! (four-rank stiffness tensor has two two-rank tensors contractions) + function get_dilatational_stiffness_tensor(cij) result(dilatational) real(kind=dp), dimension(6,6), intent(in) :: cij @@ -336,10 +357,13 @@ function get_dilatational_stiffness_tensor(cij) result(dilatational) dilatational(3,3) = cij(1,3) + cij(2,3) + cij(3,3) end function get_dilatational_stiffness_tensor - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Get voigt stiffness tensor according to Browaeys and Chevrot (2004) + function get_voigt_stiffness_tensor(cij) result(voigt) real(kind=dp), dimension(6,6), intent(in) :: cij @@ -361,10 +385,13 @@ function get_voigt_stiffness_tensor(cij) result(voigt) voigt(3,3) = cij(5,5) + cij(4,4) + cij(3,3) end function get_voigt_stiffness_tensor - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Pass triclinic elastic vector to isotropic one according to fedorov (1968) + subroutine get_isotropic_part_fedorov(triclinic,isotropic) real(kind=dp), dimension(21), intent(in) :: triclinic @@ -388,10 +415,13 @@ subroutine get_isotropic_part_fedorov(triclinic,isotropic) isotropic(7:9) = mu end subroutine get_isotropic_part_fedorov - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Define isotropic tensor with analytical formula (not efficient but still usefull) + subroutine define_isotropic_tensor_1(lambda,mu,tensor) real(kind=dp), intent(in) :: lambda, mu @@ -425,10 +455,13 @@ subroutine define_isotropic_tensor_1(lambda,mu,tensor) enddo end subroutine define_isotropic_tensor_1 - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Define hexagonal tensor with analytical fromula from cij values + subroutine define_hexagonal_tensor_1(c11,c33,c44,c66,c13,s,tensor) real(kind=dp), intent(in) :: c11, c33, c44, c66, c13 @@ -489,11 +522,13 @@ subroutine define_hexagonal_tensor_1(c11,c33,c44,c66,c13,s,tensor) enddo end subroutine define_hexagonal_tensor_1 - !-------------------------------------------------------------------------------- +! +!------------------------------------------------------------------------------------------------------------- +! - !================================================================================ ! Define hexagonal tensor with analytical fromula from thomsen parameters + subroutine define_hexagonal_tensor_2(c33,c44,eps,del,gam,s,tensor) real(kind=dp), intent(in) :: c33, c44, eps, del, gam @@ -555,11 +590,14 @@ subroutine define_hexagonal_tensor_2(c33,c44,eps,del,gam,s,tensor) enddo end subroutine define_hexagonal_tensor_2 - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Define orthorhombic tensor with analytical fromula from cij values ! a, b, and c are normal vector to the three mirror symmetry planes + subroutine define_orthorhombic_tensor1(c11,c22,c33,c23,c13,c12,c44,c55,c66,a,b,c,tensor) real(kind=dp), intent(in) :: c11, c22, c33, c44, c55, c66, c13, c23, c12 @@ -619,11 +657,13 @@ subroutine define_orthorhombic_tensor1(c11,c22,c33,c23,c13,c12,c44,c55,c66,a,b,c enddo end subroutine define_orthorhombic_tensor1 - !-------------------------------------------------------------------------------- +! +!------------------------------------------------------------------------------------------------------------- +! - !================================================================================ ! Routine for initialisation of tensor, voigt matrix and elastic vector indexing + subroutine define_indexing_vec_to_tens() implicit none @@ -715,10 +755,13 @@ subroutine define_indexing_vec_to_tens() endif end subroutine define_indexing_vec_to_tens - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Get direction cosines (th = dip angle from vertical, ph = azimuth from north (y)) + subroutine get_direction_cosines(th,ph,s) real(kind=dp), intent(in) :: th, ph @@ -734,11 +777,14 @@ subroutine get_direction_cosines(th,ph,s) s(3) = cos(thrad) end subroutine get_direction_cosines - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Get angles of symetry axis from direction cosines ! (th = dip angle from vertical, ph = azimuth from north (y)) + subroutine get_symmetry_angles(s,th,ph) real(kind=dp), dimension(3), intent(in) :: s @@ -754,10 +800,13 @@ subroutine get_symmetry_angles(s,th,ph) ph = phrad*rad2deg end subroutine get_symmetry_angles - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Small function for kronecker delta function + function delta(i,j) result(d) integer(kind=kindsi), intent(in) :: i, j @@ -770,10 +819,13 @@ function delta(i,j) result(d) endif end function delta - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Compute square norm of vector + function square_norm_vector(vector) result(norm) real(kind=dp), dimension(3), intent(in) :: vector @@ -788,10 +840,13 @@ function square_norm_vector(vector) result(norm) enddo end function square_norm_vector - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Compute square norm of a second order tensor + function square_norm_tensor_2(tensor) result(norm) real(kind=dp), dimension(3,3), intent(in) :: tensor @@ -808,10 +863,13 @@ function square_norm_tensor_2(tensor) result(norm) enddo end function square_norm_tensor_2 - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Compute square norm of a second order tensor + function square_norm_tensor_4(tensor) result(norm) real(kind=dp), dimension(3,3,3,3), intent(in) :: tensor @@ -832,10 +890,13 @@ function square_norm_tensor_4(tensor) result(norm) enddo end function square_norm_tensor_4 - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Pass a voigt matrix to elastic vector + function transform_voigt_matrix_to_vector(cij) result(vi) real(kind=dp), dimension(6,6), intent(in) :: cij @@ -881,10 +942,13 @@ function transform_voigt_matrix_to_vector(cij) result(vi) vi(21) = cij(4,5) * two_sqrt_two end function transform_voigt_matrix_to_vector - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Pass a voigt tensor to elastic vector + function transform_kelvin_tensor_to_vector(cij) result(vi) real(kind=dp), dimension(6,6), intent(in) :: cij @@ -929,12 +993,14 @@ function transform_kelvin_tensor_to_vector(cij) result(vi) vi(21) = cij(4,5) * sqrt_two end function transform_kelvin_tensor_to_vector - !-------------------------------------------------------------------------------- +! +!------------------------------------------------------------------------------------------------------------- +! - !================================================================================ ! Pass an elastic vector to a voigt tensor ! (see Browaeys and Chevrot (2004) + function transform_vector_to_voigt_tensor(vi) result(cij) real(kind=dp), dimension(21), intent(in) :: vi @@ -995,10 +1061,13 @@ function transform_vector_to_voigt_tensor(vi) result(cij) cij(6,6) = vi( 9) * 0.5_dp end function transform_vector_to_voigt_tensor - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Get Voigt m index from ij + function voigt_index(i,j) result(m) integer(kind=kindsi), intent(in) :: i, j @@ -1009,10 +1078,13 @@ function voigt_index(i,j) result(m) m = dij*i + (1-dij)*(9-i-j) end function voigt_index - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Transform fourth order tensor to second order voigt matrix + function transform_tensor_fourth_to_voigt_matrix(cijkl) result(cij) real(kind=dp), dimension(3,3,3,3), intent(in) :: cijkl @@ -1088,10 +1160,13 @@ function transform_tensor_fourth_to_voigt_matrix(cijkl) result(cij) cijkl(1,2,2,1) + cijkl(1,2,2,1)) * 0.25_dp end function transform_tensor_fourth_to_voigt_matrix - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Transform second order voigt tensor to fourth order tensor + function transform_voigt_matrix_to_fourth_order_tensor(cmn) result(cijkl) real(kind=dp), dimension(3,3), intent(in) :: cmn @@ -1122,10 +1197,13 @@ function transform_voigt_matrix_to_fourth_order_tensor(cmn) result(cijkl) enddo end function transform_voigt_matrix_to_fourth_order_tensor - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Transfrom voigt matrix 6x6 to kelvin tensor 6x6 + function voigt_matrix_to_kelvin_tensor(voigt) result(kelvin) real(kind=dp), dimension(6,6), intent(in) :: voigt @@ -1143,10 +1221,13 @@ function voigt_matrix_to_kelvin_tensor(voigt) result(kelvin) kelvin(4:6,4:6) = voigt(4:6,4:6) * 2._dp end function voigt_matrix_to_kelvin_tensor - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Transfrom voigt matrix 6x6 to kelvin tensor 6x6 + function kelvin_tensor_to_voigt_matrix(kelvin) result(voigt) real(kind=dp), dimension(6,6), intent(in) :: kelvin @@ -1164,11 +1245,14 @@ function kelvin_tensor_to_voigt_matrix(kelvin) result(voigt) voigt(4:6,4:6) = kelvin(4:6,4:6) * 0.5_dp end function kelvin_tensor_to_voigt_matrix - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Define projections of elastic tensor (Broawaeys and Cehvrot (2004)) ! (appendix A) + subroutine projection_to_higher_symmetry_class(vi,proj_type,vp,vd,dev) character(len=*), intent(in) :: proj_type @@ -1249,11 +1333,14 @@ subroutine projection_to_higher_symmetry_class(vi,proj_type,vp,vd,dev) dev = sqrt(sum(vd*vd)) end subroutine projection_to_higher_symmetry_class - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Projection of an anisotropic vector onto the isotropic subspace ! use vector + cij contractions to define viso and elastic constant K and G + subroutine projection_to_isotropic_symmetry(vi,di,vo,vp,vd,dev,k,g) real(kind=dp), dimension(3,3), intent(in) :: di, vo ! dilatational and voigt @@ -1292,11 +1379,13 @@ subroutine projection_to_isotropic_symmetry(vi,di,vo,vp,vd,dev,k,g) dev = sqrt(sum(vd*vd)) end subroutine projection_to_isotropic_symmetry - !-------------------------------------------------------------------------------- +! +!------------------------------------------------------------------------------------------------------------- +! - !================================================================================ ! Find symmetry axis of tensor + subroutine determine_tensor_symmetry_axis(cij,scc) real(kind=dp), dimension(6,6), intent(in) :: cij @@ -1418,10 +1507,13 @@ subroutine determine_tensor_symmetry_axis(cij,scc) call projection_to_higher_symmetry_class(vi,'hexagonal',vp,vd,dvm(p)) end subroutine determine_tensor_symmetry_axis - !-------------------------------------------------------------------------------- - !================================================================================ +! +!------------------------------------------------------------------------------------------------------------- +! + ! Jacobi eigenvalue decomposition Ax = lx + subroutine jacobi_eigenvalue_decomposition(a,n,tol,lambda,vector) integer(kind=kindsi), intent(in) :: n @@ -1582,7 +1674,7 @@ function threshold(a,n) result(norm) end function threshold end subroutine jacobi_eigenvalue_decomposition - !-------------------------------------------------------------------------------- + !================================================================================ ! Check tensor decomposition diff --git a/src/inverse_problem_for_model/inverse_problem_main.f90 b/src/inverse_problem_for_model/inverse_problem_main.f90 index 8db9adaa1..a60033944 100644 --- a/src/inverse_problem_for_model/inverse_problem_main.f90 +++ b/src/inverse_problem_for_model/inverse_problem_main.f90 @@ -29,12 +29,8 @@ !################################################################################################################################## !> -!! -!! !! main subroutine to launch the iterative Full Waveform Inversion. !! -!! -!! !! -- the simulation is directly based on specfem3D git devel version !! -- the MPI use directly the communicators defined in specfem !! -- additional module are used to perform an iterative full waveform inversion @@ -76,7 +72,6 @@ subroutine inverse_problem_main() integer :: ihours,iminutes,iseconds,int_tCPU double precision, external :: wtime - !!!############################################################################################################################## !!! ---------------------------------------------- INITIALIZE RUNTIME ---------------------------------------------------------- !!!############################################################################################################################## diff --git a/src/inverse_problem_for_model/inverse_problem_par.f90 b/src/inverse_problem_for_model/inverse_problem_par.f90 index 1e3a7036f..741b99f91 100644 --- a/src/inverse_problem_for_model/inverse_problem_par.f90 +++ b/src/inverse_problem_for_model/inverse_problem_par.f90 @@ -401,18 +401,18 @@ module inverse_problem_par subroutine flush_iunit(iunit) - implicit none + implicit none - integer, intent(in) :: iunit + integer, intent(in) :: iunit - ! note: Fortran2003 includes a FLUSH statement - ! which is implemented by most compilers by now - ! - ! otherwise: - ! a) comment out the line below - ! b) try to use instead: call flush(iunit) + ! note: Fortran2003 includes a FLUSH statement + ! which is implemented by most compilers by now + ! + ! otherwise: + ! a) comment out the line below + ! b) try to use instead: call flush(iunit) - flush(iunit) + flush(iunit) end subroutine flush_iunit diff --git a/src/inverse_problem_for_model/parallel_for_inverse_problem.f90 b/src/inverse_problem_for_model/parallel_for_inverse_problem.f90 index 6418a6814..a136b0fa9 100644 --- a/src/inverse_problem_for_model/parallel_for_inverse_problem.f90 +++ b/src/inverse_problem_for_model/parallel_for_inverse_problem.f90 @@ -26,7 +26,7 @@ !===================================================================== -subroutine sum_all_all_cr_for_simulatenous_runs(sendbuf, recvbuf, countval) + subroutine sum_all_all_cr_for_simulatenous_runs(sendbuf, recvbuf, countval) use my_mpi use constants, only: CUSTOM_REAL @@ -42,19 +42,19 @@ subroutine sum_all_all_cr_for_simulatenous_runs(sendbuf, recvbuf, countval) integer :: ier if (NUMBER_OF_SIMULTANEOUS_RUNS <= 1) then - recvbuf=sendbuf - return + recvbuf = sendbuf + return endif call MPI_ALLREDUCE(sendbuf, recvbuf, countval, CUSTOM_MPI_TYPE, MPI_SUM, my_local_mpi_comm_for_bcast, ier) -end subroutine sum_all_all_cr_for_simulatenous_runs + end subroutine sum_all_all_cr_for_simulatenous_runs -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -!------------------------------------------------------------------------------------------------------------- ! !------------------------------------------------------------------------------------------------------------- -subroutine max_all_all_cr_for_simulatenous_runs(sendbuf, recvbuf, countval) +! + + subroutine max_all_all_cr_for_simulatenous_runs(sendbuf, recvbuf, countval) use my_mpi use constants, only: CUSTOM_REAL @@ -70,18 +70,19 @@ subroutine max_all_all_cr_for_simulatenous_runs(sendbuf, recvbuf, countval) integer :: ier if (NUMBER_OF_SIMULTANEOUS_RUNS <= 1) then - recvbuf=sendbuf - return + recvbuf = sendbuf + return endif call MPI_ALLREDUCE(sendbuf, recvbuf, countval, CUSTOM_MPI_TYPE, MPI_MAX, my_local_mpi_comm_for_bcast, ier) -end subroutine max_all_all_cr_for_simulatenous_runs -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -!------------------------------------------------------------------------------------------------------------- + end subroutine max_all_all_cr_for_simulatenous_runs + ! !------------------------------------------------------------------------------------------------------------- -subroutine min_all_all_cr_for_simulatenous_runs(sendbuf, recvbuf, countval) +! + + subroutine min_all_all_cr_for_simulatenous_runs(sendbuf, recvbuf, countval) use my_mpi use constants, only: CUSTOM_REAL @@ -97,19 +98,19 @@ subroutine min_all_all_cr_for_simulatenous_runs(sendbuf, recvbuf, countval) integer :: ier if (NUMBER_OF_SIMULTANEOUS_RUNS <= 1) then - recvbuf=sendbuf - return + recvbuf = sendbuf + return endif call MPI_ALLREDUCE(sendbuf, recvbuf, countval, CUSTOM_MPI_TYPE, MPI_MIN, my_local_mpi_comm_for_bcast, ier) -end subroutine min_all_all_cr_for_simulatenous_runs -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -!------------------------------------------------------------------------------------------------------------- + end subroutine min_all_all_cr_for_simulatenous_runs + ! !------------------------------------------------------------------------------------------------------------- +! -subroutine synchronize_all_world() + subroutine synchronize_all_world() use my_mpi @@ -121,12 +122,13 @@ subroutine synchronize_all_world() call MPI_BARRIER(MPI_COMM_WORLD,ier) if (ier /= 0 ) stop 'Error synchronize MPI processes' -end subroutine synchronize_all_world -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -!------------------------------------------------------------------------------------------------------------- + end subroutine synchronize_all_world + ! !------------------------------------------------------------------------------------------------------------- -subroutine synchronize_for_bcast() +! + + subroutine synchronize_for_bcast() use my_mpi @@ -138,23 +140,28 @@ subroutine synchronize_for_bcast() call MPI_BARRIER(my_local_mpi_comm_for_bcast,ier) if (ier /= 0 ) stop 'Error synchronize MPI processes' -end subroutine synchronize_for_bcast + end subroutine synchronize_for_bcast + +! +!------------------------------------------------------------------------------------------------------------- +! -subroutine dummy_bcast(dummy_integer) + subroutine dummy_bcast(dummy_integer) use my_mpi + implicit none - integer dummy_integer, ier + integer :: dummy_integer, ier + call mpi_bcast(dummy_integer, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ier) if (ier /= 0 ) stop 'Error synchronize MPI processes' - -end subroutine dummy_bcast -!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% -!------------------------------------------------------------------------------------------------------------- + end subroutine dummy_bcast ! !------------------------------------------------------------------------------------------------------------- -subroutine sum_all_all_cr_array(sendbuf, recvbuf, countval) +! + + subroutine sum_all_all_cr_array(sendbuf, recvbuf, countval) use my_mpi use constants, only: CUSTOM_REAL @@ -170,4 +177,4 @@ subroutine sum_all_all_cr_array(sendbuf, recvbuf, countval) call MPI_ALLREDUCE(sendbuf, recvbuf, countval, CUSTOM_MPI_TYPE, MPI_SUM, my_local_mpi_comm_world, ier) -end subroutine sum_all_all_cr_array + end subroutine sum_all_all_cr_array