[cig-commits] [commit] : 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. (82f4c29)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Nov 14 20:16:41 PST 2013
Repository : ssh://geoshell/specfem3d
On branch :
Link : https://github.com/geodynamics/specfem2d/compare/1e201257d91c794056b990a43329e05d04f77454...0000000000000000000000000000000000000000
>---------------------------------------------------------------
commit 82f4c2901c8d8834f78f27feff3af711518c15bb
Author: Carl Tape <carltape at gi.alaska.edu>
Date: Thu Jun 3 14:12:04 2010 +0000
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.
>---------------------------------------------------------------
82f4c2901c8d8834f78f27feff3af711518c15bb
DATA/hauksson_model/README_hauksson_model | 22 ++++
.../hauksson_model/hauksson_model.f90 | 0
DATA/lin_model/{README => README_lin_model} | 47 ++++---
DATA/lin_model/hauksson_model.f90 | 13 +-
DATA/lin_model/hauksson_tomo_model.m | 144 +++++++++++++++++++--
DATA/lin_model/lin_2007/README | 7 +-
constants.h.in | 49 +++++--
create_regions_mesh.f90 | 8 +-
hauksson_model.f90 | 39 +++---
9 files changed, 243 insertions(+), 86 deletions(-)
diff --git a/DATA/hauksson_model/README_hauksson_model b/DATA/hauksson_model/README_hauksson_model
new file mode 100644
index 0000000..e57c340
--- /dev/null
+++ b/DATA/hauksson_model/README_hauksson_model
@@ -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
diff --git a/hauksson_model.f90 b/DATA/hauksson_model/hauksson_model.f90
similarity index 100%
copy from hauksson_model.f90
copy to DATA/hauksson_model/hauksson_model.f90
diff --git a/DATA/lin_model/README b/DATA/lin_model/README_lin_model
similarity index 50%
rename from DATA/lin_model/README
rename to DATA/lin_model/README_lin_model
index 8d54303..0b65a86 100644
--- a/DATA/lin_model/README
+++ b/DATA/lin_model/README_lin_model
@@ -1,10 +1,27 @@
-Carl Tape, 02-August-2007
-Tomography model of Lin-Shearer-Hauksson-Thurber (2007).
+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.
+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
@@ -12,13 +29,7 @@ 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
-----------------------------------
@@ -32,17 +43,11 @@ Lin et al. (2007) model and comment the call to the Hauksson (2000) model.
'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.
+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
+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
+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
@@ -56,3 +61,5 @@ 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
diff --git a/DATA/lin_model/hauksson_model.f90 b/DATA/lin_model/hauksson_model.f90
index 9d17e4e..5f4a456 100644
--- a/DATA/lin_model/hauksson_model.f90
+++ b/DATA/lin_model/hauksson_model.f90
@@ -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)
diff --git a/DATA/lin_model/hauksson_tomo_model.m b/DATA/lin_model/hauksson_tomo_model.m
index 7e883ff..e2ed9cc 100644
--- a/DATA/lin_model/hauksson_tomo_model.m
+++ b/DATA/lin_model/hauksson_tomo_model.m
@@ -1,10 +1,15 @@
%
% 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 @@ close all
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 = [-120.3 -114.7 32.2 36.8];
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 @@ end
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 @@ beta = 1./vpvs(:,1) .* alpha;
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 @@ npts = NX * NY;
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 @@ plot(lon(isurface),lat(isurface),'b.');
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 @@ orient tall, wysiwyg
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 @@ orient tall, wysiwyg
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);
diff --git a/DATA/lin_model/lin_2007/README b/DATA/lin_model/lin_2007/README
index befecfd..bf720cf 100644
--- a/DATA/lin_model/lin_2007/README
+++ b/DATA/lin_model/lin_2007/README
@@ -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 @@ The Lin et al. Vp and Vp/Vs models are attached. Just ignore the readme
file, it is related to something else.
Let me know how this goes.
Egill
+
+===========================================
\ No newline at end of file
diff --git a/constants.h.in b/constants.h.in
index 2abf774..6ca6ca9 100644
--- a/constants.h.in
+++ b/constants.h.in
@@ -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
diff --git a/create_regions_mesh.f90 b/create_regions_mesh.f90
index 5ae6a4a..99dd020 100644
--- a/create_regions_mesh.f90
+++ b/create_regions_mesh.f90
@@ -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/lin_model/lin_final_grid_smooth.dat')
+! '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')
open(unit=14,file=HAUKSSON_REGIONAL_MODEL_FILE,status='old',action='read')
do iy = 1,NGRID_NEW_HAUKSSON
do ix = 1,NGRID_NEW_HAUKSSON
diff --git a/hauksson_model.f90 b/hauksson_model.f90
index f295bcf..5f4a456 100644
--- a/hauksson_model.f90
+++ b/hauksson_model.f90
@@ -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