From 6b2d6aca4dc9582aa6bd5a17384c79453aa3af60 Mon Sep 17 00:00:00 2001 From: Daniel Peter Date: Wed, 10 May 2023 14:21:14 +0200 Subject: [PATCH] supports reading (optional) attenuation Q-values (Qp & Qs) in tomography file --- EXAMPLES/tomographic_ocean_model/README | 8 + setup/constants.h.in | 3 + src/shared/read_material_table.f90 | 7 +- src/specfem2D/define_external_model.f90 | 20 +- .../define_external_model_from_marmousi.f90 | 10 +- .../define_external_model_from_tomo_file.f90 | 221 ++++++++++++++++-- src/specfem2D/read_external_model.f90 | 21 +- src/specfem2D/read_materials.f90 | 15 +- src/specfem2D/save_model_files.f90 | 1 + 9 files changed, 256 insertions(+), 50 deletions(-) diff --git a/EXAMPLES/tomographic_ocean_model/README b/EXAMPLES/tomographic_ocean_model/README index dc43c7995..4d8c0ffec 100644 --- a/EXAMPLES/tomographic_ocean_model/README +++ b/EXAMPLES/tomographic_ocean_model/README @@ -7,6 +7,14 @@ with topography variations. The tomographic file will be created by script createTomographyFile.py provided in utils/ directory. +The default tomography file will provide Vp, Vs and density on a data line like: + #x #z #vp #vs #rho +Optionally, attenuation Q-values can be added to these data lines, like: + #x #z #vp #vs #rho #Qp #Qs +to specify a 2D velocity and attenuation model. +You could for example modify the script createTomographyFile.py accordingly. + + To run this example, type: ./run_this_example.sh diff --git a/setup/constants.h.in b/setup/constants.h.in index 0310ab7f2..9bd3cdfdc 100644 --- a/setup/constants.h.in +++ b/setup/constants.h.in @@ -121,6 +121,9 @@ ! for viscoacousticity logical, parameter :: USE_A_STRONG_FORMULATION_FOR_E1 = .false. +! maximum attenuation Q-value + double precision, parameter :: ATTENUATION_COMP_MAXIMUM = 9999.d0 + !!----------------------------------------------------------- !! !! noise simulations diff --git a/src/shared/read_material_table.f90 b/src/shared/read_material_table.f90 index da0fb19c5..a162515e5 100644 --- a/src/shared/read_material_table.f90 +++ b/src/shared/read_material_table.f90 @@ -35,7 +35,8 @@ subroutine read_material_table() ! reads in material definitions in DATA/Par_file - use constants, only: IMAIN,TINYVAL,ISOTROPIC_MATERIAL,ANISOTROPIC_MATERIAL,POROELASTIC_MATERIAL + use constants, only: IMAIN,TINYVAL,ISOTROPIC_MATERIAL,ANISOTROPIC_MATERIAL,POROELASTIC_MATERIAL, & + ATTENUATION_COMP_MAXIMUM use shared_parameters, only: AXISYM,nbmodels,icodemat,cp,cs, & aniso3,aniso4,aniso5,aniso6,aniso7,aniso8,aniso9,aniso10,aniso11,aniso12, & @@ -116,8 +117,8 @@ subroutine read_material_table() aniso11(:) = 0.d0 comp_g(:) = 0.0d0 - QKappa(:) = 9999.d0 - Qmu(:) = 9999.d0 + QKappa(:) = ATTENUATION_COMP_MAXIMUM + Qmu(:) = ATTENUATION_COMP_MAXIMUM rho_s_read(:) = 0.d0 rho_f_read(:) = 0.d0 diff --git a/src/specfem2D/define_external_model.f90 b/src/specfem2D/define_external_model.f90 index d00ccdc28..94c3ea6ae 100644 --- a/src/specfem2D/define_external_model.f90 +++ b/src/specfem2D/define_external_model.f90 @@ -35,7 +35,7 @@ subroutine define_external_model_dummy(coord,material_element,ibool, & rho,vp,vs,QKappa_attenuation,Qmu_attenuation, & c11,c12,c13,c15,c23,c25,c33,c35,c55) - use constants, only: CUSTOM_REAL,NGLLX,NGLLZ,NDIM,IMAIN + use constants, only: CUSTOM_REAL,NGLLX,NGLLZ,NDIM,IMAIN,ATTENUATION_COMP_MAXIMUM use specfem_par, only: myrank,nspec,nglob @@ -75,8 +75,8 @@ subroutine define_external_model_dummy(coord,material_element,ibool, & endif ! default no attenuation - QKappa_attenuation(:,:,:) = 9999.d0 - Qmu_attenuation(:,:,:) = 9999.d0 + QKappa_attenuation(:,:,:) = ATTENUATION_COMP_MAXIMUM + Qmu_attenuation(:,:,:) = ATTENUATION_COMP_MAXIMUM ! loop on all the elements of the mesh, and inside each element loop on all the GLL points do ispec = 1,nspec @@ -92,8 +92,8 @@ subroutine define_external_model_dummy(coord,material_element,ibool, & rho(i,j,ispec) = 2000.d0 + 10.d0*x - 32.17d0*z**3 vp(i,j,ispec) = 3000.d0 - 7.12d0*cos(x) vs(i,j,ispec) = vp(i,j,ispec) / sqrt(3.d0) - QKappa_attenuation(i,j,ispec) = 9999. ! this means no attenuation - Qmu_attenuation(i,j,ispec) = 9999. ! this means no attenuation + QKappa_attenuation(i,j,ispec) = ATTENUATION_COMP_MAXIMUM ! this means no attenuation + Qmu_attenuation(i,j,ispec) = ATTENUATION_COMP_MAXIMUM ! this means no attenuation c11(i,j,ispec) = 169.d9 c13(i,j,ispec) = 122.d9 c15(i,j,ispec) = 0.d0 @@ -142,7 +142,7 @@ subroutine define_external_model(coord,material_element,ibool, & rho,vp,vs,QKappa_attenuation,Qmu_attenuation, & c11,c12,c13,c15,c23,c25,c33,c35,c55) - use constants, only: CUSTOM_REAL,NGLLX,NGLLZ,NDIM,IMAIN + use constants, only: CUSTOM_REAL,NGLLX,NGLLZ,NDIM,IMAIN,ATTENUATION_COMP_MAXIMUM use specfem_par, only: myrank,nspec,nglob @@ -214,8 +214,8 @@ subroutine define_external_model(coord,material_element,ibool, & endif ! default no attenuation - QKappa_attenuation(:,:,:) = 9999.d0 - Qmu_attenuation(:,:,:) = 9999.d0 + QKappa_attenuation(:,:,:) = ATTENUATION_COMP_MAXIMUM + Qmu_attenuation(:,:,:) = ATTENUATION_COMP_MAXIMUM ! default no anisotropy c11(:,:,:) = 0.d0 @@ -1118,8 +1118,8 @@ subroutine define_external_model(coord,material_element,ibool, & ! also set fictitious attenuation to a very high value (attenuation is not used in the fluid) if (material_element(ispec) == IREGION_OUTER_CORE) then vs(i,j,ispec) = 0.d0 - Qkappa_attenuation(i,j,ispec) = 9999.d0 - Qmu_attenuation(i,j,ispec) = 9999.d0 + Qkappa_attenuation(i,j,ispec) = ATTENUATION_COMP_MAXIMUM + Qmu_attenuation(i,j,ispec) = ATTENUATION_COMP_MAXIMUM endif enddo diff --git a/src/specfem2D/define_external_model_from_marmousi.f90 b/src/specfem2D/define_external_model_from_marmousi.f90 index 716708939..36c239a85 100644 --- a/src/specfem2D/define_external_model_from_marmousi.f90 +++ b/src/specfem2D/define_external_model_from_marmousi.f90 @@ -37,7 +37,7 @@ subroutine define_external_model_from_marmousi(coord,ibool,rho,vp,vs,QKappa_atte ! reads in external model using a marmousi format which defines a compaction gradient - use constants, only: CUSTOM_REAL,NGLLX,NGLLZ,NDIM,IMAIN + use constants, only: CUSTOM_REAL,NGLLX,NGLLZ,NDIM,IMAIN,ATTENUATION_COMP_MAXIMUM use specfem_par, only: poroelastcoef,density,kmato,myrank,nspec,nglob @@ -66,8 +66,8 @@ subroutine define_external_model_from_marmousi(coord,ibool,rho,vp,vs,QKappa_atte endif ! default no attenuation - QKappa_attenuation(:,:,:) = 9999.d0 - Qmu_attenuation(:,:,:) = 9999.d0 + QKappa_attenuation(:,:,:) = ATTENUATION_COMP_MAXIMUM + Qmu_attenuation(:,:,:) = ATTENUATION_COMP_MAXIMUM ! loop on all the elements of the mesh, and inside each element loop on all the GLL points do ispec = 1,nspec @@ -93,8 +93,8 @@ subroutine define_external_model_from_marmousi(coord,ibool,rho,vp,vs,QKappa_atte ! assumes Poisson solids vs(i,j,ispec) = vp(i,j,ispec) / sqrt(3.d0) - QKappa_attenuation(i,j,ispec) = 9999. ! this means no attenuation - Qmu_attenuation(i,j,ispec) = 9999. ! this means no attenuation + QKappa_attenuation(i,j,ispec) = ATTENUATION_COMP_MAXIMUM ! this means no attenuation + Qmu_attenuation(i,j,ispec) = ATTENUATION_COMP_MAXIMUM ! this means no attenuation c11(i,j,ispec) = 0.d0 ! this means no anisotropy c13(i,j,ispec) = 0.d0 diff --git a/src/specfem2D/define_external_model_from_tomo_file.f90 b/src/specfem2D/define_external_model_from_tomo_file.f90 index 8bb0e5ea1..6a8e91d3d 100644 --- a/src/specfem2D/define_external_model_from_tomo_file.f90 +++ b/src/specfem2D/define_external_model_from_tomo_file.f90 @@ -50,6 +50,7 @@ module model_tomography_par ! models parameter records double precision, dimension(:), allocatable :: x_tomo,z_tomo double precision, dimension(:,:), allocatable :: vp_tomo,vs_tomo,rho_tomo + double precision, dimension(:,:), allocatable :: qp_tomo,qs_tomo ! models entries integer :: NX,NZ @@ -58,6 +59,9 @@ module model_tomography_par ! min/max statistics double precision :: VP_MIN,VS_MIN,RHO_MIN,VP_MAX,VS_MAX,RHO_MAX + ! (optional) Q values + logical :: has_q_values + end module model_tomography_par ! @@ -345,7 +349,7 @@ subroutine define_external_model_from_tomo_file(rhoext,vpext,vsext, & ! ---------------------------------------------------------------------------------------- use specfem_par, only: tomo_material,coord,nspec,ibool,kmato, & - QKappa_attenuationcoef,Qmu_attenuationcoef,anisotropycoef, & + Qkappa_attenuationcoef,Qmu_attenuationcoef,anisotropycoef, & poroelastcoef,density use specfem_par, only: myrank,TOMOGRAPHY_FILE @@ -353,18 +357,22 @@ subroutine define_external_model_from_tomo_file(rhoext,vpext,vsext, & use model_tomography_par use interpolation - use constants, only: NGLLX,NGLLZ,TINYVAL,IMAIN,CUSTOM_REAL + use constants, only: NGLLX,NGLLZ,TINYVAL,IMAIN,CUSTOM_REAL,ATTENUATION_COMP_MAXIMUM implicit none real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec), intent(out) :: rhoext,vpext,vsext - real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec), intent(out) :: QKappa_attenuationext,Qmu_attenuationext + real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec), intent(out) :: Qkappa_attenuationext,Qmu_attenuationext real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec), intent(out) :: c11ext,c12ext,c13ext,c15ext,c22ext,c23ext,c25ext, & c33ext,c35ext,c55ext ! local parameters integer :: i,j,ispec,iglob double precision :: xmesh,zmesh + double precision :: rho_final + double precision :: vp_final,vs_final + double precision :: qp_final,qs_final + double precision :: L_val,qmu_atten,qkappa_atten ! stats integer :: npoint_tomo,npoint_internal @@ -378,8 +386,8 @@ subroutine define_external_model_from_tomo_file(rhoext,vpext,vsext, & call read_tomo_file() ! default no attenuation - QKappa_attenuationext(:,:,:) = 9999.d0 - Qmu_attenuationext(:,:,:) = 9999.d0 + QKappa_attenuationext(:,:,:) = ATTENUATION_COMP_MAXIMUM + Qmu_attenuationext(:,:,:) = ATTENUATION_COMP_MAXIMUM ! default no anisotropy c11ext(:,:,:) = 0.d0 @@ -405,9 +413,50 @@ subroutine define_external_model_from_tomo_file(rhoext,vpext,vsext, & xmesh = coord(1,iglob) zmesh = coord(2,iglob) - rhoext(i,j,ispec) = interpolate(NX, x_tomo, NZ, z_tomo, rho_tomo, xmesh, zmesh,TINYVAL) - vpext(i,j,ispec) = interpolate(NX, x_tomo, NZ, z_tomo, vp_tomo, xmesh, zmesh,TINYVAL) - vsext(i,j,ispec) = interpolate(NX, x_tomo, NZ, z_tomo, vs_tomo, xmesh, zmesh,TINYVAL) + rho_final = interpolate(NX, x_tomo, NZ, z_tomo, rho_tomo, xmesh, zmesh, TINYVAL) + vp_final = interpolate(NX, x_tomo, NZ, z_tomo, vp_tomo, xmesh, zmesh, TINYVAL) + vs_final = interpolate(NX, x_tomo, NZ, z_tomo, vs_tomo, xmesh, zmesh, TINYVAL) + + rhoext(i,j,ispec) = rho_final + vpext(i,j,ispec) = vp_final + vsext(i,j,ispec) = vs_final + + if (has_q_values) then + qp_final = interpolate(NX, x_tomo, NZ, z_tomo, qp_tomo, xmesh, zmesh, TINYVAL) + qs_final = interpolate(NX, x_tomo, NZ, z_tomo, qs_tomo, xmesh, zmesh, TINYVAL) + + ! attenuation zero (means negligible attenuation) + if (qs_final <= TINYVAL) qs_final = ATTENUATION_COMP_MAXIMUM + if (qp_final <= TINYVAL) qp_final = ATTENUATION_COMP_MAXIMUM + + ! Anderson & Hart (1978) conversion between (Qp,Qs) and (Qkappa,Qmu) + ! factor L + L_val = 4.d0/3.d0 * (vs_final/vp_final)**2 + + ! shear attenuation + qmu_atten = qs_final + + ! converts to bulk attenuation + if (abs(qs_final - L_val * qp_final) <= TINYVAL) then + ! negligible bulk attenuation + qkappa_atten = ATTENUATION_COMP_MAXIMUM + else + qkappa_atten = (1.d0 - L_val) * qp_final * qs_final / (qs_final - L_val * qp_final) + endif + + ! attenuation zero (means negligible attenuation) + if (qmu_atten <= TINYVAL) qmu_atten = ATTENUATION_COMP_MAXIMUM + if (qkappa_atten <= TINYVAL) qkappa_atten = ATTENUATION_COMP_MAXIMUM + + ! limits Q values + if (qmu_atten < 1.0d0) qmu_atten = 1.0d0 + if (qmu_atten > ATTENUATION_COMP_MAXIMUM) qmu_atten = ATTENUATION_COMP_MAXIMUM + if (qkappa_atten < 1.0d0) qkappa_atten = 1.0d0 + if (qkappa_atten > ATTENUATION_COMP_MAXIMUM) qkappa_atten = ATTENUATION_COMP_MAXIMUM + + Qkappa_attenuationext(i,j,ispec) = qkappa_atten + Qmu_attenuationext(i,j,ispec) = qmu_atten + endif !! ABAB : The 3 following lines are important, otherwise PMLs won't work. TODO check that !! (we assign these values several times: indeed for each kmato(ispec) it can exist a lot of rhoext(i,j,ispec) ) @@ -514,6 +563,11 @@ subroutine read_tomo_file() ! ... ! x(NX) z(NZ) vp vs rho ! +! +! in case Q-values are also provided, the data lines here would have a format like: +! x(1) z(1) vp vs rho Qp Qs +! .. + ! Where : ! _x and z must be increasing ! _ORIGIN_X, END_X are, respectively, the coordinates of the initial and final tomographic @@ -538,16 +592,23 @@ subroutine read_tomo_file() use specfem_par, only: myrank,TOMOGRAPHY_FILE use model_tomography_par - use constants, only: IIN,IMAIN,HUGEVAL + use constants, only: IIN,IMAIN,HUGEVAL,MAX_STRING_LEN implicit none ! local parameters integer :: ier,irecord,i,j - character(len=150) :: string_read + character(len=MAX_STRING_LEN) :: string_read double precision, dimension(:), allocatable :: x_tomography,z_tomography,vp_tomography,vs_tomography,rho_tomography - double precision :: read_rho_min,read_rho_max,read_vp_min,read_vp_max,read_vs_min,read_vs_max,tmp_dp + double precision, dimension(:), allocatable :: qp_tomography,qs_tomography + double precision :: read_rho_min,read_rho_max + double precision :: read_vp_min,read_vp_max + double precision :: read_vs_min,read_vs_max + double precision :: read_qp_min,read_qp_max + double precision :: read_qs_min,read_qs_max + double precision :: tmp_dp + integer :: ntokens ! opens file for reading open(unit=IIN,file=trim(TOMOGRAPHY_FILE),status='old',action='read',iostat=ier) @@ -611,11 +672,60 @@ subroutine read_tomo_file() read_vs_max = - HUGEVAL read_rho_min = + HUGEVAL read_rho_max = - HUGEVAL + read_qp_min = + HUGEVAL + read_qp_max = - HUGEVAL + read_qs_min = + HUGEVAL + read_qs_max = - HUGEVAL do while (irecord < nrecord) call tomo_read_next_line(IIN,string_read) - read(string_read,*) x_tomography(irecord+1),z_tomography(irecord+1), & - vp_tomography(irecord+1),vs_tomography(irecord+1),rho_tomography(irecord+1) + + if (irecord == 0) then + ! checks number of entries of first data line + call tomo_get_number_of_tokens(string_read,ntokens) + + !debug + !print *,'tomography file: number of tokens on first data line: ',ntokens,' line: ***'//trim(string_read)//'***' + + ! data line formats: + ! #x(1) #z(1) #vp #vs #rho + ! or + ! #x(1) #z(1) #vp #vs #rho #Qp #Qs + if (ntokens /= 5 .and. ntokens /= 7) then + print *,'Error reading tomography file, data line has wrong number of entries: ',trim(string_read) + stop 'Error reading tomography file' + endif + + ! determines data format + if (ntokens == 7) then + has_q_values = .true. + else + has_q_values = .false. + endif + + if (has_q_values) then + ! allocate models parameter records + allocate(qp_tomography(nrecord),qs_tomography(nrecord),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'not enough memory to allocate tomo Q arrays') + qp_tomography(:) = 0.d0; qs_tomography(:) = 0.d0 + + allocate(qp_tomo(NX,NZ),qs_tomo(NX,NZ),stat=ier) + if (ier /= 0) call exit_MPI(myrank,'not enough memory to allocate tomo arrays') + qp_tomo(:,:) = 0.d0; qs_tomo(:,:) = 0.d0 + endif + endif + + ! reads in tomo values + if (has_q_values) then + ! #x(1) #z(1) #vp #vs #rho #Qp #Qs + read(string_read,*) x_tomography(irecord+1),z_tomography(irecord+1), & + vp_tomography(irecord+1),vs_tomography(irecord+1),rho_tomography(irecord+1), & + qp_tomography(irecord+1),qs_tomography(irecord+1) + else + ! #x(1) #z(1) #vp #vs #rho + read(string_read,*) x_tomography(irecord+1),z_tomography(irecord+1), & + vp_tomography(irecord+1),vs_tomography(irecord+1),rho_tomography(irecord+1) + endif if (irecord < NX) x_tomo(irecord+1) = x_tomography(irecord+1) @@ -627,6 +737,13 @@ subroutine read_tomo_file() if (rho_tomography(irecord+1) > read_rho_max) read_rho_max = rho_tomography(irecord+1) if (rho_tomography(irecord+1) < read_rho_min) read_rho_min = rho_tomography(irecord+1) + if (has_q_values) then + if (qp_tomography(irecord+1) > read_qp_max) read_qp_max = qp_tomography(irecord+1) + if (qp_tomography(irecord+1) < read_qp_min) read_qp_min = qp_tomography(irecord+1) + if (qs_tomography(irecord+1) > read_qs_max) read_qs_max = qs_tomography(irecord+1) + if (qs_tomography(irecord+1) < read_qs_min) read_qs_min = qs_tomography(irecord+1) + endif + ! counter irecord = irecord + 1 enddo @@ -641,6 +758,15 @@ subroutine read_tomo_file() enddo enddo + if (has_q_values) then + do i = 1,NX + do j = 1,NZ + qp_tomo(i,j) = qp_tomography(NX*(j-1)+i) + qs_tomo(i,j) = qs_tomography(NX*(j-1)+i) + enddo + enddo + endif + call synchronize_all() if (irecord /= nrecord .and. myrank == 0) then @@ -666,6 +792,19 @@ subroutine read_tomo_file() tmp_dp = read_rho_max call max_all_dp(tmp_dp,read_rho_max) + if (has_q_values) then + ! Q-model values + tmp_dp = read_qp_min + call min_all_dp(tmp_dp,read_qp_min) + tmp_dp = read_qp_max + call max_all_dp(tmp_dp,read_qp_max) + + tmp_dp = read_qs_min + call min_all_dp(tmp_dp,read_qs_min) + tmp_dp = read_qs_max + call max_all_dp(tmp_dp,read_qs_max) + endif + ! user output if (myrank == 0) then write(IMAIN,*) ' Number of grid points = NX*NZ:',nrecord @@ -674,12 +813,20 @@ subroutine read_tomo_file() write(IMAIN,*) ' vs : min/max = ',sngl(read_vs_min),' / ',sngl(read_vs_max) write(IMAIN,*) ' rho: min/max = ',sngl(read_rho_min),' / ',sngl(read_rho_max) write(IMAIN,*) + if (has_q_values) then + write(IMAIN,*) ' Qp : min/max = ',sngl(read_qp_min),' / ',sngl(read_qp_max) + write(IMAIN,*) ' Qs : min/max = ',sngl(read_qs_min),' / ',sngl(read_qs_max) + write(IMAIN,*) + endif call flush_IMAIN() endif ! closes file close(IIN) deallocate(x_tomography,z_tomography,vp_tomography,vs_tomography,rho_tomography) + if (has_q_values) then + deallocate(qp_tomography,qs_tomography) + endif end subroutine read_tomo_file @@ -691,10 +838,11 @@ subroutine tomo_read_next_line(unit_in,string_read) ! Store next line of file unit_in in string_read + use constants, only: MAX_STRING_LEN implicit none integer :: unit_in - character(len=150) :: string_read + character(len=MAX_STRING_LEN) :: string_read integer :: ier @@ -708,8 +856,10 @@ subroutine tomo_read_next_line(unit_in,string_read) ! suppress trailing carriage return (ASCII code 13) if any (e.g. if input text file coming from Windows/DOS) if (index(string_read,achar(13)) > 0) string_read = string_read(1:index(string_read,achar(13))-1) - ! exit loop when we find the first line that is not a comment or a white line + ! reads next line if empty if (len_trim(string_read) == 0) cycle + + ! exit loop when we find the first line that is not a comment or a white line if (string_read(1:1) /= '#') exit enddo @@ -725,3 +875,44 @@ subroutine tomo_read_next_line(unit_in,string_read) end subroutine tomo_read_next_line +! +!------------------------------------------------------------------------------------------------- +! + + subroutine tomo_get_number_of_tokens(string_read,ntokens) + + use constants, only: MAX_STRING_LEN + implicit none + + character(len=MAX_STRING_LEN),intent(in) :: string_read + integer,intent(out) :: ntokens + + ! local parameters + integer :: i + logical :: previous_is_delim + character(len=1), parameter :: delim_space = ' ' + character(len=1), parameter :: delim_tab = achar(9) ! tab delimiter + + ! initializes + ntokens = 0 + + ! checks if anything to do + if (len_trim(string_read) == 0) return + + ! counts tokens + ntokens = 1 + previous_is_delim = .true. + do i = 1, len_trim(string_read) + ! finds next delimiter (space or tabular) + if (string_read(i:i) == delim_space .or. string_read(i:i) == delim_tab) then + if (.not. previous_is_delim) then + ntokens = ntokens + 1 + previous_is_delim = .true. + endif + else + if (previous_is_delim) previous_is_delim = .false. + endif + enddo + + end subroutine tomo_get_number_of_tokens + diff --git a/src/specfem2D/read_external_model.f90 b/src/specfem2D/read_external_model.f90 index d2af0df9e..383227b44 100644 --- a/src/specfem2D/read_external_model.f90 +++ b/src/specfem2D/read_external_model.f90 @@ -36,7 +36,8 @@ subroutine read_external_model(rhoext,vpext,vsext,QKappa_attenuationext,Qmu_atte ! reads in external model files - use constants, only: CUSTOM_REAL,NGLLX,NGLLZ,TINYVAL,IMAIN,IN_DATA_FILES,IIN,MAX_STRING_LEN + use constants, only: CUSTOM_REAL,NGLLX,NGLLZ,TINYVAL,IMAIN,IN_DATA_FILES,IIN,MAX_STRING_LEN, & + ATTENUATION_COMP_MAXIMUM use specfem_par, only: nspec,ibool,ispec_is_elastic,ispec_is_anisotropic, & coord,kmato,MODEL,myrank,setup_with_binary_database @@ -135,8 +136,8 @@ subroutine read_external_model(rhoext,vpext,vsext,QKappa_attenuationext,Qmu_atte close(IIN) ! default no attenuation - QKappa_attenuationext(:,:,:) = 9999.d0 - Qmu_attenuationext(:,:,:) = 9999.d0 + QKappa_attenuationext(:,:,:) = ATTENUATION_COMP_MAXIMUM + Qmu_attenuationext(:,:,:) = ATTENUATION_COMP_MAXIMUM case ('ascii') ! ascii model format @@ -159,8 +160,8 @@ subroutine read_external_model(rhoext,vpext,vsext,QKappa_attenuationext,Qmu_atte enddo close(IIN) ! default no attenuation - QKappa_attenuationext(:,:,:) = 9999.d0 - Qmu_attenuationext(:,:,:) = 9999.d0 + QKappa_attenuationext(:,:,:) = ATTENUATION_COMP_MAXIMUM + Qmu_attenuationext(:,:,:) = ATTENUATION_COMP_MAXIMUM case ('binary','gll') ! binary formats @@ -205,7 +206,7 @@ subroutine read_external_model(rhoext,vpext,vsext,QKappa_attenuationext,Qmu_atte read(IIN) QKappa_attenuationext close(IIN) - Qmu_attenuationext(:,:,:) = 9999.d0 + Qmu_attenuationext(:,:,:) = ATTENUATION_COMP_MAXIMUM ! user output if (myrank == 0) then @@ -240,8 +241,8 @@ subroutine read_external_model(rhoext,vpext,vsext,QKappa_attenuationext,Qmu_atte else ! default no attenuation - QKappa_attenuationext(:,:,:) = 9999.d0 - Qmu_attenuationext(:,:,:) = 9999.d0 + QKappa_attenuationext(:,:,:) = ATTENUATION_COMP_MAXIMUM + Qmu_attenuationext(:,:,:) = ATTENUATION_COMP_MAXIMUM endif @@ -326,8 +327,8 @@ subroutine read_external_model(rhoext,vpext,vsext,QKappa_attenuationext,Qmu_atte enddo ! default no attenuation - QKappa_attenuationext(:,:,:) = 9999.d0 - Qmu_attenuationext(:,:,:) = 9999.d0 + QKappa_attenuationext(:,:,:) = ATTENUATION_COMP_MAXIMUM + Qmu_attenuationext(:,:,:) = ATTENUATION_COMP_MAXIMUM case ('external') ! generic model defined in external files diff --git a/src/specfem2D/read_materials.f90 b/src/specfem2D/read_materials.f90 index adc24ad49..2316f407d 100644 --- a/src/specfem2D/read_materials.f90 +++ b/src/specfem2D/read_materials.f90 @@ -36,7 +36,8 @@ subroutine read_materials(f0) ! reads properties of a 2D isotropic or anisotropic linear elastic element use constants, only: IIN,IMAIN,ZERO,FOUR_THIRDS,TWO_THIRDS,HALF,TINYVAL, & - ISOTROPIC_MATERIAL,ANISOTROPIC_MATERIAL,POROELASTIC_MATERIAL + ISOTROPIC_MATERIAL,ANISOTROPIC_MATERIAL,POROELASTIC_MATERIAL, & + ATTENUATION_COMP_MAXIMUM use specfem_par, only: AXISYM,density,porosity,tortuosity,anisotropycoef,permeability,poroelastcoef, & numat,myrank,QKappa_attenuationcoef,Qmu_attenuationcoef, & @@ -72,8 +73,8 @@ subroutine read_materials(f0) poroelastcoef(:,:,:) = ZERO anisotropycoef(:,:) = ZERO - QKappa_attenuationcoef(:) = 9999. - Qmu_attenuationcoef(:) = 9999. + QKappa_attenuationcoef(:) = ATTENUATION_COMP_MAXIMUM + Qmu_attenuationcoef(:) = ATTENUATION_COMP_MAXIMUM ! Index of the material that will be defined by an external tomo file if needed (TOMOGRAPHY_FILE) tomo_material = 0 @@ -97,8 +98,8 @@ subroutine read_materials(f0) c23 = ZERO c25 = ZERO c22 = ZERO - Qkappa = 9999. - Qmu = 9999. + Qkappa = ATTENUATION_COMP_MAXIMUM + Qmu = ATTENUATION_COMP_MAXIMUM ! supported model formats: ! acoustic - model_number 1 rho Vp 0 0 0 QKappa Qmu 0 0 0 0 0 0 @@ -393,8 +394,8 @@ subroutine read_materials(f0) poroelastcoef(3,1,n) = -1.0d0 poroelastcoef(4,1,n) = ZERO - QKappa_attenuationcoef(n) = 9999. - Qmu_attenuationcoef(n) = 9999. + QKappa_attenuationcoef(n) = ATTENUATION_COMP_MAXIMUM + Qmu_attenuationcoef(n) = ATTENUATION_COMP_MAXIMUM if (mu > TINYVAL) then porosity(n) = 0.d0 diff --git a/src/specfem2D/save_model_files.f90 b/src/specfem2D/save_model_files.f90 index 19d3f1f23..849c4956c 100644 --- a/src/specfem2D/save_model_files.f90 +++ b/src/specfem2D/save_model_files.f90 @@ -88,6 +88,7 @@ subroutine save_model_files() endif if (ATTENUATION_VISCOELASTIC) then + ! safety stop if (ATTENUATION_VISCOACOUSTIC) & call exit_MPI(myrank,'Not possible yet to save model with both acoustic and elastic attenuation')