[cig-commits] r16878 - in seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN: . DATA/hauksson_model DATA/lin_model DATA/lin_model/lin_2007
carltape at geodynamics.org
carltape at geodynamics.org
Thu Jun 3 07:12:05 PDT 2010
Author: carltape
Date: 2010-06-03 07:12:04 -0700 (Thu, 03 Jun 2010)
New Revision: 16878
Added:
seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/hauksson_model/README_hauksson_model
seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/hauksson_model/hauksson_model.f90
seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/lin_model/README_lin_model
Removed:
seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/lin_model/README
Modified:
seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/lin_model/hauksson_model.f90
seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/lin_model/hauksson_tomo_model.m
seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/lin_model/lin_2007/README
seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/constants.h.in
seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/create_regions_mesh.f90
seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/hauksson_model.f90
Log:
Replaced Hauksson (2000) model with Lin-Shearer-Hauksson-Thurber (2007) model. Hauksson (2000) model can be implemented with minor modifications outlined in DATA/hauksson_model/README_hauksson_model.
Added: seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/hauksson_model/README_hauksson_model
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/hauksson_model/README_hauksson_model (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/hauksson_model/README_hauksson_model 2010-06-03 14:12:04 UTC (rev 16878)
@@ -0,0 +1,22 @@
+As of June 2010, the Hauksson-2000 model has been replaced by the Lin-2007 model as the background tomographic model for southern California simulations.
+
+In order to implement the Hauksson-2000 model, the following changes should be made.
+
+1. In create_regions_mesh.f90, uncomment the call to the Hauksson (2000) model and comment the call to the Lin et al. (2007) model.
+
+ call get_value_string(HAUKSSON_REGIONAL_MODEL_FILE, &
+ 'model.HAUKSSON_REGIONAL_MODEL_FILE', &
+ 'DATA/hauksson_model/hauksson_final_grid_smooth.dat')
+! call get_value_string(HAUKSSON_REGIONAL_MODEL_FILE, &
+! 'model.HAUKSSON_REGIONAL_MODEL_FILE', &
+! 'DATA/lin_model/lin_final_grid_smooth.dat')
+
+
+2. In constants.h.in, uncomment the block of lines beginning with "Hauksson (2000)" and comment the block of lines beginning with "Lin-Shearer-Hauksson-Thurber (2007)". These minor differences are due to the differences in the input models.
+
+3. cp ./hauksson_model.f90 ../../
+ This modified file accounts for the 9 layers of the Hauksson model, rather than the 8 layers of the Lin model.
+
+See the the PDF files in DATA/lin_model/ to examine the differences between the parameterizations of the two tomographic models.
+
+===========================================
\ No newline at end of file
Added: seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/hauksson_model/hauksson_model.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/hauksson_model/hauksson_model.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/hauksson_model/hauksson_model.f90 2010-06-03 14:12:04 UTC (rev 16878)
@@ -0,0 +1,226 @@
+!=====================================================================
+!
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
+!
+! Dimitri Komatitsch and Jeroen Tromp
+! Seismological Laboratory - California Institute of Technology
+! (c) California Institute of Technology September 2006
+!
+! 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.
+!
+!=====================================================================
+
+ subroutine hauksson_model(vp,vs,utm_x_eval,utm_y_eval,z_eval,vp_final,vs_final,MOHO_MAP_LUPEI)
+
+ implicit none
+
+ include "constants.h"
+
+!! DK DK UGLY one day, we should clarify the issue of merging Hauksson's Moho
+!! DK DK UGLY with our Lupei Moho. Should not be a big issue because in
+!! DK DK UGLY principle Hauksson used Lupei's map to build his model
+
+ double precision utm_x_eval,utm_y_eval,z_eval
+ double precision vp_final,vs_final
+ logical MOHO_MAP_LUPEI
+
+ double precision, dimension(NLAYERS_HAUKSSON,NGRID_NEW_HAUKSSON,NGRID_NEW_HAUKSSON) :: vp,vs
+ double precision, dimension(NLAYERS_HAUKSSON) :: vp_interp,vs_interp
+
+ integer ilayer
+ integer icell_interp_x,icell_interp_y
+ double precision spacing_x,spacing_y
+ double precision utm_x_eval_copy,utm_y_eval_copy
+ double precision gamma_interp_x,gamma_interp_y,gamma_interp_z
+ double precision v1,v2,v3,v4
+ double precision vp_upper,vs_upper,vp_lower,vs_lower,z_upper,z_lower
+
+! copy input values
+ utm_x_eval_copy = utm_x_eval
+ utm_y_eval_copy = utm_y_eval
+
+! make sure we stay inside Hauksson's grid
+ if(utm_x_eval_copy < UTM_X_ORIG_HAUKSSON) utm_x_eval_copy = UTM_X_ORIG_HAUKSSON
+ if(utm_y_eval_copy < UTM_Y_ORIG_HAUKSSON) utm_y_eval_copy = UTM_Y_ORIG_HAUKSSON
+
+! determine spacing and cell for linear interpolation
+ spacing_x = (utm_x_eval_copy - UTM_X_ORIG_HAUKSSON) / SPACING_UTM_X_HAUKSSON
+ spacing_y = (utm_y_eval_copy - UTM_Y_ORIG_HAUKSSON) / SPACING_UTM_Y_HAUKSSON
+
+ icell_interp_x = int(spacing_x) + 1
+ icell_interp_y = int(spacing_y) + 1
+
+ gamma_interp_x = spacing_x - int(spacing_x)
+ gamma_interp_y = spacing_y - int(spacing_y)
+
+! suppress edge effects for points outside of Hauksson's model
+ if(icell_interp_x < 1) then
+ icell_interp_x = 1
+ gamma_interp_x = 0.d0
+ endif
+ if(icell_interp_x > NGRID_NEW_HAUKSSON-1) then
+ icell_interp_x = NGRID_NEW_HAUKSSON-1
+ gamma_interp_x = 1.d0
+ endif
+
+ if(icell_interp_y < 1) then
+ icell_interp_y = 1
+ gamma_interp_y = 0.d0
+ endif
+ if(icell_interp_y > NGRID_NEW_HAUKSSON-1) then
+ icell_interp_y = NGRID_NEW_HAUKSSON-1
+ gamma_interp_y = 1.d0
+ endif
+
+! make sure interpolation makes sense
+ if(gamma_interp_x < -0.001d0 .or. gamma_interp_x > 1.001d0) &
+ stop 'interpolation in x is incorrect in Hauksson'
+ if(gamma_interp_y < -0.001d0 .or. gamma_interp_y > 1.001d0) &
+ stop 'interpolation in y is incorrect in Hauksson'
+
+! interpolate Hauksson's model at right location using bilinear interpolation
+ do ilayer = 1,NLAYERS_HAUKSSON
+
+! for Vp
+ v1 = vp(ilayer,icell_interp_x,icell_interp_y)
+ v2 = vp(ilayer,icell_interp_x+1,icell_interp_y)
+ v3 = vp(ilayer,icell_interp_x+1,icell_interp_y+1)
+ v4 = vp(ilayer,icell_interp_x,icell_interp_y+1)
+
+ vp_interp(ilayer) = v1*(1.-gamma_interp_x)*(1.-gamma_interp_y) + &
+ v2*gamma_interp_x*(1.-gamma_interp_y) + &
+ v3*gamma_interp_x*gamma_interp_y + &
+ v4*(1.-gamma_interp_x)*gamma_interp_y
+
+! for Vs
+ v1 = vs(ilayer,icell_interp_x,icell_interp_y)
+ v2 = vs(ilayer,icell_interp_x+1,icell_interp_y)
+ v3 = vs(ilayer,icell_interp_x+1,icell_interp_y+1)
+ v4 = vs(ilayer,icell_interp_x,icell_interp_y+1)
+
+ vs_interp(ilayer) = v1*(1.-gamma_interp_x)*(1.-gamma_interp_y) + &
+ v2*gamma_interp_x*(1.-gamma_interp_y) + &
+ v3*gamma_interp_x*gamma_interp_y + &
+ v4*(1.-gamma_interp_x)*gamma_interp_y
+
+ enddo
+
+! choose right values depending on depth of target point
+ if(z_eval >= Z_HAUKSSON_LAYER_1) then
+ vp_final = vp_interp(1)
+ vs_final = vs_interp(1)
+ return
+
+ else if(z_eval <= Z_HAUKSSON_LAYER_9) then
+ vp_final = vp_interp(9)
+ vs_final = vs_interp(9)
+ return
+
+ else if(z_eval >= Z_HAUKSSON_LAYER_2) then
+ vp_upper = vp_interp(1)
+ vs_upper = vs_interp(1)
+ z_upper = Z_HAUKSSON_LAYER_1
+
+ vp_lower = vp_interp(2)
+ vs_lower = vs_interp(2)
+ z_lower = Z_HAUKSSON_LAYER_2
+
+ else if(z_eval >= Z_HAUKSSON_LAYER_3) then
+ vp_upper = vp_interp(2)
+ vs_upper = vs_interp(2)
+ z_upper = Z_HAUKSSON_LAYER_2
+
+ vp_lower = vp_interp(3)
+ vs_lower = vs_interp(3)
+ z_lower = Z_HAUKSSON_LAYER_3
+
+ else if(z_eval >= Z_HAUKSSON_LAYER_4) then
+ vp_upper = vp_interp(3)
+ vs_upper = vs_interp(3)
+ z_upper = Z_HAUKSSON_LAYER_3
+
+ vp_lower = vp_interp(4)
+ vs_lower = vs_interp(4)
+ z_lower = Z_HAUKSSON_LAYER_4
+
+ else if(z_eval >= Z_HAUKSSON_LAYER_5) then
+ vp_upper = vp_interp(4)
+ vs_upper = vs_interp(4)
+ z_upper = Z_HAUKSSON_LAYER_4
+
+ vp_lower = vp_interp(5)
+ vs_lower = vs_interp(5)
+ z_lower = Z_HAUKSSON_LAYER_5
+
+ else if(z_eval >= Z_HAUKSSON_LAYER_6) then
+ vp_upper = vp_interp(5)
+ vs_upper = vs_interp(5)
+ z_upper = Z_HAUKSSON_LAYER_5
+
+ vp_lower = vp_interp(6)
+ vs_lower = vs_interp(6)
+ z_lower = Z_HAUKSSON_LAYER_6
+
+ else if(z_eval >= Z_HAUKSSON_LAYER_7) then
+ vp_upper = vp_interp(6)
+ vs_upper = vs_interp(6)
+ z_upper = Z_HAUKSSON_LAYER_6
+
+ vp_lower = vp_interp(7)
+ vs_lower = vs_interp(7)
+ z_lower = Z_HAUKSSON_LAYER_7
+
+ else if(z_eval >= Z_HAUKSSON_LAYER_8) then
+ vp_upper = vp_interp(7)
+ vs_upper = vs_interp(7)
+ z_upper = Z_HAUKSSON_LAYER_7
+
+ vp_lower = vp_interp(8)
+ vs_lower = vs_interp(8)
+ z_lower = Z_HAUKSSON_LAYER_8
+
+ else
+ if(.not. MOHO_MAP_LUPEI) then
+ vp_upper = vp_interp(8)
+ vs_upper = vs_interp(8)
+ z_upper = Z_HAUKSSON_LAYER_8
+
+ vp_lower = vp_interp(9)
+ vs_lower = vs_interp(9)
+ z_lower = Z_HAUKSSON_LAYER_9
+ !!! waiting for better interpolation of Moho maps.
+ else
+ vp_upper = vp_interp(8)
+ vs_upper = vs_interp(8)
+ z_upper = Z_HAUKSSON_LAYER_8
+
+ vp_lower = vp_interp(9)
+ vs_lower = vs_interp(9)
+ z_lower = Z_HAUKSSON_LAYER_9
+ endif
+
+ endif
+
+ gamma_interp_z = (z_eval - z_lower) / (z_upper - z_lower)
+
+ if(gamma_interp_z < -0.001d0 .or. gamma_interp_z > 1.001d0) &
+ stop 'interpolation in z is incorrect in Hauksson'
+
+ vp_final = vp_upper * gamma_interp_z + vp_lower * (1.-gamma_interp_z)
+ vs_final = vs_upper * gamma_interp_z + vs_lower * (1.-gamma_interp_z)
+
+ end subroutine hauksson_model
+
Deleted: seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/lin_model/README
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/lin_model/README 2010-06-03 10:36:19 UTC (rev 16877)
+++ seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/lin_model/README 2010-06-03 14:12:04 UTC (rev 16878)
@@ -1,58 +0,0 @@
-Carl Tape, 02-August-2007
-Tomography model of Lin-Shearer-Hauksson-Thurber (2007).
-
-See lin_2007/README for details on the tomographic model.
-
-hauk-lin.pdf shows some subtle differences between the
-Hauksson (2000) model and the Lin et al. (2007) model.
-
-new_bounds.pdf shows the nodes for the Lin model.
-red dashed -- 'standard' simulation domain
-red solid -- Hauksson interpolated regular grid domain
-black dashed -- 'new' simulation domain (extended west)
-black solid -- new Hauksson interpolated regular grid domain
-
-The starting file lin_new_format.dat was generated by
-the Matlab script hauksson_tomo_model.m, which is included
-here for completeness.
-
-To generate lin_final_grid_smooth.dat, the file read in by SPECFEM3D,
-you must execute runall.csh.
-
-CHANGES TO RUN IN SPECFEM3D
------------------------------------
-
-First, in create_regions_mesh.f90, you need to uncomment the call to the
-Lin et al. (2007) model and comment the call to the Hauksson (2000) model.
-
-! call get_value_string(HAUKSSON_REGIONAL_MODEL_FILE, &
-! 'model.HAUKSSON_REGIONAL_MODEL_FILE', &
-! 'DATA/hauksson_model/hauksson_final_grid_smooth.dat')
- call get_value_string(HAUKSSON_REGIONAL_MODEL_FILE, &
- 'model.HAUKSSON_REGIONAL_MODEL_FILE', &
- 'DATA/lin_model/lin_final_grid_smooth.dat')
-
-Second, the parameters in constants.h need to be adjusted to reflect the slightly
-different parameterization used in the Lin model versus the Hauksson model.
-I have left the pertinent lines in SPECFEM3D/DATA/lin_model/constants.h; basically
-you need to comment the Hauksson model lines and un-comment the Lin model lines.
-
-Third, by having one less layer in the Lin model, I needed to make minor
-adjustments to hauksson_model.f90, which I have copied to
-SPECFEM3D/DATA/lin_model/hauksson_model.f90
-
-Fourth, make sure the extended region is specified in the Par_file,
-which I have copied to
-SPECFEM3D/DATA/lin_model/Par_file
-
-LONGITUDE_MIN = -121.6d0 ! -120.3d0
-NEX_XI = 336
-NEX_ETA = 288
-NPROC_XI = 14
-NPROC_ETA = 12
-
------------------------------------
-
-In the future, we can decide how to commit these changes. For example, perhaps having
-a LIN_MODEL option would be sensible, or eliminating the Hauksson model altogether,
-since this is presubaly an updated version of his model.
Copied: seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/lin_model/README_lin_model (from rev 16872, seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/lin_model/README)
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/lin_model/README_lin_model (rev 0)
+++ seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/lin_model/README_lin_model 2010-06-03 14:12:04 UTC (rev 16878)
@@ -0,0 +1,65 @@
+Carl Tape, 02-August-2007, 02-June-2010
+Tomography model of Lin-Shearer-Hauksson-Thurber (JGR 2007).
+
+As of June 2010, the Lin-2007 model is the default backgrodun tomographic model for the southern California simulations. This model replaces the Hauksson-2000 model, which can be implemented with minor modifications -- see DATA/hauksson_model/README_hauksson_model for details.
+
+-----------------------------------
+IMPLEMENTING THE LIN MODEL FOR SPECFEM3D
+-----------------------------------
+
+This procedure is adapted from the files in DATA/hauksson_model/
+
+1. The original model files from Egill Hauksson are in lin_2007/
+2. The file lin_new_format.dat is created in hauksson_tomo_model.m,
+ which requires additional Matlab scripts to run.
+3. runall.csh contains three command that densify the grid to
+ create lin_final_grid_smooth.dat, which is read into SPECFEM3D.
+
+-----------------------------------
+COMPARISONS BETWEEN HAUKSSON-2000 and LIN-2007 MODELS
+-----------------------------------
+
+See lin_2007/README for details on the tomographic model.
+
+hauk-lin.pdf shows some subtle differences between the Hauksson (2000) model and the Lin et al. (2007) model.
+
+new_bounds.pdf shows the nodes for the Lin model.
+red dashed -- 'standard' simulation domain
+red solid -- Hauksson interpolated regular grid domain
+black dashed -- 'new' simulation domain (extended west)
+black solid -- new Hauksson interpolated regular grid domain
+
+-----------------------------------
+CHANGES TO RUN IN SPECFEM3D
+-----------------------------------
+
+First, in create_regions_mesh.f90, you need to uncomment the call to the
+Lin et al. (2007) model and comment the call to the Hauksson (2000) model.
+
+! call get_value_string(HAUKSSON_REGIONAL_MODEL_FILE, &
+! 'model.HAUKSSON_REGIONAL_MODEL_FILE', &
+! 'DATA/hauksson_model/hauksson_final_grid_smooth.dat')
+ call get_value_string(HAUKSSON_REGIONAL_MODEL_FILE, &
+ 'model.HAUKSSON_REGIONAL_MODEL_FILE', &
+ 'DATA/lin_model/lin_final_grid_smooth.dat')
+
+Second, the parameters in constants.h need to be adjusted to reflect the slightly different parameterization used in the Lin model versus the Hauksson model. I have left the pertinent lines in SPECFEM3D/DATA/lin_model/constants.h; basically you need to comment the Hauksson model lines and un-comment the Lin model lines.
+
+Third, by having one less layer in the Lin model, I needed to make minor adjustments to hauksson_model.f90, which I have copied to SPECFEM3D/DATA/lin_model/hauksson_model.f90
+
+Fourth, make sure the extended region is specified in the Par_file, which I have copied to
+SPECFEM3D/DATA/lin_model/Par_file
+
+LONGITUDE_MIN = -121.6d0 ! -120.3d0
+NEX_XI = 336
+NEX_ETA = 288
+NPROC_XI = 14
+NPROC_ETA = 12
+
+-----------------------------------
+
+In the future, we can decide how to commit these changes. For example, perhaps having
+a LIN_MODEL option would be sensible, or eliminating the Hauksson model altogether,
+since this is presubaly an updated version of his model.
+
+=======================================================
\ No newline at end of file
Modified: seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/lin_model/hauksson_model.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/lin_model/hauksson_model.f90 2010-06-03 10:36:19 UTC (rev 16877)
+++ seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/lin_model/hauksson_model.f90 2010-06-03 14:12:04 UTC (rev 16878)
@@ -1,7 +1,7 @@
!=====================================================================
!
-! S p e c f e m 3 D B a s i n V e r s i o n 1 . 4
-! --------------------------------------------------
+! S p e c f e m 3 D V e r s i o n 1 . 4
+! ---------------------------------------
!
! Dimitri Komatitsch and Jeroen Tromp
! Seismological Laboratory - California Institute of Technology
@@ -183,15 +183,6 @@
vs_lower = vs_interp(7)
z_lower = Z_HAUKSSON_LAYER_7
-! else if(z_eval >= Z_HAUKSSON_LAYER_8) then
-! vp_upper = vp_interp(7)
-! vs_upper = vs_interp(7)
-! z_upper = Z_HAUKSSON_LAYER_7
-!
-! vp_lower = vp_interp(8)
-! vs_lower = vs_interp(8)
-! z_lower = Z_HAUKSSON_LAYER_8
-
else
if(.not. MOHO_MAP_LUPEI) then
vp_upper = vp_interp(7)
Modified: seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/lin_model/hauksson_tomo_model.m
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/lin_model/hauksson_tomo_model.m 2010-06-03 10:36:19 UTC (rev 16877)
+++ seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/lin_model/hauksson_tomo_model.m 2010-06-03 14:12:04 UTC (rev 16878)
@@ -1,11 +1,16 @@
%
% hauksson_tomo_model.m
-% CARL TAPE, 11-July-2007
-% printed xxx
+% Carl Tape, 11-July-2007
%
-% This program reads in the Hauksson (2000) and Lin-Shearer-Hauksson-Thurber
-% (2007) tomography models for southern California, and outputs some figures.
+% This program reads in the southern California tomography models of
+% Hauksson (JGR 2000)
+% Lin-Shearer-Hauksson-Thurber (JGR 2007)
+% and outputs some figures.
%
+% NOTE: THIS PROGRAM REQUIRES SOME ADDITIONAL CUSTOM SUBROUTINES.
+% If you want to run the program and generate the output
+% files and figures, email Carl Tape for the additional scripts.
+%
% calls xxx
% called by xxx
%
@@ -15,9 +20,9 @@
clc
format compact
-dir = '/net/denali/home2/carltape/gmt/tomography/';
+idir = '/home/carltape/gmt/tomography/';
-ihauk = 0;
+ihauk = 1; % =1 to also plot Hauksson (2000) model
iwrite = 0;
% bounds for 'standard' SPECFEM3D simulation
@@ -33,12 +38,81 @@
axbox_hauk = [-121 -114 32 37];
%------------------------------------------------------
+% dimensions of PoChen et al. LA Basin model (2007)
+%
+% ---------- Forwarded message ----------
+% Date: Tue, 30 Oct 2007 14:33:52 -0400
+% From: Po Chen <pseudopochen at gmail.com>
+% To: Carl Tape <carltape at gps.caltech.edu>
+% Subject: Re: BSSA paper, figure 1 bounds
+%
+% Hello Carl,
+%
+% I could not find the exact lat/lon for that box used in the paper. The
+% following numbers might be close:
+%
+% lon: -118.76 -117.22
+% lat: 33.66 34.41
+%
+% Yes, the depth extend of the inversion is 26km.
+%
+% Best,
+% Po
+if 0==1
+ ax0 = [-118.76 -117.22 33.66 34.41];
+ ax_utm = axes_utm2ll(ax0,11,0);
+ dx = diff(ax_utm(1:2));
+ dy = diff(ax_utm(3:4));
+ dz = 26*1e3;
+ V = dx*dy*dz;
+ disp(sprintf(' dx dy dz : %.2f km %.2f km %.2f km',dx*1e-3,dy*1e-3,dz*1e-3));
+ disp(sprintf(' Volume : %.4e m^3',V));
+
+ lons = [238.6100 239.0480 239.4235 238.9854];
+ lats = [35.6421 35.2860 35.5913 35.9474];
+ disp([lons'-360 lats']);
+
+ % bounds of our simulations
+ ax0 = [-121.6 -114.7 32.2 36.8];
+ ax0_utm = axes_utm2ll(ax0,11,0);
+ dx0 = diff(ax0_utm(1:2));
+ dy0 = diff(ax0_utm(3:4));
+ dz0 = 60*1e3;
+ V0 = dx0*dy0*dz0;
+ disp(sprintf(' dx dy dz : %.2f km %.2f km %.2f km',dx0*1e-3,dy0*1e-3,dz0*1e-3));
+ disp(sprintf(' Volume : %.4e m^3',V0));
+ disp(sprintf(' Volume ratio : %.2f m^3',V0/V));
+ break
+end
+%------------------------------------------------------
+
+if 0==1
+ % regrid_hauksson_regular.f90
+ % write(13,'(2i10,3e20.10)') i, j, utm_x_new(i,j), utm_y_new(i,j), distmin_all(i,j)
+ dall = load('distmin.dat');
+
+ utmx = dall(:,3);
+ utmy = dall(:,4);
+ dmin = dall(:,5)/1000; % km
+
+ figure; nr=2; nc=1;
+ subplot(nr,nc,1); hold on;
+ plot3(utmx,utmy,zeros,'r.','markersize',1);
+ plot3(utmx,utmy,dmin,'b.','markersize',1); grid on; box on;
+ xlabel('UTM-X (m)'); ylabel('UTM-Y (m)'); zlabel(' Min distance, km');
+ subplot(nr,nc,2); plot_histo(dmin,[0:1:10]); xlabel(' Min distance, km');
+ orient tall, wysiwyg
+
+ break
+end
+
+%------------------------------------------------------
% Hauksson (2000) tomo model
if ihauk==1
disp(' '); disp(' Hauksson (2000) tomography model');
- dir_hauk = [dir 'hauk_2000/'];
+ dir_hauk = [idir 'hauk_2000/'];
lines = textread([dir_hauk 'hauksson_new_format.dat'],'%s','delimiter','\n','whitespace','');
nlines = length(lines);
@@ -115,7 +189,7 @@
disp(' '); disp(' Lin-Shearer-Hauksson-Thurber (2007) tomography model');
% load Vp model
-dir_lin = [dir 'lin_2007/'];
+dir_lin = [idir 'lin_2007/'];
vp = load([dir_lin 'pa_all-1']);
alpha = vp(:,1);
lon = vp(:,3);
@@ -133,6 +207,18 @@
udep = unique(dep);
ndep = length(udep);
+% assuming a max depth of the model, compute the thickness of each layer
+dmax = 60;
+dlayer = diff(udep)/2;
+dtop = udep - [0; dlayer];
+dbot = udep + [dlayer ; dmax-udep(end)];
+dthickness = dbot-dtop;
+dall = [udep dtop dbot dthickness];
+if sum(dall(:,4)) ~= dmax, error(['thickness of all layers should equal ' num2str(dmax) ]); end
+disp('Here is how we might think of the depth layers:');
+disp(' reference top bottom thickness');
+disp(dall);
+
% find range of grid
mingY = min(gridY); maxgY = max(gridY);
mingX = min(gridX); maxgX = max(gridX);
@@ -167,6 +253,42 @@
disp(['Grid is ' num2str(NX) ' by ' num2str(NY) ' = ' num2str(npts) ' nodes']);
disp([ num2str(ndep) ' depth layers']);
+% extract the boundary lat-lon points (for GMT plots)
+% ORDERING: from top point, clockwise
+iNE_boundary = flipud( find( gridX(isurface) == mingX ) );
+iSE_boundary = find( gridY(isurface) == mingY );
+iSW_boundary = find( gridX(isurface) == maxgX );
+iNW_boundary = flipud( find( gridY(isurface) == maxgY ) );
+iboundary = isurface([iNE_boundary ; iSE_boundary(2:end) ; iSW_boundary(2:end) ; iNW_boundary]);
+%figure; plot(lon(iboundary),lat(iboundary),'b.-');
+if iwrite == 1
+ write_xy_points([dir_lin 'lin_boundary'],lon(iboundary),lat(iboundary));
+end
+
+% compute the mean 1D model -- including scaling to density from alpha
+alpha_1D = zeros(ndep,3);
+beta_1D = zeros(ndep,3);
+for ii=1:ndep
+ inds = find(dep == udep(ii));
+ alpha_1D(ii,:) = [min(alpha(inds)) mean(alpha(inds)) max(alpha(inds))] ;
+ beta_1D(ii,:) = [min(beta(inds)) mean(beta(inds)) max(beta(inds)) ] ;
+end
+rho_1D = alpha_rho(alpha_1D*1e3);
+disp(' ');
+disp('1D averaged model:');
+disp(' depth thickness vp-min vp-mean vp-max vs-min vs-mean vs-max rho-min rho-mean rho-max');
+disp([ udep dthickness alpha_1D beta_1D rho_1D/1000]);
+
+% using the 1D model and the layer thicknesses, compute the overall mean velocities
+alpha_mean = sum(alpha_1D(:,2) .* dthickness) / sum(dthickness);
+beta_mean = sum(beta_1D(:,2) .* dthickness) / sum(dthickness);
+c_mean = sqrt( alpha_mean^2 - (4/3)*beta_mean^2 );
+disp(' ');
+disp([' Overall mean velocities, assuming a bottom depth of ' num2str(dmax) ' km:'])
+disp(sprintf('%8s : %.1f (%.4f) km/s','alpha',alpha_mean,alpha_mean));
+disp(sprintf('%8s : %.1f (%.4f) km/s','beta',beta_mean,beta_mean));
+disp(sprintf('%8s : %.1f (%.4f) km/s','c',c_mean,c_mean));
+
%---------------------------------
% figures
@@ -187,7 +309,7 @@
if ihauk==1
plot(lon_hauk(isurface_hauk),lat_hauk(isurface_hauk),'r.');
end
-axis equal; ax0 = axis; axis(axes_expand(ax0,1.05));
+axis equal; ax0 = axis; axis(axes_expand(ax0,1.05,1));
plot(axbox_hauk([1 2 2 1 1]),axbox_hauk([3 3 4 4 3]),'r','linewidth',2);
plot(axbox([1 2 2 1 1]),axbox([3 3 4 4 3]),'k','linewidth',2);
if ihauk==1, legend('Lin et al. 2007','Hauksson 2000','Regular tomo in SPECFEM','SPECFEM bounds');
@@ -198,7 +320,7 @@
figure; hold on;
plot(lon(isurface),lat(isurface),'k.')
-axis equal; ax0 = axis; axis(axes_expand(ax0,1.05));
+axis equal; ax0 = axis; axis(axes_expand(ax0,1.05,1));
plot(lon(izeroY),lat(izeroY),'ro');
plot(lon(izeroX),lat(izeroX),'bo');
text(lon1,lat1,st1); text(lon2,lat2,st2); text(lon3,lat3,st3); text(lon4,lat4,st4);
@@ -209,7 +331,7 @@
figure; hold on;
plot(gridX(isurface),gridY(isurface),'k.')
-axis equal; ax0 = axis; axis(axes_expand(ax0,1.05));
+axis equal; ax0 = axis; axis(axes_expand(ax0,1.05,1));
plot(gridX(izeroY),gridY(izeroY),'ro');
plot(gridX(izeroX),gridY(izeroX),'bo');
plot(gridX(iorigin),gridY(iorigin),'kp','markersize',15);
@@ -249,4 +371,4 @@
fclose(fid);
end
-%==========================================
\ No newline at end of file
+%==========================================
Modified: seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/lin_model/lin_2007/README
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/lin_model/lin_2007/README 2010-06-03 10:36:19 UTC (rev 16877)
+++ seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/DATA/lin_model/lin_2007/README 2010-06-03 14:12:04 UTC (rev 16878)
@@ -1,17 +1,16 @@
Carl Tape, 11-July-2007
-Tomography model of Lin-Shearer-Hauksson-Thurber (2007).
+Tomography model of Lin-Shearer-Hauksson-Thurber (JGR 2007).
The files pa_all-1 and ra_all-1 were obtained from Egill Hauksson June 2007.
pa_all-1 Vp tomography model for SoCal
ra_all-1 Vp/Vs tomography model for SoCal
-In order to convert these files into the format of hauksson_new_format.dat (see haukson_model directory), I used a script in Matlab (hauksson_tomo_model.m; not included). That output file is
+In order to convert these files into the format of hauksson_new_format.dat (see haukson_model directory), I used the Matlab script lin_model/hauksson_tomo_model.m. That output file is
lin_new_format.dat
-
---------- Forwarded message ----------
Date: Fri, 15 Jun 2007 14:31:37 -0700
From: Egill Hauksson <ehauksson at gmail.com>
@@ -24,3 +23,5 @@
file, it is related to something else.
Let me know how this goes.
Egill
+
+===========================================
\ No newline at end of file
Modified: seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/constants.h.in 2010-06-03 10:36:19 UTC (rev 16877)
+++ seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/constants.h.in 2010-06-03 14:12:04 UTC (rev 16878)
@@ -302,38 +302,61 @@
double precision, parameter :: GOCAD_ST_NO_DATA_VALUE = -99999
!
-!--- larger Hauksson model for entire So-Cal, 15 km resolution
+!--- larger crustal tomographic model for entire So-Cal, 15 km resolution
+!--- default implementation is the Lin-2007 model
+!--- see DATA/hauksson_model/ to implement the Hauksson-2000 model
!
+!! Hauksson (2000)
+!! number of non-constant layers
+! integer, parameter :: NLAYERS_HAUKSSON = 9
+!! depth of layers
+! double precision, parameter :: Z_HAUKSSON_LAYER_1 = -1000.d0
+! double precision, parameter :: Z_HAUKSSON_LAYER_2 = -4000.d0
+! double precision, parameter :: Z_HAUKSSON_LAYER_3 = -6000.d0
+! double precision, parameter :: Z_HAUKSSON_LAYER_4 = -10000.d0
+! double precision, parameter :: Z_HAUKSSON_LAYER_5 = -15000.d0
+! double precision, parameter :: Z_HAUKSSON_LAYER_6 = -17000.d0
+! double precision, parameter :: Z_HAUKSSON_LAYER_7 = -22000.d0
+! double precision, parameter :: Z_HAUKSSON_LAYER_8 = -31000.d0
+! double precision, parameter :: Z_HAUKSSON_LAYER_9 = -33000.d0
+!
+! integer, parameter :: NGRID_NEW_HAUKSSON = 201
+!
+!! corners of new Hauksson interpolated grid
+! double precision, parameter :: UTM_X_ORIG_HAUKSSON = 122035.012d0
+! double precision, parameter :: UTM_X_END_HAUKSSON = 766968.628d0
+! double precision, parameter :: UTM_Y_ORIG_HAUKSSON = 3547232.986d0
+! double precision, parameter :: UTM_Y_END_HAUKSSON = 4098868.501d0
+
+! Lin-Shearer-Hauksson-Thurber (2007)
! number of non-constant layers
- integer, parameter :: NLAYERS_HAUKSSON = 9
+ integer, parameter :: NLAYERS_HAUKSSON = 8
! depth of layers
- double precision, parameter :: Z_HAUKSSON_LAYER_1 = -1000.d0
- double precision, parameter :: Z_HAUKSSON_LAYER_2 = -4000.d0
+ double precision, parameter :: Z_HAUKSSON_LAYER_1 = 0.d0
+ double precision, parameter :: Z_HAUKSSON_LAYER_2 = -3000.d0
double precision, parameter :: Z_HAUKSSON_LAYER_3 = -6000.d0
double precision, parameter :: Z_HAUKSSON_LAYER_4 = -10000.d0
double precision, parameter :: Z_HAUKSSON_LAYER_5 = -15000.d0
double precision, parameter :: Z_HAUKSSON_LAYER_6 = -17000.d0
double precision, parameter :: Z_HAUKSSON_LAYER_7 = -22000.d0
double precision, parameter :: Z_HAUKSSON_LAYER_8 = -31000.d0
- double precision, parameter :: Z_HAUKSSON_LAYER_9 = -33000.d0
- integer, parameter :: NGRID_NEW_HAUKSSON = 201
+ integer, parameter :: NGRID_NEW_HAUKSSON = 241
-! corners of new Hauksson's interpolated grid
- double precision, parameter :: UTM_X_ORIG_HAUKSSON = 122035.012d0
- double precision, parameter :: UTM_X_END_HAUKSSON = 766968.628d0
- double precision, parameter :: UTM_Y_ORIG_HAUKSSON = 3547232.986d0
- double precision, parameter :: UTM_Y_END_HAUKSSON = 4098868.501d0
+! corners of new Hauksson interpolated grid
+ double precision, parameter :: UTM_X_ORIG_HAUKSSON = 88021.568d0
+ double precision, parameter :: UTM_X_END_HAUKSSON = 861517.886d0
+ double precision, parameter :: UTM_Y_ORIG_HAUKSSON = 3404059.875d0
+ double precision, parameter :: UTM_Y_END_HAUKSSON = 4180234.582d0
double precision, parameter :: SPACING_UTM_X_HAUKSSON = (UTM_X_END_HAUKSSON - UTM_X_ORIG_HAUKSSON) / (NGRID_NEW_HAUKSSON-1.d0)
double precision, parameter :: SPACING_UTM_Y_HAUKSSON = (UTM_Y_END_HAUKSSON - UTM_Y_ORIG_HAUKSSON) / (NGRID_NEW_HAUKSSON-1.d0)
! layers in the So-Cal regional model
! DEPTH_MOHO_SOCAL = -35 km was based on Dreger and Helmberger (1990)
-! and is (July 2007) the preferred Moho depth for Dreger.
! The depth of 32 km is used in the standard processing (Wald et al., 1995)
-! of SoCal events and is the value in the original Kanamori-Hadley (1975) model.
+! of So-Cal events and is the value in the original Kanamori-Hadley (1975) model.
double precision, parameter :: DEPTH_5p5km_SOCAL = -5500.d0
double precision, parameter :: DEPTH_16km_SOCAL = -16000.d0
double precision, parameter :: DEPTH_MOHO_SOCAL = -32000.d0
Modified: seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/create_regions_mesh.f90 2010-06-03 10:36:19 UTC (rev 16877)
+++ seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/create_regions_mesh.f90 2010-06-03 14:12:04 UTC (rev 16878)
@@ -394,12 +394,12 @@
!--- read Hauksson's model
if(HAUKSSON_REGIONAL_MODEL) then
+! call get_value_string(HAUKSSON_REGIONAL_MODEL_FILE, &
+! 'model.HAUKSSON_REGIONAL_MODEL_FILE', &
+! 'DATA/hauksson_model/hauksson_final_grid_smooth.dat')
call get_value_string(HAUKSSON_REGIONAL_MODEL_FILE, &
'model.HAUKSSON_REGIONAL_MODEL_FILE', &
- 'DATA/hauksson_model/hauksson_final_grid_smooth.dat')
-! call get_value_string(HAUKSSON_REGIONAL_MODEL_FILE, &
-! 'model.HAUKSSON_REGIONAL_MODEL_FILE', &
-! 'DATA/lin_model/lin_final_grid_smooth.dat')
+ 'DATA/lin_model/lin_final_grid_smooth.dat')
open(unit=14,file=HAUKSSON_REGIONAL_MODEL_FILE,status='old',action='read')
do iy = 1,NGRID_NEW_HAUKSSON
do ix = 1,NGRID_NEW_HAUKSSON
Modified: seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/hauksson_model.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/hauksson_model.f90 2010-06-03 10:36:19 UTC (rev 16877)
+++ seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/hauksson_model.f90 2010-06-03 14:12:04 UTC (rev 16878)
@@ -124,9 +124,9 @@
vs_final = vs_interp(1)
return
- else if(z_eval <= Z_HAUKSSON_LAYER_9) then
- vp_final = vp_interp(9)
- vs_final = vs_interp(9)
+ else if(z_eval <= Z_HAUKSSON_LAYER_8) then
+ vp_final = vp_interp(8)
+ vs_final = vs_interp(8)
return
else if(z_eval >= Z_HAUKSSON_LAYER_2) then
@@ -183,33 +183,24 @@
vs_lower = vs_interp(7)
z_lower = Z_HAUKSSON_LAYER_7
- else if(z_eval >= Z_HAUKSSON_LAYER_8) then
- vp_upper = vp_interp(7)
- vs_upper = vs_interp(7)
- z_upper = Z_HAUKSSON_LAYER_7
-
- vp_lower = vp_interp(8)
- vs_lower = vs_interp(8)
- z_lower = Z_HAUKSSON_LAYER_8
-
else
if(.not. MOHO_MAP_LUPEI) then
- vp_upper = vp_interp(8)
- vs_upper = vs_interp(8)
- z_upper = Z_HAUKSSON_LAYER_8
+ vp_upper = vp_interp(7)
+ vs_upper = vs_interp(7)
+ z_upper = Z_HAUKSSON_LAYER_7
- vp_lower = vp_interp(9)
- vs_lower = vs_interp(9)
- z_lower = Z_HAUKSSON_LAYER_9
+ vp_lower = vp_interp(8)
+ vs_lower = vs_interp(8)
+ z_lower = Z_HAUKSSON_LAYER_8
!!! waiting for better interpolation of Moho maps.
else
- vp_upper = vp_interp(8)
- vs_upper = vs_interp(8)
- z_upper = Z_HAUKSSON_LAYER_8
+ vp_upper = vp_interp(7)
+ vs_upper = vs_interp(7)
+ z_upper = Z_HAUKSSON_LAYER_7
- vp_lower = vp_interp(9)
- vs_lower = vs_interp(9)
- z_lower = Z_HAUKSSON_LAYER_9
+ vp_lower = vp_interp(8)
+ vs_lower = vs_interp(8)
+ z_lower = Z_HAUKSSON_LAYER_8
endif
endif
More information about the CIG-COMMITS
mailing list