[cig-commits] r16307 - seismo/3D/SPECFEM3D_GLOBE/trunk
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Mon Feb 22 11:44:40 PST 2010
Author: danielpeter
Date: 2010-02-22 11:44:40 -0800 (Mon, 22 Feb 2010)
New Revision: 16307
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/compute_element_properties.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/get_model.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/get_model_parameters.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D_models.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/moho_stretching.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90
Log:
adding feature to use a 1D crust together with 3D mantle models: simply append _1Dcrust to the model name in the DATA/Par_file, e.g. s20rts_1Dcrust, in order to have the 1D crustal model from the corresponding 1D reference model; making the distinction between elements in crust and mantle more explicit such that there is a sharp transition on all GLL points within purely crustal elements to mantle elements when assigning the velocities in case of stretched element layers
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/compute_element_properties.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/compute_element_properties.f90 2010-02-22 18:52:21 UTC (rev 16306)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/compute_element_properties.f90 2010-02-22 19:44:40 UTC (rev 16307)
@@ -109,11 +109,12 @@
double precision:: zigll(NGLLZ)
! Parameter used to decide whether this element is in the crust or not
- logical:: elem_in_crust
+ logical:: elem_in_crust,elem_in_mantle
! add topography of the Moho *before* adding the 3D crustal velocity model so that the streched
! mesh gets assigned the right model values
elem_in_crust = .false.
+ elem_in_mantle = .false.
if( iregion_code == IREGION_CRUST_MANTLE ) then
if( CRUSTAL .and. CASE_3D ) then
if( idoubling(ispec) == IFLAG_CRUST &
@@ -121,7 +122,7 @@
.or. idoubling(ispec) == IFLAG_80_MOHO ) then
! Stretch mesh to honor smoothed moho thickness from crust2.0
call moho_stretching_honor_crust(myrank,xelm,yelm,zelm,RMOHO_FICTITIOUS_IN_MESHER,&
- R220,RMIDDLE_CRUST,elem_in_crust)
+ R220,RMIDDLE_CRUST,elem_in_crust,elem_in_mantle)
endif
endif
endif
@@ -144,7 +145,7 @@
R771,R400,R120,R80,RMIDDLE_CRUST,ROCEAN, &
tau_s,tau_e_store,Qmu_store,T_c_source, &
size(tau_e_store,2),size(tau_e_store,3),size(tau_e_store,4),size(tau_e_store,5), &
- ABSORBING_CONDITIONS,elem_in_crust )
+ ABSORBING_CONDITIONS,elem_in_crust,elem_in_mantle)
! either use GLL points or anchor points to capture TOPOGRAPHY and ELLIPTICITY
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/get_model.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/get_model.f90 2010-02-22 18:52:21 UTC (rev 16306)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/get_model.f90 2010-02-22 19:44:40 UTC (rev 16307)
@@ -36,7 +36,7 @@
rmin,rmax,RCMB,RICB,R670,RMOHO,RTOPDDOUBLEPRIME,R600,R220, &
R771,R400,R120,R80,RMIDDLE_CRUST,ROCEAN, &
tau_s,tau_e_store,Qmu_store,T_c_source,vx,vy,vz,vnspec, &
- ABSORBING_CONDITIONS,elem_in_crust)
+ ABSORBING_CONDITIONS,elem_in_crust,elem_in_mantle)
use meshfem3D_models_par
@@ -75,7 +75,7 @@
double precision T_c_source
logical ABSORBING_CONDITIONS
- logical elem_in_crust
+ logical elem_in_crust,elem_in_mantle
! local parameters
double precision xmesh,ymesh,zmesh
@@ -163,13 +163,15 @@
c33,c34,c35,c36,c44,c45,c46,c55,c56,c66)
! gets the 3-D crustal model
- if(CRUSTAL) &
- call meshfem3D_models_get3Dcrust_val(iregion_code,xmesh,ymesh,zmesh,r, &
+ if( CRUSTAL ) then
+ if( .not. elem_in_mantle) &
+ call meshfem3D_models_get3Dcrust_val(iregion_code,xmesh,ymesh,zmesh,r, &
vpv,vph,vsv,vsh,rho,eta_aniso,dvp, &
c11,c12,c13,c14,c15,c16,c22,c23,c24,c25, &
c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66, &
elem_in_crust,moho)
-
+ endif
+
! overwrites with tomographic model values (from iteration step) here, given at all GLL points
call meshfem3D_models_impose_val(vpv,vph,vsv,vsh,rho,dvp,eta_aniso,&
myrank,iregion_code,ispec,i,j,k)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/get_model_parameters.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/get_model_parameters.f90 2010-02-22 18:52:21 UTC (rev 16306)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/get_model_parameters.f90 2010-02-22 19:44:40 UTC (rev 16307)
@@ -100,9 +100,12 @@
ISOTROPIC_3D_MANTLE,ONE_CRUST,TRANSVERSE_ISOTROPY
logical OCEANS,TOPOGRAPHY
+ ! local parameters
character(len=4) ending
- character(len=150) MODEL_ROOT
+ character(len=8) ending_1Dcrust
+ character(len=150) MODEL_ROOT
+ logical :: impose_1Dcrust
! defaults:
!
@@ -135,6 +138,17 @@
MODEL_ROOT = MODEL
endif
+ ! 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( ending_1Dcrust == '_1Dcrust' ) then
+ impose_1Dcrust = .true.
+ ! in case it has an ending for the inner core, remove it from the name
+ MODEL_ROOT = MODEL_ROOT(1: len_trim(MODEL)-8)
+ endif
+
+
!---
!
! ADD YOUR MODEL HERE
@@ -379,6 +393,19 @@
TOPOGRAPHY = .false.
endif
+ ! additional option for 3D mantle models:
+ ! this takes crust from reference 1D model rather than a 3D crust;
+ if( impose_1Dcrust ) then
+ ! no 3D crust
+ CRUSTAL = .false.
+ ! no crustal moho stretching
+ CASE_3D = .false.
+ ! mesh honors the 1D moho depth
+ HONOR_1D_SPHERICAL_MOHO = .true.
+ ! 2 element layers in top crust region rather than just one
+ ONE_CRUST = .false.
+ endif
+
! checks flag consistency for crust
if( HONOR_1D_SPHERICAL_MOHO .and. CRUSTAL ) &
stop 'honor 1D spherical moho excludes having 3D crustal structure'
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D_models.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D_models.f90 2010-02-22 18:52:21 UTC (rev 16306)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D_models.f90 2010-02-22 19:44:40 UTC (rev 16307)
@@ -727,7 +727,7 @@
double precision :: r_used,r_dummy,theta,phi
double precision :: dvs,drho,vp,vs
real(kind=4) :: xcolat,xlon,xrad,dvpv,dvph,dvsv,dvsh
- logical :: found_crust
+ logical :: found_crust,suppress_mantle_extension
! initializes perturbation values
dvs = ZERO
@@ -738,7 +738,8 @@
dvsv = 0.
dvsh = 0.
r_used = ZERO
-
+ suppress_mantle_extension = .false.
+
! gets point's theta/phi
call xyz_2_rthetaphi_dble(xmesh,ymesh,zmesh,r_dummy,theta,phi)
call reduce(theta,phi)
@@ -748,6 +749,11 @@
! ADD YOUR MODEL HERE
!
!---
+
+ ! sets flag when mantle should not be extended to surface
+ if(r_prem >= RMOHO/R_EARTH .and. .not. CRUSTAL) then
+ suppress_mantle_extension = .true.
+ endif
! gets parameters for isotropic 3D mantle model
!
@@ -755,18 +761,20 @@
! like kappav,kappah,muv,muh and eta_aniso are used for these simulations
!
! note: in general, models here make use of perturbation values with respect to their
- ! corresponding 1-D reference models
- if( ISOTROPIC_3D_MANTLE .and. r_prem > RCMB/R_EARTH ) then
+ ! corresponding 1-D reference models
+ if( ISOTROPIC_3D_MANTLE .and. r_prem > RCMB/R_EARTH .and. .not. suppress_mantle_extension) then
! extend 3-D mantle model above the Moho to the surface before adding the crust
if(r_prem > RCMB/R_EARTH .and. r_prem < RMOHO/R_EARTH) then
! GLL point is in mantle region, takes exact location
r_used = r
else ! else if(r_prem >= RMOHO/R_EARTH) then
- ! GLL point is above moho
- ! takes radius slightly below moho radius, this will then "extend the mantle up to the surface";
- ! crustal values will be superimposed later on
- r_used = 0.999999d0*RMOHO/R_EARTH
+ if( CRUSTAL ) then
+ ! GLL point is above moho
+ ! takes radius slightly below moho radius, this will then "extend the mantle up to the surface";
+ ! crustal values will be superimposed later on
+ r_used = 0.999999d0*RMOHO/R_EARTH
+ endif
endif
! gets model parameters
@@ -866,7 +874,7 @@
endif ! ISOTROPIC_3D_MANTLE
! heterogen model
- if( HETEROGEN_3D_MANTLE ) then
+ if( HETEROGEN_3D_MANTLE .and. .not. suppress_mantle_extension ) then
call xyz_2_rthetaphi_dble(xmesh,ymesh,zmesh,r_used,theta,phi)
call reduce(theta,phi)
call model_heterogen_mantle(r_used,theta,phi,dvs,dvp,drho,HMM)
@@ -884,20 +892,21 @@
if(ANISOTROPIC_3D_MANTLE .and. iregion_code == IREGION_CRUST_MANTLE) then
! anisotropic model between the Moho and 670 km (change to CMB if desired)
- if( r_prem > R670/R_EARTH) then
+ if( r_prem > R670/R_EARTH .and. .not. suppress_mantle_extension ) then
! extend 3-D mantle model above the Moho to the surface before adding the crust
if( r_prem < RMOHO/R_EARTH) then
r_used = r_prem
else
- ! fills 3-D mantle model above the Moho with the values at moho depth
- r_used = RMOHO/R_EARTH
+ if( CRUSTAL ) then
+ ! fills 3-D mantle model above the Moho with the values at moho depth
+ r_used = RMOHO/R_EARTH
+ endif
endif
call model_aniso_mantle(r_used,theta,phi,rho,c11,c12,c13,c14,c15,c16, &
c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66,AMM_V)
else
-
! fills the rest of the mantle with the isotropic model
c11 = rho*vpv*vpv
c12 = rho*(vpv*vpv-2.*vsv*vsv)
@@ -953,7 +962,7 @@
c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66, &
elem_in_crust,moho)
-! returns velocities and density for points in crustal region
+! returns velocities and density for points in 3D crustal region
use meshfem3D_models_par
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/moho_stretching.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/moho_stretching.f90 2010-02-22 18:52:21 UTC (rev 16306)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/moho_stretching.f90 2010-02-22 19:44:40 UTC (rev 16307)
@@ -27,7 +27,7 @@
subroutine moho_stretching_honor_crust(myrank,xelm,yelm,zelm,RMOHO_FICTITIOUS_IN_MESHER,&
- R220,RMIDDLE_CRUST,elem_in_crust)
+ R220,RMIDDLE_CRUST,elem_in_crust,elem_in_mantle)
! stretching the moho according to the crust 2.0
! input: myrank, xelm, yelm, zelm, RMOHO_FICTITIOUS_IN_MESHER R220,RMIDDLE_CRUST, CM_V
@@ -43,10 +43,10 @@
double precision R220,RMIDDLE_CRUST
double precision RMOHO_FICTITIOUS_IN_MESHER
integer :: myrank
- logical :: elem_in_crust
+ logical :: elem_in_crust,elem_in_mantle
! local parameters
- integer:: ia,count
+ integer:: ia,count_crust,count_mantle
double precision:: r,theta,phi,lat,lon
double precision:: vpc,vsc,rhoc,moho,elevation,gamma
logical:: found_crust
@@ -62,7 +62,8 @@
R_middlecrust = RMIDDLE_CRUST/R_EARTH
! loops over element's anchor points
- count = 0
+ count_crust = 0
+ count_mantle = 0
do ia = 1,NGNOD
x = xelm(ia)
y = yelm(ia)
@@ -166,17 +167,28 @@
end if
! counts corners in above moho
+ ! note: uses a small tolerance
if ( r >= 0.9999d0*moho ) then
- count = count + 1
- end if
+ count_crust = count_crust + 1
+ endif
+ ! counts corners below moho
+ ! again within a small tolerance
+ if ( r <= 1.0001d0*moho ) then
+ count_mantle = count_mantle + 1
+ endif
end do
! sets flag when all corners are above moho
- if ( count == NGNOD) then
+ if( count_crust == NGNOD) then
elem_in_crust = .true.
end if
+ ! sets flag when all corners are below moho
+ if( count_mantle == NGNOD) then
+ elem_in_mantle = .true.
+ end if
+ ! small stretch check: stretching should affect only points above R220
if( r*R_EARTH < R220 ) then
print*,'error moho stretching: ',r*R_EARTH,R220,moho*R_EARTH
call exit_mpi(myrank,'incorrect moho stretching')
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90 2010-02-22 18:52:21 UTC (rev 16306)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90 2010-02-22 19:44:40 UTC (rev 16307)
@@ -1910,6 +1910,19 @@
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
ibool_outer_core)
+! ! uses deville et al. (2002) routine
+! call compute_forces_outer_core_Dev(time,deltat,two_omega_earth, &
+! A_array_rotation,B_array_rotation,d_ln_density_dr_table, &
+! minus_rho_g_over_kappa_fluid,displ_outer_core,accel_outer_core,div_displ_outer_core, &
+! xstore_outer_core,ystore_outer_core,zstore_outer_core, &
+! xix_outer_core,xiy_outer_core,xiz_outer_core, &
+! etax_outer_core,etay_outer_core,etaz_outer_core, &
+! gammax_outer_core,gammay_outer_core,gammaz_outer_core, &
+! hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
+! wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+! ibool_outer_core)
+
+
if (SIMULATION_TYPE == 3) then
call compute_forces_outer_core(time,b_deltat,b_two_omega_earth, &
b_A_array_rotation,b_B_array_rotation,d_ln_density_dr_table, &
@@ -1921,6 +1934,20 @@
hprime_xx,hprime_yy,hprime_zz,hprimewgll_xx,hprimewgll_yy,hprimewgll_zz, &
wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
ibool_outer_core)
+
+! ! uses deville et al. (2002) routine
+! call compute_forces_outer_core_Dev(time,b_deltat,b_two_omega_earth, &
+! b_A_array_rotation,b_B_array_rotation,d_ln_density_dr_table, &
+! minus_rho_g_over_kappa_fluid,b_displ_outer_core,b_accel_outer_core,b_div_displ_outer_core, &
+! xstore_outer_core,ystore_outer_core,zstore_outer_core, &
+! xix_outer_core,xiy_outer_core,xiz_outer_core, &
+! etax_outer_core,etay_outer_core,etaz_outer_core, &
+! gammax_outer_core,gammay_outer_core,gammaz_outer_core, &
+! hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT, &
+! wgllwgll_xy,wgllwgll_xz,wgllwgll_yz,wgll_cube, &
+! ibool_outer_core)
+
+
endif
! Stacey absorbing boundaries
More information about the CIG-COMMITS
mailing list