[cig-commits] r11780 - seismo/3D/SPECFEM3D_GLOBE/trunk
vala at geodynamics.org
vala at geodynamics.org
Tue Apr 8 22:30:07 PDT 2008
Author: vala
Date: 2008-04-08 22:30:06 -0700 (Tue, 08 Apr 2008)
New Revision: 11780
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/attenuation_model.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/get_model.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/model_ref.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/write_AVS_DX_global_chunks_data.f90
Log:
Fixed a problem that caused the mantle not to be filled in properly for model_ref when the crust was thinner than the PREM crust. Also set
eta_aniso to always be 1.0 within the crust, and for model s362iso.
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/attenuation_model.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/attenuation_model.f90 2008-04-09 05:16:07 UTC (rev 11779)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/attenuation_model.f90 2008-04-09 05:30:06 UTC (rev 11780)
@@ -209,7 +209,7 @@
call define_model_1066a(.FALSE., M1066a_V)
AM_V%Qn = NR_1066A
else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_REF) then
- call define_model_ref(.FALSE., Mref_V)
+ call define_model_ref(Mref_V)
AM_V%Qn = NR_REF
else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_JP1D) then
AM_V%Qn = 12
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90 2008-04-09 05:16:07 UTC (rev 11779)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90 2008-04-09 05:30:06 UTC (rev 11780)
@@ -697,7 +697,7 @@
elseif(REFERENCE_1D_MODEL == REFERENCE_MODEL_AK135) then
call define_model_ak135(CRUSTAL, Mak135_V)
elseif(REFERENCE_1D_MODEL == REFERENCE_MODEL_REF) then
- call define_model_ref(CRUSTAL, Mref_V)
+ call define_model_ref(Mref_V)
elseif(REFERENCE_1D_MODEL == REFERENCE_MODEL_SEA1D) then
call define_model_sea1d(CRUSTAL, SEA1DM_V)
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/get_model.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/get_model.f90 2008-04-09 05:16:07 UTC (rev 11779)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/get_model.f90 2008-04-09 05:30:06 UTC (rev 11780)
@@ -410,7 +410,7 @@
R600,R670,R220,R771,R400,R80,RMOHO,RMIDDLE_CRUST,ROCEAN)
else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_REF) then
- call model_ref(r_prem,rho,vpv,vph,vsv,vsh,eta_aniso,Qkappa,Qmu,iregion_code,Mref_V)
+ call model_ref(r_prem,rho,vpv,vph,vsv,vsh,eta_aniso,Qkappa,Qmu,iregion_code,CRUSTAL,Mref_V)
else
stop 'unknown 1D transversely isotropic reference Earth model in get_model'
@@ -434,7 +434,7 @@
call model_ak135(r_prem,rho,vp,vs,Qkappa,Qmu,iregion_code,Mak135_V)
else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_REF) then
- call model_ref(r_prem,rho,vpv,vph,vsv,vsh,eta_aniso,Qkappa,Qmu,iregion_code,Mref_V)
+ call model_ref(r_prem,rho,vpv,vph,vsv,vsh,eta_aniso,Qkappa,Qmu,iregion_code,CRUSTAL,Mref_V)
if(.not. ISOTROPIC_3D_MANTLE) then
vp = sqrt(((8.d0+4.d0*eta_aniso)*vph*vph + 3.d0*vpv*vpv + (8.d0 - 8.d0*eta_aniso)*vsv*vsv)/15.d0)
vs = sqrt(((1.d0-2.d0*eta_aniso)*vph*vph + vpv*vpv + 5.d0*vsh*vsh + (6.d0+4.d0*eta_aniso)*vsv*vsv)/15.d0)
@@ -564,6 +564,7 @@
vph=vp
vsv=vs
vsh=vs
+ eta_aniso=1.0d0
endif
else
stop 'unknown 3D Earth model in get_model'
@@ -660,6 +661,7 @@
vph=vp
vsv=vs
vsh=vs
+ eta_aniso=1.0d0
endif
else
stop 'unknown 3D Earth model in get_model'
@@ -775,6 +777,7 @@
vsv=vsc
vsh=vsc
rho=rhoc
+ eta_aniso=1.0d0
if(ANISOTROPIC_3D_MANTLE .and. iregion_code == IREGION_CRUST_MANTLE) then
c11 = rho*vpv*vpv
c12 = rho*(vpv*vpv-2.*vsv*vsv)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/model_ref.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/model_ref.f90 2008-04-09 05:16:07 UTC (rev 11779)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/model_ref.f90 2008-04-09 05:30:06 UTC (rev 11780)
@@ -26,7 +26,7 @@
!=====================================================================
- subroutine model_ref(x,rho,vpv,vph,vsv,vsh,eta,Qkappa,Qmu,iregion_code,Mref_V)
+ subroutine model_ref(x,rho,vpv,vph,vsv,vsh,eta,Qkappa,Qmu,iregion_code,CRUSTAL,Mref_V)
implicit none
@@ -69,6 +69,7 @@
integer i
double precision r,frac,scaleval
+ logical CRUSTAL
! compute real physical radius in meters
r = x * R_EARTH
@@ -85,7 +86,9 @@
if(iregion_code == IREGION_OUTER_CORE .and. i > 358) i = 358
if(iregion_code == IREGION_CRUST_MANTLE .and. i < 360) i = 360
+ if(CRUSTAL .and. i > 717) i = 717
+
if(i == 1) then
rho = Mref_V%density_ref(i)
vpv = Mref_V%vpv_ref(i)
@@ -133,7 +136,7 @@
!-------------------
- subroutine define_model_ref(USE_EXTERNAL_CRUSTAL_MODEL,Mref_V)
+ subroutine define_model_ref(Mref_V)
implicit none
include "constants.h"
@@ -155,10 +158,7 @@
type (model_ref_variables) Mref_V
! model_ref_variables
- logical USE_EXTERNAL_CRUSTAL_MODEL
- integer i
-
! define the 1D REF model of Kustowski et al. (2007)
Mref_V%radius_ref( 1 : 30 ) = (/ &
@@ -7369,19 +7369,6 @@
Mref_V%vsh_ref(718:750) = Mref_V%vsh_ref(717)
endif
-! strip the crust and replace it by mantle if we use an external crustal model
- if(USE_EXTERNAL_CRUSTAL_MODEL) then
- do i=NR_REF-10,NR_REF
- Mref_V%density_ref(i) = Mref_V%density_ref(NR_REF-11)
- Mref_V%vpv_ref(i) = Mref_V%vpv_ref(NR_REF-11)
- Mref_V%vph_ref(i) = Mref_V%vph_ref(NR_REF-11)
- Mref_V%vsv_ref(i) = Mref_V%vsv_ref(NR_REF-11)
- Mref_V%vsh_ref(i) = Mref_V%vsh_ref(NR_REF-11)
- Mref_V%eta_ref(i) = Mref_V%eta_ref(NR_REF-11)
- Mref_V%Qkappa_ref(i) = Mref_V%Qkappa_ref(NR_REF-11)
- Mref_V%Qmu_ref(i) = Mref_V%Qmu_ref(NR_REF-11)
- enddo
- endif
end subroutine define_model_ref
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/write_AVS_DX_global_chunks_data.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/write_AVS_DX_global_chunks_data.f90 2008-04-09 05:16:07 UTC (rev 11779)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/write_AVS_DX_global_chunks_data.f90 2008-04-09 05:30:06 UTC (rev 11780)
@@ -629,7 +629,7 @@
call model_ak135(r,rho,vp,vs,Qkappa,Qmu,idoubling(ispec),Mak135_V)
else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_REF) then
- call model_ref(r,rho,vpv,vph,vsv,vsh,eta_aniso,Qkappa,Qmu,idoubling(ispec),Mref_V)
+ call model_ref(r,rho,vpv,vph,vsv,vsh,eta_aniso,Qkappa,Qmu,idoubling(ispec),CRUSTAL,Mref_V)
vp = vpv
vs = vsv
else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_JP1D) then
More information about the cig-commits
mailing list