[cig-commits] r22747 - in seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src: shared specfem3D
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Fri Aug 30 08:40:46 PDT 2013
Author: danielpeter
Date: 2013-08-30 08:40:46 -0700 (Fri, 30 Aug 2013)
New Revision: 22747
Modified:
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/get_model_parameters.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/save_header_file.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_receivers.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/save_kernels.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/save_regular_kernels.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/specfem3D_par.F90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_seismograms.f90
Log:
updates locate receivers routine; updates model ak135 modifications; removes some unused variables from save_header_file output
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/get_model_parameters.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/get_model_parameters.f90 2013-08-30 12:45:44 UTC (rev 22746)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/get_model_parameters.f90 2013-08-30 15:40:46 UTC (rev 22747)
@@ -192,7 +192,7 @@
TRANSVERSE_ISOTROPY = .true.
else if(MODEL_ROOT == '1D_iasp91' .or. MODEL_ROOT == '1D_1066a' .or. &
- MODEL_ROOT == '1D_ak135' .or. MODEL_ROOT == '1D_jp3d' .or. &
+ MODEL_ROOT == '1D_ak135f_no_mud' .or. MODEL_ROOT == '1D_jp3d' .or. &
MODEL_ROOT == '1D_sea99') then
HONOR_1D_SPHERICAL_MOHO = .true.
if(MODEL_ROOT == '1D_iasp91') then
@@ -228,7 +228,7 @@
ONE_CRUST = .true.
else if(MODEL_ROOT == '1D_iasp91_onecrust' .or. MODEL_ROOT == '1D_1066a_onecrust' &
- .or. MODEL_ROOT == '1D_ak135_onecrust') then
+ .or. MODEL_ROOT == '1D_ak135f_no_mud_onecrust') then
HONOR_1D_SPHERICAL_MOHO = .true.
ONE_CRUST = .true.
if(MODEL_ROOT == '1D_iasp91_onecrust') then
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/save_header_file.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/save_header_file.F90 2013-08-30 12:45:44 UTC (rev 22746)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/shared/save_header_file.F90 2013-08-30 15:40:46 UTC (rev 22747)
@@ -45,7 +45,15 @@
NSPEC_CRUST_MANTLE_STACEY,NSPEC_OUTER_CORE_STACEY, &
NGLOB_CRUST_MANTLE_OCEANS,NSPEC_OUTER_CORE_ROTATION )
+! daniel note: the comment below is wrong, since e.g. NSPEC_CRUST_MANTLE_ADJOINT is either 1 (dummy value)
+! for SIMULATION_TYPE == 1 or equal to NSPEC_CRUST_MANTLE for SIMULATION_TYPE == 3 or SAVE_FORWARD set to .true.
+! the value is determined in routine memory_eval() and passed here as argument.
!
+! thus, in order to be able to run a kernel simulation, users still have to re-compile the binaries with the right
+! parameters set in Par_file otherwise the kernel runs will crash.
+!
+! there is hardly another way of this usage when trying to take advantage of static compilation for the solver.
+! please ignore the following remark:
! ****************************************************************************************************
! IMPORTANT: this routine must *NOT* use flag SIMULATION_TYPE (nor SAVE_FORWARD), i.e. none of the parameters it computes
! should depend on SIMULATION_TYPE, because most users do not recompile the code nor rerun the mesher
@@ -65,9 +73,9 @@
CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH, &
NSTEP,NEX_XI,NEX_ETA, &
NPROC_XI,NPROC_ETA, &
- SAVE_REGULAR_KL,NOISE_TOMOGRAPHY, &
- APPROXIMATE_HESS_KL,ANISOTROPIC_KL,PARTIAL_PHYS_DISPERSION_ONLY, &
- SAVE_SOURCE_MASK,ABSORBING_CONDITIONS,USE_LDDRK,EXACT_MASS_MATRIX_FOR_ROTATION, &
+ SAVE_REGULAR_KL, &
+ PARTIAL_PHYS_DISPERSION_ONLY, &
+ ABSORBING_CONDITIONS,USE_LDDRK,EXACT_MASS_MATRIX_FOR_ROTATION, &
ATTENUATION_1D_WITH_3D_STORAGE, &
ATT1,ATT2,ATT3,ATT4,ATT5, &
MOVIE_VOLUME
@@ -345,23 +353,26 @@
write(IOUT,*) 'integer, parameter :: NSPEC_OUTER_CORE_ADJOINT = ',NSPEC_OUTER_CORE_ADJOINT
write(IOUT,*) 'integer, parameter :: NSPEC_INNER_CORE_ADJOINT = ',NSPEC_INNER_CORE_ADJOINT
- if(ANISOTROPIC_KL) then
- write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_ADJOINT_ANISO_KL = ',NSPEC_CRUST_MANTLE_ADJOINT
- else
- write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_ADJOINT_ANISO_KL = ',1
- endif
+ ! unused... (dynamic allocation used)
+ !if(ANISOTROPIC_KL) then
+ ! write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_ADJOINT_ANISO_KL = ',NSPEC_CRUST_MANTLE_ADJOINT
+ !else
+ ! write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_ADJOINT_ANISO_KL = ',1
+ !endif
- if(APPROXIMATE_HESS_KL) then
- write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_ADJOINT_HESS = ',NSPEC_CRUST_MANTLE_ADJOINT
- else
- write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_ADJOINT_HESS = ',1
- endif
+ ! unused... (dynamic allocation used)
+ !if(APPROXIMATE_HESS_KL) then
+ ! write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_ADJOINT_HESS = ',NSPEC_CRUST_MANTLE_ADJOINT
+ !else
+ ! write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_ADJOINT_HESS = ',1
+ !endif
- if(NOISE_TOMOGRAPHY > 0) then
- write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_ADJOINT_NOISE = ',NSPEC_CRUST_MANTLE_ADJOINT
- else
- write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_ADJOINT_NOISE = ',1
- endif
+ ! unused... (dynamic allocation used)
+ !if(NOISE_TOMOGRAPHY > 0) then
+ ! write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_ADJOINT_NOISE = ',NSPEC_CRUST_MANTLE_ADJOINT
+ !else
+ ! write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_ADJOINT_NOISE = ',1
+ !endif
write(IOUT,*) 'integer, parameter :: NGLOB_CRUST_MANTLE_ADJOINT = ',NGLOB_CRUST_MANTLE_ADJOINT
write(IOUT,*) 'integer, parameter :: NGLOB_OUTER_CORE_ADJOINT = ',NGLOB_OUTER_CORE_ADJOINT
@@ -532,12 +543,13 @@
endif
write(IOUT,*)
- if (SAVE_SOURCE_MASK) then
- write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_MASK_SOURCE = NSPEC_CRUST_MANTLE'
- else
- write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_MASK_SOURCE = 1'
- endif
- write(IOUT,*)
+ ! unused... (dynamic allocation used)
+ !if (SAVE_SOURCE_MASK) then
+ ! write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_MASK_SOURCE = NSPEC_CRUST_MANTLE'
+ !else
+ ! write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_MASK_SOURCE = 1'
+ !endif
+ !write(IOUT,*)
! in the case of Stacey boundary conditions, add C*delta/2 contribution to the mass matrix
! on the Stacey edges for the crust_mantle and outer_core regions but not for the inner_core region
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_receivers.f90 2013-08-30 12:45:44 UTC (rev 22746)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/locate_receivers.f90 2013-08-30 15:40:46 UTC (rev 22747)
@@ -110,12 +110,12 @@
double precision, dimension(:), allocatable :: final_distance
double precision, dimension(:,:), allocatable :: final_distance_all
- double precision distmin,final_distance_max
+ double precision :: distmin,final_distance_max
! receiver information
! timing information for the stations
! station information for writing the seismograms
- integer nsamp
+ integer :: nsamp
integer, dimension(nrec) :: islice_selected_rec_found,ispec_selected_rec_found
double precision, dimension(nrec) :: xi_receiver_found,eta_receiver_found,gamma_receiver_found
@@ -123,19 +123,25 @@
character(len=MAX_LENGTH_STATION_NAME), dimension(nrec) :: station_name_found
character(len=MAX_LENGTH_NETWORK_NAME), dimension(nrec) :: network_name_found
double precision, dimension(nrec) :: stlat_found,stlon_found,stele_found,stbur_found,epidist_found
- character(len=150) STATIONS
+ character(len=150) :: STATIONS
integer, allocatable, dimension(:,:) :: ispec_selected_rec_all
double precision, allocatable, dimension(:,:) :: xi_receiver_all,eta_receiver_all,gamma_receiver_all
+ double precision, allocatable, dimension(:) :: xi_receiver_subset,eta_receiver_subset,gamma_receiver_subset
- double precision x_target_rec,y_target_rec,z_target_rec
+ double precision :: x_target_rec,y_target_rec,z_target_rec
- double precision typical_size
- logical located_target
+ double precision :: typical_size
+ logical :: located_target
- character(len=150) OUTPUT_FILES
- character(len=2) bic
+ character(len=150) :: OUTPUT_FILES
+ character(len=2) :: bic
+ integer :: nrec_SUBSET_current_size,irec_in_this_subset,irec_already_done
+ double precision, allocatable, dimension(:) :: x_found_subset,y_found_subset,z_found_subset
+ double precision, allocatable, dimension(:) :: final_distance_subset
+ integer, allocatable, dimension(:) :: ispec_selected_rec_subset
+
integer, allocatable, dimension(:) :: station_duplet
! timing
@@ -186,19 +192,11 @@
x_found(nrec), &
y_found(nrec), &
z_found(nrec), &
- final_distance(nrec), &
- ispec_selected_rec_all(nrec,0:NPROCTOT_VAL-1), &
- xi_receiver_all(nrec,0:NPROCTOT_VAL-1), &
- eta_receiver_all(nrec,0:NPROCTOT_VAL-1), &
- gamma_receiver_all(nrec,0:NPROCTOT_VAL-1), &
- x_found_all(nrec,0:NPROCTOT_VAL-1), &
- y_found_all(nrec,0:NPROCTOT_VAL-1), &
- z_found_all(nrec,0:NPROCTOT_VAL-1), &
- final_distance_all(nrec,0:NPROCTOT_VAL-1),stat=ier)
+ final_distance(nrec),stat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'error allocating temporary receiver arrays')
+
! initializes
final_distance(:) = HUGEVAL
- final_distance_all(:,:) = HUGEVAL
! read that STATIONS file on the master
if(myrank == 0) then
@@ -468,173 +466,265 @@
! find the best (xi,eta) for each receiver
! ****************************************
- ! loop on all the receivers to iterate in that slice
- do irec = 1,nrec
+ ! loop on all the receivers
+ ! gather source information in subsets to reduce memory requirements
+ islice_selected_rec(:) = -1
- ispec_iterate = ispec_selected_rec(irec)
+ ! loop over subsets of receivers
+ do irec_already_done = 0, nrec, nrec_SUBSET_MAX
- ! define coordinates of the control points of the element
- do ia=1,NGNOD
+ ! the size of the subset can be the maximum size, or less (if we are in the last subset,
+ ! or if there are fewer sources than the maximum size of a subset)
+ nrec_SUBSET_current_size = min(nrec_SUBSET_MAX, nrec - irec_already_done)
- iax = 0
- if(iaddx(ia) == 0) then
- iax = 1
- else if(iaddx(ia) == 1) then
- iax = MIDX
- else if(iaddx(ia) == 2) then
- iax = NGLLX
- else
- call exit_MPI(myrank,'incorrect value of iaddx')
- endif
+ ! allocate arrays specific to each subset
+ allocate(ispec_selected_rec_subset(nrec_SUBSET_current_size), &
+ xi_receiver_subset(nrec_SUBSET_current_size), &
+ eta_receiver_subset(nrec_SUBSET_current_size), &
+ gamma_receiver_subset(nrec_SUBSET_current_size), &
+ x_found_subset(nrec_SUBSET_current_size), &
+ y_found_subset(nrec_SUBSET_current_size), &
+ z_found_subset(nrec_SUBSET_current_size), &
+ final_distance_subset(nrec_SUBSET_current_size),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating temporary receiver arrays')
- iay = 0
- if(iaddy(ia) == 0) then
- iay = 1
- else if(iaddy(ia) == 1) then
- iay = MIDY
- else if(iaddy(ia) == 2) then
- iay = NGLLY
- else
- call exit_MPI(myrank,'incorrect value of iaddy')
- endif
+ ! gather arrays
+ if( myrank == 0 ) then
+ ! only master process needs full arrays allocated
+ allocate(ispec_selected_rec_all(nrec_SUBSET_current_size,0:NPROCTOT_VAL-1), &
+ xi_receiver_all(nrec_SUBSET_current_size,0:NPROCTOT_VAL-1), &
+ eta_receiver_all(nrec_SUBSET_current_size,0:NPROCTOT_VAL-1), &
+ gamma_receiver_all(nrec_SUBSET_current_size,0:NPROCTOT_VAL-1), &
+ x_found_all(nrec_SUBSET_current_size,0:NPROCTOT_VAL-1), &
+ y_found_all(nrec_SUBSET_current_size,0:NPROCTOT_VAL-1), &
+ z_found_all(nrec_SUBSET_current_size,0:NPROCTOT_VAL-1), &
+ final_distance_all(nrec_SUBSET_current_size,0:NPROCTOT_VAL-1),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating temporary gather receiver arrays')
+ else
+ ! dummy arrays
+ allocate(ispec_selected_rec_all(1,1), &
+ xi_receiver_all(1,1), &
+ eta_receiver_all(1,1), &
+ gamma_receiver_all(1,1), &
+ x_found_all(1,1), &
+ y_found_all(1,1), &
+ z_found_all(1,1), &
+ final_distance_all(1,1),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating temporary dummy receiver arrays')
+ endif
+ final_distance_all(:,:) = HUGEVAL
- iaz = 0
- if(iaddr(ia) == 0) then
- iaz = 1
- else if(iaddr(ia) == 1) then
- iaz = MIDZ
- else if(iaddr(ia) == 2) then
- iaz = NGLLZ
- else
- call exit_MPI(myrank,'incorrect value of iaddr')
- endif
+ ! loop over the stations within this subset
+ do irec_in_this_subset = 1,nrec_SUBSET_current_size
- iglob = ibool(iax,iay,iaz,ispec_iterate)
- xelm(ia) = dble(xstore(iglob))
- yelm(ia) = dble(ystore(iglob))
- zelm(ia) = dble(zstore(iglob))
+ ! mapping from station number in current subset to real station number in all the subsets
+ irec = irec_in_this_subset + irec_already_done
- enddo
+ ispec_selected_rec_subset(irec_in_this_subset) = ispec_selected_rec(irec)
+ ispec_iterate = ispec_selected_rec(irec)
- ! use initial guess in xi and eta
- xi = xigll(ix_initial_guess(irec))
- eta = yigll(iy_initial_guess(irec))
- gamma = zigll(iz_initial_guess(irec))
+ ! define coordinates of the control points of the element
+ do ia=1,NGNOD
- x_target_rec = x_target(irec)
- y_target_rec = y_target(irec)
- z_target_rec = z_target(irec)
+ iax = 0
+ if(iaddx(ia) == 0) then
+ iax = 1
+ else if(iaddx(ia) == 1) then
+ iax = MIDX
+ else if(iaddx(ia) == 2) then
+ iax = NGLLX
+ else
+ call exit_MPI(myrank,'incorrect value of iaddx')
+ endif
- ! iterate to solve the non linear system
- do iter_loop = 1,NUM_ITER
+ iay = 0
+ if(iaddy(ia) == 0) then
+ iay = 1
+ else if(iaddy(ia) == 1) then
+ iay = MIDY
+ else if(iaddy(ia) == 2) then
+ iay = NGLLY
+ else
+ call exit_MPI(myrank,'incorrect value of iaddy')
+ endif
- ! impose receiver exactly at the surface
+ iaz = 0
+ if(iaddr(ia) == 0) then
+ iaz = 1
+ else if(iaddr(ia) == 1) then
+ iaz = MIDZ
+ else if(iaddr(ia) == 2) then
+ iaz = NGLLZ
+ else
+ call exit_MPI(myrank,'incorrect value of iaddr')
+ endif
+
+ iglob = ibool(iax,iay,iaz,ispec_iterate)
+ xelm(ia) = dble(xstore(iglob))
+ yelm(ia) = dble(ystore(iglob))
+ zelm(ia) = dble(zstore(iglob))
+
+ enddo
+
+ ! use initial guess in xi and eta
+ xi = xigll(ix_initial_guess(irec))
+ eta = yigll(iy_initial_guess(irec))
+ gamma = zigll(iz_initial_guess(irec))
+
+ x_target_rec = x_target(irec)
+ y_target_rec = y_target(irec)
+ z_target_rec = z_target(irec)
+
+ ! iterate to solve the non linear system
+ do iter_loop = 1,NUM_ITER
+
+ ! impose receiver exactly at the surface
+ if(.not. RECEIVERS_CAN_BE_BURIED) gamma = 1.d0
+
+ ! recompute jacobian for the new point
+ call recompute_jacobian(xelm,yelm,zelm,xi,eta,gamma,x,y,z, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
+
+ ! compute distance to target location
+ dx = - (x - x_target_rec)
+ dy = - (y - y_target_rec)
+ dz = - (z - z_target_rec)
+
+ ! compute increments
+ ! gamma does not change since we know the receiver is exactly on the surface
+ dxi = xix*dx + xiy*dy + xiz*dz
+ deta = etax*dx + etay*dy + etaz*dz
+
+ ! update values
+ xi = xi + dxi
+ eta = eta + deta
+
+ ! buried receivers vary in z depth
+ if(RECEIVERS_CAN_BE_BURIED) then
+ dgamma = gammax*dx + gammay*dy + gammaz*dz
+ gamma = gamma + dgamma
+ endif
+
+ ! impose that we stay in that element
+ ! (useful if user gives a receiver outside the mesh for instance)
+ ! we can go slightly outside the [1,1] segment since with finite elements
+ ! the polynomial solution is defined everywhere
+ ! can be useful for convergence of iterative scheme with distorted elements
+ if (xi > 1.10d0) xi = 1.10d0
+ if (xi < -1.10d0) xi = -1.10d0
+ if (eta > 1.10d0) eta = 1.10d0
+ if (eta < -1.10d0) eta = -1.10d0
+ if (gamma > 1.10d0) gamma = 1.10d0
+ if (gamma < -1.10d0) gamma = -1.10d0
+
+ ! end of non linear iterations
+ enddo
+
+ ! impose receiver exactly at the surface after final iteration
if(.not. RECEIVERS_CAN_BE_BURIED) gamma = 1.d0
- ! recompute jacobian for the new point
+ ! compute final coordinates of point found
call recompute_jacobian(xelm,yelm,zelm,xi,eta,gamma,x,y,z, &
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
- ! compute distance to target location
- dx = - (x - x_target_rec)
- dy = - (y - y_target_rec)
- dz = - (z - z_target_rec)
+ ! store xi,eta and x,y,z of point found
+ xi_receiver_subset(irec_in_this_subset) = xi
+ eta_receiver_subset(irec_in_this_subset) = eta
+ gamma_receiver_subset(irec_in_this_subset) = gamma
- ! compute increments
- ! gamma does not change since we know the receiver is exactly on the surface
- dxi = xix*dx + xiy*dy + xiz*dz
- deta = etax*dx + etay*dy + etaz*dz
+ x_found(irec) = x
+ y_found(irec) = y
+ z_found(irec) = z
- ! update values
- xi = xi + dxi
- eta = eta + deta
+ x_found_subset(irec_in_this_subset) = x_found(irec)
+ y_found_subset(irec_in_this_subset) = y_found(irec)
+ z_found_subset(irec_in_this_subset) = z_found(irec)
- ! buried receivers vary in z depth
- if(RECEIVERS_CAN_BE_BURIED) then
- dgamma = gammax*dx + gammay*dy + gammaz*dz
- gamma = gamma + dgamma
- endif
+ ! compute final distance between asked and found (converted to km)
+ final_distance(irec) = dsqrt((x_target_rec-x)**2 + &
+ (y_target_rec-y)**2 + &
+ (z_target_rec-z)**2)*R_EARTH/1000.d0
- ! impose that we stay in that element
- ! (useful if user gives a receiver outside the mesh for instance)
- ! we can go slightly outside the [1,1] segment since with finite elements
- ! the polynomial solution is defined everywhere
- ! can be useful for convergence of iterative scheme with distorted elements
- if (xi > 1.10d0) xi = 1.10d0
- if (xi < -1.10d0) xi = -1.10d0
- if (eta > 1.10d0) eta = 1.10d0
- if (eta < -1.10d0) eta = -1.10d0
- if (gamma > 1.10d0) gamma = 1.10d0
- if (gamma < -1.10d0) gamma = -1.10d0
+ final_distance_subset(irec_in_this_subset) = final_distance(irec)
- ! end of non linear iterations
- enddo
+ enddo ! end of loop on all stations within current subset
- ! impose receiver exactly at the surface after final iteration
- if(.not. RECEIVERS_CAN_BE_BURIED) gamma = 1.d0
+ ! for MPI version, gather information from all the nodes
+ ispec_selected_rec_all(:,:) = -1
- ! compute final coordinates of point found
- call recompute_jacobian(xelm,yelm,zelm,xi,eta,gamma,x,y,z, &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
+ call gather_all_i(ispec_selected_rec_subset,nrec_SUBSET_current_size, &
+ ispec_selected_rec_all,nrec_SUBSET_current_size,NPROCTOT_VAL)
- ! store xi,eta and x,y,z of point found
- xi_receiver(irec) = xi
- eta_receiver(irec) = eta
- gamma_receiver(irec) = gamma
- x_found(irec) = x
- y_found(irec) = y
- z_found(irec) = z
+ ! this is executed by main process only
+ if(myrank == 0) then
+ ! check that the gather operation went well
+ if(any(ispec_selected_rec_all(:,:) == -1)) call exit_MPI(myrank,'gather operation failed for receivers')
+ endif
- ! compute final distance between asked and found (converted to km)
- final_distance(irec) = dsqrt((x_target_rec-x)**2 + &
- (y_target_rec-y)**2 + &
- (z_target_rec-z)**2)*R_EARTH/1000.d0
+ call gather_all_dp(xi_receiver_subset,nrec_SUBSET_current_size,xi_receiver_all,nrec_SUBSET_current_size,NPROCTOT_VAL)
+ call gather_all_dp(eta_receiver_subset,nrec_SUBSET_current_size,eta_receiver_all,nrec_SUBSET_current_size,NPROCTOT_VAL)
+ call gather_all_dp(gamma_receiver_subset,nrec_SUBSET_current_size,gamma_receiver_all,nrec_SUBSET_current_size,NPROCTOT_VAL)
+ call gather_all_dp(final_distance_subset,nrec_SUBSET_current_size,final_distance_all,nrec_SUBSET_current_size,NPROCTOT_VAL)
+ call gather_all_dp(x_found_subset,nrec_SUBSET_current_size,x_found_all,nrec_SUBSET_current_size,NPROCTOT_VAL)
+ call gather_all_dp(y_found_subset,nrec_SUBSET_current_size,y_found_all,nrec_SUBSET_current_size,NPROCTOT_VAL)
+ call gather_all_dp(z_found_subset,nrec_SUBSET_current_size,z_found_all,nrec_SUBSET_current_size,NPROCTOT_VAL)
- enddo
+ ! this is executed by main process only
+ if(myrank == 0) then
- ! for MPI version, gather information from all the nodes
- ispec_selected_rec_all(:,:) = -1
- call gather_all_i(ispec_selected_rec,nrec,ispec_selected_rec_all,nrec,NPROCTOT_VAL)
+ ! MPI loop on all the results to determine the best slice
+ islice_selected_rec(:) = -1
+ do irec_in_this_subset = 1,nrec_SUBSET_current_size
- call gather_all_dp(xi_receiver,nrec,xi_receiver_all,nrec,NPROCTOT_VAL)
- call gather_all_dp(eta_receiver,nrec,eta_receiver_all,nrec,NPROCTOT_VAL)
- call gather_all_dp(gamma_receiver,nrec,gamma_receiver_all,nrec,NPROCTOT_VAL)
- call gather_all_dp(final_distance,nrec,final_distance_all,nrec,NPROCTOT_VAL)
- call gather_all_dp(x_found,nrec,x_found_all,nrec,NPROCTOT_VAL)
- call gather_all_dp(y_found,nrec,y_found_all,nrec,NPROCTOT_VAL)
- call gather_all_dp(z_found,nrec,z_found_all,nrec,NPROCTOT_VAL)
+ ! mapping from station number in current subset to real station number in all the subsets
+ irec = irec_in_this_subset + irec_already_done
- ! this is executed by main process only
- if(myrank == 0) then
+ distmin = HUGEVAL
+ do iprocloop = 0,NPROCTOT_VAL-1
+ if(final_distance_all(irec_in_this_subset,iprocloop) < distmin) then
+ distmin = final_distance_all(irec_in_this_subset,iprocloop)
+ islice_selected_rec(irec) = iprocloop
+ ispec_selected_rec(irec) = ispec_selected_rec_all(irec_in_this_subset,iprocloop)
+ xi_receiver(irec) = xi_receiver_all(irec_in_this_subset,iprocloop)
+ eta_receiver(irec) = eta_receiver_all(irec_in_this_subset,iprocloop)
+ gamma_receiver(irec) = gamma_receiver_all(irec_in_this_subset,iprocloop)
+ x_found(irec) = x_found_all(irec_in_this_subset,iprocloop)
+ y_found(irec) = y_found_all(irec_in_this_subset,iprocloop)
+ z_found(irec) = z_found_all(irec_in_this_subset,iprocloop)
+ endif
+ enddo
+ final_distance(irec) = distmin
+ enddo
+ endif ! end of section executed by main process only
- ! check that the gather operation went well
- if(any(ispec_selected_rec_all(:,:) == -1)) call exit_MPI(myrank,'gather operation failed for receivers')
+ deallocate(ispec_selected_rec_subset)
+ deallocate(ispec_selected_rec_all)
+ deallocate(xi_receiver_subset)
+ deallocate(eta_receiver_subset)
+ deallocate(gamma_receiver_subset)
+ deallocate(xi_receiver_all)
+ deallocate(eta_receiver_all)
+ deallocate(gamma_receiver_all)
+ deallocate(x_found_subset)
+ deallocate(y_found_subset)
+ deallocate(z_found_subset)
+ deallocate(x_found_all)
+ deallocate(y_found_all)
+ deallocate(z_found_all)
+ deallocate(final_distance_subset)
+ deallocate(final_distance_all)
- ! MPI loop on all the results to determine the best slice
- islice_selected_rec(:) = -1
- do irec = 1,nrec
- distmin = HUGEVAL
- do iprocloop = 0,NPROCTOT_VAL-1
- if(final_distance_all(irec,iprocloop) < distmin) then
- distmin = final_distance_all(irec,iprocloop)
- islice_selected_rec(irec) = iprocloop
- ispec_selected_rec(irec) = ispec_selected_rec_all(irec,iprocloop)
- xi_receiver(irec) = xi_receiver_all(irec,iprocloop)
- eta_receiver(irec) = eta_receiver_all(irec,iprocloop)
- gamma_receiver(irec) = gamma_receiver_all(irec,iprocloop)
- x_found(irec) = x_found_all(irec,iprocloop)
- y_found(irec) = y_found_all(irec,iprocloop)
- z_found(irec) = z_found_all(irec,iprocloop)
- endif
- enddo
- final_distance(irec) = distmin
- enddo
+ enddo ! end of loop over all station subsets
+ ! this is executed by the main process only
+ if(myrank == 0) then
+
! appends receiver locations to sr.vtk file
open(IOVTK,file=trim(OUTPUT_FILES)//'/sr_tmp.vtk', &
position='append',status='old',iostat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'Error opening and appending receivers to file sr_tmp.vtk')
+ ! chooses best receivers locations
islice_selected_rec_found(:) = -1
nrec_found = 0
do irec=1,nrec
@@ -762,10 +852,6 @@
deallocate(x_target,y_target,z_target)
deallocate(x_found,y_found,z_found)
deallocate(final_distance)
- deallocate(ispec_selected_rec_all)
- deallocate(xi_receiver_all,eta_receiver_all,gamma_receiver_all)
- deallocate(x_found_all,y_found_all,z_found_all)
- deallocate(final_distance_all)
! main process broadcasts the results to all the slices
call bcast_all_i(nrec,1)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90 2013-08-30 12:45:44 UTC (rev 22746)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/prepare_timerun.f90 2013-08-30 15:40:46 UTC (rev 22747)
@@ -1055,8 +1055,14 @@
rho_kl_crust_mantle(:,:,:,:) = 0._CUSTOM_REAL
beta_kl_crust_mantle(:,:,:,:) = 0._CUSTOM_REAL
alpha_kl_crust_mantle(:,:,:,:) = 0._CUSTOM_REAL
- if (NOISE_TOMOGRAPHY == 3) Sigma_kl_crust_mantle(:,:,:,:) = 0._CUSTOM_REAL
+ ! noise strength kernel
+ if (NOISE_TOMOGRAPHY == 3) then
+ allocate( Sigma_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating noise sigma kernel')
+ Sigma_kl_crust_mantle(:,:,:,:) = 0._CUSTOM_REAL
+ endif
+
! approximate hessian
if( APPROXIMATE_HESS_KL ) then
allocate( hess_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT),stat=ier)
@@ -1065,7 +1071,11 @@
endif
! For anisotropic kernels (in crust_mantle only)
- cijkl_kl_crust_mantle(:,:,:,:,:) = 0._CUSTOM_REAL
+ if( ANISOTROPIC_KL ) then
+ allocate( cijkl_kl_crust_mantle(21,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT),stat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error allocating full cijkl kernel in crust_mantle')
+ cijkl_kl_crust_mantle(:,:,:,:,:) = 0._CUSTOM_REAL
+ endif
rho_kl_outer_core(:,:,:,:) = 0._CUSTOM_REAL
alpha_kl_outer_core(:,:,:,:) = 0._CUSTOM_REAL
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/save_kernels.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/save_kernels.f90 2013-08-30 12:45:44 UTC (rev 22746)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/save_kernels.f90 2013-08-30 15:40:46 UTC (rev 22747)
@@ -322,96 +322,96 @@
if( SAVE_TRANSVERSE_KL_ONLY ) then
! transverse isotropic kernels
! (alpha_v, alpha_h, beta_v, beta_h, eta, rho ) parameterization
- open(unit=27,file=trim(prname)//'alphav_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) alphav_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'alphah_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) alphah_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'betav_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) betav_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'betah_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) betah_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'eta_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) eta_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'rho_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) rho_kl_crust_mantle
- close(27)
+ open(unit=IOUT,file=trim(prname)//'alphav_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) alphav_kl_crust_mantle
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'alphah_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) alphah_kl_crust_mantle
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'betav_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) betav_kl_crust_mantle
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'betah_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) betah_kl_crust_mantle
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'eta_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) eta_kl_crust_mantle
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'rho_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) rho_kl_crust_mantle
+ close(IOUT)
! in case one is interested in primary kernel K_rho
- !open(unit=27,file=trim(prname)//'rhonotprime_kernel.bin',status='unknown',form='unformatted',action='write')
- !write(27) rhonotprime_kl_crust_mantle
- !close(27)
+ !open(unit=IOUT,file=trim(prname)//'rhonotprime_kernel.bin',status='unknown',form='unformatted',action='write')
+ !write(IOUT) rhonotprime_kl_crust_mantle
+ !close(IOUT)
! (bulk, beta_v, beta_h, eta, rho ) parameterization: K_eta and K_rho same as above
- open(unit=27,file=trim(prname)//'bulk_c_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) bulk_c_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'bulk_betav_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) bulk_betav_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'bulk_betah_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) bulk_betah_kl_crust_mantle
- close(27)
+ open(unit=IOUT,file=trim(prname)//'bulk_c_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) bulk_c_kl_crust_mantle
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'bulk_betav_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) bulk_betav_kl_crust_mantle
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'bulk_betah_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) bulk_betah_kl_crust_mantle
+ close(IOUT)
! to check: isotropic kernels
- open(unit=27,file=trim(prname)//'alpha_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) alpha_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'beta_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) beta_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'bulk_beta_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) bulk_beta_kl_crust_mantle
- close(27)
+ open(unit=IOUT,file=trim(prname)//'alpha_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) alpha_kl_crust_mantle
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'beta_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) beta_kl_crust_mantle
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'bulk_beta_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) bulk_beta_kl_crust_mantle
+ close(IOUT)
else
! fully anisotropic kernels
! note: the C_ij and density kernels are not for relative perturbations (delta ln( m_i) = delta m_i / m_i),
! but absolute perturbations (delta m_i = m_i - m_0)
- open(unit=27,file=trim(prname)//'rho_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) - rho_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'cijkl_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) - cijkl_kl_crust_mantle
- close(27)
+ open(unit=IOUT,file=trim(prname)//'rho_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) - rho_kl_crust_mantle
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'cijkl_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) - cijkl_kl_crust_mantle
+ close(IOUT)
endif
else
! primary kernels: (rho,kappa,mu) parameterization
- open(unit=27,file=trim(prname)//'rhonotprime_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) rhonotprime_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'kappa_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) kappa_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'mu_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) mu_kl_crust_mantle
- close(27)
+ open(unit=IOUT,file=trim(prname)//'rhonotprime_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) rhonotprime_kl_crust_mantle
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'kappa_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) kappa_kl_crust_mantle
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'mu_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) mu_kl_crust_mantle
+ close(IOUT)
! (rho, alpha, beta ) parameterization
- open(unit=27,file=trim(prname)//'rho_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) rho_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'alpha_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) alpha_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'beta_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) beta_kl_crust_mantle
- close(27)
+ open(unit=IOUT,file=trim(prname)//'rho_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) rho_kl_crust_mantle
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'alpha_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) alpha_kl_crust_mantle
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'beta_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) beta_kl_crust_mantle
+ close(IOUT)
! (rho, bulk, beta ) parameterization, K_rho same as above
- open(unit=27,file=trim(prname)//'bulk_c_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) bulk_c_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'bulk_beta_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) bulk_beta_kl_crust_mantle
- close(27)
+ open(unit=IOUT,file=trim(prname)//'bulk_c_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) bulk_c_kl_crust_mantle
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'bulk_beta_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) bulk_beta_kl_crust_mantle
+ close(IOUT)
endif
@@ -467,12 +467,12 @@
call create_name_database(prname,myrank,IREGION_OUTER_CORE,LOCAL_TMP_PATH)
- open(unit=27,file=trim(prname)//'rho_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) rho_kl_outer_core
- close(27)
- open(unit=27,file=trim(prname)//'alpha_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) alpha_kl_outer_core
- close(27)
+ open(unit=IOUT,file=trim(prname)//'rho_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) rho_kl_outer_core
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'alpha_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) alpha_kl_outer_core
+ close(IOUT)
end subroutine save_kernels_outer_core
@@ -517,15 +517,15 @@
call create_name_database(prname,myrank,IREGION_INNER_CORE,LOCAL_TMP_PATH)
- open(unit=27,file=trim(prname)//'rho_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) rho_kl_inner_core
- close(27)
- open(unit=27,file=trim(prname)//'alpha_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) alpha_kl_inner_core
- close(27)
- open(unit=27,file=trim(prname)//'beta_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) beta_kl_inner_core
- close(27)
+ open(unit=IOUT,file=trim(prname)//'rho_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) rho_kl_inner_core
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'alpha_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) alpha_kl_inner_core
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'beta_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) beta_kl_inner_core
+ close(IOUT)
end subroutine save_kernels_inner_core
@@ -557,28 +557,28 @@
call create_name_database(prname,myrank,IREGION_CRUST_MANTLE,LOCAL_TMP_PATH)
if (.not. SUPPRESS_CRUSTAL_MESH .and. HONOR_1D_SPHERICAL_MOHO) then
- open(unit=27,file=trim(prname)//'moho_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) moho_kl
- close(27)
+ open(unit=IOUT,file=trim(prname)//'moho_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) moho_kl
+ close(IOUT)
endif
- open(unit=27,file=trim(prname)//'d400_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) d400_kl
- close(27)
+ open(unit=IOUT,file=trim(prname)//'d400_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) d400_kl
+ close(IOUT)
- open(unit=27,file=trim(prname)//'d670_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) d670_kl
- close(27)
+ open(unit=IOUT,file=trim(prname)//'d670_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) d670_kl
+ close(IOUT)
- open(unit=27,file=trim(prname)//'CMB_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) cmb_kl
- close(27)
+ open(unit=IOUT,file=trim(prname)//'CMB_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) cmb_kl
+ close(IOUT)
call create_name_database(prname,myrank,IREGION_OUTER_CORE,LOCAL_PATH)
- open(unit=27,file=trim(prname)//'ICB_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) icb_kl
- close(27)
+ open(unit=IOUT,file=trim(prname)//'ICB_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) icb_kl
+ close(IOUT)
end subroutine save_kernels_boundary_kl
@@ -613,7 +613,7 @@
shdur_der(irec_local) = shdur_der(irec_local) * scale_displ**2
write(outputname,'(a,i5.5)') 'OUTPUT_FILES/src_frechet.',number_receiver_global(irec_local)
- open(unit=27,file=trim(outputname),status='unknown',action='write')
+ open(unit=IOUT,file=trim(outputname),status='unknown',action='write')
!
! r -> z, theta -> -n, phi -> e, plus factor 2 for Mrt,Mrp,Mtp, and 1e-7 to dyne.cm
! Mrr = Mzz
@@ -625,18 +625,18 @@
! for consistency, location derivatives are in the order of [Xr,Xt,Xp]
! minus sign for sloc_der(3,irec_local) to get derivative for depth instead of radius
- write(27,'(g16.5)') moment_der(3,3,irec_local) * 1e-7
- write(27,'(g16.5)') moment_der(1,1,irec_local) * 1e-7
- write(27,'(g16.5)') moment_der(2,2,irec_local) * 1e-7
- write(27,'(g16.5)') -2*moment_der(1,3,irec_local) * 1e-7
- write(27,'(g16.5)') 2*moment_der(2,3,irec_local) * 1e-7
- write(27,'(g16.5)') -2*moment_der(1,2,irec_local) * 1e-7
- write(27,'(g16.5)') sloc_der(2,irec_local)
- write(27,'(g16.5)') sloc_der(1,irec_local)
- write(27,'(g16.5)') -sloc_der(3,irec_local)
- write(27,'(g16.5)') stshift_der(irec_local)
- write(27,'(g16.5)') shdur_der(irec_local)
- close(27)
+ write(IOUT,'(g16.5)') moment_der(3,3,irec_local) * 1e-7
+ write(IOUT,'(g16.5)') moment_der(1,1,irec_local) * 1e-7
+ write(IOUT,'(g16.5)') moment_der(2,2,irec_local) * 1e-7
+ write(IOUT,'(g16.5)') -2*moment_der(1,3,irec_local) * 1e-7
+ write(IOUT,'(g16.5)') 2*moment_der(2,3,irec_local) * 1e-7
+ write(IOUT,'(g16.5)') -2*moment_der(1,2,irec_local) * 1e-7
+ write(IOUT,'(g16.5)') sloc_der(2,irec_local)
+ write(IOUT,'(g16.5)') sloc_der(1,irec_local)
+ write(IOUT,'(g16.5)') -sloc_der(3,irec_local)
+ write(IOUT,'(g16.5)') stshift_der(irec_local)
+ write(IOUT,'(g16.5)') shdur_der(irec_local)
+ close(IOUT)
enddo
end subroutine save_kernels_source_derivatives
@@ -664,9 +664,9 @@
! stores into file
call create_name_database(prname,myrank,IREGION_CRUST_MANTLE,LOCAL_TMP_PATH)
- open(unit=27,file=trim(prname)//'hess_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) hess_kl_crust_mantle
- close(27)
+ open(unit=IOUT,file=trim(prname)//'hess_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) hess_kl_crust_mantle
+ close(IOUT)
end subroutine save_kernels_hessian
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/save_regular_kernels.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/save_regular_kernels.f90 2013-08-30 12:45:44 UTC (rev 22746)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/save_regular_kernels.f90 2013-08-30 15:40:46 UTC (rev 22747)
@@ -394,96 +394,96 @@
if (SAVE_TRANSVERSE_KL_ONLY) then
! transverse isotropic kernels
! (alpha_v, alpha_h, beta_v, beta_h, eta, rho ) parameterization
- open(unit=27,file=trim(prname)//'alphav_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) alphav_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'alphah_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) alphah_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'betav_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) betav_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'betah_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) betah_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'eta_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) eta_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'rho_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) rho_kl_crust_mantle_reg
- close(27)
+ open(unit=IOUT,file=trim(prname)//'alphav_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) alphav_kl_crust_mantle
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'alphah_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) alphah_kl_crust_mantle
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'betav_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) betav_kl_crust_mantle
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'betah_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) betah_kl_crust_mantle
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'eta_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) eta_kl_crust_mantle
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'rho_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) rho_kl_crust_mantle_reg
+ close(IOUT)
! in case one is interested in primary kernel K_rho
- !open(unit=27,file=trim(prname)//'rhonotprime_kernel.bin',status='unknown',form='unformatted',action='write')
- !write(27) rhonotprime_kl_crust_mantle
- !close(27)
+ !open(unit=IOUT,file=trim(prname)//'rhonotprime_kernel.bin',status='unknown',form='unformatted',action='write')
+ !write(IOUT) rhonotprime_kl_crust_mantle
+ !close(IOUT)
! (bulk, beta_v, beta_h, eta, rho ) parameterization: K_eta and K_rho same as above
- open(unit=27,file=trim(prname)//'bulk_c_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) bulk_c_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'bulk_betav_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) bulk_betav_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'bulk_betah_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) bulk_betah_kl_crust_mantle
- close(27)
+ open(unit=IOUT,file=trim(prname)//'bulk_c_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) bulk_c_kl_crust_mantle
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'bulk_betav_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) bulk_betav_kl_crust_mantle
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'bulk_betah_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) bulk_betah_kl_crust_mantle
+ close(IOUT)
! to check: isotropic kernels
- open(unit=27,file=trim(prname)//'alpha_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) alpha_kl_crust_mantle_reg
- close(27)
- open(unit=27,file=trim(prname)//'beta_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) beta_kl_crust_mantle_reg
- close(27)
- open(unit=27,file=trim(prname)//'bulk_beta_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) bulk_beta_kl_crust_mantle
- close(27)
+ open(unit=IOUT,file=trim(prname)//'alpha_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) alpha_kl_crust_mantle_reg
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'beta_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) beta_kl_crust_mantle_reg
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'bulk_beta_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) bulk_beta_kl_crust_mantle
+ close(IOUT)
else
! fully anisotropic kernels
! note: the C_ij and density kernels are not for relative perturbations (delta ln( m_i) = delta m_i / m_i),
! but absolute perturbations (delta m_i = m_i - m_0)
- open(unit=27,file=trim(prname)//'rho_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) - rho_kl_crust_mantle_reg
- close(27)
- open(unit=27,file=trim(prname)//'cijkl_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) - cijkl_kl_crust_mantle_reg
- close(27)
+ open(unit=IOUT,file=trim(prname)//'rho_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) - rho_kl_crust_mantle_reg
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'cijkl_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) - cijkl_kl_crust_mantle_reg
+ close(IOUT)
endif
else
! primary kernels: (rho,kappa,mu) parameterization
- open(unit=27,file=trim(prname)//'rhonotprime_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) rhonotprime_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'kappa_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) kappa_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'mu_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) mu_kl_crust_mantle
- close(27)
+ open(unit=IOUT,file=trim(prname)//'rhonotprime_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) rhonotprime_kl_crust_mantle
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'kappa_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) kappa_kl_crust_mantle
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'mu_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) mu_kl_crust_mantle
+ close(IOUT)
! (rho, alpha, beta ) parameterization
- open(unit=27,file=trim(prname)//'rho_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) rho_kl_crust_mantle_reg
- close(27)
- open(unit=27,file=trim(prname)//'alpha_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) alpha_kl_crust_mantle_reg
- close(27)
- open(unit=27,file=trim(prname)//'beta_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) beta_kl_crust_mantle_reg
- close(27)
+ open(unit=IOUT,file=trim(prname)//'rho_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) rho_kl_crust_mantle_reg
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'alpha_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) alpha_kl_crust_mantle_reg
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'beta_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) beta_kl_crust_mantle_reg
+ close(IOUT)
! (rho, bulk, beta ) parameterization, K_rho same as above
- open(unit=27,file=trim(prname)//'bulk_c_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) bulk_c_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'bulk_beta_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) bulk_beta_kl_crust_mantle
- close(27)
+ open(unit=IOUT,file=trim(prname)//'bulk_c_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) bulk_c_kl_crust_mantle
+ close(IOUT)
+ open(unit=IOUT,file=trim(prname)//'bulk_beta_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(IOUT) bulk_beta_kl_crust_mantle
+ close(IOUT)
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/specfem3D_par.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/specfem3D_par.F90 2013-08-30 12:45:44 UTC (rev 22746)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/specfem3D_par.F90 2013-08-30 15:40:46 UTC (rev 22747)
@@ -440,10 +440,12 @@
! kernels
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT) :: &
- rho_kl_crust_mantle,beta_kl_crust_mantle, alpha_kl_crust_mantle, Sigma_kl_crust_mantle
+ rho_kl_crust_mantle,beta_kl_crust_mantle,alpha_kl_crust_mantle
+
+ ! noise strength kernel
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:),allocatable :: Sigma_kl_crust_mantle
! For anisotropic kernels (see compute_kernels.f90 for a definition of the array)
- real(kind=CUSTOM_REAL), dimension(21,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT) :: &
- cijkl_kl_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:,:),allocatable :: cijkl_kl_crust_mantle
! approximate hessian
real(kind=CUSTOM_REAL), dimension(:,:,:,:),allocatable :: hess_kl_crust_mantle
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_seismograms.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_seismograms.f90 2013-08-30 12:45:44 UTC (rev 22746)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/specfem3D/write_seismograms.f90 2013-08-30 15:40:46 UTC (rev 22747)
@@ -144,6 +144,7 @@
integer :: iproc,sender,irec_local,irec,ier,receiver
integer :: nrec_local_received
+ integer :: nrec_tot_found
integer :: total_seismos,total_seismos_local
integer,dimension(:),allocatable:: islice_num_rec_local
character(len=256) :: sisname
@@ -154,6 +155,11 @@
allocate(one_seismogram(NDIM,NTSTEP_BETWEEN_OUTPUT_SEISMOS),stat=ier)
if(ier /= 0) call exit_mpi(myrank,'error while allocating one temporary seismogram')
+ ! check that the sum of the number of receivers in each slice is nrec
+ call sum_all_i(nrec_local,nrec_tot_found)
+ if(myrank == 0 .and. nrec_tot_found /= nrec) &
+ call exit_MPI(myrank,'total number of receivers is incorrect')
+
! get the base pathname for output files
call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
More information about the CIG-COMMITS
mailing list