[cig-commits] r16121 - seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN
pieyre at geodynamics.org
pieyre at geodynamics.org
Wed Dec 30 06:46:29 PST 2009
Author: pieyre
Date: 2009-12-30 06:46:28 -0800 (Wed, 30 Dec 2009)
New Revision: 16121
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/interpolate_gocad_block_HR.f90
seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/interpolate_gocad_block_MR.f90
Log:
fixed the problem in the smoothing of the transition between HR and MR blocks
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 2009-12-29 23:03:41 UTC (rev 16120)
+++ seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/create_regions_mesh.f90 2009-12-30 14:46:28 UTC (rev 16121)
@@ -205,6 +205,7 @@
double precision c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
! for the Harvard 3-D basin model
+ logical is_second_call
double precision vp_block_gocad_MR(0:NX_GOCAD_MR-1,0:NY_GOCAD_MR-1,0:NZ_GOCAD_MR-1)
double precision vp_block_gocad_HR(0:NX_GOCAD_HR-1,0:NY_GOCAD_HR-1,0:NZ_GOCAD_HR-1)
integer irecord,nrecord,i_vp
@@ -414,6 +415,7 @@
!--- read the Harvard 3-D basin model
if(HARVARD_3D_GOCAD_MODEL) then
+ is_second_call = .false.
! read medium-resolution model
! initialize array to undefined values everywhere
@@ -709,7 +711,7 @@
VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM, &
IMPOSE_MINIMUM_VP_GOCAD,THICKNESS_TAPER_BLOCK_MR, &
vp_hauksson,vs_hauksson,doubling_index,HAUKSSON_REGIONAL_MODEL,&
- MOHO_MAP_LUPEI)
+ MOHO_MAP_LUPEI,is_second_call)
! then superimpose high-resolution model
if(xmesh >= ORIG_X_GOCAD_HR &
Modified: seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/interpolate_gocad_block_HR.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/interpolate_gocad_block_HR.f90 2009-12-29 23:03:41 UTC (rev 16120)
+++ seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/interpolate_gocad_block_HR.f90 2009-12-30 14:46:28 UTC (rev 16121)
@@ -9,6 +9,8 @@
include "constants.h"
+ logical is_second_call
+
double precision vp_block_gocad_HR(0:NX_GOCAD_HR-1,0:NY_GOCAD_HR-1,0:NZ_GOCAD_HR-1)
double precision vp_block_gocad_MR(0:NX_GOCAD_MR-1,0:NY_GOCAD_MR-1,0:NZ_GOCAD_MR-1)
@@ -77,10 +79,10 @@
v7 < 19000. .and. v8 < 19000.) then
! set flag indicating whether point is in the sediments
- point_is_in_sediments = .true.
+ point_is_in_sediments = .true.
! use trilinear interpolation in cell to define Vp
- vp_final = &
+ vp_final = &
v1*(1.-gamma_interp_x)*(1.-gamma_interp_y)*(1.-gamma_interp_z) + &
v2*gamma_interp_x*(1.-gamma_interp_y)*(1.-gamma_interp_z) + &
v3*gamma_interp_x*gamma_interp_y*(1.-gamma_interp_z) + &
@@ -91,82 +93,83 @@
v8*(1.-gamma_interp_x)*gamma_interp_y*gamma_interp_z
! impose minimum velocity if needed
- if(IMPOSE_MINIMUM_VP_GOCAD .and. vp_final < VP_MIN_GOCAD) vp_final = VP_MIN_GOCAD
+ if(IMPOSE_MINIMUM_VP_GOCAD .and. vp_final < VP_MIN_GOCAD) vp_final = VP_MIN_GOCAD
! taper edges to make smooth transition between MR and HR blocks
! get value from edge of medium-resolution block
! then use linear interpolation from edge of the model
- if(TAPER_GOCAD_TRANSITIONS) then
+ if(TAPER_GOCAD_TRANSITIONS) then
+ is_second_call = .true.
! x = xmin
- if(utm_x_eval < ORIG_X_GOCAD_HR + THICKNESS_TAPER_BLOCK_HR) then
- gamma_interp_x = (utm_x_eval - ORIG_X_GOCAD_HR) / THICKNESS_TAPER_BLOCK_HR
- call interpolate_gocad_block_MR(vp_block_gocad_MR, &
- ORIG_X_GOCAD_HR,utm_y_eval,z_eval,rho_ref_MR,vp_ref_MR,vs_ref_MR,dummy_flag, &
- VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM, &
- IMPOSE_MINIMUM_VP_GOCAD,THICKNESS_TAPER_BLOCK_HR, &
- vp_hauksson,vs_hauksson,doubling_index,HAUKSSON_REGIONAL_MODEL,MOHO_MAP_LUPEI)
+ if(utm_x_eval < ORIG_X_GOCAD_HR + THICKNESS_TAPER_BLOCK_HR) then
+ gamma_interp_x = (utm_x_eval - ORIG_X_GOCAD_HR) / THICKNESS_TAPER_BLOCK_HR
+ call interpolate_gocad_block_MR(vp_block_gocad_MR, &
+ ORIG_X_GOCAD_HR,utm_y_eval,z_eval,rho_ref_MR,vp_ref_MR,vs_ref_MR,dummy_flag, &
+ VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM, &
+ IMPOSE_MINIMUM_VP_GOCAD,THICKNESS_TAPER_BLOCK_HR, &
+ vp_hauksson,vs_hauksson,doubling_index,HAUKSSON_REGIONAL_MODEL,MOHO_MAP_LUPEI,is_second_call)
! the vp_final we start from is the one we computed by interpolation in the HR block,
! and as we get closer to the transition with the MR block we mix it with the value on the edge of MR;
! it is not a linear interpolation between two fixed values but rather a linear variation of the amount of mixing
- vp_final = vp_ref_MR * (1. - gamma_interp_x) + vp_final * gamma_interp_x
+ vp_final = vp_ref_MR * (1. - gamma_interp_x) + vp_final * gamma_interp_x
! x = xmax
- else if(utm_x_eval > END_X_GOCAD_HR - THICKNESS_TAPER_BLOCK_HR) then
- gamma_interp_x = (utm_x_eval - (END_X_GOCAD_HR - THICKNESS_TAPER_BLOCK_HR)) / THICKNESS_TAPER_BLOCK_HR
- call interpolate_gocad_block_MR(vp_block_gocad_MR, &
- END_X_GOCAD_HR,utm_y_eval,z_eval,rho_ref_MR,vp_ref_MR,vs_ref_MR,dummy_flag, &
- VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM, &
- IMPOSE_MINIMUM_VP_GOCAD,THICKNESS_TAPER_BLOCK_HR, &
- vp_hauksson,vs_hauksson,doubling_index,HAUKSSON_REGIONAL_MODEL, MOHO_MAP_LUPEI)
+ else if(utm_x_eval > END_X_GOCAD_HR - THICKNESS_TAPER_BLOCK_HR) then
+ gamma_interp_x = (utm_x_eval - (END_X_GOCAD_HR - THICKNESS_TAPER_BLOCK_HR)) / THICKNESS_TAPER_BLOCK_HR
+ call interpolate_gocad_block_MR(vp_block_gocad_MR, &
+ END_X_GOCAD_HR,utm_y_eval,z_eval,rho_ref_MR,vp_ref_MR,vs_ref_MR,dummy_flag, &
+ VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM, &
+ IMPOSE_MINIMUM_VP_GOCAD,THICKNESS_TAPER_BLOCK_HR, &
+ vp_hauksson,vs_hauksson,doubling_index,HAUKSSON_REGIONAL_MODEL, MOHO_MAP_LUPEI,is_second_call)
! the vp_final we start from is the one we computed by interpolation in the HR block,
! and as we get closer to the transition with the MR block we mix it with the value on the edge of MR;
! it is not a linear interpolation between two fixed values but rather a linear variation of the amount of mixing
- vp_final = vp_ref_MR * gamma_interp_x + vp_final * (1. - gamma_interp_x)
+ vp_final = vp_ref_MR * gamma_interp_x + vp_final * (1. - gamma_interp_x)
! y = ymin
- else if(utm_y_eval < ORIG_Y_GOCAD_HR + THICKNESS_TAPER_BLOCK_HR) then
- gamma_interp_y = (utm_y_eval - ORIG_Y_GOCAD_HR) / THICKNESS_TAPER_BLOCK_HR
- call interpolate_gocad_block_MR(vp_block_gocad_MR, &
- utm_x_eval,ORIG_Y_GOCAD_HR,z_eval,rho_ref_MR,vp_ref_MR,vs_ref_MR,dummy_flag, &
- VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM, &
- IMPOSE_MINIMUM_VP_GOCAD,THICKNESS_TAPER_BLOCK_HR, &
- vp_hauksson,vs_hauksson,doubling_index,HAUKSSON_REGIONAL_MODEL, MOHO_MAP_LUPEI)
+ else if(utm_y_eval < ORIG_Y_GOCAD_HR + THICKNESS_TAPER_BLOCK_HR) then
+ gamma_interp_y = (utm_y_eval - ORIG_Y_GOCAD_HR) / THICKNESS_TAPER_BLOCK_HR
+ call interpolate_gocad_block_MR(vp_block_gocad_MR, &
+ utm_x_eval,ORIG_Y_GOCAD_HR,z_eval,rho_ref_MR,vp_ref_MR,vs_ref_MR,dummy_flag, &
+ VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM, &
+ IMPOSE_MINIMUM_VP_GOCAD,THICKNESS_TAPER_BLOCK_HR, &
+ vp_hauksson,vs_hauksson,doubling_index,HAUKSSON_REGIONAL_MODEL, MOHO_MAP_LUPEI,is_second_call)
! the vp_final we start from is the one we computed by interpolation in the HR block,
! and as we get closer to the transition with the MR block we mix it with the value on the edge of MR;
! it is not a linear interpolation between two fixed values but rather a linear variation of the amount of mixing
- vp_final = vp_ref_MR * (1. - gamma_interp_y) + vp_final * gamma_interp_y
+ vp_final = vp_ref_MR * (1. - gamma_interp_y) + vp_final * gamma_interp_y
! y = ymax
- else if(utm_y_eval > END_Y_GOCAD_HR - THICKNESS_TAPER_BLOCK_HR) then
- gamma_interp_y = (utm_y_eval - (END_Y_GOCAD_HR - THICKNESS_TAPER_BLOCK_HR)) / THICKNESS_TAPER_BLOCK_HR
- call interpolate_gocad_block_MR(vp_block_gocad_MR, &
- utm_x_eval,END_Y_GOCAD_HR,z_eval,rho_ref_MR,vp_ref_MR,vs_ref_MR,dummy_flag, &
- VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM, &
- IMPOSE_MINIMUM_VP_GOCAD,THICKNESS_TAPER_BLOCK_HR, &
- vp_hauksson,vs_hauksson,doubling_index,HAUKSSON_REGIONAL_MODEL, MOHO_MAP_LUPEI)
+ else if(utm_y_eval > END_Y_GOCAD_HR - THICKNESS_TAPER_BLOCK_HR) then
+ gamma_interp_y = (utm_y_eval - (END_Y_GOCAD_HR - THICKNESS_TAPER_BLOCK_HR)) / THICKNESS_TAPER_BLOCK_HR
+ call interpolate_gocad_block_MR(vp_block_gocad_MR, &
+ utm_x_eval,END_Y_GOCAD_HR,z_eval,rho_ref_MR,vp_ref_MR,vs_ref_MR,dummy_flag, &
+ VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM, &
+ IMPOSE_MINIMUM_VP_GOCAD,THICKNESS_TAPER_BLOCK_HR, &
+ vp_hauksson,vs_hauksson,doubling_index,HAUKSSON_REGIONAL_MODEL, MOHO_MAP_LUPEI,is_second_call)
! the vp_final we start from is the one we computed by interpolation in the HR block,
! and as we get closer to the transition with the MR block we mix it with the value on the edge of MR;
! it is not a linear interpolation between two fixed values but rather a linear variation of the amount of mixing
- vp_final = vp_ref_MR * gamma_interp_y + vp_final * (1. - gamma_interp_y)
+ vp_final = vp_ref_MR * gamma_interp_y + vp_final * (1. - gamma_interp_y)
- endif
+ endif
+
+ endif
- endif
-
! use linear variation of vp/vs ratio with depth, between 0. and 8.5 km
- vp_vs_ratio = VP_VS_RATIO_GOCAD_BOTTOM + &
+ vp_vs_ratio = VP_VS_RATIO_GOCAD_BOTTOM + &
(VP_VS_RATIO_GOCAD_TOP - VP_VS_RATIO_GOCAD_BOTTOM) * &
(z_eval - (-8500.d0)) / (0.d0 - (-8500.d0))
! make sure ratio remains in interval
- if(vp_vs_ratio < VP_VS_RATIO_GOCAD_BOTTOM) vp_vs_ratio = VP_VS_RATIO_GOCAD_BOTTOM
- if(vp_vs_ratio > VP_VS_RATIO_GOCAD_TOP) vp_vs_ratio = VP_VS_RATIO_GOCAD_TOP
+ if(vp_vs_ratio < VP_VS_RATIO_GOCAD_BOTTOM) vp_vs_ratio = VP_VS_RATIO_GOCAD_BOTTOM
+ if(vp_vs_ratio > VP_VS_RATIO_GOCAD_TOP) vp_vs_ratio = VP_VS_RATIO_GOCAD_TOP
+
+ vs_final = vp_final / vp_vs_ratio
+ call compute_rho_estimate(rho_final,vp_final)
+
+ endif
- vs_final = vp_final / vp_vs_ratio
- call compute_rho_estimate(rho_final,vp_final)
+ end subroutine interpolate_gocad_block_HR
- endif
-
- end subroutine interpolate_gocad_block_HR
-
Modified: seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/interpolate_gocad_block_MR.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/interpolate_gocad_block_MR.f90 2009-12-29 23:03:41 UTC (rev 16120)
+++ seismo/3D/SPECFEM3D_SESAME/tags/v1.4.4_last_BASIN/interpolate_gocad_block_MR.f90 2009-12-30 14:46:28 UTC (rev 16121)
@@ -3,7 +3,7 @@
utm_x_eval,utm_y_eval,z_eval,rho_final,vp_final,vs_final,point_is_in_sediments, &
VP_MIN_GOCAD,VP_VS_RATIO_GOCAD_TOP,VP_VS_RATIO_GOCAD_BOTTOM, &
IMPOSE_MINIMUM_VP_GOCAD,THICKNESS_TAPER_BLOCK_MR, &
- vp_hauksson,vs_hauksson,doubling_index,HAUKSSON_REGIONAL_MODEL,MOHO_MAP_LUPEI)
+ vp_hauksson,vs_hauksson,doubling_index,HAUKSSON_REGIONAL_MODEL,MOHO_MAP_LUPEI,is_second_call)
implicit none
@@ -22,6 +22,7 @@
double precision xmesh,ymesh,zmesh,vs_dummy,rho_dummy
logical point_is_in_sediments,IMPOSE_MINIMUM_VP_GOCAD
+ logical is_second_call
! for Hauksson's model
integer doubling_index
@@ -90,10 +91,10 @@
v7 < 19000. .and. v8 < 19000.) then
! set flag indicating whether point is in the sediments
- point_is_in_sediments = .true.
+ point_is_in_sediments = .true.
! use trilinear interpolation in cell to define Vp
- vp_final = &
+ vp_final = &
v1*(1.-gamma_interp_x)*(1.-gamma_interp_y)*(1.-gamma_interp_z) + &
v2*gamma_interp_x*(1.-gamma_interp_y)*(1.-gamma_interp_z) + &
v3*gamma_interp_x*gamma_interp_y*(1.-gamma_interp_z) + &
@@ -104,82 +105,99 @@
v8*(1.-gamma_interp_x)*gamma_interp_y*gamma_interp_z
! impose minimum velocity if needed
- if(IMPOSE_MINIMUM_VP_GOCAD .and. vp_final < VP_MIN_GOCAD) vp_final = VP_MIN_GOCAD
+ if(IMPOSE_MINIMUM_VP_GOCAD .and. vp_final < VP_MIN_GOCAD) vp_final = VP_MIN_GOCAD
! taper edges to make smooth transition between Hauksson and MR blocks
! get value from edge of medium-resolution block
! then use linear interpolation from edge of the model
- if(TAPER_GOCAD_TRANSITIONS) then
+ if(TAPER_GOCAD_TRANSITIONS) then
! x = xmin
- if(utm_x_eval < ORIG_X_GOCAD_MR + THICKNESS_TAPER_BLOCK_MR) then
- xmesh = ORIG_X_GOCAD_MR
- ymesh = utm_y_eval
- zmesh = z_eval
- if(HAUKSSON_REGIONAL_MODEL) then
- call hauksson_model(vp_hauksson,vs_hauksson,xmesh,ymesh,zmesh,vp_ref_hauksson,vs_dummy, MOHO_MAP_LUPEI)
- else
- call socal_model(doubling_index,rho_dummy,vp_ref_hauksson,vs_dummy)
- endif
- gamma_interp_x = (utm_x_eval - ORIG_X_GOCAD_MR) / THICKNESS_TAPER_BLOCK_MR
- vp_final = vp_ref_hauksson * (1. - gamma_interp_x) + vp_final * gamma_interp_x
+ if(utm_x_eval < ORIG_X_GOCAD_MR + THICKNESS_TAPER_BLOCK_MR) then
+ xmesh = ORIG_X_GOCAD_MR
+ ymesh = utm_y_eval
+ zmesh = z_eval
+ if(HAUKSSON_REGIONAL_MODEL) then
+ call hauksson_model(vp_hauksson,vs_hauksson,xmesh,ymesh,zmesh,vp_ref_hauksson,vs_dummy, MOHO_MAP_LUPEI)
+ else
+ call socal_model(doubling_index,rho_dummy,vp_ref_hauksson,vs_dummy)
+ endif
+ gamma_interp_x = (utm_x_eval - ORIG_X_GOCAD_MR) / THICKNESS_TAPER_BLOCK_MR
+ vp_final = vp_ref_hauksson * (1. - gamma_interp_x) + vp_final * gamma_interp_x
! x = xmax
- else if(utm_x_eval > END_X_GOCAD_MR - THICKNESS_TAPER_BLOCK_MR) then
- xmesh = END_X_GOCAD_MR
- ymesh = utm_y_eval
- zmesh = z_eval
- if(HAUKSSON_REGIONAL_MODEL) then
- call hauksson_model(vp_hauksson,vs_hauksson,xmesh,ymesh,zmesh,vp_ref_hauksson,vs_dummy, MOHO_MAP_LUPEI)
- else
- call socal_model(doubling_index,rho_dummy,vp_ref_hauksson,vs_dummy)
- endif
- gamma_interp_x = (utm_x_eval - (END_X_GOCAD_MR - THICKNESS_TAPER_BLOCK_MR)) / THICKNESS_TAPER_BLOCK_MR
- vp_final = vp_ref_hauksson * gamma_interp_x + vp_final * (1. - gamma_interp_x)
+ else if(utm_x_eval > END_X_GOCAD_MR - THICKNESS_TAPER_BLOCK_MR) then
+ xmesh = END_X_GOCAD_MR
+ ymesh = utm_y_eval
+ zmesh = z_eval
+ if(HAUKSSON_REGIONAL_MODEL) then
+ call hauksson_model(vp_hauksson,vs_hauksson,xmesh,ymesh,zmesh,vp_ref_hauksson,vs_dummy, MOHO_MAP_LUPEI)
+ else
+ call socal_model(doubling_index,rho_dummy,vp_ref_hauksson,vs_dummy)
+ endif
+ gamma_interp_x = (utm_x_eval - (END_X_GOCAD_MR - THICKNESS_TAPER_BLOCK_MR)) / THICKNESS_TAPER_BLOCK_MR
+ vp_final = vp_ref_hauksson * gamma_interp_x + vp_final * (1. - gamma_interp_x)
! y = ymin
- else if(utm_y_eval < ORIG_Y_GOCAD_MR + THICKNESS_TAPER_BLOCK_MR) then
- xmesh = utm_x_eval
- ymesh = ORIG_Y_GOCAD_MR
- zmesh = z_eval
- if(HAUKSSON_REGIONAL_MODEL) then
- call hauksson_model(vp_hauksson,vs_hauksson,xmesh,ymesh,zmesh,vp_ref_hauksson,vs_dummy, MOHO_MAP_LUPEI)
- else
- call socal_model(doubling_index,rho_dummy,vp_ref_hauksson,vs_dummy)
- endif
- gamma_interp_y = (utm_y_eval - ORIG_Y_GOCAD_MR) / THICKNESS_TAPER_BLOCK_MR
- vp_final = vp_ref_hauksson * (1. - gamma_interp_y) + vp_final * gamma_interp_y
+ else if(utm_y_eval < ORIG_Y_GOCAD_MR + THICKNESS_TAPER_BLOCK_MR) then
+ xmesh = utm_x_eval
+ ymesh = ORIG_Y_GOCAD_MR
+ zmesh = z_eval
+ if(HAUKSSON_REGIONAL_MODEL) then
+ call hauksson_model(vp_hauksson,vs_hauksson,xmesh,ymesh,zmesh,vp_ref_hauksson,vs_dummy, MOHO_MAP_LUPEI)
+ else
+ call socal_model(doubling_index,rho_dummy,vp_ref_hauksson,vs_dummy)
+ endif
+ gamma_interp_y = (utm_y_eval - ORIG_Y_GOCAD_MR) / THICKNESS_TAPER_BLOCK_MR
+ vp_final = vp_ref_hauksson * (1. - gamma_interp_y) + vp_final * gamma_interp_y
! y = ymax
- else if(utm_y_eval > END_Y_GOCAD_MR - THICKNESS_TAPER_BLOCK_MR) then
- xmesh = utm_x_eval
- ymesh = END_Y_GOCAD_MR
- zmesh = z_eval
- if(HAUKSSON_REGIONAL_MODEL) then
- call hauksson_model(vp_hauksson,vs_hauksson,xmesh,ymesh,zmesh,vp_ref_hauksson,vs_dummy, MOHO_MAP_LUPEI)
- else
- call socal_model(doubling_index,rho_dummy,vp_ref_hauksson,vs_dummy)
- endif
- gamma_interp_y = (utm_y_eval - (END_Y_GOCAD_MR - THICKNESS_TAPER_BLOCK_MR)) / THICKNESS_TAPER_BLOCK_MR
- vp_final = vp_ref_hauksson * gamma_interp_y + vp_final * (1. - gamma_interp_y)
+ else if(utm_y_eval > END_Y_GOCAD_MR - THICKNESS_TAPER_BLOCK_MR) then
+ xmesh = utm_x_eval
+ ymesh = END_Y_GOCAD_MR
+ zmesh = z_eval
+ if(HAUKSSON_REGIONAL_MODEL) then
+ call hauksson_model(vp_hauksson,vs_hauksson,xmesh,ymesh,zmesh,vp_ref_hauksson,vs_dummy, MOHO_MAP_LUPEI)
+ else
+ call socal_model(doubling_index,rho_dummy,vp_ref_hauksson,vs_dummy)
+ endif
+ gamma_interp_y = (utm_y_eval - (END_Y_GOCAD_MR - THICKNESS_TAPER_BLOCK_MR)) / THICKNESS_TAPER_BLOCK_MR
+ vp_final = vp_ref_hauksson * gamma_interp_y + vp_final * (1. - gamma_interp_y)
+
+ endif
+
+ endif
- endif
+! use linear variation of vp/vs ratio with depth, between 0. and 8.5 km
+ vp_vs_ratio = VP_VS_RATIO_GOCAD_BOTTOM + &
+ (VP_VS_RATIO_GOCAD_TOP - VP_VS_RATIO_GOCAD_BOTTOM) * &
+ (z_eval - (-8500.d0)) / (0.d0 - (-8500.d0))
- endif
+! make sure ratio remains in interval
+ if(vp_vs_ratio < VP_VS_RATIO_GOCAD_BOTTOM) vp_vs_ratio = VP_VS_RATIO_GOCAD_BOTTOM
+ if(vp_vs_ratio > VP_VS_RATIO_GOCAD_TOP) vp_vs_ratio = VP_VS_RATIO_GOCAD_TOP
+ vs_final = vp_final / vp_vs_ratio
+ call compute_rho_estimate(rho_final,vp_final)
+
+! fix for the smooth transition between HR and MR blocks : search vp_final in Hauksson's model if MR element is undefined
+ elseif(is_second_call .eqv. .true.) then
+
+ call hauksson_model(vp_hauksson,vs_hauksson,utm_x_eval,utm_y_eval,z_eval,vp_final,vs_final,MOHO_MAP_LUPEI)
+
! use linear variation of vp/vs ratio with depth, between 0. and 8.5 km
- vp_vs_ratio = VP_VS_RATIO_GOCAD_BOTTOM + &
+ vp_vs_ratio = VP_VS_RATIO_GOCAD_BOTTOM + &
(VP_VS_RATIO_GOCAD_TOP - VP_VS_RATIO_GOCAD_BOTTOM) * &
(z_eval - (-8500.d0)) / (0.d0 - (-8500.d0))
! make sure ratio remains in interval
- if(vp_vs_ratio < VP_VS_RATIO_GOCAD_BOTTOM) vp_vs_ratio = VP_VS_RATIO_GOCAD_BOTTOM
- if(vp_vs_ratio > VP_VS_RATIO_GOCAD_TOP) vp_vs_ratio = VP_VS_RATIO_GOCAD_TOP
+ if(vp_vs_ratio < VP_VS_RATIO_GOCAD_BOTTOM) vp_vs_ratio = VP_VS_RATIO_GOCAD_BOTTOM
+ if(vp_vs_ratio > VP_VS_RATIO_GOCAD_TOP) vp_vs_ratio = VP_VS_RATIO_GOCAD_TOP
+
+ vs_final = vp_final / vp_vs_ratio
+ call compute_rho_estimate(rho_final,vp_final)
- vs_final = vp_final / vp_vs_ratio
- call compute_rho_estimate(rho_final,vp_final)
+ endif
+
+ end subroutine interpolate_gocad_block_MR
- endif
-
- end subroutine interpolate_gocad_block_MR
-
More information about the CIG-COMMITS
mailing list