[cig-commits] r17259 - seismo/3D/SPECFEM3D_GLOBE/trunk
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Mon Oct 11 14:16:31 PDT 2010
Author: danielpeter
Date: 2010-10-11 14:16:31 -0700 (Mon, 11 Oct 2010)
New Revision: 17259
Added:
seismo/3D/SPECFEM3D_GLOBE/trunk/model_gapp2.f90
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/Makefile.in
seismo/3D/SPECFEM3D_GLOBE/trunk/check_simulation_stability.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in
seismo/3D/SPECFEM3D_GLOBE/trunk/create_chunk_buffers.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/get_model_parameters.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/initialize_simulation.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/locate_receivers.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D_models.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/model_gll.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/model_ppm.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/model_s362ani.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/prepare_timerun.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/save_arrays_solver.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/save_header_file.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90
Log:
adds 3D mantle model from Masayuki Obayashi; minor modifications in models GLL and PPM; minor code bug fix in model s362ani
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/Makefile.in 2010-10-11 09:41:46 UTC (rev 17258)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/Makefile.in 2010-10-11 21:16:31 UTC (rev 17259)
@@ -129,6 +129,7 @@
$O/model_jp1d.o \
$O/model_jp3d.o \
$O/model_ppm.o \
+ $O/model_gapp2.o \
$O/model_prem.o \
$O/model_1dref.o \
$O/model_s20rts.o \
@@ -809,6 +810,9 @@
$O/model_ppm.o: constants.h $S/model_ppm.f90
${MPIFCCOMPILE_CHECK} -c -o $O/model_ppm.o ${FCFLAGS_f90} $S/model_ppm.f90
+$O/model_gapp2.o: constants.h $S/model_gapp2.f90
+ ${MPIFCCOMPILE_CHECK} -c -o $O/model_gapp2.o ${FCFLAGS_f90} $S/model_gapp2.f90
+
$O/model_s20rts.o: constants.h $S/model_s20rts.f90
${MPIFCCOMPILE_CHECK} -c -o $O/model_s20rts.o ${FCFLAGS_f90} $S/model_s20rts.f90
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/check_simulation_stability.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/check_simulation_stability.f90 2010-10-11 09:41:46 UTC (rev 17258)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/check_simulation_stability.f90 2010-10-11 21:16:31 UTC (rev 17259)
@@ -87,7 +87,7 @@
timestamp_remote,year_remote,mon_remote,day_remote,hr_remote,minutes_remote,day_of_week_remote
integer :: ier
integer, external :: idaywk
-
+
double precision,parameter :: scale_displ = R_EARTH
@@ -98,7 +98,7 @@
maxval(sqrt(displ_inner_core(1,:)**2 + displ_inner_core(2,:)**2 + displ_inner_core(3,:)**2)))
Ufluidnorm = maxval(abs(displ_outer_core))
-
+
! compute the maximum of the maxima for all the slices using an MPI reduction
call MPI_REDUCE(Usolidnorm,Usolidnorm_all,1,CUSTOM_MPI_TYPE,MPI_MAX,0, &
MPI_COMM_WORLD,ier)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in 2010-10-11 09:41:46 UTC (rev 17258)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in 2010-10-11 21:16:31 UTC (rev 17259)
@@ -352,6 +352,7 @@
integer, parameter :: THREE_D_MODEL_PPM = 9 ! format for point profile models
integer, parameter :: THREE_D_MODEL_GLL = 10 ! format for iterations with GLL mesh
integer, parameter :: THREE_D_MODEL_S40RTS = 11
+ integer, parameter :: THREE_D_MODEL_GAPP2 = 12
! define flag for regions of the global Earth for attenuation
integer, parameter :: NUM_REGIONS_ATTENUATION = 5
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/create_chunk_buffers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/create_chunk_buffers.f90 2010-10-11 09:41:46 UTC (rev 17258)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/create_chunk_buffers.f90 2010-10-11 21:16:31 UTC (rev 17259)
@@ -161,6 +161,19 @@
write(IMAIN,*)
endif
+ ! initializes counters
+ NUM_FACES = 0
+ NUM_MSG_TYPES = 0
+ iproc_xi_send = 0
+ iproc_xi_receive = 0
+ iproc_eta_send = 0
+ iproc_eta_receive = 0
+ iproc_edge_send = 0
+ iproc_edge_receive = 0
+ iedge = 0
+ ichunk_receive = 0
+ ichunk_send = 0
+
! number of corners and faces shared between chunks and number of message types
if(NCHUNKS == 1 .or. NCHUNKS == 2) then
NCORNERSCHUNKS = 1
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/get_model_parameters.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/get_model_parameters.f90 2010-10-11 09:41:46 UTC (rev 17258)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/get_model_parameters.f90 2010-10-11 21:16:31 UTC (rev 17259)
@@ -141,7 +141,8 @@
! checks with '_1Dcrust' option
impose_1Dcrust = .false.
ending_1Dcrust = ' '
- if( len_trim(MODEL_ROOT) > 8 ) ending_1Dcrust = MODEL_ROOT(len_trim(MODEL_ROOT)-7:len_trim(MODEL_ROOT))
+ if( len_trim(MODEL_ROOT) > 8 ) &
+ ending_1Dcrust = MODEL_ROOT(len_trim(MODEL_ROOT)-7:len_trim(MODEL_ROOT))
if( ending_1Dcrust == '_1Dcrust' ) then
impose_1Dcrust = .true.
! in case it has an ending for the inner core, remove it from the name
@@ -387,6 +388,15 @@
! model, based upon the s29ea model, but putting mgll_v as parameter to this
! routine involves too many changes. )
+ else if(MODEL == 'gapp2') then
+ CASE_3D = .true.
+ CRUSTAL = .true.
+ ISOTROPIC_3D_MANTLE = .true.
+ ONE_CRUST = .true.
+ REFERENCE_1D_MODEL = REFERENCE_MODEL_PREM
+ THREE_D_MODEL = THREE_D_MODEL_GAPP2
+ TRANSVERSE_ISOTROPY = .true.
+
else
print*
print*,'error model: ',trim(MODEL)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/initialize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/initialize_simulation.f90 2010-10-11 09:41:46 UTC (rev 17258)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/initialize_simulation.f90 2010-10-11 21:16:31 UTC (rev 17259)
@@ -411,18 +411,18 @@
if (SIMULATION_TYPE /= 1 .and. NSOURCES > 999999) &
call exit_MPI(myrank, 'for adjoint simulations, NSOURCES <= 999999, if you need more change i6.6 in write_seismograms.f90')
-
- if((SIMULATION_TYPE == 1 .and. SAVE_FORWARD) .or. SIMULATION_TYPE == 3) then
+
+ if((SIMULATION_TYPE == 1 .and. SAVE_FORWARD) .or. SIMULATION_TYPE == 3) then
if ( ATTENUATION_VAL) then
! checks mimic flag:
! attenuation for adjoint simulations must have USE_ATTENUATION_MIMIC set by xcreate_header_file
if( USE_ATTENUATION_MIMIC .eqv. .false. ) &
- call exit_MPI(myrank,'error in compiled attenuation parameters, please recompile solver 17')
+ call exit_MPI(myrank,'error in compiled attenuation parameters, please recompile solver 17b')
! user output
if( myrank == 0 ) write(IMAIN,*) 'incorporates ATTENUATION for time-reversed simulation'
- endif
-
+ endif
+
! checks adjoint array dimensions
if(NSPEC_CRUST_MANTLE_ADJOINT /= NSPEC_CRUST_MANTLE &
.or. NSPEC_OUTER_CORE_ADJOINT /= NSPEC_OUTER_CORE &
@@ -438,7 +438,7 @@
if (NSPEC_CRUST_MANTLE_ATTENUAT /= NSPEC_CRUST_MANTLE) &
call exit_MPI(myrank, 'NSPEC_CRUST_MANTLE_ATTENUAT /= NSPEC_CRUST_MANTLE, exit')
if (NSPEC_INNER_CORE_ATTENUATION /= NSPEC_INNER_CORE) &
- call exit_MPI(myrank, 'NSPEC_INNER_CORE_ATTENUATION /= NSPEC_INNER_CORE, exit')
+ call exit_MPI(myrank, 'NSPEC_INNER_CORE_ATTENUATION /= NSPEC_INNER_CORE, exit')
endif
! checks strain storage
@@ -448,7 +448,7 @@
call exit_MPI(myrank, 'error in compiled compute_and_store_strain parameter, please recompile solver 19')
else
if( COMPUTE_AND_STORE_STRAIN .neqv. .false. ) &
- call exit_MPI(myrank, 'error in compiled compute_and_store_strain parameter, please recompile solver 20')
+ call exit_MPI(myrank, 'error in compiled compute_and_store_strain parameter, please recompile solver 20')
endif
if (SIMULATION_TYPE == 3 .and. (ANISOTROPIC_3D_MANTLE_VAL .or. ANISOTROPIC_INNER_CORE_VAL)) &
@@ -461,7 +461,7 @@
endif
if( SIMULATION_TYPE == 3 ) then
if( .not. ANISOTROPIC_KL ) then
- call exit_mpi(myrank,'error SAVE_TRANSVERSE_KL: needs anisotropic kernel calculations')
+ call exit_mpi(myrank,'error SAVE_TRANSVERSE_KL: needs anisotropic kernel calculations')
endif
endif
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/locate_receivers.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/locate_receivers.f90 2010-10-11 09:41:46 UTC (rev 17258)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/locate_receivers.f90 2010-10-11 21:16:31 UTC (rev 17259)
@@ -198,10 +198,16 @@
! read that STATIONS file on the master
if(myrank == 0) then
call get_value_string(STATIONS, 'solver.STATIONS', rec_filename)
- open(unit=1,file=STATIONS,status='old',action='read')
+ open(unit=1,file=STATIONS,status='old',action='read',iostat=ier)
+ if( ier /= 0 ) call exit_MPI(myrank,'error opening STATIONS file')
+
! loop on all the stations to read station information
do irec = 1,nrec
- read(1,*) station_name(irec),network_name(irec),stlat(irec),stlon(irec),stele(irec),stbur(irec)
+ read(1,*,iostat=ier) station_name(irec),network_name(irec),stlat(irec),stlon(irec),stele(irec),stbur(irec)
+ if( ier /= 0 ) then
+ write(IMAIN,*) 'error reading in station ',irec
+ call exit_MPI(myrank,'error reading in station in STATIONS file')
+ endif
enddo
! close receiver file
close(1)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D_models.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D_models.f90 2010-10-11 09:41:46 UTC (rev 17258)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D_models.f90 2010-10-11 21:16:31 UTC (rev 17259)
@@ -501,6 +501,10 @@
! could use EUcrust07 Vp crustal structure
call model_eucrust_broadcast(myrank,EUCM_V)
+ case(THREE_D_MODEL_GAPP2)
+ ! GAP model
+ call model_gapp2_broadcast(myrank)
+
case default
call exit_MPI(myrank,'3D model not defined')
@@ -904,6 +908,15 @@
vsh=vsh*(1.0d0+dvs)
rho=rho*(1.0d0+drho)
+ case(THREE_D_MODEL_GAPP2 )
+ ! 3D GAP model (Obayashi)
+ call mantle_gapmodel(r_used,theta,phi,dvs,dvp,drho)
+ vpv=vpv*(1.0d0+dvp)
+ vph=vph*(1.0d0+dvp)
+ vsv=vsv*(1.0d0+dvs)
+ vsh=vsh*(1.0d0+dvs)
+ rho=rho*(1.0d0+drho)
+
case default
stop 'unknown 3D Earth model in meshfem3D_models_get3Dmntl_val() '
@@ -1313,6 +1326,7 @@
endif
! takes stored gll values from file
+ ! ( note that these values are non-dimensionalized)
if(CUSTOM_REAL == SIZE_REAL) then
vp = dble( MGLL_V%vp_new(i,j,k,ispec) )
vs = dble( MGLL_V%vs_new(i,j,k,ispec) )
@@ -1322,17 +1336,12 @@
vs = MGLL_V%vs_new(i,j,k,ispec)
rho = MGLL_V%rho_new(i,j,k,ispec)
endif
- ! non-dimensionalize
- vp = vp * MGLL_V%scale_velocity
- vs = vs * MGLL_V%scale_velocity
- rho = rho * MGLL_V%scale_density
! isotropic model
vpv = vp
vph = vp
vsv = vs
vsh = vs
rho = rho
- dvp = 0.0d0
eta_aniso = 1.0d0
! transverse isotropic model
@@ -1346,28 +1355,22 @@
! takes stored gll values from file
if(CUSTOM_REAL == SIZE_REAL) then
vph = dble( MGLL_V%vph_new(i,j,k,ispec) )
+ vpv = dble( MGLL_V%vpv_new(i,j,k,ispec) )
vsh = dble( MGLL_V%vsh_new(i,j,k,ispec) )
- vpv = dble( MGLL_V%vpv_new(i,j,k,ispec) )
vsv = dble( MGLL_V%vsv_new(i,j,k,ispec) )
rho = dble( MGLL_V%rho_new(i,j,k,ispec) )
eta_aniso = dble( MGLL_V%eta_new(i,j,k,ispec) )
else
vph = MGLL_V%vph_new(i,j,k,ispec)
+ vpv = MGLL_V%vpv_new(i,j,k,ispec)
vsh = MGLL_V%vsh_new(i,j,k,ispec)
- vpv = MGLL_V%vpv_new(i,j,k,ispec)
vsv = MGLL_V%vsv_new(i,j,k,ispec)
rho = MGLL_V%rho_new(i,j,k,ispec)
eta_aniso = MGLL_V%eta_new(i,j,k,ispec)
endif
- ! non-dimensionalize
- ! transverse isotropic model
- vpv = vpv * MGLL_V%scale_velocity
- vph = vph * MGLL_V%scale_velocity
- vsv = vsv * MGLL_V%scale_velocity
- vsh = vsh * MGLL_V%scale_velocity
- rho = rho * MGLL_V%scale_density
- dvp = 0.0d0
endif
+ ! no mantle vp perturbation
+ dvp = 0.0d0
endif ! MODEL_GLL
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/model_gapp2.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/model_gapp2.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/model_gapp2.f90 2010-10-11 21:16:31 UTC (rev 17259)
@@ -0,0 +1,224 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 5 . 0
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Princeton University, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) Princeton University / California Institute of Technology and University of Pau / CNRS / INRIA
+! March 2010
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+!--------------------------------------------------------------------------------------------------
+!
+! GAP P2 model - Global automatic parameterization model
+!
+! 3D Vp mantle model (version P2) from Masayuki Obayashi
+!
+!--------------------------------------------------------------------------------------------------
+
+
+ module gapp2_mantle_model_constants
+ ! data file resolution
+ integer, parameter :: ma=128,mo=256,mr=33,mr1=66
+ integer no,na,nnr,nr1
+ real dela,delo
+ ! allocatable model arrays
+ real,dimension(:),allocatable :: dep,dep1,vp1
+ real,dimension(:,:,:),allocatable :: vp3
+ end module gapp2_mantle_model_constants
+
+!
+!--------------------------------------------------------------------------------------------------
+!
+
+ subroutine model_gapp2_broadcast(myrank)
+
+! standard routine to setup model
+
+ use gapp2_mantle_model_constants
+
+ implicit none
+ include "constants.h"
+ ! standard include of the MPI library
+ include 'mpif.h'
+ integer :: myrank
+ integer :: ier
+
+ ! allocates arrays only when called and needed
+ allocate(dep(0:mr),dep1(0:mr1),vp1(0:mr1),vp3(ma,mo,mr), &
+ stat=ier)
+ if( ier /= 0 ) then
+ call exit_mpi(myrank,'error allocation GAP model')
+ endif
+
+ ! the variables read are declared in the module
+ if(myrank == 0) call read_mantle_gapmodel()
+
+ ! master process broadcasts data to all processes
+ call MPI_BCAST( dep,mr+1,MPI_REAL,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST(dep1,mr1+1,MPI_REAL,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST( vp1,mr1+1,MPI_REAL,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST( vp3,ma*mo*mr,MPI_REAL,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST( nnr,1,MPI_INTEGER,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST( nr1,1,MPI_INTEGER,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST( no,1,MPI_INTEGER,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST( na,1,MPI_INTEGER,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST( dela,1,MPI_REAL,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST( delo,1,MPI_REAL,0,MPI_COMM_WORLD,ier)
+
+ end subroutine model_gapp2_broadcast
+
+!
+!--------------------------------------------------------------------------------------------------
+!
+
+ subroutine read_mantle_gapmodel()
+
+ use gapp2_mantle_model_constants
+
+ implicit none
+ include "constants.h"
+ integer i,ir,ia,io
+ character(len=150) GAPP2
+
+!...........................................input data
+
+ ! default model: 3dvpGAP_P2
+ call get_value_string(GAPP2, 'model.GAPP2', 'DATA/3dvpGAP_P2')
+
+ ! reads in GAP-P2 model from Obayashi
+ open(unit=10,file=GAPP2,status='old',action='read')
+
+ read(10,'(3i4,2f10.6)') no,na,nnr,dela,delo
+ read(10,'(34f8.2)') (dep(i),i=0,nnr)
+ read(10,*) nr1
+ read(10,'(67f8.2)') (dep1(i),i=0,nr1)
+ read(10,'(67f8.3)') (vp1(i),i=0,nr1)
+ do ir=1,nnr
+ do ia=1,na
+ read(10,'(256f7.3)') (vp3(ia,io,ir),io=1,no)
+ enddo
+ enddo
+ write(6,*) vp3(1,1,1),vp3(na,no,nnr)
+ close(10)
+
+ end subroutine read_mantle_gapmodel
+
+!
+!--------------------------------------------------------------------------------------------------
+!
+
+ subroutine mantle_gapmodel(radius,theta,phi,dvs,dvp,drho)
+
+ use gapp2_mantle_model_constants
+
+ implicit none
+ include "constants.h"
+ integer id,ia,io,icon
+ real d,dtheta,dphi
+
+ double precision radius,theta,phi,dvs,dvp,drho
+
+! factor to convert perturbations in shear speed to perturbations in density
+ double precision, parameter :: SCALE_VS = 1.40d0
+ double precision, parameter :: SCALE_RHO = 0.0d0
+
+ double precision, parameter :: R_EARTH_ = 6371.d0
+ double precision, parameter :: ZERO_ = 0.d0
+
+!.....................................
+
+ dvs = ZERO_
+ dvp = ZERO_
+ drho = ZERO_
+
+ ! increments in latitude/longitude (in rad)
+ dtheta = dela * PI / 180.0
+ dphi = delo * PI / 180.0
+
+ ! depth given in km
+ d=R_EARTH_-radius*R_EARTH_
+
+ call d2id(d,nnr,dep,id,icon)
+ if(icon.ne.0) then
+ write(6,*)icon
+ write(6,*) radius,theta,phi,dvp,dvs,drho
+ endif
+
+ ! latitude
+ if(theta.ge.PI) then
+ ia = na
+ else
+ ia = theta / dtheta + 1
+ endif
+ ! longitude
+ if(phi .lt. 0.0d0) phi = phi + 2.*PI
+ io=phi / dphi + 1
+ if(io.gt.no) io=io-no
+
+ ! velocity and density perturbations
+ dvp = vp3(ia,io,id)/100.d0
+ dvs = SCALE_VS*dvp
+ drho = SCALE_RHO*dvs
+
+ end subroutine mantle_gapmodel
+
+!
+!--------------------------------------------------------------------------------------------------
+!
+
+ subroutine d2id(d,mr,di,id,icon)
+!.................................................................
+! radial section index for a given depth d
+!.................................................................
+! d i depth(km)
+! mr i number of radial division
+! di i depth table
+! id o depth section index for d
+! shallow .... di(id-1) <= d < di(id) .... deep
+! icon o condition code
+! 0:normal, -99:above the surface, 99:below the cmb
+!.................................................................
+ integer i, mr, id, icon
+ real d,dmax,dmin
+ real di(0:mr)
+ icon=0
+ dmax=di(mr)
+ dmin=di(0)
+ if(d.gt.dmax) then
+ icon=99
+ else if(d.lt.dmin) then
+ icon=-99
+ else if(d.eq.dmax) then
+ id=mr+1
+ else
+ do i = 0, mr
+ if(d.lt.di(i)) then
+ id=i
+ goto 900
+ endif
+ enddo
+ end if
+900 continue
+
+!..................................................................
+ return
+
+ end subroutine d2id
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/model_gll.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/model_gll.f90 2010-10-11 09:41:46 UTC (rev 17258)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/model_gll.f90 2010-10-11 21:16:31 UTC (rev 17259)
@@ -45,6 +45,7 @@
include "constants.h"
! standard include of the MPI library
include 'mpif.h'
+ include "precision.h"
! GLL model_variables
type model_gll_variables
@@ -65,9 +66,11 @@
integer :: myrank
! local parameters
- double precision :: min_dvs,max_dvs,min_dvs_all,max_dvs_all,scaleval
+ double precision :: scaleval
+ real(kind=CUSTOM_REAL) :: min,max,min_all,max_all
integer :: ier
+ ! allocates arrays
! differs for isotropic model or transverse isotropic models
if( .not. TRANSVERSE_ISOTROPY ) then
! isotropic model
@@ -76,37 +79,132 @@
else
! transverse isotropic model
allocate( MGLL_V%vpv_new(NGLLX,NGLLY,NGLLZ,NSPEC(IREGION_CRUST_MANTLE)) )
+ allocate( MGLL_V%vph_new(NGLLX,NGLLY,NGLLZ,NSPEC(IREGION_CRUST_MANTLE)) )
allocate( MGLL_V%vsv_new(NGLLX,NGLLY,NGLLZ,NSPEC(IREGION_CRUST_MANTLE)) )
- allocate( MGLL_V%vph_new(NGLLX,NGLLY,NGLLZ,NSPEC(IREGION_CRUST_MANTLE)) )
allocate( MGLL_V%vsh_new(NGLLX,NGLLY,NGLLZ,NSPEC(IREGION_CRUST_MANTLE)) )
allocate( MGLL_V%eta_new(NGLLX,NGLLY,NGLLZ,NSPEC(IREGION_CRUST_MANTLE)) )
endif
- allocate( MGLL_V%rho_new(NGLLX,NGLLY,NGLLZ,NSPEC(IREGION_CRUST_MANTLE)) )
-
- ! non-dimensionalize scaling values
- ! (model velocities must be given as km/s)
- scaleval = dsqrt(PI*GRAV*RHOAV)
- MGLL_V%scale_velocity = 1000.0d0/(R_EARTH*scaleval)
- MGLL_V%scale_density = 1000.0d0/RHOAV
+ allocate( MGLL_V%rho_new(NGLLX,NGLLY,NGLLZ,NSPEC(IREGION_CRUST_MANTLE)) )
+ ! reads in model files for each process
call read_gll_model(myrank,MGLL_V,NSPEC)
! checks velocity range
if( .not. TRANSVERSE_ISOTROPY ) then
- max_dvs = maxval( MGLL_V%vs_new )
- min_dvs = minval( MGLL_V%vs_new )
+
+ ! isotropic model
+ if( myrank == 0 ) then
+ write(IMAIN,*)'model GLL: isotropic'
+ endif
+
+ ! Vs
+ max = maxval( MGLL_V%vs_new )
+ min = minval( MGLL_V%vs_new )
+ call mpi_reduce(max, max_all, 1, CUSTOM_MPI_TYPE, MPI_MAX, 0, MPI_COMM_WORLD,ier)
+ call mpi_reduce(min, min_all, 1, CUSTOM_MPI_TYPE, MPI_MIN, 0, MPI_COMM_WORLD,ier)
+ if( myrank == 0 ) then
+ write(IMAIN,*) ' vs new min/max: ',min_all,max_all
+ endif
+ ! Vp
+ max = maxval( MGLL_V%vp_new )
+ min = minval( MGLL_V%vp_new )
+ call mpi_reduce(max, max_all, 1, CUSTOM_MPI_TYPE, MPI_MAX, 0, MPI_COMM_WORLD,ier)
+ call mpi_reduce(min, min_all, 1, CUSTOM_MPI_TYPE, MPI_MIN, 0, MPI_COMM_WORLD,ier)
+ if( myrank == 0 ) then
+ write(IMAIN,*) ' vp new min/max: ',min_all,max_all
+ endif
+ ! density
+ max = maxval( MGLL_V%rho_new )
+ min = minval( MGLL_V%rho_new )
+ call mpi_reduce(max, max_all, 1, CUSTOM_MPI_TYPE, MPI_MAX, 0, MPI_COMM_WORLD,ier)
+ call mpi_reduce(min, min_all, 1, CUSTOM_MPI_TYPE, MPI_MIN, 0, MPI_COMM_WORLD,ier)
+ if( myrank == 0 ) then
+ write(IMAIN,*) ' rho new min/max: ',min_all,max_all
+ write(IMAIN,*)
+ endif
+
else
- max_dvs = maxval( MGLL_V%vsv_new + MGLL_V%vsh_new )
- min_dvs = minval( MGLL_V%vsv_new + MGLL_V%vsh_new )
+
+ ! transverse isotropic model
+ if( myrank == 0 ) then
+ write(IMAIN,*)'model GLL: transverse isotropic'
+ endif
+
+ ! Vsv
+ max = maxval( MGLL_V%vsv_new )
+ min = minval( MGLL_V%vsv_new )
+ call mpi_reduce(max, max_all, 1, CUSTOM_MPI_TYPE, MPI_MAX, 0, MPI_COMM_WORLD,ier)
+ call mpi_reduce(min, min_all, 1, CUSTOM_MPI_TYPE, MPI_MIN, 0, MPI_COMM_WORLD,ier)
+ if( myrank == 0 ) then
+ write(IMAIN,*) ' vsv new min/max: ',min_all,max_all
+ endif
+ ! Vsh
+ max = maxval( MGLL_V%vsh_new )
+ min = minval( MGLL_V%vsh_new )
+ call mpi_reduce(max, max_all, 1, CUSTOM_MPI_TYPE, MPI_MAX, 0, MPI_COMM_WORLD,ier)
+ call mpi_reduce(min, min_all, 1, CUSTOM_MPI_TYPE, MPI_MIN, 0, MPI_COMM_WORLD,ier)
+ if( myrank == 0 ) then
+ write(IMAIN,*) ' vsh new min/max: ',min_all,max_all
+ endif
+ ! Vpv
+ max = maxval( MGLL_V%vpv_new )
+ min = minval( MGLL_V%vpv_new )
+ call mpi_reduce(max, max_all, 1, CUSTOM_MPI_TYPE, MPI_MAX, 0, MPI_COMM_WORLD,ier)
+ call mpi_reduce(min, min_all, 1, CUSTOM_MPI_TYPE, MPI_MIN, 0, MPI_COMM_WORLD,ier)
+ if( myrank == 0 ) then
+ write(IMAIN,*) ' vpv new min/max: ',min_all,max_all
+ endif
+ ! Vph
+ max = maxval( MGLL_V%vph_new )
+ min = minval( MGLL_V%vph_new )
+ call mpi_reduce(max, max_all, 1, CUSTOM_MPI_TYPE, MPI_MAX, 0, MPI_COMM_WORLD,ier)
+ call mpi_reduce(min, min_all, 1, CUSTOM_MPI_TYPE, MPI_MIN, 0, MPI_COMM_WORLD,ier)
+ if( myrank == 0 ) then
+ write(IMAIN,*) ' vph new min/max: ',min_all,max_all
+ endif
+ ! density
+ max = maxval( MGLL_V%rho_new )
+ min = minval( MGLL_V%rho_new )
+ call mpi_reduce(max, max_all, 1, CUSTOM_MPI_TYPE, MPI_MAX, 0, MPI_COMM_WORLD,ier)
+ call mpi_reduce(min, min_all, 1, CUSTOM_MPI_TYPE, MPI_MIN, 0, MPI_COMM_WORLD,ier)
+ if( myrank == 0 ) then
+ write(IMAIN,*) ' rho new min/max: ',min_all,max_all
+ endif
+ ! eta
+ max = maxval( MGLL_V%eta_new )
+ min = minval( MGLL_V%eta_new )
+ call mpi_reduce(max, max_all, 1, CUSTOM_MPI_TYPE, MPI_MAX, 0, MPI_COMM_WORLD,ier)
+ call mpi_reduce(min, min_all, 1, CUSTOM_MPI_TYPE, MPI_MIN, 0, MPI_COMM_WORLD,ier)
+ if( myrank == 0 ) then
+ write(IMAIN,*) ' eta new min/max: ',min_all,max_all
+ write(IMAIN,*)
+ endif
+
endif
- call mpi_reduce(max_dvs, max_dvs_all, 1, MPI_DOUBLE_PRECISION, MPI_MAX, 0, MPI_COMM_WORLD,ier)
- call mpi_reduce(min_dvs, min_dvs_all, 1, MPI_DOUBLE_PRECISION, MPI_MIN, 0, MPI_COMM_WORLD,ier)
- if( myrank == 0 ) then
- write(IMAIN,*)'model GLL:'
- write(IMAIN,*) ' vs new min/max: ',min_dvs_all,max_dvs_all
- write(IMAIN,*)
- endif
+ ! non-dimensionalizes model values
+ ! (SPECFEM3D_GLOBE uses non-dimensionalized values in subsequent computations)
+ ! scaling values
+ ! (model velocities must be given as km/s)
+ scaleval = dsqrt(PI*GRAV*RHOAV)
+ MGLL_V%scale_velocity = 1000.0d0/(R_EARTH*scaleval)
+ MGLL_V%scale_density = 1000.0d0/RHOAV
+ if( .not. TRANSVERSE_ISOTROPY ) then
+ ! non-dimensionalize isotropic values
+ MGLL_V%vp_new = MGLL_V%vp_new * MGLL_V%scale_velocity
+ MGLL_V%vs_new = MGLL_V%vs_new * MGLL_V%scale_velocity
+ MGLL_V%rho_new = MGLL_V%rho_new * MGLL_V%scale_density
+ else
+ ! non-dimensionalize
+ ! transverse isotropic model
+ MGLL_V%vpv_new = MGLL_V%vpv_new * MGLL_V%scale_velocity
+ MGLL_V%vph_new = MGLL_V%vph_new * MGLL_V%scale_velocity
+ MGLL_V%vsv_new = MGLL_V%vsv_new * MGLL_V%scale_velocity
+ MGLL_V%vsh_new = MGLL_V%vsh_new * MGLL_V%scale_velocity
+ MGLL_V%rho_new = MGLL_V%rho_new * MGLL_V%scale_density
+ ! eta is already non-dimensional
+ endif
+
end subroutine model_gll_broadcast
!
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/model_ppm.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/model_ppm.f90 2010-10-11 09:41:46 UTC (rev 17258)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/model_ppm.f90 2010-10-11 21:16:31 UTC (rev 17259)
@@ -638,13 +638,14 @@
real(kind=CUSTOM_REAL) maxlat,maxlon,maxdepth
real(kind=CUSTOM_REAL) minlat,minlon,mindepth
real(kind=CUSTOM_REAL) radius,theta,phi,lat,lon,r_depth,margin_v,margin_h
+ real(kind=CUSTOM_REAL) dist_h,dist_v
!----------------------------------------------------------------------------------------------------
! smoothing parameters
logical,parameter:: GAUSS_SMOOTHING = .false. ! set to true to use this smoothing routine
sigma_h = 100.0 ! km, horizontal
- sigma_v = 50.0 ! km, vertical
+ sigma_v = 100.0 ! km, vertical
! check if smoothing applies
if( .not. GAUSS_SMOOTHING ) return
@@ -940,10 +941,19 @@
! --- only double loop over the elements in the search radius ---
do ispec2 = 1, nspec
- ! checks distance between centers of elements
- if ( (cx(ispec2)-cx0(ispec))**2 + (cy(ispec2)-cy0(ispec))** 2 > sigma_h3 &
- .or. (cz(ispec2)-cz0(ispec))** 2 > sigma_v3 ) cycle
+ ! calculates horizontal and vertical distance between two element centers
+
+ ! vector approximation
+ call get_distance_vec(dist_h,dist_v,cx0(ispec),cy0(ispec),cz0(ispec),&
+ cx(ispec2),cy(ispec2),cz(ispec2))
+
+ ! note: distances and sigmah, sigmav are normalized by R_EARTH
+
+ ! checks distance between centers of elements
+ if ( dist_h > sigma_h3 .or. abs(dist_v) > sigma_v3 ) cycle
+
+
factor(:,:,:) = jacobian(:,:,:,ispec2) * wgll_cube(:,:,:) ! integration factors
! loop over GLL points of the elements in current slice (ispec)
@@ -951,15 +961,19 @@
do j = 1, NGLLY
do i = 1, NGLLX
- x0 = xl(i,j,k,ispec)
- y0 = yl(i,j,k,ispec)
- z0 = zl(i,j,k,ispec) ! current point (i,j,k,ispec)
+ ! current point (i,j,k,ispec) location, cartesian coordinates
+ x0 = xl(i,j,k,ispec)
+ y0 = yl(i,j,k,ispec)
+ z0 = zl(i,j,k,ispec)
- ! gaussian function
- exp_val(:,:,:) = exp( -(xx(:,:,:,ispec2)-x0)**2/(2.0*sigma_h2) &
- -(yy(:,:,:,ispec2)-y0)**2/(2.0*sigma_h2) &
- -(zz(:,:,:,ispec2)-z0)**2/(2.0*sigma_v2) ) * factor(:,:,:)
+ ! calculate weights based on gaussian smoothing
+ call smoothing_weights_vec(x0,y0,z0,ispec2,sigma_h2,sigma_v2,exp_val,&
+ xx(:,:,:,ispec2),yy(:,:,:,ispec2),zz(:,:,:,ispec2))
+
+ ! adds GLL integration weights
+ exp_val(:,:,:) = exp_val(:,:,:) * factor(:,:,:)
+
! smoothed kernel values
tk_rho(i,j,k,ispec) = tk_rho(i,j,k,ispec) + sum(exp_val(:,:,:) * ks_rho(:,:,:,ispec2))
tk_kv(i,j,k,ispec) = tk_kv(i,j,k,ispec) + sum(exp_val(:,:,:) * ks_kv(:,:,:,ispec2))
@@ -996,7 +1010,7 @@
margin_v = sigma_v*R_EARTH/1000.0 ! in km
margin_h = sigma_h*R_EARTH/1000.0 * 180.0/(R_EARTH_KM*PI) ! in degree
- ! compute the smoothed kernel values
+ ! computes the smoothed values
do ispec = 1, nspec
! depth of given radius (in km)
@@ -1088,8 +1102,104 @@
end subroutine
+!
+! -----------------------------------------------------------------------------
+!
+ subroutine smoothing_weights_vec(x0,y0,z0,ispec2,sigma_h2,sigma_v2,exp_val,&
+ xx_elem,yy_elem,zz_elem)
+ implicit none
+ include "constants.h"
+
+ real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ),intent(out) :: exp_val
+ real(kind=CUSTOM_REAL),dimension(NGLLX,NGLLY,NGLLZ),intent(in) :: xx_elem, yy_elem, zz_elem
+ real(kind=CUSTOM_REAL),intent(in) :: x0,y0,z0,sigma_h2,sigma_v2
+ integer,intent(in) :: ispec2
+
+ ! local parameters
+ integer :: ii,jj,kk
+ real(kind=CUSTOM_REAL) :: dist_h,dist_v
+ !real(kind=CUSTOM_REAL) :: r0,r1,theta1
+
+ ! >>>>>
+ ! uniform sigma
+ ! just to avoid compiler warning
+ ii = ispec2
+ !exp_val(:,:,:) = exp( -((xx(:,:,:,ispec2)-x0)**2+(yy(:,:,:,ispec2)-y0)**2 &
+ ! +(zz(:,:,:,ispec2)-z0)**2 )/(2*sigma2) )*factor(:,:,:)
+
+ ! from basin code smoothing:
+ ! gaussian function
+ !exp_val(:,:,:) = exp( -(xx(:,:,:,ispec2)-x0)**2/(sigma_h2) &
+ ! -(yy(:,:,:,ispec2)-y0)**2/(sigma_h2) &
+ ! -(zz(:,:,:,ispec2)-z0)**2/(sigma_v2) ) * factor(:,:,:)
+ ! >>>>>
+
+ do kk = 1, NGLLZ
+ do jj = 1, NGLLY
+ do ii = 1, NGLLX
+ ! point in second slice
+
+ ! vector approximation:
+ call get_distance_vec(dist_h,dist_v,x0,y0,z0, &
+ xx_elem(ii,jj,kk),yy_elem(ii,jj,kk),zz_elem(ii,jj,kk))
+
+ ! gaussian function
+ exp_val(ii,jj,kk) = exp( - dist_h*dist_h/sigma_h2 &
+ - dist_v*dist_v/sigma_v2 ) ! * factor(ii,jj,kk)
+
+ enddo
+ enddo
+ enddo
+
+ end subroutine smoothing_weights_vec
+
+
!
+! -----------------------------------------------------------------------------
+!
+
+ subroutine get_distance_vec(dist_h,dist_v,x0,y0,z0,x1,y1,z1)
+
+! returns vector lengths as distances in radial and horizontal direction
+
+ implicit none
+ include "constants.h"
+
+ real(kind=CUSTOM_REAL),intent(out) :: dist_h,dist_v
+ real(kind=CUSTOM_REAL),intent(in) :: x0,y0,z0,x1,y1,z1
+
+ ! local parameters
+ real(kind=CUSTOM_REAL) :: r0,r1,alpha
+ real(kind=CUSTOM_REAL) :: vx,vy,vz
+
+ ! vertical distance
+ r0 = sqrt( x0*x0 + y0*y0 + z0*z0 ) ! length of first position vector
+ r1 = sqrt( x1*x1 + y1*y1 + z1*z1 )
+ dist_v = r1 - r0
+ ! only for flat earth with z in depth: dist_v = sqrt( (cz(ispec2)-cz0(ispec))** 2)
+
+ ! horizontal distance
+ ! length of vector from point 0 to point 1
+ ! assuming small earth curvature (since only for neighboring elements)
+
+ ! scales r0 to have same length as r1
+ alpha = r1 / r0
+ vx = alpha * x0
+ vy = alpha * y0
+ vz = alpha * z0
+
+ ! vector in horizontal between new r0 and r1
+ vx = x1 - vx
+ vy = y1 - vy
+ vz = z1 - vz
+
+ ! distance is vector length
+ dist_h = sqrt( vx*vx + vy*vy + vz*vz )
+
+ end subroutine get_distance_vec
+
+!
!--------------------------------------------------------------------------------------------------
!
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/model_s362ani.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/model_s362ani.f90 2010-10-11 09:41:46 UTC (rev 17258)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/model_s362ani.f90 2010-10-11 21:16:31 UTC (rev 17259)
@@ -1100,11 +1100,22 @@
implicit none
+ integer :: ncon,nver
+
+!daniel: original
+! integer icon(1)
+!
+! real(kind=4) verlat(1)
+! real(kind=4) verlon(1)
+! real(kind=4) verrad(1)
+! real(kind=4) con(1)
+
+!daniel: avoiding out-of-bounds errors
+ real(kind=4) verlat(nver)
+ real(kind=4) verlon(nver)
+ real(kind=4) verrad(nver)
+
integer icon(1)
-
- real(kind=4) verlat(1)
- real(kind=4) verlon(1)
- real(kind=4) verrad(1)
real(kind=4) con(1)
double precision dd
@@ -1114,8 +1125,7 @@
double precision ver8
double precision xla8
- integer :: ncon,iver,nver
-
+ integer :: iver
real(kind=4) :: xlat,xlon
xrad=3.14159265358979/180.d0
@@ -1243,9 +1253,15 @@
call ylm(y,x,lmax,ylmcof(1,ihpa),wk1,wk2,wk3)
else if(itypehpa(ihpa) == 2) then
numcof=numcoe(ihpa)
- call splcon(y,x,numcof,xlaspl(1,ihpa), &
- xlospl(1,ihpa),radspl(1,ihpa), &
+!daniel
+! call splcon(y,x,numcof,xlaspl(1,ihpa), &
+! xlospl(1,ihpa),radspl(1,ihpa), &
+! nconpt(ihpa),iconpt(1,ihpa),conpt(1,ihpa))
+
+ call splcon(y,x,numcof,xlaspl(1:numcof,ihpa), &
+ xlospl(1:numcof,ihpa),radspl(1:numcof,ihpa), &
nconpt(ihpa),iconpt(1,ihpa),conpt(1,ihpa))
+
else
write(6,"('problem 1')")
endif
@@ -1397,9 +1413,17 @@
call ylm(y,x,lmax,ylmcof(1,ihpa),wk1,wk2,wk3)
else if(itypehpa(ihpa) == 2) then
numcof=numcoe(ihpa)
- call splcon(y,x,numcof,xlaspl(1,ihpa), &
- xlospl(1,ihpa),radspl(1,ihpa), &
+
+!daniel
+! call splcon(y,x,numcof,xlaspl(1,ihpa), &
+! xlospl(1,ihpa),radspl(1,ihpa), &
+! nconpt(ihpa),iconpt(1,ihpa),conpt(1,ihpa))
+
+ call splcon(y,x,numcof,xlaspl(1:numcof,ihpa), &
+ xlospl(1:numcof,ihpa),radspl(1:numcof,ihpa), &
nconpt(ihpa),iconpt(1,ihpa),conpt(1,ihpa))
+
+
else
write(6,"('problem 1')")
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/prepare_timerun.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/prepare_timerun.f90 2010-10-11 09:41:46 UTC (rev 17258)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/prepare_timerun.f90 2010-10-11 21:16:31 UTC (rev 17259)
@@ -364,7 +364,7 @@
two_omega_earth = - 2.d0 * TWO_PI / (HOURS_PER_DAY * 3600.d0 * scale_t_inv)
endif
endif
-
+
A_array_rotation = 0._CUSTOM_REAL
B_array_rotation = 0._CUSTOM_REAL
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/save_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/save_arrays_solver.f90 2010-10-11 09:41:46 UTC (rev 17258)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/save_arrays_solver.f90 2010-10-11 21:16:31 UTC (rev 17259)
@@ -51,24 +51,24 @@
include "constants.h"
! model_attenuation_variables
- type model_attenuation_variables
- sequence
- double precision min_period, max_period
- double precision :: QT_c_source ! Source Frequency
- double precision, dimension(:), pointer :: Qtau_s ! tau_sigma
- double precision, dimension(:), pointer :: QrDisc ! Discontinutitues Defined
- double precision, dimension(:), pointer :: Qr ! Radius
- double precision, dimension(:), pointer :: Qmu ! Shear Attenuation
- double precision, dimension(:,:), pointer :: Qtau_e ! tau_epsilon
- double precision, dimension(:), pointer :: Qomsb, Qomsb2 ! one_minus_sum_beta
- double precision, dimension(:,:), pointer :: Qfc, Qfc2 ! factor_common
- double precision, dimension(:), pointer :: Qsf, Qsf2 ! scale_factor
- integer, dimension(:), pointer :: Qrmin ! Max and Mins of idoubling
- integer, dimension(:), pointer :: Qrmax ! Max and Mins of idoubling
- integer, dimension(:), pointer :: interval_Q ! Steps
- integer :: Qn ! Number of points
- integer dummy_pad ! padding 4 bytes to align the structure
- end type model_attenuation_variables
+! type model_attenuation_variables
+! sequence
+! double precision min_period, max_period
+! double precision :: QT_c_source ! Source Frequency
+! double precision, dimension(:), pointer :: Qtau_s ! tau_sigma
+! double precision, dimension(:), pointer :: QrDisc ! Discontinutitues Defined
+! double precision, dimension(:), pointer :: Qr ! Radius
+! double precision, dimension(:), pointer :: Qmu ! Shear Attenuation
+! double precision, dimension(:,:), pointer :: Qtau_e ! tau_epsilon
+! double precision, dimension(:), pointer :: Qomsb, Qomsb2 ! one_minus_sum_beta
+! double precision, dimension(:,:), pointer :: Qfc, Qfc2 ! factor_common
+! double precision, dimension(:), pointer :: Qsf, Qsf2 ! scale_factor
+! integer, dimension(:), pointer :: Qrmin ! Max and Mins of idoubling
+! integer, dimension(:), pointer :: Qrmax ! Max and Mins of idoubling
+! integer, dimension(:), pointer :: interval_Q ! Steps
+! integer :: Qn ! Number of points
+! integer dummy_pad ! padding 4 bytes to align the structure
+! end type model_attenuation_variables
logical ATTENUATION
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/save_header_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/save_header_file.f90 2010-10-11 09:41:46 UTC (rev 17258)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/save_header_file.f90 2010-10-11 21:16:31 UTC (rev 17259)
@@ -105,8 +105,8 @@
integer :: SIMULATION_TYPE
logical :: SAVE_FORWARD,MOVIE_VOLUME
-
-
+
+
! copy number of elements and points in an include file for the solver
call get_value_string(HEADER_FILE, 'solver.HEADER_FILE', 'OUTPUT_FILES/values_from_mesher.h')
open(unit=IOUT,file=HEADER_FILE,status='unknown')
@@ -494,9 +494,10 @@
write(IOUT,*) 'logical, parameter :: USE_DEVILLE_PRODUCTS_VAL = .false.'
endif
- ! backward/reconstruction of forward wavefield:
- ! can only mimic attenuation effects on velocity at this point, since no full wavefield snapshots are stored
+ ! backward/reconstruction of forward wavefield:
+ ! can only mimic attenuation effects on velocity at this point, since no full wavefield snapshots are stored
if((SIMULATION_TYPE == 1 .and. SAVE_FORWARD) .or. SIMULATION_TYPE == 3) then
+
! attenuation mimic:
! mimicking effect of attenuation on apparent velocities, not amplitudes. that is,
! phase shifts should be correctly accounted for, but amplitudes will differ in adjoint simulations
@@ -505,9 +506,11 @@
else
write(IOUT,*) 'logical, parameter :: USE_ATTENUATION_MIMIC = .false.'
endif
-
+
else
- write(IOUT,*) 'logical, parameter :: USE_ATTENUATION_MIMIC = .false.'
+
+ ! calculates full attenuation (phase & amplitude effects) if used
+ write(IOUT,*) 'logical, parameter :: USE_ATTENUATION_MIMIC = .false.'
endif
! attenuation and/or adjoint simulations
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90 2010-10-11 09:41:46 UTC (rev 17258)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90 2010-10-11 21:16:31 UTC (rev 17259)
@@ -1156,7 +1156,7 @@
! define local to global receiver numbering mapping
allocate(number_receiver_global(nrec_local))
! define and store Lagrange interpolators at all the receivers
- if (SIMULATION_TYPE == 2 ) then
+ if (SIMULATION_TYPE == 2) then
nadj_hprec_local = nrec_local
else
nadj_hprec_local = 1
@@ -1589,10 +1589,10 @@
normal_y_noise(:) = 0._CUSTOM_REAL
normal_z_noise(:) = 0._CUSTOM_REAL
mask_noise(:) = 0._CUSTOM_REAL
-
+
call read_parameters_noise(myrank,nrec,NSTEP,nmovie_points, &
islice_selected_rec,xi_receiver,eta_receiver,gamma_receiver,nu, &
- noise_sourcearray,xigll,yigll,zigll,NSPEC2D_TOP(IREGION_CRUST_MANTLE), &
+ noise_sourcearray,xigll,yigll,zigll,NSPEC2D_TOP(IREGION_CRUST_MANTLE), &
NIT, ibool_crust_mantle, ibelm_top_crust_mantle, &
xstore_crust_mantle,ystore_crust_mantle,zstore_crust_mantle, &
irec_master_noise,normal_x_noise,normal_y_noise,normal_z_noise,mask_noise)
@@ -2122,8 +2122,6 @@
veloc_outer_core(i+3) = veloc_outer_core(i+3) + deltatover2*accel_outer_core(i+3)
enddo
-
-
if (SIMULATION_TYPE == 3) then
call assemble_MPI_scalar(myrank,b_accel_outer_core,NGLOB_OUTER_CORE, &
iproc_xi,iproc_eta,ichunk,addressing, &
@@ -2165,7 +2163,6 @@
endif
-
! ****************************************************
! big loop over all spectral elements in the solid
! ****************************************************
@@ -2226,7 +2223,7 @@
size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5) )
endif
- if (SIMULATION_TYPE == 3) then
+ if (SIMULATION_TYPE == 3 ) then
! for anisotropy and gravity, x y and z contain r theta and phi
if( USE_DEVILLE_PRODUCTS_VAL ) then
call compute_forces_crust_mantle_Dev(minus_gravity_table,density_table,minus_deriv_gravity_table, &
@@ -2282,7 +2279,6 @@
endif
endif
-
! Stacey
if(NCHUNKS_VAL /= 6 .and. ABSORBING_CONDITIONS) then
call compute_stacey_crust_mantle(ichunk,SIMULATION_TYPE, &
@@ -2433,7 +2429,7 @@
normal_x_noise,normal_y_noise,normal_z_noise,mask_noise, &
store_val_ux,store_val_uy,store_val_uz, &
ibelm_top_crust_mantle,ibool_crust_mantle,NSPEC2D_TOP(IREGION_CRUST_MANTLE), &
- NIT,NSTEP-it+1,LOCAL_PATH,jacobian2D_top_crust_mantle,wgllwgll_xy)
+ NIT,NSTEP-it+1,LOCAL_PATH,jacobian2D_top_crust_mantle,wgllwgll_xy)
! be careful, since ensemble forward sources are reversals of generating wavefield "eta"
! hence the "NSTEP-it+1", i.e., start reading from the last timestep
! note the ensemble forward sources are generally distributed on the surface of the earth
@@ -2448,7 +2444,7 @@
normal_x_noise,normal_y_noise,normal_z_noise,mask_noise, &
store_val_ux,store_val_uy,store_val_uz, &
ibelm_top_crust_mantle,ibool_crust_mantle,NSPEC2D_TOP(IREGION_CRUST_MANTLE), &
- NIT,it,LOCAL_PATH,jacobian2D_top_crust_mantle,wgllwgll_xy)
+ NIT,it,LOCAL_PATH,jacobian2D_top_crust_mantle,wgllwgll_xy)
endif
!>YANGL
@@ -2570,9 +2566,6 @@
accel_crust_mantle(3,i+3) = accel_crust_mantle(3,i+3)*rmass_crust_mantle(i+3)
enddo
-
-
-
if (SIMULATION_TYPE == 3) then
! assemble all the contributions between slices using MPI
@@ -2662,10 +2655,8 @@
b_accel_crust_mantle(3,i+3) = b_accel_crust_mantle(3,i+3)*rmass_crust_mantle(i+3)
enddo
-
endif
-
! couples ocean with crust mantle
if(OCEANS_VAL) &
call compute_coupling_ocean(accel_crust_mantle,b_accel_crust_mantle, &
@@ -2743,8 +2734,6 @@
veloc_inner_core(:,i+2) = veloc_inner_core(:,i+2) + deltatover2*accel_inner_core(:,i+2)
enddo
-
-
if (SIMULATION_TYPE == 3) then
! way 1:
! do i=1,NGLOB_CRUST_MANTLE
@@ -2834,7 +2823,7 @@
seismo_current,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
seismograms)
- else if (SIMULATION_TYPE == 2 ) then
+ else if (SIMULATION_TYPE == 2) then
call compute_seismograms_adjoint(NSOURCES,nrec_local,displ_crust_mantle, &
eps_trace_over_3_crust_mantle,epsilondev_crust_mantle, &
nu_source,Mxx,Myy,Mzz,Mxy,Mxz,Myz, &
@@ -2850,7 +2839,7 @@
ibool_crust_mantle,ispec_selected_source,number_receiver_global, &
NSTEP,it,nit_written)
- else if (SIMULATION_TYPE == 3 ) then
+ else if (SIMULATION_TYPE == 3) then
call compute_seismograms_backward(nrec_local,nrec,b_displ_crust_mantle, &
nu,hxir_store,hetar_store,hgammar_store, &
scale_displ,ibool_crust_mantle, &
@@ -3133,11 +3122,11 @@
k_bot,ibelm_top_inner_core,normal_bottom_outer_core, &
icb_kl_bot,fluid_solid_boundary,NSPEC2D_ICB)
- icb_kl = icb_kl + (icb_kl_top - icb_kl_bot) * deltat
- endif
+ icb_kl = icb_kl + (icb_kl_top - icb_kl_bot) * deltat
+ endif
- endif ! end computing kernels
-
+ endif ! end computing kernels
+
!
!-------------------------------------------------------------------------------------------------
!-------------------------------------------------------------------------------------------------
@@ -3323,7 +3312,7 @@
rhostore_crust_mantle,muvstore_crust_mantle, &
kappavstore_crust_mantle,ibool_crust_mantle, &
kappahstore_crust_mantle,muhstore_crust_mantle, &
- eta_anisostore_crust_mantle,idoubling_crust_mantle, &
+ eta_anisostore_crust_mantle,idoubling_crust_mantle, &
LOCAL_PATH)
!<YANGL
More information about the CIG-COMMITS
mailing list