[cig-commits] [commit] devel: Modif coupling (186efbb)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Wed Oct 22 17:14:35 PDT 2014
Repository : https://github.com/geodynamics/specfem3d
On branch : devel
Link : https://github.com/geodynamics/specfem3d/compare/0479069f470e39347345ab3d4aff453320ac453e...57f2934cebe2e3d82d758afd2da6bf77d96c901d
>---------------------------------------------------------------
commit 186efbb9fd0392354d1106519aade5304221f670
Author: Clément Durochat <c.durochat at gmail.com>
Date: Fri Oct 17 17:13:19 2014 +0200
Modif coupling
>---------------------------------------------------------------
186efbb9fd0392354d1106519aade5304221f670
.../example_simple_small/DATA/Par_file | 6 ++--
src/generate_databases/get_model.f90 | 38 +++++++++++++++++-----
src/generate_databases/model_external_values.f90 | 31 ++++--------------
3 files changed, 40 insertions(+), 35 deletions(-)
diff --git a/EXAMPLES/DSM_FOR_SPECFEM3D/example_simple_small/DATA/Par_file b/EXAMPLES/DSM_FOR_SPECFEM3D/example_simple_small/DATA/Par_file
index be0f26c..f2f4d06 100644
--- a/EXAMPLES/DSM_FOR_SPECFEM3D/example_simple_small/DATA/Par_file
+++ b/EXAMPLES/DSM_FOR_SPECFEM3D/example_simple_small/DATA/Par_file
@@ -35,7 +35,7 @@ NGNOD = 8
# 1d_prem,1d_socal,1d_cascadia
# 3D models available are:
# aniso,external,gll,salton_trough,tomo
-MODEL = 1d_prem
+MODEL = external
# if you are using a SEP model (oil-industry format)
SEP_MODEL_DIRECTORY = ./DATA/my_SEP_model/
@@ -80,9 +80,9 @@ STACEY_INSTEAD_OF_FREE_SURFACE = .false.
CREATE_SHAKEMAP = .false.
MOVIE_SURFACE = .false.
MOVIE_TYPE = 1
-MOVIE_VOLUME = .true.
+MOVIE_VOLUME = .false.
SAVE_DISPLACEMENT = .false.
-USE_HIGHRES_FOR_MOVIES = .true.
+USE_HIGHRES_FOR_MOVIES = .false.
NTSTEP_BETWEEN_FRAMES = 80
HDUR_MOVIE = 0.0
diff --git a/src/generate_databases/get_model.f90 b/src/generate_databases/get_model.f90
index fa0483b..cbbeb1b 100644
--- a/src/generate_databases/get_model.f90
+++ b/src/generate_databases/get_model.f90
@@ -89,18 +89,17 @@
!! find the # layer where the middle of the element is located
if (COUPLE_WITH_EXTERNAL_CODE) then
- if (NGLL == 5) then
+ if( (NGLLX == 5) .and. (NGLLY == 5) .and. (NGLLZ == 5) ) then
! gets xyz coordinates of GLL point
iglob = ibool(3,3,3,1)
xmesh = xstore_dummy(iglob)
ymesh = ystore_dummy(iglob)
zmesh = zstore_dummy(iglob)
- else
+ call FindLayer(xmesh,ymesh,zmesh)
+ else
stop 'bad number of GLL points for coupling with DSM'
- endif
-
- call FindLayer(xmesh,ymesh,zmesh)
- end if
+ end if
+ endif
! ! Piero, read bedrock file
! in case, see file model_interface_bedrock.f90:
@@ -173,9 +172,9 @@
!! find the # layer where the middle of the element is located
if (COUPLE_WITH_EXTERNAL_CODE) then
- if (NGLL == 5) then
+ if( (NGLLX == 5) .and. (NGLLY == 5) .and. (NGLLZ == 5) ) then
if (i==3 .and. j==3 .and. k==3) call FindLayer(xmesh,ymesh,zmesh)
- else
+ else
stop 'bad number of GLL points for coupling with DSM'
endif
@@ -594,3 +593,26 @@
end select
end subroutine get_model_binaries
+
+!----------------------------------------------------------------
+
+ subroutine FindLayer(x,y,z)
+ use external_model
+ implicit none
+ integer il
+ double precision radius
+ double precision :: x,y,z
+ radius = dsqrt(x**2 + y**2 + (z+zref)**2) / 1000.d0
+
+ !write(124,*) 'RADIUS ',radius,x,y,z,z+zref,zref
+ il = 1
+ do while (radius .gt. zlayer(il).and.il.lt.nlayer)
+ il = il + 1
+ end do
+ il = il - 1
+ ilayer = il
+
+ !write(124,*) 'r, i : ',z,zref,radius, ilayer
+ end subroutine FindLayer
+
+!----------------------------------------------------------------
diff --git a/src/generate_databases/model_external_values.f90 b/src/generate_databases/model_external_values.f90
index 6dac64f..5cf6482 100644
--- a/src/generate_databases/model_external_values.f90
+++ b/src/generate_databases/model_external_values.f90
@@ -202,9 +202,8 @@
integer :: idomain_id
! local parameters
- double precision :: radius
- real(kind=CUSTOM_REAL) :: x,y,z
- real(kind=CUSTOM_REAL) :: xmin,xmax,ymin,ymax,zmin,zmax
+ double precision :: radius, x,y,z
+ real(kind=CUSTOM_REAL) :: xmin,xmax,ymin,ymax,zmin,zmax,x_target,y_target
real(kind=CUSTOM_REAL) :: depth
real(kind=CUSTOM_REAL) :: elevation,distmin
@@ -232,9 +231,11 @@
ymax = 134000._CUSTOM_REAL ! maxval(ystore_dummy)
zmin = 0._CUSTOM_REAL ! minval(zstore_dummy)
zmax = 60000._CUSTOM_REAL ! maxval(zstore_dummy)
+ x_target = x
+ y_target = y
! get approximate topography elevation at target coordinates from free surface
- call get_topo_elevation_free_closest(x,y,elevation,distmin, &
+ call get_topo_elevation_free_closest(x_target,y_target,elevation,distmin, &
nspec,nglob_dummy,ibool,xstore_dummy,ystore_dummy,zstore_dummy, &
num_free_surface_faces,free_surface_ispec,free_surface_ijk)
@@ -283,32 +284,14 @@
!! !! ================= VM VM CUSTOM SUBROUTINE FOR DSM COUPLING
!----------------------------------------------------------------
- subroutine FindLayer(x,y,z)
- use external_model
- implicit none
- integer il
- double precision radius,x,y,z
- radius = dsqrt(x**2 + y**2 + (z+zref)**2) / 1000.d0
-
- !write(124,*) 'RADIUS ',radius,x,y,z,z+zref,zref
- il = 1
- do while (radius .gt. zlayer(il).and.il.lt.nlayer)
- il = il + 1
- end do
- il = il - 1
- ilayer = il
-
- !write(124,*) 'r, i : ',z,zref,radius, ilayer
- end subroutine FindLayer
-
-!----------------------------------------------------------------
subroutine model_1D(x_eval,y_eval,z_eval, &
rho_final,vp_final,vs_final,r1)
use external_model
implicit none
- double precision r1,radius,x_eval,y_eval,z_eval
+ double precision r1,radius
double precision rho,vp,vs
+ double precision x_eval,y_eval,z_eval
real(kind=CUSTOM_REAL) rho_final,vp_final,vs_final
double precision Interpol,Xtol
More information about the CIG-COMMITS
mailing list