[cig-commits] r15877 - in seismo/3D/SPECFEM3D_SESAME/trunk: . EXAMPLES/homogeneous_halfspace EXAMPLES/layered_halfspace decompose_mesh_SCOTCH/OUTPUT_FILES
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Sun Oct 25 16:54:44 PDT 2009
Author: danielpeter
Date: 2009-10-25 16:54:42 -0700 (Sun, 25 Oct 2009)
New Revision: 15877
Modified:
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/block_mesh.py
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/cubit2specfem3d.py
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8-nodoubling.py
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8.py
seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/cubit2specfem3d.py
seismo/3D/SPECFEM3D_SESAME/trunk/Makefile.in
seismo/3D/SPECFEM3D_SESAME/trunk/aniso_model.f90
seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_elastic.f90
seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_no_Deville.f90
seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_with_Deville.f90
seismo/3D/SPECFEM3D_SESAME/trunk/constants.h.in
seismo/3D/SPECFEM3D_SESAME/trunk/create_movie_shakemap_AVS_DX_GMT.f90
seismo/3D/SPECFEM3D_SESAME/trunk/create_regions_mesh.f90
seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/OUTPUT_FILES/cubit2specfem3d.py
seismo/3D/SPECFEM3D_SESAME/trunk/detect_mesh_surfaces.f90
seismo/3D/SPECFEM3D_SESAME/trunk/generate_databases.f90
seismo/3D/SPECFEM3D_SESAME/trunk/initialize_simulation.f90
seismo/3D/SPECFEM3D_SESAME/trunk/iterate_time.f90
seismo/3D/SPECFEM3D_SESAME/trunk/prepare_timerun.f90
seismo/3D/SPECFEM3D_SESAME/trunk/read_mesh_databases.f90
seismo/3D/SPECFEM3D_SESAME/trunk/save_arrays_solver.f90
seismo/3D/SPECFEM3D_SESAME/trunk/save_header_file.f90
seismo/3D/SPECFEM3D_SESAME/trunk/setup_movie_meshes.f90
seismo/3D/SPECFEM3D_SESAME/trunk/specfem3D_par.f90
Log:
putting back anisotropy; updated cubit2specfem3d.py scripts; new velocity options for movies created with create_movie_shakemap_AVS_DX_GMT.f90
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/block_mesh.py
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/block_mesh.py 2009-10-25 23:06:34 UTC (rev 15876)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/block_mesh.py 2009-10-25 23:54:42 UTC (rev 15877)
@@ -28,12 +28,13 @@
#### Define material properties for the 3 volumes ################
cubit.cmd('#### DEFINE MATERIAL PROPERTIES #######################')
-cubit.cmd('block 1 attribute count 5')
+cubit.cmd('block 1 attribute count 6')
cubit.cmd('block 1 attribute index 1 1') # flag for material properties: 1 for 1. volume
-cubit.cmd('block 1 attribute index 2 2800') # vp
-cubit.cmd('block 1 attribute index 3 1500') # vs
-cubit.cmd('block 1 attribute index 4 2300') # rho
-cubit.cmd('block 1 attribute index 5 13') # Q flag (see constants.h: IATTENUATION_ ... )
+cubit.cmd('block 1 attribute index 2 2800') # vp
+cubit.cmd('block 1 attribute index 3 1500') # vs
+cubit.cmd('block 1 attribute index 4 2300') # rho
+cubit.cmd('block 1 attribute index 5 13') # Q flag (see constants.h: IATTENUATION_ ... )
+cubit.cmd('block 1 attribute index 6 0 ') # anisotropy_flag
cubit.cmd('export mesh "top.e" dimension 3 overwrite')
cubit.cmd('save as "meshing.cub" overwrite')
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/cubit2specfem3d.py
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/cubit2specfem3d.py 2009-10-25 23:06:34 UTC (rev 15876)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/homogeneous_halfspace/cubit2specfem3d.py 2009-10-25 23:54:42 UTC (rev 15877)
@@ -34,11 +34,13 @@
#or
# manually following the convention:
# - each material should have a block defined by name,flag of the material (integer),p velocity
-# (or the full description: name, flag, vp, vs, rho, Q ... if not present these last 3 parameters will be interpolated by module mat_parameter)
+# (or the full description: name, flag, vp, vs, rho, Q ... if not present these last 3 parameters will be
+# interpolated by module mat_parameter)
# - each mesh should have the block definition for the face on the free_surface (topography),
# the name of this block must be 'face_topo' or you can change the default name in mesh.topo defined in profile.
# - each mesh should have the block definition for the faces on the absorbing boundaries,
-# one block for each surface with x=Xmin,x=Xmax,y=Ymin,y=Ymax and z=bottom. The names of the blocks should contain the strings "xmin,xmax,ymin,ymax,bottom"
+# one block for each surface with x=Xmin,x=Xmax,y=Ymin,y=Ymax and z=bottom. The names of
+# the blocks should contain the strings "xmin,xmax,ymin,ymax,bottom"
#
#############################################################################
#RUN
@@ -79,14 +81,20 @@
# .....
# flag 'tomography' file_name #for interpolation with tomography
# .....
-# flag 'interface' file_name flag_for_the_gll_below_the_interface flag_for_the_gll_above_the_interface #for interpolation with interface
+# flag 'interface' file_name flag_for_the_gll_below_the_interface
+# flag_for_the_gll_above_the_interface #for interpolation with interface
#__________________________________________________________________________________________
##absname='absorbing_surface_file' -> this file contains all the face in all the absorbing boundaries
-##absname_local='absorbing_surface_file'+'_xmin' -> this file contains all the face in the absorbing boundary defined by x=Xmin
-##absname_local='absorbing_surface_file'+'_xmax' -> this file contains all the face in the absorbing boundary defined by x=Xmax
-##absname_local='absorbing_surface_file'+'_ymin' -> this file contains all the face in the absorbing boundary defined by y=Ymin
-##absname_local='absorbing_surface_file'+'_ymax' -> this file contains all the face in the absorbing boundary defined by y=Ymax
-##absname_local='absorbing_surface_file'+'_bottom' -> this file contains all the face in the absorbing boundary defined by z=bottom
+##absname_local='absorbing_surface_file'+'_xmin' -> this file contains all the face in the
+# absorbing boundary defined by x=Xmin
+##absname_local='absorbing_surface_file'+'_xmax' -> this file contains all the face in the
+# absorbing boundary defined by x=Xmax
+##absname_local='absorbing_surface_file'+'_ymin' -> this file contains all the face in the
+# absorbing boundary defined by y=Ymin
+##absname_local='absorbing_surface_file'+'_ymax' -> this file contains all the face in the
+# absorbing boundary defined by y=Ymax
+##absname_local='absorbing_surface_file'+'_bottom' -> this file contains all the face in the
+# absorbing boundary defined by z=bottom
# format:
# number of faces
# id_(element containg the face) id_node1_face id_node2_face id_node3_face id_node4_face
@@ -99,7 +107,8 @@
# id_(element containg the face) id_node1_face id_node2_face id_node3_face id_node4_face
#
#__________________________________________________________________________________________
-# it is possible save only one (or more) file singularly: for example if you want only the nodecoord_file call the module mesh.nodescoord_write(full path name)
+# it is possible save only one (or more) file singularly: for example if you want only the nodecoord_file
+# call the module mesh.nodescoord_write(full path name)
#
#############################################################################
@@ -120,7 +129,8 @@
def __repr__(self):
txt='Meshing for frequency up to '+str(self.frequency)+'Hz\n'
for surf,vp in zip(self.list_surf,self.list_vp):
- txt=txt+'surface '+str(surf)+', vp ='+str(vp)+' -> size '+str(self.freq2meshsize(vp)[0])+' -> dt '+str(self.freq2meshsize(vp)[0])+'\n'
+ txt=txt+'surface '+str(surf)+', vp ='+str(vp)+' -> size '+str(self.freq2meshsize(vp)[0])\
+ +' -> dt '+str(self.freq2meshsize(vp)[0])+'\n'
return txt
def freq2meshsize(self,vp):
velocity=vp*.5
@@ -182,8 +192,10 @@
"""Tools for the mesh
#########
dt,edge_dt,freq,edge_freq=seismic_resolution(edges,velocity,bins_d=None,bins_u=None,sidelist=None,ngll=5,np=8)
- Given the velocity of a list of edges, seismic_resolution provides the minimum Dt required for the stability condition (and the corrisponding edge).
- Furthermore, given the number of gll point in the element (ngll) and the number of GLL point for wavelength, it provide the maximum resolved frequency.
+ Given the velocity of a list of edges, seismic_resolution provides the minimum Dt
+ required for the stability condition (and the corrisponding edge).
+ Furthermore, given the number of gll point in the element (ngll) and the number
+ of GLL point for wavelength, it provide the maximum resolved frequency.
#########
length=edge_length(edge)
return the length of a edge
@@ -197,8 +209,10 @@
def seismic_resolution(self,edges,velocity,bins_d=None,bins_u=None,sidelist=None):
"""
dt,edge_dt,freq,edge_freq=seismic_resolution(edges,velocity,bins_d=None,bins_u=None,sidelist=None,ngll=5,np=8)
- Given the velocity of a list of edges, seismic_resolution provides the minimum Dt required for the stability condition (and the corrisponding edge).
- Furthermore, given the number of gll point in the element (ngll) and the number of GLL point for wavelength, it provide the maximum resolved frequency.
+ Given the velocity of a list of edges, seismic_resolution provides the minimum Dt
+ required for the stability condition (and the corrisponding edge).
+ Furthermore, given the number of gll point in the element (ngll) and the number
+ of GLL point for wavelength, it provide the maximum resolved frequency.
"""
ratiostore=1e10
dtstore=1e10
@@ -307,7 +321,8 @@
es=cubit.get_sub_elements("face", face, 1)
edges=edges+list(es)
edges=Set(edges)
- dtstore,edgedtstore,ratiostore,edgeratiostore=self.seismic_resolution(edges,velocity,bins_d,bins_u,sidelist)
+ dtstore,edgedtstore,ratiostore,edgeratiostore=self.seismic_resolution(edges,\
+ velocity,bins_d,bins_u,sidelist)
dt.append(dtstore)
ed_dt.append(edgedtstore)
r.append(ratiostore)
@@ -364,24 +379,28 @@
vs=None
rho=None
q=None
+ ani=None
if nattrib != 0:
flag=int(cubit.get_block_attribute_value(block,0))
if flag > 0 and nattrib >= 2:
- vel=cubit.get_block_attribute_value(block,1)
- if nattrib >= 3:
- vs=cubit.get_block_attribute_value(block,2)
- if nattrib >= 4:
- #density
- rho=cubit.get_block_attribute_value(block,3)
- if nattrib == 5:
- #Q_flag
- q=cubit.get_block_attribute_value(block,4)
- # for q to be valid: it must be an integer flag between 1 and 13
- # (see constants.h for IATTENUATION_SEDIMENT_40, etc. )
- if q < 0 or q > 13:
- print 'error, q flag invalid:', q
- print ' check with constants.h for IATTENUATION flags'
- break
+ vel=cubit.get_block_attribute_value(block,1)
+ if nattrib >= 3:
+ vs=cubit.get_block_attribute_value(block,2)
+ if nattrib >= 4:
+ #density
+ rho=cubit.get_block_attribute_value(block,3)
+ if nattrib >= 5:
+ #Q_flag
+ q=cubit.get_block_attribute_value(block,4)
+ # for q to be valid: it must be an integer flag between 1 and 13
+ # (see constants.h for IATTENUATION_SEDIMENT_40, etc. )
+ if q < 0 or q > 13:
+ print 'error, q flag invalid:', q
+ print ' check with constants.h for IATTENUATION flags'
+ break
+ if nattrib == 6:
+ #anisotropy_flag
+ ani=cubit.get_block_attribute_value(block,5)
elif flag < 0:
vel=name
attrib=cubit.get_block_attribute_value(block,1)
@@ -393,11 +412,11 @@
kind='tomography'
else:
flag=block
- vel,vs,rho,q=(name,0,0,0)
+ vel,vs,rho,q,ani=(name,0,0,0,0)
block_flag.append(int(flag))
block_mat.append(block)
if flag > 0:
- par=tuple([flag,vel,vs,rho,q])
+ par=tuple([flag,vel,vs,rho,q,ani])
elif flag < 0:
if kind=='interface':
par=tuple([flag,kind,name,flag_down,flag_up])
@@ -441,7 +460,8 @@
print bc
print topography
def mat_parameter(self,properties):
- #TODO: attenuation q .... where?
+ #TODO: material property acoustic/elastic/poroelastic ? .... where?
+ print "#material properties:"
print properties
flag=properties[0]
if flag > 0:
@@ -455,15 +475,17 @@
rho=(1.6612*vp-0.472*vp**2+0.0671*vp**3-0.0043*vp**4+0.000106*vp**4)*m2km
txt='%3i %20f %20f %20f %1i %1i\n' % (properties[0],rho,vel,vel/(3**.5),0,0)
elif type(vel) != str:
- #format nummaterials file: #material_id #rho #vp #vs #Q_flag #0
- txt='%3i %20f %20f %20f %2i %1i\n' % (properties[0],properties[3],properties[1],properties[2],properties[4],0)
+ #format nummaterials file: #material_id #rho #vp #vs #Q_flag #anisotropy_flag
+ txt='%3i %20f %20f %20f %2i %2i\n' % (properties[0],properties[3], \
+ properties[1],properties[2],properties[4],properties[5])
else:
txt='%3i %s \n' % (properties[0],properties[1])
elif flag < 0:
if properties[1] == 'tomography':
txt='%3i %s %s\n' % (properties[0],properties[1],properties[2])
elif properties[1] == 'interface':
- txt='%3i %s %s %1i %1i\n' % (properties[0],properties[1],properties[2],properties[3],properties[4])
+ txt='%3i %s %s %1i %1i\n' % (properties[0],properties[1],properties[2],\
+ properties[3],properties[4])
return txt
def nummaterial_write(self,nummaterial_name):
print 'Writing '+nummaterial_name+'.....'
@@ -542,7 +564,8 @@
#print f
nodes=cubit.get_connectivity('Face',f)
nodes_ok=self.normal_check(nodes,normal)
- txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],nodes_ok[1],nodes_ok[2],nodes_ok[3])
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],\
+ nodes_ok[1],nodes_ok[2],nodes_ok[3])
freehex.write(txt)
freehex.close()
print 'Ok'
@@ -606,9 +629,11 @@
nodes=cubit.get_connectivity('Face',f)
if not absflag:
nodes_ok=self.normal_check(nodes,normal)
- txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],nodes_ok[1],nodes_ok[2],nodes_ok[3])
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],\
+ nodes_ok[1],nodes_ok[2],nodes_ok[3])
else:
- txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],nodes[1],nodes[2],nodes[3])
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
+ nodes[1],nodes[2],nodes[3])
abshex_local.write(txt)
abshex_local.close()
print 'Ok'
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8-nodoubling.py
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8-nodoubling.py 2009-10-25 23:06:34 UTC (rev 15876)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8-nodoubling.py 2009-10-25 23:54:42 UTC (rev 15877)
@@ -50,7 +50,7 @@
cubit.cmd('imprint all')
# Meshing the volumes
-#elementsize = 1196.4
+#elementsize = 1196.4 #hi-resolution
#elementsize = 1500.0 # mid-resolution
elementsize = 3000.0 # low-resolution
@@ -73,26 +73,29 @@
#### Define material properties for the 3 volumes ################
cubit.cmd('#### DEFINE MATERIAL PROPERTIES #######################')
-cubit.cmd('block 1 attribute count 5')
+cubit.cmd('block 1 attribute count 6')
cubit.cmd('block 1 attribute index 1 1 ') # volume 1
cubit.cmd('block 1 attribute index 2 2800 ') # vp
cubit.cmd('block 1 attribute index 3 1500 ') # vs
cubit.cmd('block 1 attribute index 4 2300 ') # rho
cubit.cmd('block 1 attribute index 5 6 ') # Q_flag
+cubit.cmd('block 1 attribute index 6 1 ') # anisotropy_flag
-cubit.cmd('block 2 attribute count 5')
+cubit.cmd('block 2 attribute count 6')
cubit.cmd('block 2 attribute index 1 2 ') # volume 2
-cubit.cmd('block 2 attribute index 2 7500 ')
-cubit.cmd('block 2 attribute index 3 4300 ')
-cubit.cmd('block 2 attribute index 4 3200 ')
-cubit.cmd('block 2 attribute index 5 6')
+cubit.cmd('block 2 attribute index 2 7500 ') # vp
+cubit.cmd('block 2 attribute index 3 4300 ') # vs
+cubit.cmd('block 2 attribute index 4 3200 ') # rho
+cubit.cmd('block 2 attribute index 5 6') # Q_flag
+cubit.cmd('block 2 attribute index 6 0 ') # anisotropy_flag
-cubit.cmd('block 3 attribute count 5')
+cubit.cmd('block 3 attribute count 6')
cubit.cmd('block 3 attribute index 1 2 ') # same material properties as for volume 2
cubit.cmd('block 3 attribute index 2 7500 ')
cubit.cmd('block 3 attribute index 3 4300 ')
cubit.cmd('block 3 attribute index 4 3200 ')
cubit.cmd('block 3 attribute index 5 6')
+cubit.cmd('block 3 attribute index 6 0')
cubit.cmd('export mesh "top.e" dimension 3 overwrite')
cubit.cmd('save as "meshing.cub" overwrite')
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8.py
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8.py 2009-10-25 23:06:34 UTC (rev 15876)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/2lay_mesh_boundary_fig8.py 2009-10-25 23:54:42 UTC (rev 15877)
@@ -70,26 +70,29 @@
#### Define material properties for the 3 volumes ################
cubit.cmd('#### DEFINE MATERIAL PROPERTIES #######################')
-cubit.cmd('block 1 attribute count 5')
+cubit.cmd('block 1 attribute count 6')
cubit.cmd('block 1 attribute index 1 1 ') # volume 1
cubit.cmd('block 1 attribute index 2 2800 ') # vp
cubit.cmd('block 1 attribute index 3 1500 ') # vs
cubit.cmd('block 1 attribute index 4 2300 ') # rho
cubit.cmd('block 1 attribute index 5 6 ') # Q_flag
+cubit.cmd('block 1 attribute index 6 0 ') # anisotropy_flag
-cubit.cmd('block 2 attribute count 5')
+cubit.cmd('block 2 attribute count 6')
cubit.cmd('block 2 attribute index 1 2 ') # volume 2
cubit.cmd('block 2 attribute index 2 7500 ')
cubit.cmd('block 2 attribute index 3 4300 ')
cubit.cmd('block 2 attribute index 4 3200 ')
cubit.cmd('block 2 attribute index 5 6 ')
+cubit.cmd('block 2 attribute index 6 0 ') # anisotropy_flag
-cubit.cmd('block 3 attribute count 5')
+cubit.cmd('block 3 attribute count 6')
cubit.cmd('block 3 attribute index 1 2 ') # same properties as for volume 2
cubit.cmd('block 3 attribute index 2 7500 ')
cubit.cmd('block 3 attribute index 3 4300 ')
cubit.cmd('block 3 attribute index 4 3200 ')
cubit.cmd('block 3 attribute index 5 6 ')
+cubit.cmd('block 3 attribute index 6 0 ') # anisotropy_flag
cubit.cmd('export mesh "top.e" dimension 3 overwrite')
cubit.cmd('save as "meshing.cub" overwrite')
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/cubit2specfem3d.py
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/cubit2specfem3d.py 2009-10-25 23:06:34 UTC (rev 15876)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/EXAMPLES/layered_halfspace/cubit2specfem3d.py 2009-10-25 23:54:42 UTC (rev 15877)
@@ -34,11 +34,13 @@
#or
# manually following the convention:
# - each material should have a block defined by name,flag of the material (integer),p velocity
-# (or the full description: name, flag, vp, vs, rho, Q ... if not present these last 3 parameters will be interpolated by module mat_parameter)
+# (or the full description: name, flag, vp, vs, rho, Q ... if not present these last 3 parameters will be
+# interpolated by module mat_parameter)
# - each mesh should have the block definition for the face on the free_surface (topography),
# the name of this block must be 'face_topo' or you can change the default name in mesh.topo defined in profile.
# - each mesh should have the block definition for the faces on the absorbing boundaries,
-# one block for each surface with x=Xmin,x=Xmax,y=Ymin,y=Ymax and z=bottom. The names of the blocks should contain the strings "xmin,xmax,ymin,ymax,bottom"
+# one block for each surface with x=Xmin,x=Xmax,y=Ymin,y=Ymax and z=bottom. The names of
+# the blocks should contain the strings "xmin,xmax,ymin,ymax,bottom"
#
#############################################################################
#RUN
@@ -79,14 +81,20 @@
# .....
# flag 'tomography' file_name #for interpolation with tomography
# .....
-# flag 'interface' file_name flag_for_the_gll_below_the_interface flag_for_the_gll_above_the_interface #for interpolation with interface
+# flag 'interface' file_name flag_for_the_gll_below_the_interface
+# flag_for_the_gll_above_the_interface #for interpolation with interface
#__________________________________________________________________________________________
##absname='absorbing_surface_file' -> this file contains all the face in all the absorbing boundaries
-##absname_local='absorbing_surface_file'+'_xmin' -> this file contains all the face in the absorbing boundary defined by x=Xmin
-##absname_local='absorbing_surface_file'+'_xmax' -> this file contains all the face in the absorbing boundary defined by x=Xmax
-##absname_local='absorbing_surface_file'+'_ymin' -> this file contains all the face in the absorbing boundary defined by y=Ymin
-##absname_local='absorbing_surface_file'+'_ymax' -> this file contains all the face in the absorbing boundary defined by y=Ymax
-##absname_local='absorbing_surface_file'+'_bottom' -> this file contains all the face in the absorbing boundary defined by z=bottom
+##absname_local='absorbing_surface_file'+'_xmin' -> this file contains all the face in the
+# absorbing boundary defined by x=Xmin
+##absname_local='absorbing_surface_file'+'_xmax' -> this file contains all the face in the
+# absorbing boundary defined by x=Xmax
+##absname_local='absorbing_surface_file'+'_ymin' -> this file contains all the face in the
+# absorbing boundary defined by y=Ymin
+##absname_local='absorbing_surface_file'+'_ymax' -> this file contains all the face in the
+# absorbing boundary defined by y=Ymax
+##absname_local='absorbing_surface_file'+'_bottom' -> this file contains all the face in the
+# absorbing boundary defined by z=bottom
# format:
# number of faces
# id_(element containg the face) id_node1_face id_node2_face id_node3_face id_node4_face
@@ -99,7 +107,8 @@
# id_(element containg the face) id_node1_face id_node2_face id_node3_face id_node4_face
#
#__________________________________________________________________________________________
-# it is possible save only one (or more) file singularly: for example if you want only the nodecoord_file call the module mesh.nodescoord_write(full path name)
+# it is possible save only one (or more) file singularly: for example if you want only the nodecoord_file
+# call the module mesh.nodescoord_write(full path name)
#
#############################################################################
@@ -120,7 +129,8 @@
def __repr__(self):
txt='Meshing for frequency up to '+str(self.frequency)+'Hz\n'
for surf,vp in zip(self.list_surf,self.list_vp):
- txt=txt+'surface '+str(surf)+', vp ='+str(vp)+' -> size '+str(self.freq2meshsize(vp)[0])+' -> dt '+str(self.freq2meshsize(vp)[0])+'\n'
+ txt=txt+'surface '+str(surf)+', vp ='+str(vp)+' -> size '+str(self.freq2meshsize(vp)[0])\
+ +' -> dt '+str(self.freq2meshsize(vp)[0])+'\n'
return txt
def freq2meshsize(self,vp):
velocity=vp*.5
@@ -182,8 +192,10 @@
"""Tools for the mesh
#########
dt,edge_dt,freq,edge_freq=seismic_resolution(edges,velocity,bins_d=None,bins_u=None,sidelist=None,ngll=5,np=8)
- Given the velocity of a list of edges, seismic_resolution provides the minimum Dt required for the stability condition (and the corrisponding edge).
- Furthermore, given the number of gll point in the element (ngll) and the number of GLL point for wavelength, it provide the maximum resolved frequency.
+ Given the velocity of a list of edges, seismic_resolution provides the minimum Dt
+ required for the stability condition (and the corrisponding edge).
+ Furthermore, given the number of gll point in the element (ngll) and the number
+ of GLL point for wavelength, it provide the maximum resolved frequency.
#########
length=edge_length(edge)
return the length of a edge
@@ -197,8 +209,10 @@
def seismic_resolution(self,edges,velocity,bins_d=None,bins_u=None,sidelist=None):
"""
dt,edge_dt,freq,edge_freq=seismic_resolution(edges,velocity,bins_d=None,bins_u=None,sidelist=None,ngll=5,np=8)
- Given the velocity of a list of edges, seismic_resolution provides the minimum Dt required for the stability condition (and the corrisponding edge).
- Furthermore, given the number of gll point in the element (ngll) and the number of GLL point for wavelength, it provide the maximum resolved frequency.
+ Given the velocity of a list of edges, seismic_resolution provides the minimum Dt
+ required for the stability condition (and the corrisponding edge).
+ Furthermore, given the number of gll point in the element (ngll) and the number
+ of GLL point for wavelength, it provide the maximum resolved frequency.
"""
ratiostore=1e10
dtstore=1e10
@@ -307,7 +321,8 @@
es=cubit.get_sub_elements("face", face, 1)
edges=edges+list(es)
edges=Set(edges)
- dtstore,edgedtstore,ratiostore,edgeratiostore=self.seismic_resolution(edges,velocity,bins_d,bins_u,sidelist)
+ dtstore,edgedtstore,ratiostore,edgeratiostore=self.seismic_resolution(edges,\
+ velocity,bins_d,bins_u,sidelist)
dt.append(dtstore)
ed_dt.append(edgedtstore)
r.append(ratiostore)
@@ -364,24 +379,28 @@
vs=None
rho=None
q=None
+ ani=None
if nattrib != 0:
flag=int(cubit.get_block_attribute_value(block,0))
if flag > 0 and nattrib >= 2:
- vel=cubit.get_block_attribute_value(block,1)
- if nattrib >= 3:
- vs=cubit.get_block_attribute_value(block,2)
- if nattrib >= 4:
- #density
- rho=cubit.get_block_attribute_value(block,3)
- if nattrib == 5:
- #Q_flag
- q=cubit.get_block_attribute_value(block,4)
- # for q to be valid: it must be an integer flag between 1 and 13
- # (see constants.h for IATTENUATION_SEDIMENT_40, etc. )
- if q < 0 or q > 13:
- print 'error, q flag invalid:', q
- print ' check with constants.h for IATTENUATION flags'
- break
+ vel=cubit.get_block_attribute_value(block,1)
+ if nattrib >= 3:
+ vs=cubit.get_block_attribute_value(block,2)
+ if nattrib >= 4:
+ #density
+ rho=cubit.get_block_attribute_value(block,3)
+ if nattrib >= 5:
+ #Q_flag
+ q=cubit.get_block_attribute_value(block,4)
+ # for q to be valid: it must be an integer flag between 1 and 13
+ # (see constants.h for IATTENUATION_SEDIMENT_40, etc. )
+ if q < 0 or q > 13:
+ print 'error, q flag invalid:', q
+ print ' check with constants.h for IATTENUATION flags'
+ break
+ if nattrib == 6:
+ #anisotropy_flag
+ ani=cubit.get_block_attribute_value(block,5)
elif flag < 0:
vel=name
attrib=cubit.get_block_attribute_value(block,1)
@@ -393,11 +412,11 @@
kind='tomography'
else:
flag=block
- vel,vs,rho,q=(name,0,0,0)
+ vel,vs,rho,q,ani=(name,0,0,0,0)
block_flag.append(int(flag))
block_mat.append(block)
if flag > 0:
- par=tuple([flag,vel,vs,rho,q])
+ par=tuple([flag,vel,vs,rho,q,ani])
elif flag < 0:
if kind=='interface':
par=tuple([flag,kind,name,flag_down,flag_up])
@@ -441,7 +460,8 @@
print bc
print topography
def mat_parameter(self,properties):
- #TODO: attenuation q .... where?
+ #TODO: material property acoustic/elastic/poroelastic ? .... where?
+ print "#material properties:"
print properties
flag=properties[0]
if flag > 0:
@@ -455,15 +475,17 @@
rho=(1.6612*vp-0.472*vp**2+0.0671*vp**3-0.0043*vp**4+0.000106*vp**4)*m2km
txt='%3i %20f %20f %20f %1i %1i\n' % (properties[0],rho,vel,vel/(3**.5),0,0)
elif type(vel) != str:
- #format nummaterials file: #material_id #rho #vp #vs #Q_flag #0
- txt='%3i %20f %20f %20f %2i %1i\n' % (properties[0],properties[3],properties[1],properties[2],properties[4],0)
+ #format nummaterials file: #material_id #rho #vp #vs #Q_flag #anisotropy_flag
+ txt='%3i %20f %20f %20f %2i %2i\n' % (properties[0],properties[3], \
+ properties[1],properties[2],properties[4],properties[5])
else:
txt='%3i %s \n' % (properties[0],properties[1])
elif flag < 0:
if properties[1] == 'tomography':
txt='%3i %s %s\n' % (properties[0],properties[1],properties[2])
elif properties[1] == 'interface':
- txt='%3i %s %s %1i %1i\n' % (properties[0],properties[1],properties[2],properties[3],properties[4])
+ txt='%3i %s %s %1i %1i\n' % (properties[0],properties[1],properties[2],\
+ properties[3],properties[4])
return txt
def nummaterial_write(self,nummaterial_name):
print 'Writing '+nummaterial_name+'.....'
@@ -542,7 +564,8 @@
#print f
nodes=cubit.get_connectivity('Face',f)
nodes_ok=self.normal_check(nodes,normal)
- txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],nodes_ok[1],nodes_ok[2],nodes_ok[3])
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],\
+ nodes_ok[1],nodes_ok[2],nodes_ok[3])
freehex.write(txt)
freehex.close()
print 'Ok'
@@ -606,9 +629,11 @@
nodes=cubit.get_connectivity('Face',f)
if not absflag:
nodes_ok=self.normal_check(nodes,normal)
- txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],nodes_ok[1],nodes_ok[2],nodes_ok[3])
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],\
+ nodes_ok[1],nodes_ok[2],nodes_ok[3])
else:
- txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],nodes[1],nodes[2],nodes[3])
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
+ nodes[1],nodes[2],nodes[3])
abshex_local.write(txt)
abshex_local.close()
print 'Ok'
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/Makefile.in 2009-10-25 23:06:34 UTC (rev 15876)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/Makefile.in 2009-10-25 23:54:42 UTC (rev 15877)
@@ -57,6 +57,7 @@
O = obj
libspecfem_a_OBJECTS = \
+ $O/aniso_model.o \
$O/assemble_MPI_scalar.o \
$O/calc_jacobian.o \
$O/check_mesh_resolution.o \
@@ -108,7 +109,6 @@
###
###--obsolete routines
###
-# $O/aniso_model.o \
# $O/compute_rho_estimate.o \
# $O/define_subregions.o \
# $O/define_subregions_heuristic.o \
@@ -308,6 +308,9 @@
### serial compilation without optimization
###
+$O/aniso_model.o: constants.h aniso_model.f90
+ ${FCCOMPILE_CHECK} -c -o $O/aniso_model.o aniso_model.f90
+
$O/serial.o: constants.h exit_mpi.f90
${FCCOMPILE_CHECK} -c -o $O/serial.o serial.f90
@@ -523,10 +526,6 @@
# ${FCCOMPILE_CHECK} -c -o $O/read_moho_map.o read_moho_map.f90
#--obsolete
-#$O/aniso_model.o: constants.h aniso_model.f90
-# ${FCCOMPILE_CHECK} -c -o $O/aniso_model.o aniso_model.f90
-
-#--obsolete
#$O/compute_rho_estimate.o: constants.h compute_rho_estimate.f90
# ${FCCOMPILE_CHECK} -c -o $O/compute_rho_estimate.o compute_rho_estimate.f90
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/aniso_model.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/aniso_model.f90 2009-10-25 23:06:34 UTC (rev 15876)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/aniso_model.f90 2009-10-25 23:54:42 UTC (rev 15877)
@@ -29,109 +29,128 @@
! anisotropic models.
!=====================================================================
- subroutine aniso_model(idoubling,zmesh,rho,vp,vs,c11,c12,c13,c14,c15,c16, &
+ subroutine aniso_model(iflag_aniso,rho,vp,vs,c11,c12,c13,c14,c15,c16, &
c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66)
implicit none
include "constants.h"
- include "constants_gocad.h"
+
+! include "constants_gocad.h"
!------------------------------------------------------------------------------
! for anisotropy simulations in a halfspace model
! only related to body waves
! one-zeta term
- double precision, parameter :: FACTOR_CS1p_A = 0.01d0
- double precision, parameter :: FACTOR_CS1sv_A = 0.0d0
- double precision, parameter :: FACTOR_CS1sh_N = 0.d0
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_CS1p_A = 0.2_CUSTOM_REAL
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_CS1sv_A = 0.0_CUSTOM_REAL
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_CS1sh_N = 0._CUSTOM_REAL
! three-zeta term
- double precision, parameter :: FACTOR_CS3_L = 0.0d0
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_CS3_L = 0.0_CUSTOM_REAL
! Relative to Love wave
! four-zeta term
- double precision, parameter :: FACTOR_N = 0.d0
- double precision, parameter :: FACTOR_E_N = 0.d0
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_N = 0._CUSTOM_REAL
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_E_N = 0._CUSTOM_REAL
! Relative to Rayleigh wave
! two-zeta term
- double precision, parameter :: FACTOR_A = 0.d0
- double precision, parameter :: FACTOR_C = 0.d0
- double precision, parameter :: FACTOR_F = 0.d0
- double precision, parameter :: FACTOR_H_F = 0.d0
- double precision, parameter :: FACTOR_B_A = 0.d0
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_A = 0._CUSTOM_REAL
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_C = 0._CUSTOM_REAL
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_F = 0._CUSTOM_REAL
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_H_F = 0._CUSTOM_REAL
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_B_A = 0._CUSTOM_REAL
! Relative to both Love wave and Rayleigh wave
! two-zeta term
- double precision, parameter :: FACTOR_L = 0.d0
- double precision, parameter :: FACTOR_G_L = 0.d0
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_L = 0._CUSTOM_REAL
+ real(kind=CUSTOM_REAL), parameter :: FACTOR_G_L = 0._CUSTOM_REAL
!------------------------------------------------------------------------------
- integer idoubling
-
- double precision zmesh
- double precision rho,vp,vs
- double precision vpv,vph,vsv,vsh,eta_aniso
- double precision aa,cc,nn,ll,ff
- double precision A,C,F,AL,AN,Bc,Bs,Gc,Gs,Hc,Hs,Ec,Es,C1p,C1sv,C1sh,C3,S1p,S1sv,S1sh,S3
- double precision c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,c33,c34,c35,c36, &
+ !integer idoubling
+ integer iflag_aniso
+
+ !real(kind=CUSTOM_REAL) zmesh
+ real(kind=CUSTOM_REAL) rho,vp,vs
+ real(kind=CUSTOM_REAL) c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,c33,c34,c35,c36, &
c44,c45,c46,c55,c56,c66
- double precision d11,d12,d13,d14,d15,d16,d22,d23,d24,d25,d26,d33,d34,d35,d36, &
+
+! local parameters
+ real(kind=CUSTOM_REAL) vpv,vph,vsv,vsh,eta_aniso
+ real(kind=CUSTOM_REAL) aa,cc,nn,ll,ff
+ real(kind=CUSTOM_REAL) A,C,F,AL,AN,Bc,Bs,Gc,Gs,Hc,Hs,Ec,Es,C1p,C1sv,C1sh,C3,S1p,S1sv,S1sh,S3
+ real(kind=CUSTOM_REAL) d11,d12,d13,d14,d15,d16,d22,d23,d24,d25,d26,d33,d34,d35,d36, &
d44,d45,d46,d55,d56,d66
-! implement the background model
- if(idoubling == IFLAG_HALFSPACE_MOHO) then
- vp=7.8d0
- vs=4.5d0
- rho=3.0d0
- vph = vp
- vpv = vp
- vsh = vs
- vsv = vs
- eta_aniso = 1.0d0
+! not anymore, takes vp,vs rho from given isotropic input model...
+!
+!! implement the background model
+! if(iflag_aniso == IANISOTROPY_HALFSPACE_MOHO) then
+! vp=7.8_CUSTOM_REAL
+! vs=4.5_CUSTOM_REAL
+! rho=3.0_CUSTOM_REAL
+! vph = vp
+! vpv = vp
+! vsh = vs
+! vsv = vs
+! eta_aniso = 1.0_CUSTOM_REAL
+!
+! else if(iflag_aniso == IANISOTROPY_MOHO_16km) then
+! vp=7.8_CUSTOM_REAL
+! vs=4.5_CUSTOM_REAL
+! rho=3.0_CUSTOM_REAL
+! vph = vp
+! vpv = vp
+! vsh = vs
+! vsv = vs
+! eta_aniso = 1.0_CUSTOM_REAL
+!
+! else if(zmesh >= DEPTH_5p5km_SOCAL) then
+! vp=7.8_CUSTOM_REAL
+! vs=4.5_CUSTOM_REAL
+! rho=3.0_CUSTOM_REAL
+! vph = vp
+! vpv = vp
+! vsh = vs
+! vsv = vs
+! eta_aniso = 1.0_CUSTOM_REAL
+!
+! else
+! vp=7.8_CUSTOM_REAL
+! vs=4.5_CUSTOM_REAL
+! rho=3.0_CUSTOM_REAL
+! vph = vp
+! vpv = vp
+! vsh = vs
+! vsv = vs
+! eta_aniso = 1.0_CUSTOM_REAL
+!
+! endif
- else if(idoubling == IFLAG_MOHO_16km) then
- vp=7.8d0
- vs=4.5d0
- rho=3.0d0
- vph = vp
- vpv = vp
- vsh = vs
- vsv = vs
- eta_aniso = 1.0d0
+! scale to standard units
+! vp = vp * 1000._CUSTOM_REAL
+! vs = vs * 1000._CUSTOM_REAL
+! vph = vph * 1000._CUSTOM_REAL
+! vpv = vpv * 1000._CUSTOM_REAL
+! vsh = vsh * 1000._CUSTOM_REAL
+! vsv = vsv * 1000._CUSTOM_REAL
+! rho = rho * 1000._CUSTOM_REAL
- else if(zmesh >= DEPTH_5p5km_SOCAL) then
- vp=7.8d0
- vs=4.5d0
- rho=3.0d0
- vph = vp
- vpv = vp
- vsh = vs
- vsv = vs
- eta_aniso = 1.0d0
+! assumes vp,vs given in m/s, rho in kg/m**3
+ vph = vp
+ vpv = vp
+ vsh = vs
+ vsv = vs
+ eta_aniso = 1.0_CUSTOM_REAL
- else
- vp=7.8d0
- vs=4.5d0
- rho=3.0d0
- vph = vp
- vpv = vp
- vsh = vs
- vsv = vs
- eta_aniso = 1.0d0
+! see for example:
+!
+! M. Chen & J. Tromp, 2006. Theoretical & numerical investigations
+! of global and regional seismic wave propagation in weakly anisotropic earth models,
+! GJI, 168, 1130-1152.
- endif
-
-! scale to standard units
- vp = vp * 1000.d0
- vs = vs * 1000.d0
- vph = vph * 1000.d0
- vpv = vpv * 1000.d0
- vsh = vsh * 1000.d0
- vsv = vsv * 1000.d0
- rho = rho * 1000.d0
-
aa = rho*vph*vph
cc = rho*vpv*vpv
nn = rho*vsh*vsh
@@ -140,28 +159,67 @@
! Add anisotropic perturbation in the whole halfspace
! You can also add different perturbations to different layers
- A = aa*(1.0d0 + FACTOR_A)
- C = cc*(1.0d0 + FACTOR_C)
- AN = nn*(1.0d0 + FACTOR_N)
- AL = ll*(1.0d0 + FACTOR_L)
- F = ff*(1.0d0 + FACTOR_F)
- C1p = FACTOR_CS1p_A*aa
- S1p = 0.d0
- C1sv = FACTOR_CS1sv_A*aa
- S1sv = 0.d0
- C1sh = FACTOR_CS1sh_N*nn
- S1sh = 0.d0
- Gc = FACTOR_G_L*ll
- Gs = 0.d0
- Bc = FACTOR_B_A*aa
- Bs = 0.d0
- Hc = FACTOR_H_F*ff
- Hs = 0.d0
- C3 = FACTOR_CS3_L*ll
- S3 = 0.d0
- Ec = FACTOR_E_N*nn
- Es = 0.d0
+! no anisotropic perturbation
+ if( iflag_aniso <= 0 ) then
+ A = aa
+ C = cc
+ AN = nn
+ AL = ll
+ F = ff
+ C1p = 0.0
+ C1sv = 0.0
+ C1sh = 0.0
+ Gc = 0.0
+ Bc = 0.0
+ Hc = 0.0
+ C3 = 0.0
+ Ec = 0.0
+ endif
+
+! perturbation model 1
+ if( iflag_aniso == IANISOTROPY_MODEL1 ) then
+ A = aa*(1.0_CUSTOM_REAL + FACTOR_A)
+ C = cc*(1.0_CUSTOM_REAL + FACTOR_C)
+ AN = nn*(1.0_CUSTOM_REAL + FACTOR_N)
+ AL = ll*(1.0_CUSTOM_REAL + FACTOR_L)
+ F = ff*(1.0_CUSTOM_REAL + FACTOR_F)
+ C1p = FACTOR_CS1p_A*aa
+ C1sv = FACTOR_CS1sv_A*aa
+ C1sh = FACTOR_CS1sh_N*nn
+ Gc = FACTOR_G_L*ll
+ Bc = FACTOR_B_A*aa
+ Hc = FACTOR_H_F*ff
+ C3 = FACTOR_CS3_L*ll
+ Ec = FACTOR_E_N*nn
+ endif
+
+! perturbation model 2
+ if( iflag_aniso == IANISOTROPY_MODEL2 ) then
+ A = aa*(1.0_CUSTOM_REAL + FACTOR_A + 0.1)
+ C = cc*(1.0_CUSTOM_REAL + FACTOR_C + 0.1)
+ AN = nn*(1.0_CUSTOM_REAL + FACTOR_N + 0.1)
+ AL = ll*(1.0_CUSTOM_REAL + FACTOR_L + 0.1)
+ F = ff*(1.0_CUSTOM_REAL + FACTOR_F + 0.1)
+ C1p = FACTOR_CS1p_A*aa
+ C1sv = FACTOR_CS1sv_A*aa
+ C1sh = FACTOR_CS1sh_N*nn
+ Gc = FACTOR_G_L*ll
+ Bc = FACTOR_B_A*aa
+ Hc = FACTOR_H_F*ff
+ C3 = FACTOR_CS3_L*ll
+ Ec = FACTOR_E_N*nn
+ endif
+
+ S1p = 0._CUSTOM_REAL
+ S1sv = 0._CUSTOM_REAL
+ S1sh = 0._CUSTOM_REAL
+ Gs = 0._CUSTOM_REAL
+ Bs = 0._CUSTOM_REAL
+ Hs = 0._CUSTOM_REAL
+ S3 = 0._CUSTOM_REAL
+ Es = 0._CUSTOM_REAL
+
! The mapping from the elastic coefficients to the elastic tensor elements
! in the local Cartesian coordinate system (classical geographic) used in the
! global code (1---South, 2---East, 3---up)
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_elastic.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_elastic.f90 2009-10-25 23:06:34 UTC (rev 15876)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_elastic.f90 2009-10-25 23:54:42 UTC (rev 15877)
@@ -59,8 +59,14 @@
NSPEC_ATTENUATION_AB,R_xx,R_yy,R_xy,R_xz,R_yz, &
epsilondev_xx,epsilondev_yy,epsilondev_xy, &
epsilondev_xz,epsilondev_yz,iflag_attenuation_store, &
- rho_vs )
+ rho_vs, &
+ ANISOTROPY,NSPEC_ANISO, &
+ c11store,c12store,c13store,c14store,c15store,c16store,&
+ c22store,c23store,c24store,c25store,c26store,c33store,&
+ c34store,c35store,c36store,c44store,c45store,c46store,&
+ c55store,c56store,c66store )
+
!call compute_forces_with_Deville( phase_is_inner ,NSPEC_AB,NGLOB_AB,&
! ATTENUATION,USE_OLSEN_ATTENUATION,displ,accel,&
! xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
@@ -90,7 +96,12 @@
NSPEC_ATTENUATION_AB,R_xx,R_yy,R_xy,R_xz,R_yz, &
epsilondev_xx,epsilondev_yy,epsilondev_xy,&
epsilondev_xz,epsilondev_yz,iflag_attenuation_store,&
- rho_vs)
+ rho_vs, &
+ ANISOTROPY,NSPEC_ANISO, &
+ c11store,c12store,c13store,c14store,c15store,c16store,&
+ c22store,c23store,c24store,c25store,c26store,c33store,&
+ c34store,c35store,c36store,c44store,c45store,c46store,&
+ c55store,c56store,c66store)
endif
! adds elastic absorbing boundary term to acceleration (Stacey conditions)
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_no_Deville.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_no_Deville.f90 2009-10-25 23:06:34 UTC (rev 15876)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_no_Deville.f90 2009-10-25 23:54:42 UTC (rev 15877)
@@ -36,7 +36,12 @@
NSPEC_ATTENUATION_AB,R_xx,R_yy,R_xy,R_xz,R_yz, &
epsilondev_xx,epsilondev_yy,epsilondev_xy,&
epsilondev_xz,epsilondev_yz,iflag_attenuation_store,&
- rho_vs)
+ rho_vs,&
+ ANISOTROPY,NSPEC_ANISO, &
+ c11store,c12store,c13store,c14store,c15store,c16store,&
+ c22store,c23store,c24store,c25store,c26store,c33store,&
+ c34store,c35store,c36store,c44store,c45store,c46store,&
+ c55store,c56store,c66store)
! NSOURCES,myrank,islice_selected_source,&
! ispec_selected_source,xi_source,eta_source,&
@@ -84,6 +89,15 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: rho_vs
+! anisotropy
+ logical :: ANISOTROPY
+ integer :: NSPEC_ANISO
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO) :: &
+ c11store,c12store,c13store,c14store,c15store,c16store, &
+ c22store,c23store,c24store,c25store,c26store,c33store, &
+ c34store,c35store,c36store,c44store,c45store,c46store, &
+ c55store,c56store,c66store
+
! source
! integer :: NSOURCES,myrank,it
! integer, dimension(NSOURCES) :: islice_selected_source,ispec_selected_source
@@ -119,6 +133,10 @@
real(kind=CUSTOM_REAL) lambdal,mul,lambdalplus2mul
real(kind=CUSTOM_REAL) kappal
+! local anisotropy parameters
+ real(kind=CUSTOM_REAL) c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,&
+ c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
+
! local attenuation parameters
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: epsilondev_xx_loc, &
epsilondev_yy_loc, epsilondev_xy_loc, epsilondev_xz_loc, epsilondev_yz_loc
@@ -241,18 +259,96 @@
endif
- lambdalplus2mul = kappal + FOUR_THIRDS * mul
- lambdal = lambdalplus2mul - 2.*mul
+! full anisotropic case, stress calculations
+ if(ANISOTROPY) then
+ c11 = c11store(i,j,k,ispec)
+ c12 = c12store(i,j,k,ispec)
+ c13 = c13store(i,j,k,ispec)
+ c14 = c14store(i,j,k,ispec)
+ c15 = c15store(i,j,k,ispec)
+ c16 = c16store(i,j,k,ispec)
+ c22 = c22store(i,j,k,ispec)
+ c23 = c23store(i,j,k,ispec)
+ c24 = c24store(i,j,k,ispec)
+ c25 = c25store(i,j,k,ispec)
+ c26 = c26store(i,j,k,ispec)
+ c33 = c33store(i,j,k,ispec)
+ c34 = c34store(i,j,k,ispec)
+ c35 = c35store(i,j,k,ispec)
+ c36 = c36store(i,j,k,ispec)
+ c44 = c44store(i,j,k,ispec)
+ c45 = c45store(i,j,k,ispec)
+ c46 = c46store(i,j,k,ispec)
+ c55 = c55store(i,j,k,ispec)
+ c56 = c56store(i,j,k,ispec)
+ c66 = c66store(i,j,k,ispec)
+ !if(ATTENUATION .and. not_fully_in_bedrock(ispec)) then
+ ! mul = c44
+ ! c11 = c11 + FOUR_THIRDS * minus_sum_beta * mul
+ ! c12 = c12 - TWO_THIRDS * minus_sum_beta * mul
+ ! c13 = c13 - TWO_THIRDS * minus_sum_beta * mul
+ ! c22 = c22 + FOUR_THIRDS * minus_sum_beta * mul
+ ! c23 = c23 - TWO_THIRDS * minus_sum_beta * mul
+ ! c33 = c33 + FOUR_THIRDS * minus_sum_beta * mul
+ ! c44 = c44 + minus_sum_beta * mul
+ ! c55 = c55 + minus_sum_beta * mul
+ ! c66 = c66 + minus_sum_beta * mul
+ !endif
- ! compute stress sigma
- sigma_xx = lambdalplus2mul*duxdxl + lambdal*duydyl_plus_duzdzl
- sigma_yy = lambdalplus2mul*duydyl + lambdal*duxdxl_plus_duzdzl
- sigma_zz = lambdalplus2mul*duzdzl + lambdal*duxdxl_plus_duydyl
+ sigma_xx = c11*duxdxl + c16*duxdyl_plus_duydxl + c12*duydyl + &
+ c15*duzdxl_plus_duxdzl + c14*duzdyl_plus_duydzl + c13*duzdzl
+ sigma_yy = c12*duxdxl + c26*duxdyl_plus_duydxl + c22*duydyl + &
+ c25*duzdxl_plus_duxdzl + c24*duzdyl_plus_duydzl + c23*duzdzl
+ sigma_zz = c13*duxdxl + c36*duxdyl_plus_duydxl + c23*duydyl + &
+ c35*duzdxl_plus_duxdzl + c34*duzdyl_plus_duydzl + c33*duzdzl
+ sigma_xy = c16*duxdxl + c66*duxdyl_plus_duydxl + c26*duydyl + &
+ c56*duzdxl_plus_duxdzl + c46*duzdyl_plus_duydzl + c36*duzdzl
+ sigma_xz = c15*duxdxl + c56*duxdyl_plus_duydxl + c25*duydyl + &
+ c55*duzdxl_plus_duxdzl + c45*duzdyl_plus_duydzl + c35*duzdzl
+ sigma_yz = c14*duxdxl + c46*duxdyl_plus_duydxl + c24*duydyl + &
+ c45*duzdxl_plus_duxdzl + c44*duzdyl_plus_duydzl + c34*duzdzl
- sigma_xy = mul*duxdyl_plus_duydxl
- sigma_xz = mul*duzdxl_plus_duxdzl
- sigma_yz = mul*duzdyl_plus_duydzl
+ !if (SIMULATION_TYPE == 3) then
+ ! b_sigma_xx = c11*b_duxdxl + c16*b_duxdyl_plus_duydxl + c12*b_duydyl + &
+ ! c15*b_duzdxl_plus_duxdzl + c14*b_duzdyl_plus_duydzl + c13*b_duzdzl
+ ! b_sigma_yy = c12*b_duxdxl + c26*b_duxdyl_plus_duydxl + c22*b_duydyl + &
+ ! c25*b_duzdxl_plus_duxdzl + c24*b_duzdyl_plus_duydzl + c23*b_duzdzl
+ ! b_sigma_zz = c13*b_duxdxl + c36*b_duxdyl_plus_duydxl + c23*b_duydyl + &
+ ! c35*b_duzdxl_plus_duxdzl + c34*b_duzdyl_plus_duydzl + c33*b_duzdzl
+ ! b_sigma_xy = c16*b_duxdxl + c66*b_duxdyl_plus_duydxl + c26*b_duydyl + &
+ ! c56*b_duzdxl_plus_duxdzl + c46*b_duzdyl_plus_duydzl + c36*b_duzdzl
+ ! b_sigma_xz = c15*b_duxdxl + c56*b_duxdyl_plus_duydxl + c25*b_duydyl + &
+ ! c55*b_duzdxl_plus_duxdzl + c45*b_duzdyl_plus_duydzl + c35*b_duzdzl
+ ! b_sigma_yz = c14*b_duxdxl + c46*b_duxdyl_plus_duydxl + c24*b_duydyl + &
+ ! c45*b_duzdxl_plus_duxdzl + c44*b_duzdyl_plus_duydzl + c34*b_duzdzl
+ !endif
+ else
+! isotropic case
+ lambdalplus2mul = kappal + FOUR_THIRDS * mul
+ lambdal = lambdalplus2mul - 2.*mul
+
+ ! compute stress sigma
+ sigma_xx = lambdalplus2mul*duxdxl + lambdal*duydyl_plus_duzdzl
+ sigma_yy = lambdalplus2mul*duydyl + lambdal*duxdxl_plus_duzdzl
+ sigma_zz = lambdalplus2mul*duzdzl + lambdal*duxdxl_plus_duydyl
+
+ sigma_xy = mul*duxdyl_plus_duydxl
+ sigma_xz = mul*duzdxl_plus_duxdzl
+ sigma_yz = mul*duzdyl_plus_duydzl
+
+ !if (SIMULATION_TYPE == 3) then
+ ! b_sigma_xx = lambdalplus2mul*b_duxdxl + lambdal*b_duydyl_plus_duzdzl
+ ! b_sigma_yy = lambdalplus2mul*b_duydyl + lambdal*b_duxdxl_plus_duzdzl
+ ! b_sigma_zz = lambdalplus2mul*b_duzdzl + lambdal*b_duxdxl_plus_duydyl
+ !
+ ! b_sigma_xy = mul*b_duxdyl_plus_duydxl
+ ! b_sigma_xz = mul*b_duzdxl_plus_duxdzl
+ ! b_sigma_yz = mul*b_duzdyl_plus_duydzl
+ !endif
+
+ endif ! ANISOTROPY
+
! subtract memory variables if attenuation
if(ATTENUATION) then
do i_sls = 1,N_SLS
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_with_Deville.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_with_Deville.f90 2009-10-25 23:06:34 UTC (rev 15876)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/compute_forces_with_Deville.f90 2009-10-25 23:54:42 UTC (rev 15877)
@@ -37,7 +37,12 @@
NSPEC_ATTENUATION_AB,R_xx,R_yy,R_xy,R_xz,R_yz, &
epsilondev_xx,epsilondev_yy,epsilondev_xy, &
epsilondev_xz,epsilondev_yz,iflag_attenuation_store, &
- rho_vs )
+ rho_vs, &
+ ANISOTROPY,NSPEC_ANISO, &
+ c11store,c12store,c13store,c14store,c15store,c16store,&
+ c22store,c23store,c24store,c25store,c26store,c33store,&
+ c34store,c35store,c36store,c44store,c45store,c46store,&
+ c55store,c56store,c66store )
! computes elastic tensor term
@@ -82,6 +87,15 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: rho_vs
+! anisotropy
+ logical :: ANISOTROPY
+ integer :: NSPEC_ANISO
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO) :: &
+ c11store,c12store,c13store,c14store,c15store,c16store, &
+ c22store,c23store,c24store,c25store,c26store,c33store, &
+ c34store,c35store,c36store,c44store,c45store,c46store, &
+ c55store,c56store,c66store
+
! local parameters
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dummyx_loc,dummyy_loc,dummyz_loc, &
newtempx1,newtempx2,newtempx3,newtempy1,newtempy2,newtempy3,newtempz1,newtempz2,newtempz3
@@ -138,6 +152,10 @@
real(kind=CUSTOM_REAL) lambdal,mul,lambdalplus2mul
real(kind=CUSTOM_REAL) kappal
+
+! local anisotropy parameters
+ real(kind=CUSTOM_REAL) c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,&
+ c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
integer i_SLS,iselected
@@ -264,7 +282,7 @@
duzdyl = xiyl*tempz1(i,j,k) + etayl*tempz2(i,j,k) + gammayl*tempz3(i,j,k)
duzdzl = xizl*tempz1(i,j,k) + etazl*tempz2(i,j,k) + gammazl*tempz3(i,j,k)
- ! precompute some sums to save CPU time
+! precompute some sums to save CPU time
duxdxl_plus_duydyl = duxdxl + duydyl
duxdxl_plus_duzdzl = duxdxl + duzdzl
duydyl_plus_duzdzl = duydyl + duzdzl
@@ -274,7 +292,8 @@
kappal = kappastore(i,j,k,ispec)
mul = mustore(i,j,k,ispec)
-
+
+! attenuation
if(ATTENUATION) then
! compute deviatoric strain
epsilon_trace_over_3 = ONE_THIRD * (duxdxl + duydyl + duzdzl)
@@ -307,30 +326,119 @@
endif
- lambdalplus2mul = kappal + FOUR_THIRDS * mul
- lambdal = lambdalplus2mul - 2.*mul
+! full anisotropic case, stress calculations
+ if(ANISOTROPY) then
+ c11 = c11store(i,j,k,ispec)
+ c12 = c12store(i,j,k,ispec)
+ c13 = c13store(i,j,k,ispec)
+ c14 = c14store(i,j,k,ispec)
+ c15 = c15store(i,j,k,ispec)
+ c16 = c16store(i,j,k,ispec)
+ c22 = c22store(i,j,k,ispec)
+ c23 = c23store(i,j,k,ispec)
+ c24 = c24store(i,j,k,ispec)
+ c25 = c25store(i,j,k,ispec)
+ c26 = c26store(i,j,k,ispec)
+ c33 = c33store(i,j,k,ispec)
+ c34 = c34store(i,j,k,ispec)
+ c35 = c35store(i,j,k,ispec)
+ c36 = c36store(i,j,k,ispec)
+ c44 = c44store(i,j,k,ispec)
+ c45 = c45store(i,j,k,ispec)
+ c46 = c46store(i,j,k,ispec)
+ c55 = c55store(i,j,k,ispec)
+ c56 = c56store(i,j,k,ispec)
+ c66 = c66store(i,j,k,ispec)
+ !if(ATTENUATION .and. not_fully_in_bedrock(ispec)) then
+ ! mul = c44
+ ! c11 = c11 + FOUR_THIRDS * minus_sum_beta * mul
+ ! c12 = c12 - TWO_THIRDS * minus_sum_beta * mul
+ ! c13 = c13 - TWO_THIRDS * minus_sum_beta * mul
+ ! c22 = c22 + FOUR_THIRDS * minus_sum_beta * mul
+ ! c23 = c23 - TWO_THIRDS * minus_sum_beta * mul
+ ! c33 = c33 + FOUR_THIRDS * minus_sum_beta * mul
+ ! c44 = c44 + minus_sum_beta * mul
+ ! c55 = c55 + minus_sum_beta * mul
+ ! c66 = c66 + minus_sum_beta * mul
+ !endif
- ! compute stress sigma
- sigma_xx = lambdalplus2mul*duxdxl + lambdal*duydyl_plus_duzdzl
- sigma_yy = lambdalplus2mul*duydyl + lambdal*duxdxl_plus_duzdzl
- sigma_zz = lambdalplus2mul*duzdzl + lambdal*duxdxl_plus_duydyl
+ sigma_xx = c11*duxdxl + c16*duxdyl_plus_duydxl + c12*duydyl + &
+ c15*duzdxl_plus_duxdzl + c14*duzdyl_plus_duydzl + c13*duzdzl
+ sigma_yy = c12*duxdxl + c26*duxdyl_plus_duydxl + c22*duydyl + &
+ c25*duzdxl_plus_duxdzl + c24*duzdyl_plus_duydzl + c23*duzdzl
+ sigma_zz = c13*duxdxl + c36*duxdyl_plus_duydxl + c23*duydyl + &
+ c35*duzdxl_plus_duxdzl + c34*duzdyl_plus_duydzl + c33*duzdzl
+ sigma_xy = c16*duxdxl + c66*duxdyl_plus_duydxl + c26*duydyl + &
+ c56*duzdxl_plus_duxdzl + c46*duzdyl_plus_duydzl + c36*duzdzl
+ sigma_xz = c15*duxdxl + c56*duxdyl_plus_duydxl + c25*duydyl + &
+ c55*duzdxl_plus_duxdzl + c45*duzdyl_plus_duydzl + c35*duzdzl
+ sigma_yz = c14*duxdxl + c46*duxdyl_plus_duydxl + c24*duydyl + &
+ c45*duzdxl_plus_duxdzl + c44*duzdyl_plus_duydzl + c34*duzdzl
- sigma_xy = mul*duxdyl_plus_duydxl
- sigma_xz = mul*duzdxl_plus_duxdzl
- sigma_yz = mul*duzdyl_plus_duydzl
+ !if (SIMULATION_TYPE == 3) then
+ ! b_sigma_xx = c11*b_duxdxl + c16*b_duxdyl_plus_duydxl + c12*b_duydyl + &
+ ! c15*b_duzdxl_plus_duxdzl + c14*b_duzdyl_plus_duydzl + c13*b_duzdzl
+ ! b_sigma_yy = c12*b_duxdxl + c26*b_duxdyl_plus_duydxl + c22*b_duydyl + &
+ ! c25*b_duzdxl_plus_duxdzl + c24*b_duzdyl_plus_duydzl + c23*b_duzdzl
+ ! b_sigma_zz = c13*b_duxdxl + c36*b_duxdyl_plus_duydxl + c23*b_duydyl + &
+ ! c35*b_duzdxl_plus_duxdzl + c34*b_duzdyl_plus_duydzl + c33*b_duzdzl
+ ! b_sigma_xy = c16*b_duxdxl + c66*b_duxdyl_plus_duydxl + c26*b_duydyl + &
+ ! c56*b_duzdxl_plus_duxdzl + c46*b_duzdyl_plus_duydzl + c36*b_duzdzl
+ ! b_sigma_xz = c15*b_duxdxl + c56*b_duxdyl_plus_duydxl + c25*b_duydyl + &
+ ! c55*b_duzdxl_plus_duxdzl + c45*b_duzdyl_plus_duydzl + c35*b_duzdzl
+ ! b_sigma_yz = c14*b_duxdxl + c46*b_duxdyl_plus_duydxl + c24*b_duydyl + &
+ ! c45*b_duzdxl_plus_duxdzl + c44*b_duzdyl_plus_duydzl + c34*b_duzdzl
+ !endif
+ else
+! isotropic case
+ lambdalplus2mul = kappal + FOUR_THIRDS * mul
+ lambdal = lambdalplus2mul - 2.*mul
+
+ ! compute stress sigma
+ sigma_xx = lambdalplus2mul*duxdxl + lambdal*duydyl_plus_duzdzl
+ sigma_yy = lambdalplus2mul*duydyl + lambdal*duxdxl_plus_duzdzl
+ sigma_zz = lambdalplus2mul*duzdzl + lambdal*duxdxl_plus_duydyl
+
+ sigma_xy = mul*duxdyl_plus_duydxl
+ sigma_xz = mul*duzdxl_plus_duxdzl
+ sigma_yz = mul*duzdyl_plus_duydzl
+
+ !if (SIMULATION_TYPE == 3) then
+ ! b_sigma_xx = lambdalplus2mul*b_duxdxl + lambdal*b_duydyl_plus_duzdzl
+ ! b_sigma_yy = lambdalplus2mul*b_duydyl + lambdal*b_duxdxl_plus_duzdzl
+ ! b_sigma_zz = lambdalplus2mul*b_duzdzl + lambdal*b_duxdxl_plus_duydyl
+ !
+ ! b_sigma_xy = mul*b_duxdyl_plus_duydxl
+ ! b_sigma_xz = mul*b_duzdxl_plus_duxdzl
+ ! b_sigma_yz = mul*b_duzdyl_plus_duydzl
+ !endif
+
+ endif ! ANISOTROPY
+
! subtract memory variables if attenuation
if(ATTENUATION) then
- do i_sls = 1,N_SLS
- R_xx_val = R_xx(i,j,k,ispec,i_sls)
- R_yy_val = R_yy(i,j,k,ispec,i_sls)
- sigma_xx = sigma_xx - R_xx_val
- sigma_yy = sigma_yy - R_yy_val
- sigma_zz = sigma_zz + R_xx_val + R_yy_val
- sigma_xy = sigma_xy - R_xy(i,j,k,ispec,i_sls)
- sigma_xz = sigma_xz - R_xz(i,j,k,ispec,i_sls)
- sigma_yz = sigma_yz - R_yz(i,j,k,ispec,i_sls)
- enddo
+ do i_sls = 1,N_SLS
+ R_xx_val = R_xx(i,j,k,ispec,i_sls)
+ R_yy_val = R_yy(i,j,k,ispec,i_sls)
+ sigma_xx = sigma_xx - R_xx_val
+ sigma_yy = sigma_yy - R_yy_val
+ sigma_zz = sigma_zz + R_xx_val + R_yy_val
+ sigma_xy = sigma_xy - R_xy(i,j,k,ispec,i_sls)
+ sigma_xz = sigma_xz - R_xz(i,j,k,ispec,i_sls)
+ sigma_yz = sigma_yz - R_yz(i,j,k,ispec,i_sls)
+ enddo
+
+ !if (SIMULATION_TYPE == 3) then
+ ! b_R_xx_val = b_R_xx(i,j,k,ispec,i_sls)
+ ! b_R_yy_val = b_R_yy(i,j,k,ispec,i_sls)
+ ! b_sigma_xx = b_sigma_xx - b_R_xx_val
+ ! b_sigma_yy = b_sigma_yy - b_R_yy_val
+ ! b_sigma_zz = b_sigma_zz + b_R_xx_val + b_R_yy_val
+ ! b_sigma_xy = b_sigma_xy - b_R_xy(i,j,k,ispec,i_sls)
+ ! b_sigma_xz = b_sigma_xz - b_R_xz(i,j,k,ispec,i_sls)
+ ! b_sigma_yz = b_sigma_yz - b_R_yz(i,j,k,ispec,i_sls)
+ !endif
endif
! form dot product with test vector, symmetric form
@@ -556,8 +664,568 @@
end subroutine compute_forces_with_Deville
+!
+!-------------------------------------------------------------------------------------------------
+!
+!subroutine compute_forces_with_Deville_noanisotropy( phase_is_inner ,NSPEC_AB,NGLOB_AB, &
+! displ,accel, &
+! xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz, &
+! hprime_xx,hprime_xxT, &
+! hprimewgll_xx,hprimewgll_xxT, &
+! wgllwgll_xy,wgllwgll_xz,wgllwgll_yz, &
+! kappastore,mustore,jacobian,ibool, &
+! ispec_is_inner, &
+! ATTENUATION,USE_OLSEN_ATTENUATION, &
+! one_minus_sum_beta,factor_common,alphaval,betaval,gammaval,&
+! NSPEC_ATTENUATION_AB,R_xx,R_yy,R_xy,R_xz,R_yz, &
+! epsilondev_xx,epsilondev_yy,epsilondev_xy, &
+! epsilondev_xz,epsilondev_yz,iflag_attenuation_store, &
+! rho_vs )
+!
+!! computes elastic tensor term
+!
+! implicit none
+!
+! include "constants.h"
+!! include values created by the mesher
+!! include "OUTPUT_FILES/values_from_mesher.h"
+!
+! integer :: NSPEC_AB,NGLOB_AB
+!
+!! displacement and acceleration
+! real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_AB) :: displ,accel
+!
+!! arrays with mesh parameters per slice
+! integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: ibool
+! real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+! xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+! real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: &
+! kappastore,mustore,jacobian
+!
+!! array with derivatives of Lagrange polynomials and precalculated products
+! real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprime_xxT,hprimewgll_xx,hprimewgll_xxT
+! real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY) :: wgllwgll_xy
+! real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ) :: wgllwgll_xz
+! real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLZ) :: wgllwgll_yz
+!
+!! communication overlap
+! logical, dimension(NSPEC_AB) :: ispec_is_inner
+! logical :: phase_is_inner
+!
+!! memory variables and standard linear solids for attenuation
+! logical :: ATTENUATION,USE_OLSEN_ATTENUATION
+! integer :: NSPEC_ATTENUATION_AB
+! integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: iflag_attenuation_store
+! real(kind=CUSTOM_REAL), dimension(NUM_REGIONS_ATTENUATION) :: one_minus_sum_beta
+! real(kind=CUSTOM_REAL), dimension(NUM_REGIONS_ATTENUATION,N_SLS) :: factor_common, alphaval,betaval,gammaval
+! real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB,N_SLS) :: &
+! R_xx,R_yy,R_xy,R_xz,R_yz
+! real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ATTENUATION_AB) :: &
+! epsilondev_xx,epsilondev_yy,epsilondev_xy,epsilondev_xz,epsilondev_yz
+!
+! real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_AB) :: rho_vs
+!
+!! local parameters
+! real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dummyx_loc,dummyy_loc,dummyz_loc, &
+! newtempx1,newtempx2,newtempx3,newtempy1,newtempy2,newtempy3,newtempz1,newtempz2,newtempz3
+!
+! real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
+! tempx1,tempx2,tempx3,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3
+!
+!! manually inline the calls to the Deville et al. (2002) routines
+! real(kind=4), dimension(NGLLX,m2) :: B1_m1_m2_5points,B2_m1_m2_5points,B3_m1_m2_5points
+! real(kind=4), dimension(m1,m2) :: C1_m1_m2_5points,C2_m1_m2_5points,C3_m1_m2_5points
+! real(kind=4), dimension(m1,m2) :: E1_m1_m2_5points,E2_m1_m2_5points,E3_m1_m2_5points
+!
+! equivalence(dummyx_loc,B1_m1_m2_5points)
+! equivalence(dummyy_loc,B2_m1_m2_5points)
+! equivalence(dummyz_loc,B3_m1_m2_5points)
+! equivalence(tempx1,C1_m1_m2_5points)
+! equivalence(tempy1,C2_m1_m2_5points)
+! equivalence(tempz1,C3_m1_m2_5points)
+! equivalence(newtempx1,E1_m1_m2_5points)
+! equivalence(newtempy1,E2_m1_m2_5points)
+! equivalence(newtempz1,E3_m1_m2_5points)
+!
+! real(kind=4), dimension(m2,NGLLX) :: A1_mxm_m2_m1_5points,A2_mxm_m2_m1_5points,A3_mxm_m2_m1_5points
+! real(kind=4), dimension(m2,m1) :: C1_mxm_m2_m1_5points,C2_mxm_m2_m1_5points,C3_mxm_m2_m1_5points
+! real(kind=4), dimension(m2,m1) :: E1_mxm_m2_m1_5points,E2_mxm_m2_m1_5points,E3_mxm_m2_m1_5points
+!
+! equivalence(dummyx_loc,A1_mxm_m2_m1_5points)
+! equivalence(dummyy_loc,A2_mxm_m2_m1_5points)
+! equivalence(dummyz_loc,A3_mxm_m2_m1_5points)
+! equivalence(tempx3,C1_mxm_m2_m1_5points)
+! equivalence(tempy3,C2_mxm_m2_m1_5points)
+! equivalence(tempz3,C3_mxm_m2_m1_5points)
+! equivalence(newtempx3,E1_mxm_m2_m1_5points)
+! equivalence(newtempy3,E2_mxm_m2_m1_5points)
+! equivalence(newtempz3,E3_mxm_m2_m1_5points)
+!
+!! local attenuation parameters
+! real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: epsilondev_xx_loc, &
+! epsilondev_yy_loc, epsilondev_xy_loc, epsilondev_xz_loc, epsilondev_yz_loc
+! real(kind=CUSTOM_REAL) R_xx_val,R_yy_val
+! real(kind=CUSTOM_REAL) factor_loc,alphaval_loc,betaval_loc,gammaval_loc,Sn,Snp1
+! real(kind=CUSTOM_REAL) epsilon_trace_over_3
+! real(kind=CUSTOM_REAL) vs_val
+!
+! real(kind=CUSTOM_REAL) xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl,jacobianl
+! real(kind=CUSTOM_REAL) duxdxl,duxdyl,duxdzl,duydxl,duydyl,duydzl,duzdxl,duzdyl,duzdzl
+!
+! real(kind=CUSTOM_REAL) duxdxl_plus_duydyl,duxdxl_plus_duzdzl,duydyl_plus_duzdzl
+! real(kind=CUSTOM_REAL) duxdyl_plus_duydxl,duzdxl_plus_duxdzl,duzdyl_plus_duydzl
+!
+! real(kind=CUSTOM_REAL) sigma_xx,sigma_yy,sigma_zz,sigma_xy,sigma_xz,sigma_yz
+!
+! real(kind=CUSTOM_REAL) fac1,fac2,fac3
+!
+! real(kind=CUSTOM_REAL) lambdal,mul,lambdalplus2mul
+! real(kind=CUSTOM_REAL) kappal
+!
+! integer i_SLS,iselected
+!
+! integer ispec,iglob
+! integer i,j,k
+!
+!! loops over all elements
+! do ispec = 1,NSPEC_AB
+!
+! if (ispec_is_inner(ispec) .eqv. phase_is_inner) then
+!
+! do k=1,NGLLZ
+! do j=1,NGLLY
+! do i=1,NGLLX
+! iglob = ibool(i,j,k,ispec)
+! dummyx_loc(i,j,k) = displ(1,iglob)
+! dummyy_loc(i,j,k) = displ(2,iglob)
+! dummyz_loc(i,j,k) = displ(3,iglob)
+! enddo
+! enddo
+! enddo
+!
+! ! subroutines adapted from Deville, Fischer and Mund, High-order methods
+! ! for incompressible fluid flow, Cambridge University Press (2002),
+! ! pages 386 and 389 and Figure 8.3.1
+! ! call mxm_m1_m2_5points(hprime_xx,dummyx_loc,dummyy_loc,dummyz_loc,tempx1,tempy1,tempz1)
+! do j=1,m2
+! do i=1,m1
+! C1_m1_m2_5points(i,j) = hprime_xx(i,1)*B1_m1_m2_5points(1,j) + &
+! hprime_xx(i,2)*B1_m1_m2_5points(2,j) + &
+! hprime_xx(i,3)*B1_m1_m2_5points(3,j) + &
+! hprime_xx(i,4)*B1_m1_m2_5points(4,j) + &
+! hprime_xx(i,5)*B1_m1_m2_5points(5,j)
+!
+! C2_m1_m2_5points(i,j) = hprime_xx(i,1)*B2_m1_m2_5points(1,j) + &
+! hprime_xx(i,2)*B2_m1_m2_5points(2,j) + &
+! hprime_xx(i,3)*B2_m1_m2_5points(3,j) + &
+! hprime_xx(i,4)*B2_m1_m2_5points(4,j) + &
+! hprime_xx(i,5)*B2_m1_m2_5points(5,j)
+!
+! C3_m1_m2_5points(i,j) = hprime_xx(i,1)*B3_m1_m2_5points(1,j) + &
+! hprime_xx(i,2)*B3_m1_m2_5points(2,j) + &
+! hprime_xx(i,3)*B3_m1_m2_5points(3,j) + &
+! hprime_xx(i,4)*B3_m1_m2_5points(4,j) + &
+! hprime_xx(i,5)*B3_m1_m2_5points(5,j)
+! enddo
+! enddo
+!
+! ! call mxm_m1_m1_5points(dummyx_loc(1,1,k),dummyy_loc(1,1,k),dummyz_loc(1,1,k), &
+! ! hprime_xxT,tempx2(1,1,k),tempy2(1,1,k),tempz2(1,1,k))
+! do j=1,m1
+! do i=1,m1
+! ! for efficiency it is better to leave this loop on k inside, it leads to slightly faster code
+! do k = 1,NGLLX
+! tempx2(i,j,k) = dummyx_loc(i,1,k)*hprime_xxT(1,j) + &
+! dummyx_loc(i,2,k)*hprime_xxT(2,j) + &
+! dummyx_loc(i,3,k)*hprime_xxT(3,j) + &
+! dummyx_loc(i,4,k)*hprime_xxT(4,j) + &
+! dummyx_loc(i,5,k)*hprime_xxT(5,j)
+!
+! tempy2(i,j,k) = dummyy_loc(i,1,k)*hprime_xxT(1,j) + &
+! dummyy_loc(i,2,k)*hprime_xxT(2,j) + &
+! dummyy_loc(i,3,k)*hprime_xxT(3,j) + &
+! dummyy_loc(i,4,k)*hprime_xxT(4,j) + &
+! dummyy_loc(i,5,k)*hprime_xxT(5,j)
+!
+! tempz2(i,j,k) = dummyz_loc(i,1,k)*hprime_xxT(1,j) + &
+! dummyz_loc(i,2,k)*hprime_xxT(2,j) + &
+! dummyz_loc(i,3,k)*hprime_xxT(3,j) + &
+! dummyz_loc(i,4,k)*hprime_xxT(4,j) + &
+! dummyz_loc(i,5,k)*hprime_xxT(5,j)
+! enddo
+! enddo
+! enddo
+!
+! ! call mxm_m2_m1_5points(dummyx_loc,dummyy_loc,dummyz_loc,tempx3,tempy3,tempz3)
+! do j=1,m1
+! do i=1,m2
+! C1_mxm_m2_m1_5points(i,j) = A1_mxm_m2_m1_5points(i,1)*hprime_xxT(1,j) + &
+! A1_mxm_m2_m1_5points(i,2)*hprime_xxT(2,j) + &
+! A1_mxm_m2_m1_5points(i,3)*hprime_xxT(3,j) + &
+! A1_mxm_m2_m1_5points(i,4)*hprime_xxT(4,j) + &
+! A1_mxm_m2_m1_5points(i,5)*hprime_xxT(5,j)
+!
+! C2_mxm_m2_m1_5points(i,j) = A2_mxm_m2_m1_5points(i,1)*hprime_xxT(1,j) + &
+! A2_mxm_m2_m1_5points(i,2)*hprime_xxT(2,j) + &
+! A2_mxm_m2_m1_5points(i,3)*hprime_xxT(3,j) + &
+! A2_mxm_m2_m1_5points(i,4)*hprime_xxT(4,j) + &
+! A2_mxm_m2_m1_5points(i,5)*hprime_xxT(5,j)
+!
+! C3_mxm_m2_m1_5points(i,j) = A3_mxm_m2_m1_5points(i,1)*hprime_xxT(1,j) + &
+! A3_mxm_m2_m1_5points(i,2)*hprime_xxT(2,j) + &
+! A3_mxm_m2_m1_5points(i,3)*hprime_xxT(3,j) + &
+! A3_mxm_m2_m1_5points(i,4)*hprime_xxT(4,j) + &
+! A3_mxm_m2_m1_5points(i,5)*hprime_xxT(5,j)
+! enddo
+! enddo
+!
+! do k=1,NGLLZ
+! do j=1,NGLLY
+! do i=1,NGLLX
+!
+! ! get derivatives of ux, uy and uz with respect to x, y and z
+! xixl = xix(i,j,k,ispec)
+! xiyl = xiy(i,j,k,ispec)
+! xizl = xiz(i,j,k,ispec)
+! etaxl = etax(i,j,k,ispec)
+! etayl = etay(i,j,k,ispec)
+! etazl = etaz(i,j,k,ispec)
+! gammaxl = gammax(i,j,k,ispec)
+! gammayl = gammay(i,j,k,ispec)
+! gammazl = gammaz(i,j,k,ispec)
+! jacobianl = jacobian(i,j,k,ispec)
+!
+! duxdxl = xixl*tempx1(i,j,k) + etaxl*tempx2(i,j,k) + gammaxl*tempx3(i,j,k)
+! duxdyl = xiyl*tempx1(i,j,k) + etayl*tempx2(i,j,k) + gammayl*tempx3(i,j,k)
+! duxdzl = xizl*tempx1(i,j,k) + etazl*tempx2(i,j,k) + gammazl*tempx3(i,j,k)
+!
+! duydxl = xixl*tempy1(i,j,k) + etaxl*tempy2(i,j,k) + gammaxl*tempy3(i,j,k)
+! duydyl = xiyl*tempy1(i,j,k) + etayl*tempy2(i,j,k) + gammayl*tempy3(i,j,k)
+! duydzl = xizl*tempy1(i,j,k) + etazl*tempy2(i,j,k) + gammazl*tempy3(i,j,k)
+!
+! duzdxl = xixl*tempz1(i,j,k) + etaxl*tempz2(i,j,k) + gammaxl*tempz3(i,j,k)
+! duzdyl = xiyl*tempz1(i,j,k) + etayl*tempz2(i,j,k) + gammayl*tempz3(i,j,k)
+! duzdzl = xizl*tempz1(i,j,k) + etazl*tempz2(i,j,k) + gammazl*tempz3(i,j,k)
+!
+!! precompute some sums to save CPU time
+! duxdxl_plus_duydyl = duxdxl + duydyl
+! duxdxl_plus_duzdzl = duxdxl + duzdzl
+! duydyl_plus_duzdzl = duydyl + duzdzl
+! duxdyl_plus_duydxl = duxdyl + duydxl
+! duzdxl_plus_duxdzl = duzdxl + duxdzl
+! duzdyl_plus_duydzl = duzdyl + duydzl
+!
+! kappal = kappastore(i,j,k,ispec)
+! mul = mustore(i,j,k,ispec)
+!
+!! attenuation
+! if(ATTENUATION) then
+! ! compute deviatoric strain
+! epsilon_trace_over_3 = ONE_THIRD * (duxdxl + duydyl + duzdzl)
+! epsilondev_xx_loc(i,j,k) = duxdxl - epsilon_trace_over_3
+! epsilondev_yy_loc(i,j,k) = duydyl - epsilon_trace_over_3
+! epsilondev_xy_loc(i,j,k) = 0.5 * duxdyl_plus_duydxl
+! epsilondev_xz_loc(i,j,k) = 0.5 * duzdxl_plus_duxdzl
+! epsilondev_yz_loc(i,j,k) = 0.5 * duzdyl_plus_duydzl
+!
+! !if (SIMULATION_TYPE == 3) then
+! ! b_epsilon_trace_over_3 = ONE_THIRD * (b_duxdxl + b_duydyl + b_duzdzl)
+! ! b_epsilondev_xx_loc(i,j,k) = b_duxdxl - b_epsilon_trace_over_3
+! ! b_epsilondev_yy_loc(i,j,k) = b_duydyl - b_epsilon_trace_over_3
+! ! b_epsilondev_xy_loc(i,j,k) = 0.5 * b_duxdyl_plus_duydxl
+! ! b_epsilondev_xz_loc(i,j,k) = 0.5 * b_duzdxl_plus_duxdzl
+! ! b_epsilondev_yz_loc(i,j,k) = 0.5 * b_duzdyl_plus_duydzl
+! !endif
+!
+! ! uses scaling rule similar to Olsen et al. (2003) or mesh flag
+! if(USE_OLSEN_ATTENUATION) then
+! vs_val = mustore(i,j,k,ispec) / rho_vs(i,j,k,ispec)
+! call get_attenuation_model_Olsen_sediment( vs_val, iselected )
+! else
+! ! iflag from (CUBIT) mesh
+! iselected = iflag_attenuation_store(i,j,k,ispec)
+! endif
+!
+! ! use unrelaxed parameters if attenuation
+! mul = mul * one_minus_sum_beta(iselected)
+!
+! endif
+!
+!! isotropic case
+! lambdalplus2mul = kappal + FOUR_THIRDS * mul
+! lambdal = lambdalplus2mul - 2.*mul
+!
+! ! compute stress sigma
+! sigma_xx = lambdalplus2mul*duxdxl + lambdal*duydyl_plus_duzdzl
+! sigma_yy = lambdalplus2mul*duydyl + lambdal*duxdxl_plus_duzdzl
+! sigma_zz = lambdalplus2mul*duzdzl + lambdal*duxdxl_plus_duydyl
+!
+! sigma_xy = mul*duxdyl_plus_duydxl
+! sigma_xz = mul*duzdxl_plus_duxdzl
+! sigma_yz = mul*duzdyl_plus_duydzl
+!
+! !if (SIMULATION_TYPE == 3) then
+! ! b_sigma_xx = lambdalplus2mul*b_duxdxl + lambdal*b_duydyl_plus_duzdzl
+! ! b_sigma_yy = lambdalplus2mul*b_duydyl + lambdal*b_duxdxl_plus_duzdzl
+! ! b_sigma_zz = lambdalplus2mul*b_duzdzl + lambdal*b_duxdxl_plus_duydyl
+! !
+! ! b_sigma_xy = mul*b_duxdyl_plus_duydxl
+! ! b_sigma_xz = mul*b_duzdxl_plus_duxdzl
+! ! b_sigma_yz = mul*b_duzdyl_plus_duydzl
+! !endif
+!
+! ! subtract memory variables if attenuation
+! if(ATTENUATION) then
+! do i_sls = 1,N_SLS
+! R_xx_val = R_xx(i,j,k,ispec,i_sls)
+! R_yy_val = R_yy(i,j,k,ispec,i_sls)
+! sigma_xx = sigma_xx - R_xx_val
+! sigma_yy = sigma_yy - R_yy_val
+! sigma_zz = sigma_zz + R_xx_val + R_yy_val
+! sigma_xy = sigma_xy - R_xy(i,j,k,ispec,i_sls)
+! sigma_xz = sigma_xz - R_xz(i,j,k,ispec,i_sls)
+! sigma_yz = sigma_yz - R_yz(i,j,k,ispec,i_sls)
+! enddo
+!
+! !if (SIMULATION_TYPE == 3) then
+! ! b_R_xx_val = b_R_xx(i,j,k,ispec,i_sls)
+! ! b_R_yy_val = b_R_yy(i,j,k,ispec,i_sls)
+! ! b_sigma_xx = b_sigma_xx - b_R_xx_val
+! ! b_sigma_yy = b_sigma_yy - b_R_yy_val
+! ! b_sigma_zz = b_sigma_zz + b_R_xx_val + b_R_yy_val
+! ! b_sigma_xy = b_sigma_xy - b_R_xy(i,j,k,ispec,i_sls)
+! ! b_sigma_xz = b_sigma_xz - b_R_xz(i,j,k,ispec,i_sls)
+! ! b_sigma_yz = b_sigma_yz - b_R_yz(i,j,k,ispec,i_sls)
+! !endif
+! endif
+!
+! ! form dot product with test vector, symmetric form
+! tempx1(i,j,k) = jacobianl * (sigma_xx*xixl + sigma_xy*xiyl + sigma_xz*xizl)
+! tempy1(i,j,k) = jacobianl * (sigma_xy*xixl + sigma_yy*xiyl + sigma_yz*xizl)
+! tempz1(i,j,k) = jacobianl * (sigma_xz*xixl + sigma_yz*xiyl + sigma_zz*xizl)
+!
+! tempx2(i,j,k) = jacobianl * (sigma_xx*etaxl + sigma_xy*etayl + sigma_xz*etazl)
+! tempy2(i,j,k) = jacobianl * (sigma_xy*etaxl + sigma_yy*etayl + sigma_yz*etazl)
+! tempz2(i,j,k) = jacobianl * (sigma_xz*etaxl + sigma_yz*etayl + sigma_zz*etazl)
+!
+! tempx3(i,j,k) = jacobianl * (sigma_xx*gammaxl + sigma_xy*gammayl + sigma_xz*gammazl)
+! tempy3(i,j,k) = jacobianl * (sigma_xy*gammaxl + sigma_yy*gammayl + sigma_yz*gammazl)
+! tempz3(i,j,k) = jacobianl * (sigma_xz*gammaxl + sigma_yz*gammayl + sigma_zz*gammazl)
+!
+! enddo
+! enddo
+! enddo
+!
+! ! subroutines adapted from Deville, Fischer and Mund, High-order methods
+! ! for incompressible fluid flow, Cambridge University Press (2002),
+! ! pages 386 and 389 and Figure 8.3.1
+! ! call mxm_m1_m2_5points(hprimewgll_xxT,tempx1,tempy1,tempz1,newtempx1,newtempy1,newtempz1)
+! do j=1,m2
+! do i=1,m1
+! E1_m1_m2_5points(i,j) = hprimewgll_xxT(i,1)*C1_m1_m2_5points(1,j) + &
+! hprimewgll_xxT(i,2)*C1_m1_m2_5points(2,j) + &
+! hprimewgll_xxT(i,3)*C1_m1_m2_5points(3,j) + &
+! hprimewgll_xxT(i,4)*C1_m1_m2_5points(4,j) + &
+! hprimewgll_xxT(i,5)*C1_m1_m2_5points(5,j)
+!
+! E2_m1_m2_5points(i,j) = hprimewgll_xxT(i,1)*C2_m1_m2_5points(1,j) + &
+! hprimewgll_xxT(i,2)*C2_m1_m2_5points(2,j) + &
+! hprimewgll_xxT(i,3)*C2_m1_m2_5points(3,j) + &
+! hprimewgll_xxT(i,4)*C2_m1_m2_5points(4,j) + &
+! hprimewgll_xxT(i,5)*C2_m1_m2_5points(5,j)
+!
+! E3_m1_m2_5points(i,j) = hprimewgll_xxT(i,1)*C3_m1_m2_5points(1,j) + &
+! hprimewgll_xxT(i,2)*C3_m1_m2_5points(2,j) + &
+! hprimewgll_xxT(i,3)*C3_m1_m2_5points(3,j) + &
+! hprimewgll_xxT(i,4)*C3_m1_m2_5points(4,j) + &
+! hprimewgll_xxT(i,5)*C3_m1_m2_5points(5,j)
+! enddo
+! enddo
+!
+! ! call mxm_m1_m1_5points(tempx2(1,1,k),tempy2(1,1,k),tempz2(1,1,k), &
+! ! hprimewgll_xx,newtempx2(1,1,k),newtempy2(1,1,k),newtempz2(1,1,k))
+! do i=1,m1
+! do j=1,m1
+! ! for efficiency it is better to leave this loop on k inside, it leads to slightly faster code
+! do k = 1,NGLLX
+! newtempx2(i,j,k) = tempx2(i,1,k)*hprimewgll_xx(1,j) + &
+! tempx2(i,2,k)*hprimewgll_xx(2,j) + &
+! tempx2(i,3,k)*hprimewgll_xx(3,j) + &
+! tempx2(i,4,k)*hprimewgll_xx(4,j) + &
+! tempx2(i,5,k)*hprimewgll_xx(5,j)
+!
+! newtempy2(i,j,k) = tempy2(i,1,k)*hprimewgll_xx(1,j) + &
+! tempy2(i,2,k)*hprimewgll_xx(2,j) + &
+! tempy2(i,3,k)*hprimewgll_xx(3,j) + &
+! tempy2(i,4,k)*hprimewgll_xx(4,j) + &
+! tempy2(i,5,k)*hprimewgll_xx(5,j)
+!
+! newtempz2(i,j,k) = tempz2(i,1,k)*hprimewgll_xx(1,j) + &
+! tempz2(i,2,k)*hprimewgll_xx(2,j) + &
+! tempz2(i,3,k)*hprimewgll_xx(3,j) + &
+! tempz2(i,4,k)*hprimewgll_xx(4,j) + &
+! tempz2(i,5,k)*hprimewgll_xx(5,j)
+! enddo
+! enddo
+! enddo
+!
+! ! call mxm_m2_m1_5points(tempx3,tempy3,tempz3,hprimewgll_xx,newtempx3,newtempy3,newtempz3)
+! do j=1,m1
+! do i=1,m2
+! E1_mxm_m2_m1_5points(i,j) = C1_mxm_m2_m1_5points(i,1)*hprimewgll_xx(1,j) + &
+! C1_mxm_m2_m1_5points(i,2)*hprimewgll_xx(2,j) + &
+! C1_mxm_m2_m1_5points(i,3)*hprimewgll_xx(3,j) + &
+! C1_mxm_m2_m1_5points(i,4)*hprimewgll_xx(4,j) + &
+! C1_mxm_m2_m1_5points(i,5)*hprimewgll_xx(5,j)
+!
+! E2_mxm_m2_m1_5points(i,j) = C2_mxm_m2_m1_5points(i,1)*hprimewgll_xx(1,j) + &
+! C2_mxm_m2_m1_5points(i,2)*hprimewgll_xx(2,j) + &
+! C2_mxm_m2_m1_5points(i,3)*hprimewgll_xx(3,j) + &
+! C2_mxm_m2_m1_5points(i,4)*hprimewgll_xx(4,j) + &
+! C2_mxm_m2_m1_5points(i,5)*hprimewgll_xx(5,j)
+!
+! E3_mxm_m2_m1_5points(i,j) = C3_mxm_m2_m1_5points(i,1)*hprimewgll_xx(1,j) + &
+! C3_mxm_m2_m1_5points(i,2)*hprimewgll_xx(2,j) + &
+! C3_mxm_m2_m1_5points(i,3)*hprimewgll_xx(3,j) + &
+! C3_mxm_m2_m1_5points(i,4)*hprimewgll_xx(4,j) + &
+! C3_mxm_m2_m1_5points(i,5)*hprimewgll_xx(5,j)
+! enddo
+! enddo
+!
+! do k=1,NGLLZ
+! do j=1,NGLLY
+! do i=1,NGLLX
+!
+! fac1 = wgllwgll_yz(j,k)
+! fac2 = wgllwgll_xz(i,k)
+! fac3 = wgllwgll_xy(i,j)
+!
+! ! sum contributions from each element to the global mesh using indirect addressing
+! iglob = ibool(i,j,k,ispec)
+! accel(1,iglob) = accel(1,iglob) - fac1*newtempx1(i,j,k) - &
+! fac2*newtempx2(i,j,k) - fac3*newtempx3(i,j,k)
+! accel(2,iglob) = accel(2,iglob) - fac1*newtempy1(i,j,k) - &
+! fac2*newtempy2(i,j,k) - fac3*newtempy3(i,j,k)
+! accel(3,iglob) = accel(3,iglob) - fac1*newtempz1(i,j,k) - &
+! fac2*newtempz2(i,j,k) - fac3*newtempz3(i,j,k)
+!
+! ! update memory variables based upon the Runge-Kutta scheme
+! if(ATTENUATION) then
+!
+! ! use Runge-Kutta scheme to march in time
+! do i_sls = 1,N_SLS
+!
+! ! get coefficients for that standard linear solid
+! if( USE_OLSEN_ATTENUATION ) then
+! vs_val = mustore(i,j,k,ispec) / rho_vs(i,j,k,ispec)
+! call get_attenuation_model_Olsen_sediment( vs_val, iselected )
+! else
+! iselected = iflag_attenuation_store(i,j,k,ispec)
+! endif
+!
+! factor_loc = mustore(i,j,k,ispec) * factor_common(iselected,i_sls)
+! alphaval_loc = alphaval(iselected,i_sls)
+! betaval_loc = betaval(iselected,i_sls)
+! gammaval_loc = gammaval(iselected,i_sls)
+!
+! ! term in xx
+! Sn = factor_loc * epsilondev_xx(i,j,k,ispec)
+! Snp1 = factor_loc * epsilondev_xx_loc(i,j,k)
+! R_xx(i,j,k,ispec,i_sls) = alphaval_loc * R_xx(i,j,k,ispec,i_sls) + &
+! betaval_loc * Sn + gammaval_loc * Snp1
+!
+! ! term in yy
+! Sn = factor_loc * epsilondev_yy(i,j,k,ispec)
+! Snp1 = factor_loc * epsilondev_yy_loc(i,j,k)
+! R_yy(i,j,k,ispec,i_sls) = alphaval_loc * R_yy(i,j,k,ispec,i_sls) + &
+! betaval_loc * Sn + gammaval_loc * Snp1
+!
+! ! term in zz not computed since zero trace
+!
+! ! term in xy
+! Sn = factor_loc * epsilondev_xy(i,j,k,ispec)
+! Snp1 = factor_loc * epsilondev_xy_loc(i,j,k)
+! R_xy(i,j,k,ispec,i_sls) = alphaval_loc * R_xy(i,j,k,ispec,i_sls) + &
+! betaval_loc * Sn + gammaval_loc * Snp1
+!
+! ! term in xz
+! Sn = factor_loc * epsilondev_xz(i,j,k,ispec)
+! Snp1 = factor_loc * epsilondev_xz_loc(i,j,k)
+! R_xz(i,j,k,ispec,i_sls) = alphaval_loc * R_xz(i,j,k,ispec,i_sls) + &
+! betaval_loc * Sn + gammaval_loc * Snp1
+!
+! ! term in yz
+! Sn = factor_loc * epsilondev_yz(i,j,k,ispec)
+! Snp1 = factor_loc * epsilondev_yz_loc(i,j,k)
+! R_yz(i,j,k,ispec,i_sls) = alphaval_loc * R_yz(i,j,k,ispec,i_sls) + &
+! betaval_loc * Sn + gammaval_loc * Snp1
+!
+! !if (SIMULATION_TYPE == 3) then
+! ! b_alphaval_loc = b_alphaval(iselected,i_sls)
+! ! b_betaval_loc = b_betaval(iselected,i_sls)
+! ! b_gammaval_loc = b_gammaval(iselected,i_sls)
+! ! ! term in xx
+! ! b_Sn = factor_loc * b_epsilondev_xx(i,j,k,ispec)
+! ! b_Snp1 = factor_loc * b_epsilondev_xx_loc(i,j,k)
+! ! b_R_xx(i,j,k,ispec,i_sls) = b_alphaval_loc * b_R_xx(i,j,k,ispec,i_sls) + &
+! ! b_betaval_loc * b_Sn + b_gammaval_loc * b_Snp1
+! ! ! term in yy
+! ! b_Sn = factor_loc * b_epsilondev_yy(i,j,k,ispec)
+! ! b_Snp1 = factor_loc * b_epsilondev_yy_loc(i,j,k)
+! ! b_R_yy(i,j,k,ispec,i_sls) = b_alphaval_loc * b_R_yy(i,j,k,ispec,i_sls) + &
+! ! b_betaval_loc * b_Sn + b_gammaval_loc * b_Snp1
+! ! ! term in zz not computed since zero trace
+! ! ! term in xy
+! ! b_Sn = factor_loc * b_epsilondev_xy(i,j,k,ispec)
+! ! b_Snp1 = factor_loc * b_epsilondev_xy_loc(i,j,k)
+! ! b_R_xy(i,j,k,ispec,i_sls) = b_alphaval_loc * b_R_xy(i,j,k,ispec,i_sls) + &
+! ! b_betaval_loc * b_Sn + b_gammaval_loc * b_Snp1
+! ! ! term in xz
+! ! b_Sn = factor_loc * b_epsilondev_xz(i,j,k,ispec)
+! ! b_Snp1 = factor_loc * b_epsilondev_xz_loc(i,j,k)
+! ! b_R_xz(i,j,k,ispec,i_sls) = b_alphaval_loc * b_R_xz(i,j,k,ispec,i_sls) + &
+! ! b_betaval_loc * b_Sn + b_gammaval_loc * b_Snp1
+! ! ! term in yz
+! ! b_Sn = factor_loc * b_epsilondev_yz(i,j,k,ispec)
+! ! b_Snp1 = factor_loc * b_epsilondev_yz_loc(i,j,k)
+! ! b_R_yz(i,j,k,ispec,i_sls) = b_alphaval_loc * b_R_yz(i,j,k,ispec,i_sls) + &
+! ! b_betaval_loc * b_Sn + b_gammaval_loc * b_Snp1
+! !endif
+!
+! enddo ! end of loop on memory variables
+!
+! endif ! end attenuation
+!
+! enddo
+! enddo
+! enddo
+!
+! ! save deviatoric strain for Runge-Kutta scheme
+! if(ATTENUATION) then
+! epsilondev_xx(:,:,:,ispec) = epsilondev_xx_loc(:,:,:)
+! epsilondev_yy(:,:,:,ispec) = epsilondev_yy_loc(:,:,:)
+! epsilondev_xy(:,:,:,ispec) = epsilondev_xy_loc(:,:,:)
+! epsilondev_xz(:,:,:,ispec) = epsilondev_xz_loc(:,:,:)
+! epsilondev_yz(:,:,:,ispec) = epsilondev_yz_loc(:,:,:)
+! !if (SIMULATION_TYPE == 3) then
+! ! b_epsilondev_xx(:,:,:,ispec) = b_epsilondev_xx_loc(:,:,:)
+! ! b_epsilondev_yy(:,:,:,ispec) = b_epsilondev_yy_loc(:,:,:)
+! ! b_epsilondev_xy(:,:,:,ispec) = b_epsilondev_xy_loc(:,:,:)
+! ! b_epsilondev_xz(:,:,:,ispec) = b_epsilondev_xz_loc(:,:,:)
+! ! b_epsilondev_yz(:,:,:,ispec) = b_epsilondev_yz_loc(:,:,:)
+! !endif
+! endif
+!
+! endif ! if (ispec_is_inner(ispec) .eqv. phase_is_inner)
+!
+! enddo ! spectral element loop
+!
+!end subroutine compute_forces_with_Deville_noanisotropy
+
+
+
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/constants.h.in 2009-10-25 23:06:34 UTC (rev 15876)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/constants.h.in 2009-10-25 23:54:42 UTC (rev 15877)
@@ -209,6 +209,10 @@
! number of standard linear solids in parallel for attenuation
integer, parameter :: N_SLS = 3
+! define flag for regions for anisotropy
+ integer, parameter :: IANISOTROPY_MODEL1 = 1
+ integer, parameter :: IANISOTROPY_MODEL2 = 2
+
! flag for projection from latitude/longitude to UTM, and back
integer, parameter :: ILONGLAT2UTM = 0, IUTM2LONGLAT = 1
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/create_movie_shakemap_AVS_DX_GMT.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/create_movie_shakemap_AVS_DX_GMT.f90 2009-10-25 23:06:34 UTC (rev 15876)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/create_movie_shakemap_AVS_DX_GMT.f90 2009-10-25 23:54:42 UTC (rev 15877)
@@ -35,17 +35,21 @@
include "constants.h"
include "OUTPUT_FILES/surface_from_mesher.h"
-! number of points in each AVS or OpenDX quadrangular cell for movies
- integer, parameter :: NGNOD2D_AVS_DX = 4
-
+!-------------------------------------------------------------------------------------------------
+! user parameters
! threshold in percent of the maximum below which we cut the amplitude
- logical, parameter :: APPLY_THRESHOLD = .true.
+ logical, parameter :: APPLY_THRESHOLD = .false.
real(kind=CUSTOM_REAL), parameter :: THRESHOLD = 1._CUSTOM_REAL / 100._CUSTOM_REAL
! coefficient of power law used for non linear scaling
- logical, parameter :: NONLINEAR_SCALING = .true.
+ logical, parameter :: NONLINEAR_SCALING = .false.
real(kind=CUSTOM_REAL), parameter :: POWER_SCALING = 0.13_CUSTOM_REAL
+!-------------------------------------------------------------------------------------------------
+
+! number of points in each AVS or OpenDX quadrangular cell for movies
+ integer, parameter :: NGNOD2D_AVS_DX = 4
+
integer it,it1,it2,ivalue,nspectot_AVS_max,ispec
integer iformat,nframes,iframe,inumber,inorm,iscaling_shake
integer ibool_number,ibool_number1,ibool_number2,ibool_number3,ibool_number4
@@ -296,7 +300,11 @@
else
print *
print *,'movie data:'
- print *,' norm of velocity vector will be displayed'
+ print *,'1= norm of velocity 2=velocity x-comp 3=velocity y-comp 4=velocity z-comp'
+ print *
+ read(5,*) inorm
+ if(inorm < 1 .or. inorm > 4) stop 'incorrect value of inorm'
+ !print *,' norm of velocity vector will be displayed'
endif
! define the total number of elements at the surface
@@ -341,450 +349,487 @@
! loop on all the time steps in the range entered
do it = it1,it2
-! check if time step corresponds to a movie frame
- if(mod(it,NTSTEP_BETWEEN_FRAMES) == 0 .or. plot_shaking_map) then
+ ! check if time step corresponds to a movie frame
+ if(mod(it,NTSTEP_BETWEEN_FRAMES) == 0 .or. plot_shaking_map) then
- iframe = iframe + 1
+ iframe = iframe + 1
- print *
- if(plot_shaking_map) then
- print *,'reading shaking map snapshot'
- else
- print *,'reading snapshot time step ',it,' out of ',NSTEP
- endif
- print *
+ print *
+ if(plot_shaking_map) then
+ print *,'reading shaking map snapshot'
+ else
+ print *,'reading snapshot time step ',it,' out of ',NSTEP
+ endif
+ print *
-! read all the elements from the same file
- if(plot_shaking_map) then
- write(outputname,"('/shakingdata')")
- else
- write(outputname,"('/moviedata',i6.6)") it
- endif
- open(unit=IOUT,file=trim(OUTPUT_FILES)//outputname,status='old',action='read',form='unformatted')
- read(IOUT) store_val_x
- read(IOUT) store_val_y
- read(IOUT) store_val_z
- read(IOUT) store_val_ux
- read(IOUT) store_val_uy
- read(IOUT) store_val_uz
- close(IOUT)
+ ! read all the elements from the same file
+ if(plot_shaking_map) then
+ write(outputname,"('/shakingdata')")
+ else
+ write(outputname,"('/moviedata',i6.6)") it
+ endif
+ open(unit=IOUT,file=trim(OUTPUT_FILES)//outputname,status='old',action='read',form='unformatted')
+ read(IOUT) store_val_x
+ read(IOUT) store_val_y
+ read(IOUT) store_val_z
+ read(IOUT) store_val_ux
+ read(IOUT) store_val_uy
+ read(IOUT) store_val_uz
+ close(IOUT)
-! clear number of elements kept
- ispec = 0
+ ! clear number of elements kept
+ ispec = 0
-! read points for all the slices
- do iproc = 0,NPROC-1
+ ! read points for all the slices
+ do iproc = 0,NPROC-1
-! reset point number
- ipoin = 0
+ ! reset point number
+ ipoin = 0
- !do ispecloc = 1,NEX_PER_PROC_XI*NEX_PER_PROC_ETA
- do ispecloc = 1, NSPEC_SURFACE_EXT_MESH
+ !do ispecloc = 1,NEX_PER_PROC_XI*NEX_PER_PROC_ETA
+ do ispecloc = 1, NSPEC_SURFACE_EXT_MESH
- if(USE_HIGHRES_FOR_MOVIES) then
-! assign the OpenDX "elements"
+ if(USE_HIGHRES_FOR_MOVIES) then
+ ! assign the OpenDX "elements"
- do j = 1,NGLLY
- do i = 1,NGLLX
+ do j = 1,NGLLY
+ do i = 1,NGLLX
- ipoin = ipoin + 1
+ ipoin = ipoin + 1
- xcoord = store_val_x(ipoin,iproc)
- ycoord = store_val_y(ipoin,iproc)
- zcoord = store_val_z(ipoin,iproc)
+ ! x,y,z coordinates
+ xcoord = store_val_x(ipoin,iproc)
+ ycoord = store_val_y(ipoin,iproc)
+ zcoord = store_val_z(ipoin,iproc)
- vectorx = store_val_ux(ipoin,iproc)
- vectory = store_val_uy(ipoin,iproc)
- vectorz = store_val_uz(ipoin,iproc)
+ ! note:
+ ! for shakemaps: ux = norm displacement, uy = norm velocity, uz = norm acceleration
+ ! for movies: ux = velocity x-component, uy = velocity y-component, uz = velocity z-component
+ vectorx = store_val_ux(ipoin,iproc)
+ vectory = store_val_uy(ipoin,iproc)
+ vectorz = store_val_uz(ipoin,iproc)
- x(i,j) = xcoord
- y(i,j) = ycoord
- z(i,j) = zcoord
+ x(i,j) = xcoord
+ y(i,j) = ycoord
+ z(i,j) = zcoord
- if(plot_shaking_map) then
-!!!! NL NL mute value near source
- if ( (sqrt(((x(i,j) - (X_SOURCE_EXT_MESH))**2 + &
- (y(i,j) - (Y_SOURCE_EXT_MESH))**2 + &
- (z(i,j) - (Z_SOURCE_EXT_MESH))**2)) < RADIUS_TO_MUTE) &
- .and. MUTE_SOURCE) then
+ ! shakemap
+ if(plot_shaking_map) then
+ !!!! NL NL mute value near source
+ if ( (sqrt(((x(i,j) - (X_SOURCE_EXT_MESH))**2 + &
+ (y(i,j) - (Y_SOURCE_EXT_MESH))**2 + &
+ (z(i,j) - (Z_SOURCE_EXT_MESH))**2)) < RADIUS_TO_MUTE) &
+ .and. MUTE_SOURCE) then
- display(i,j) = 0.
- else
-
- if(inorm == 1) then
- display(i,j) = vectorx
- else if(inorm == 2) then
- display(i,j) = vectory
+ display(i,j) = 0.
+ else
+ ! chooses norm
+ if(inorm == 1) then
+ ! norm displacement
+ display(i,j) = vectorx
+ else if(inorm == 2) then
+ ! norm velocity
+ display(i,j) = vectory
+ else
+ ! norm acceleration
+ display(i,j) = vectorz
+ endif
+ endif
else
- display(i,j) = vectorz
+ ! movie
+ if(inorm == 1) then
+ ! norm of velocity
+ display(i,j) = sqrt(vectorz**2+vectory**2+vectorx**2)
+ else if( inorm == 2 ) then
+ ! velocity x-component
+ display(i,j) = vectorx
+ else if( inorm == 3 ) then
+ ! velocity y-component
+ display(i,j) = vectory
+ else if( inorm == 4 ) then
+ ! velocity z-component
+ display(i,j) = vectorz
+ endif
endif
- endif
- else
- display(i,j) = sqrt(vectorz**2+vectory**2+vectorx**2)
- endif
- enddo
- enddo
+ enddo
+ enddo
-! assign the values of the corners of the OpenDX "elements"
- ispec = ispec + 1
- ielm = (NGLLX-1)*(NGLLY-1)*(ispec-1)
+ ! assign the values of the corners of the OpenDX "elements"
+ ispec = ispec + 1
+ ielm = (NGLLX-1)*(NGLLY-1)*(ispec-1)
- do j = 1,NGLLY-1
- do i = 1,NGLLX-1
- ieoff = NGNOD2D_AVS_DX*(ielm+(i-1)+(j-1)*(NGLLX-1))
- do ilocnum = 1,NGNOD2D_AVS_DX
-! do k = 1,NGNOD2D_AVS_DX
+ do j = 1,NGLLY-1
+ do i = 1,NGLLX-1
+ ieoff = NGNOD2D_AVS_DX*(ielm+(i-1)+(j-1)*(NGLLX-1))
+ do ilocnum = 1,NGNOD2D_AVS_DX
+ ! do k = 1,NGNOD2D_AVS_DX
- if(ilocnum == 1) then
- xp(ieoff+ilocnum) = dble(x(i,j))
- yp(ieoff+ilocnum) = dble(y(i,j))
- zp(ieoff+ilocnum) = dble(z(i,j))
- field_display(ieoff+ilocnum) = dble(display(i,j))
- elseif(ilocnum == 2) then
+ if(ilocnum == 1) then
+ xp(ieoff+ilocnum) = dble(x(i,j))
+ yp(ieoff+ilocnum) = dble(y(i,j))
+ zp(ieoff+ilocnum) = dble(z(i,j))
+ field_display(ieoff+ilocnum) = dble(display(i,j))
+ elseif(ilocnum == 2) then
-! accounts for different ordering of square points
- xp(ieoff+ilocnum) = dble(x(i+1,j+1))
- yp(ieoff+ilocnum) = dble(y(i+1,j+1))
- zp(ieoff+ilocnum) = dble(z(i+1,j+1))
- field_display(ieoff+ilocnum) = dble(display(i+1,j+1))
+ ! accounts for different ordering of square points
+ xp(ieoff+ilocnum) = dble(x(i+1,j+1))
+ yp(ieoff+ilocnum) = dble(y(i+1,j+1))
+ zp(ieoff+ilocnum) = dble(z(i+1,j+1))
+ field_display(ieoff+ilocnum) = dble(display(i+1,j+1))
-! xp(ieoff+ilocnum) = dble(x(i+1,j))
-! yp(ieoff+ilocnum) = dble(y(i+1,j))
-! zp(ieoff+ilocnum) = dble(z(i+1,j))
-! field_display(ieoff+ilocnum) = dble(display(i+1,j))
+ ! xp(ieoff+ilocnum) = dble(x(i+1,j))
+ ! yp(ieoff+ilocnum) = dble(y(i+1,j))
+ ! zp(ieoff+ilocnum) = dble(z(i+1,j))
+ ! field_display(ieoff+ilocnum) = dble(display(i+1,j))
- elseif(ilocnum == 3) then
+ elseif(ilocnum == 3) then
-! accounts for different ordering of square points
- xp(ieoff+ilocnum) = dble(x(i+1,j))
- yp(ieoff+ilocnum) = dble(y(i+1,j))
- zp(ieoff+ilocnum) = dble(z(i+1,j))
- field_display(ieoff+ilocnum) = dble(display(i+1,j))
+ ! accounts for different ordering of square points
+ xp(ieoff+ilocnum) = dble(x(i+1,j))
+ yp(ieoff+ilocnum) = dble(y(i+1,j))
+ zp(ieoff+ilocnum) = dble(z(i+1,j))
+ field_display(ieoff+ilocnum) = dble(display(i+1,j))
-! xp(ieoff+ilocnum) = dble(x(i+1,j+1))
-! yp(ieoff+ilocnum) = dble(y(i+1,j+1))
-! zp(ieoff+ilocnum) = dble(z(i+1,j+1))
-! field_display(ieoff+ilocnum) = dble(display(i+1,j+1))
- else
- xp(ieoff+ilocnum) = dble(x(i,j+1))
- yp(ieoff+ilocnum) = dble(y(i,j+1))
- zp(ieoff+ilocnum) = dble(z(i,j+1))
- field_display(ieoff+ilocnum) = dble(display(i,j+1))
- endif
+ ! xp(ieoff+ilocnum) = dble(x(i+1,j+1))
+ ! yp(ieoff+ilocnum) = dble(y(i+1,j+1))
+ ! zp(ieoff+ilocnum) = dble(z(i+1,j+1))
+ ! field_display(ieoff+ilocnum) = dble(display(i+1,j+1))
+ else
+ xp(ieoff+ilocnum) = dble(x(i,j+1))
+ yp(ieoff+ilocnum) = dble(y(i,j+1))
+ zp(ieoff+ilocnum) = dble(z(i,j+1))
+ field_display(ieoff+ilocnum) = dble(display(i,j+1))
+ endif
+ enddo
+
+ !if( j==1 .and. ispec==1) then
+ !print*,'p1',xp(ieoff+1),yp(ieoff+1),zp(ieoff+1)
+ !print*,'p2',xp(ieoff+2),yp(ieoff+2),zp(ieoff+2)
+ !print*,'p3',xp(ieoff+3),yp(ieoff+3),zp(ieoff+3)
+ !print*,'p4',xp(ieoff+4),yp(ieoff+4),zp(ieoff+4)
+ !endif
+
+ enddo
enddo
-
- !if( j==1 .and. ispec==1) then
- !print*,'p1',xp(ieoff+1),yp(ieoff+1),zp(ieoff+1)
- !print*,'p2',xp(ieoff+2),yp(ieoff+2),zp(ieoff+2)
- !print*,'p3',xp(ieoff+3),yp(ieoff+3),zp(ieoff+3)
- !print*,'p4',xp(ieoff+4),yp(ieoff+4),zp(ieoff+4)
- !endif
-
- enddo
- enddo
- else
+ else
+ ! low-resolution (only spectral element corners)
+ ispec = ispec + 1
+ ieoff = NGNOD2D_AVS_DX*(ispec-1)
- ispec = ispec + 1
- ieoff = NGNOD2D_AVS_DX*(ispec-1)
+ ! four points for each element
+ do i = 1,NGNOD2D_AVS_DX
-! four points for each element
- do i = 1,NGNOD2D_AVS_DX
+ ! accounts for different ordering of square points
+ ilocnum = iorder(i)
+
+ ipoin = ipoin + 1
- ! accounts for different ordering of square points
- ilocnum = iorder(i)
-
- ipoin = ipoin + 1
+ xcoord = store_val_x(ipoin,iproc)
+ ycoord = store_val_y(ipoin,iproc)
+ zcoord = store_val_z(ipoin,iproc)
- xcoord = store_val_x(ipoin,iproc)
- ycoord = store_val_y(ipoin,iproc)
- zcoord = store_val_z(ipoin,iproc)
+ vectorx = store_val_ux(ipoin,iproc)
+ vectory = store_val_uy(ipoin,iproc)
+ vectorz = store_val_uz(ipoin,iproc)
- vectorx = store_val_ux(ipoin,iproc)
- vectory = store_val_uy(ipoin,iproc)
- vectorz = store_val_uz(ipoin,iproc)
+ xp(ilocnum+ieoff) = dble(xcoord)
+ yp(ilocnum+ieoff) = dble(ycoord)
+ zp(ilocnum+ieoff) = dble(zcoord)
- xp(ilocnum+ieoff) = dble(xcoord)
- yp(ilocnum+ieoff) = dble(ycoord)
- zp(ilocnum+ieoff) = dble(zcoord)
-
-! show vertical component of displacement or velocity in the movie
-! or show norm of vector if shaking map
-! for shaking map, norm of U stored in ux, V in uy and A in uz
- if(plot_shaking_map) then
-!!!! NL NL mute value near source
- if ( (sqrt(((dble(xcoord) - (X_SOURCE_EXT_MESH))**2 + &
- (dble(ycoord) - (Y_SOURCE_EXT_MESH))**2 + &
- (dble(zcoord) - (Z_SOURCE_EXT_MESH))**2)) < RADIUS_TO_MUTE) &
- .and. MUTE_SOURCE) then
- field_display(ilocnum+ieoff) = 0.
- else
- if(inorm == 1) then
- field_display(ilocnum+ieoff) = dble(vectorx)
- else if(inorm == 2) then
- field_display(ilocnum+ieoff) = dble(vectory)
+ ! shakemap
+ if(plot_shaking_map) then
+ !!!! NL NL mute value near source
+ if ( (sqrt(((dble(xcoord) - (X_SOURCE_EXT_MESH))**2 + &
+ (dble(ycoord) - (Y_SOURCE_EXT_MESH))**2 + &
+ (dble(zcoord) - (Z_SOURCE_EXT_MESH))**2)) < RADIUS_TO_MUTE) &
+ .and. MUTE_SOURCE) then
+ field_display(ilocnum+ieoff) = 0.
+ else
+ if(inorm == 1) then
+ ! norm of displacement
+ field_display(ilocnum+ieoff) = dble(vectorx)
+ else if(inorm == 2) then
+ ! norm of velocity
+ field_display(ilocnum+ieoff) = dble(vectory)
+ else
+ ! norm of acceleration
+ field_display(ilocnum+ieoff) = dble(vectorz)
+ endif
+ endif
else
- field_display(ilocnum+ieoff) = dble(vectorz)
+ ! movie
+ if(inorm == 1) then
+ ! norm of velocity
+ field_display(ilocnum+ieoff) = sqrt(vectorz**2+vectory**2+vectorx**2)
+ else if( inorm == 2 ) then
+ ! velocity x-component
+ field_display(ilocnum+ieoff) = vectorx
+ else if( inorm == 3 ) then
+ ! velocity y-component
+ field_display(ilocnum+ieoff) = vectory
+ else
+ ! velocity z-component
+ field_display(ilocnum+ieoff) = vectorz
+ endif
+ ! takes norm of velocity vector
+ !field_display(ilocnum+ieoff) =sqrt(vectorz**2+vectory**2+vectorx**2)
endif
- endif
- else
- ! takes norm of velocity vector
- field_display(ilocnum+ieoff) =sqrt(vectorz**2+vectory**2+vectorx**2)
- endif
- enddo
+ enddo
- endif
+ endif ! USE_HIGHRES_FOR_MOVIES
- enddo
- enddo
+ enddo ! NSPEC_SURFACE_EXT_MESH
+ enddo ! NPROC
-! copy coordinate arrays since the sorting routine does not preserve them
- xp_save(:) = xp(:)
- yp_save(:) = yp(:)
- zp_save(:) = zp(:)
+ ! copy coordinate arrays since the sorting routine does not preserve them
+ xp_save(:) = xp(:)
+ yp_save(:) = yp(:)
+ zp_save(:) = zp(:)
-!--- sort the list based upon coordinates to get rid of multiples
- print *,'sorting list of points'
- call get_global_AVS(nspectot_AVS_max,xp,yp,zp,iglob,loc,ifseg,nglob,npointot, &
- dble(minval(store_val_x(:,0))),dble(maxval(store_val_x(:,0))))
+ !--- sort the list based upon coordinates to get rid of multiples
+ print *,'sorting list of points'
+ call get_global_AVS(nspectot_AVS_max,xp,yp,zp,iglob,loc,ifseg,nglob,npointot, &
+ dble(minval(store_val_x(:,0))),dble(maxval(store_val_x(:,0))))
-!--- print total number of points found
- print *
- print *,'found a total of ',nglob,' points'
- print *,'initial number of points (with multiples) was ',npointot
+ !--- print total number of points found
+ print *
+ print *,'found a total of ',nglob,' points'
+ print *,'initial number of points (with multiples) was ',npointot
-!--- normalize and scale vector field
+ !--- normalize and scale vector field
-! compute min and max of data value to normalize
- min_field_current = minval(field_display(:))
- max_field_current = maxval(field_display(:))
+ ! compute min and max of data value to normalize
+ min_field_current = minval(field_display(:))
+ max_field_current = maxval(field_display(:))
-! print minimum and maximum amplitude in current snapshot
- print *
- print *,'minimum amplitude in current snapshot = ',min_field_current
- print *,'maximum amplitude in current snapshot = ',max_field_current
- print *
+ ! print minimum and maximum amplitude in current snapshot
+ print *
+ print *,'minimum amplitude in current snapshot = ',min_field_current
+ print *,'maximum amplitude in current snapshot = ',max_field_current
+ print *
-!-----------------------------------------
+ !-----------------------------------------
- if(plot_shaking_map) then
+ if(plot_shaking_map) then
-! compute min and max of data value to normalize
- min_field_current = minval(field_display(:))
- max_field_current = maxval(field_display(:))
+ ! compute min and max of data value to normalize
+ min_field_current = minval(field_display(:))
+ max_field_current = maxval(field_display(:))
-! print minimum and maximum amplitude in current snapshot
- print *
- print *,'minimum amplitude in current snapshot after removal = ',min_field_current
- print *,'maximum amplitude in current snapshot after removal = ',max_field_current
- print *
+ ! print minimum and maximum amplitude in current snapshot
+ print *
+ print *,'minimum amplitude in current snapshot after removal = ',min_field_current
+ print *,'maximum amplitude in current snapshot after removal = ',max_field_current
+ print *
- endif
+ endif
-!-----------------------------------------
+ !-----------------------------------------
-! apply scaling in all cases for movies
- if(.not. plot_shaking_map) then
+ ! apply scaling in all cases for movies
+ if(.not. plot_shaking_map) then
-! make sure range is always symmetric and center is in zero
-! this assumption works only for fields that can be negative
-! would not work for norm of vector for instance
-! (we would lose half of the color palette if no negative values)
- max_absol = max(abs(min_field_current),abs(max_field_current))
- min_field_current = - max_absol
- max_field_current = + max_absol
+ ! make sure range is always symmetric and center is in zero
+ ! this assumption works only for fields that can be negative
+ ! would not work for norm of vector for instance
+ ! (we would lose half of the color palette if no negative values)
+ max_absol = max(abs(min_field_current),abs(max_field_current))
+ min_field_current = - max_absol
+ max_field_current = + max_absol
-! normalize field to [0:1]
- field_display(:) = (field_display(:) - min_field_current) / (max_field_current - min_field_current)
+ ! normalize field to [0:1]
+ field_display(:) = (field_display(:) - min_field_current) / (max_field_current - min_field_current)
-! rescale to [-1,1]
- field_display(:) = 2.*field_display(:) - 1.
+ ! rescale to [-1,1]
+ field_display(:) = 2.*field_display(:) - 1.
-! apply threshold to normalized field
- if(APPLY_THRESHOLD) &
- where(abs(field_display(:)) <= THRESHOLD) field_display = 0.
+ ! apply threshold to normalized field
+ if(APPLY_THRESHOLD) &
+ where(abs(field_display(:)) <= THRESHOLD) field_display = 0.
-! apply non linear scaling to normalized field if needed
- if(NONLINEAR_SCALING) then
- where(field_display(:) >= 0.)
- field_display = field_display ** POWER_SCALING
- elsewhere
- field_display = - abs(field_display) ** POWER_SCALING
- endwhere
- endif
+ ! apply non linear scaling to normalized field if needed
+ if(NONLINEAR_SCALING) then
+ where(field_display(:) >= 0.)
+ field_display = field_display ** POWER_SCALING
+ elsewhere
+ field_display = - abs(field_display) ** POWER_SCALING
+ endwhere
+ endif
-! map back to [0,1]
- field_display(:) = (field_display(:) + 1.) / 2.
+ ! map back to [0,1]
+ field_display(:) = (field_display(:) + 1.) / 2.
-! map field to [0:255] for AVS color scale
- field_display(:) = 255. * field_display(:)
+ ! map field to [0:255] for AVS color scale
+ field_display(:) = 255. * field_display(:)
-! apply scaling only if selected for shaking map
+ ! apply scaling only if selected for shaking map
- else if(NONLINEAR_SCALING .and. iscaling_shake == 1) then
+ else if(NONLINEAR_SCALING .and. iscaling_shake == 1) then
-! normalize field to [0:1]
- field_display(:) = field_display(:) / max_field_current
+ ! normalize field to [0:1]
+ field_display(:) = field_display(:) / max_field_current
-! apply non linear scaling to normalized field
- field_display = field_display ** POWER_SCALING
+ ! apply non linear scaling to normalized field
+ field_display = field_display ** POWER_SCALING
-! map field to [0:255] for AVS color scale
- field_display(:) = 255. * field_display(:)
+ ! map field to [0:255] for AVS color scale
+ field_display(:) = 255. * field_display(:)
- endif
+ endif
-!--- ****** create AVS file using sorted list ******
+ !--- ****** create AVS file using sorted list ******
- if(.not. plot_shaking_map) then
- if(inumber == 1) then
- ivalue = iframe
- else
- ivalue = it
- endif
- endif
+ if(.not. plot_shaking_map) then
+ if(inumber == 1) then
+ ivalue = iframe
+ else
+ ivalue = it
+ endif
+ endif
-! create file name and open file
- if(plot_shaking_map) then
+ ! create file name and open file
+ if(plot_shaking_map) then
- if(USE_OPENDX) then
- write(outputname,"('/DX_shaking_map.dx')")
- open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
- write(11,*) 'object 1 class array type float rank 1 shape 3 items ',nglob,' data follows'
- else if(USE_AVS) then
- write(outputname,"('/AVS_shaking_map.inp')")
- open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
- write(11,*) nglob,' ',nspectot_AVS_max,' 1 0 0'
- else
- stop 'wrong output format selected'
- endif
+ if(USE_OPENDX) then
+ write(outputname,"('/DX_shaking_map.dx')")
+ open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
+ write(11,*) 'object 1 class array type float rank 1 shape 3 items ',nglob,' data follows'
+ else if(USE_AVS) then
+ write(outputname,"('/AVS_shaking_map.inp')")
+ open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
+ write(11,*) nglob,' ',nspectot_AVS_max,' 1 0 0'
+ else
+ stop 'wrong output format selected'
+ endif
- else
+ else
- if(USE_OPENDX) then
- write(outputname,"('/DX_movie_',i6.6,'.dx')") ivalue
- open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
- write(11,*) 'object 1 class array type float rank 1 shape 3 items ',nglob,' data follows'
- else if(USE_AVS) then
- write(outputname,"('/AVS_movie_',i6.6,'.inp')") ivalue
- open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
- write(11,*) nglob,' ',nspectot_AVS_max,' 1 0 0'
- else
- stop 'wrong output format selected'
- endif
+ if(USE_OPENDX) then
+ write(outputname,"('/DX_movie_',i6.6,'.dx')") ivalue
+ open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
+ write(11,*) 'object 1 class array type float rank 1 shape 3 items ',nglob,' data follows'
+ else if(USE_AVS) then
+ write(outputname,"('/AVS_movie_',i6.6,'.inp')") ivalue
+ open(unit=11,file=trim(OUTPUT_FILES)//outputname,status='unknown')
+ write(11,*) nglob,' ',nspectot_AVS_max,' 1 0 0'
+ else
+ stop 'wrong output format selected'
+ endif
- endif
+ endif
- if(.false.) then
- else
-
-! output list of points
- mask_point = .false.
- ipoin = 0
- do ispec=1,nspectot_AVS_max
- ieoff = NGNOD2D_AVS_DX*(ispec-1)
-! four points for each element
- do ilocnum = 1,NGNOD2D_AVS_DX
- ibool_number = iglob(ilocnum+ieoff)
- if(.not. mask_point(ibool_number)) then
- ipoin = ipoin + 1
- ireorder(ibool_number) = ipoin
- if(USE_OPENDX) then
- write(11,*) xp_save(ilocnum+ieoff),yp_save(ilocnum+ieoff),zp_save(ilocnum+ieoff)
- else if(USE_AVS) then
- write(11,'(i,3f)') ireorder(ibool_number),xp_save(ilocnum+ieoff), &
- yp_save(ilocnum+ieoff),zp_save(ilocnum+ieoff)
- endif
- endif
- mask_point(ibool_number) = .true.
- enddo
- enddo
-
- if(USE_OPENDX) &
- write(11,*) 'object 2 class array type int rank 1 shape 4 items ',nspectot_AVS_max,' data follows'
-
-! output list of elements
- do ispec=1,nspectot_AVS_max
- ieoff = NGNOD2D_AVS_DX*(ispec-1)
-! four points for each element
- ibool_number1 = iglob(ieoff + 1)
- ibool_number2 = iglob(ieoff + 2)
- ibool_number3 = iglob(ieoff + 3)
- ibool_number4 = iglob(ieoff + 4)
- if(USE_OPENDX) then
-! point order in OpenDX is 1,4,2,3 *not* 1,2,3,4 as in AVS
- write(11,"(i10,1x,i10,1x,i10,1x,i10)") ireorder(ibool_number1)-1, &
- ireorder(ibool_number4)-1,ireorder(ibool_number2)-1,ireorder(ibool_number3)-1
+ if(.false.) then
+ ! GMT format not implemented yet
else
- write(11,"(i10,' 1 quad ',i10,1x,i10,1x,i10,1x,i10)") ispec,ireorder(ibool_number1), &
- ireorder(ibool_number4),ireorder(ibool_number2),ireorder(ibool_number3)
- endif
- enddo
- if(USE_OPENDX) then
- write(11,*) 'attribute "element type" string "quads"'
- write(11,*) 'attribute "ref" string "positions"'
- write(11,*) 'object 3 class array type float rank 0 items ',nglob,' data follows'
- else
-! dummy text for labels
- write(11,*) '1 1'
- write(11,*) 'a, b'
- endif
+ ! output list of points
+ mask_point = .false.
+ ipoin = 0
+ do ispec=1,nspectot_AVS_max
+ ieoff = NGNOD2D_AVS_DX*(ispec-1)
+ ! four points for each element
+ do ilocnum = 1,NGNOD2D_AVS_DX
+ ibool_number = iglob(ilocnum+ieoff)
+ if(.not. mask_point(ibool_number)) then
+ ipoin = ipoin + 1
+ ireorder(ibool_number) = ipoin
+ if(USE_OPENDX) then
+ write(11,*) xp_save(ilocnum+ieoff),yp_save(ilocnum+ieoff),zp_save(ilocnum+ieoff)
+ else if(USE_AVS) then
+ write(11,'(i,3f)') ireorder(ibool_number),xp_save(ilocnum+ieoff), &
+ yp_save(ilocnum+ieoff),zp_save(ilocnum+ieoff)
+ endif
+ endif
+ mask_point(ibool_number) = .true.
+ enddo
+ enddo
-! output data values
- mask_point = .false.
+ if(USE_OPENDX) &
+ write(11,*) 'object 2 class array type int rank 1 shape 4 items ',nspectot_AVS_max,' data follows'
-! output point data
- do ispec=1,nspectot_AVS_max
- ieoff = NGNOD2D_AVS_DX*(ispec-1)
-! four points for each element
- do ilocnum = 1,NGNOD2D_AVS_DX
- ibool_number = iglob(ilocnum+ieoff)
- if(.not. mask_point(ibool_number)) then
+ ! output list of elements
+ do ispec=1,nspectot_AVS_max
+ ieoff = NGNOD2D_AVS_DX*(ispec-1)
+ ! four points for each element
+ ibool_number1 = iglob(ieoff + 1)
+ ibool_number2 = iglob(ieoff + 2)
+ ibool_number3 = iglob(ieoff + 3)
+ ibool_number4 = iglob(ieoff + 4)
if(USE_OPENDX) then
- if(plot_shaking_map) then
- write(11,*) sngl(field_display(ilocnum+ieoff))
- else
- write(11,"(f7.2)") field_display(ilocnum+ieoff)
- endif
+ ! point order in OpenDX is 1,4,2,3 *not* 1,2,3,4 as in AVS
+ write(11,"(i10,1x,i10,1x,i10,1x,i10)") ireorder(ibool_number1)-1, &
+ ireorder(ibool_number4)-1,ireorder(ibool_number2)-1,ireorder(ibool_number3)-1
else
- if(plot_shaking_map) then
- write(11,*) ireorder(ibool_number),field_display(ilocnum+ieoff)
- else
- write(11,"(i10,1x,f7.2)") ireorder(ibool_number),field_display(ilocnum+ieoff)
- endif
+ write(11,"(i10,' 1 quad ',i10,1x,i10,1x,i10,1x,i10)") ispec,ireorder(ibool_number1), &
+ ireorder(ibool_number4),ireorder(ibool_number2),ireorder(ibool_number3)
endif
+ enddo
+
+ if(USE_OPENDX) then
+ write(11,*) 'attribute "element type" string "quads"'
+ write(11,*) 'attribute "ref" string "positions"'
+ write(11,*) 'object 3 class array type float rank 0 items ',nglob,' data follows'
+ else
+ ! dummy text for labels
+ write(11,*) '1 1'
+ write(11,*) 'a, b'
endif
- mask_point(ibool_number) = .true.
- enddo
- enddo
-! define OpenDX field
- if(USE_OPENDX) then
- write(11,*) 'attribute "dep" string "positions"'
- write(11,*) 'object "irregular positions irregular connections" class field'
- write(11,*) 'component "positions" value 1'
- write(11,*) 'component "connections" value 2'
- write(11,*) 'component "data" value 3'
- write(11,*) 'end'
- endif
+ ! output data values
+ mask_point = .false.
-! end of test for GMT format
- endif
+ ! output point data
+ do ispec=1,nspectot_AVS_max
+ ieoff = NGNOD2D_AVS_DX*(ispec-1)
+ ! four points for each element
+ do ilocnum = 1,NGNOD2D_AVS_DX
+ ibool_number = iglob(ilocnum+ieoff)
+ if(.not. mask_point(ibool_number)) then
+ if(USE_OPENDX) then
+ if(plot_shaking_map) then
+ write(11,*) sngl(field_display(ilocnum+ieoff))
+ else
+ write(11,"(f7.2)") field_display(ilocnum+ieoff)
+ endif
+ else
+ if(plot_shaking_map) then
+ write(11,*) ireorder(ibool_number),field_display(ilocnum+ieoff)
+ else
+ write(11,"(i10,1x,f7.2)") ireorder(ibool_number),field_display(ilocnum+ieoff)
+ endif
+ endif
+ endif
+ mask_point(ibool_number) = .true.
+ enddo
+ enddo
- close(11)
+ ! define OpenDX field
+ if(USE_OPENDX) then
+ write(11,*) 'attribute "dep" string "positions"'
+ write(11,*) 'object "irregular positions irregular connections" class field'
+ write(11,*) 'component "positions" value 1'
+ write(11,*) 'component "connections" value 2'
+ write(11,*) 'component "data" value 3'
+ write(11,*) 'end'
+ endif
-! end of loop and test on all the time steps for all the movie images
- endif
+ ! end of test for GMT format
+ endif
+
+ close(11)
+
+ ! end of loop and test on all the time steps for all the movie images
+ endif
enddo ! it
print *
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/create_regions_mesh.f90 2009-10-25 23:06:34 UTC (rev 15876)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/create_regions_mesh.f90 2009-10-25 23:54:42 UTC (rev 15877)
@@ -39,7 +39,8 @@
ibelm_xmin, ibelm_xmax, ibelm_ymin, ibelm_ymax, ibelm_bottom, ibelm_top, &
nodes_ibelm_xmin,nodes_ibelm_xmax,nodes_ibelm_ymin,nodes_ibelm_ymax,&
nodes_ibelm_bottom,nodes_ibelm_top, &
- SAVE_MESH_FILES,nglob)
+ SAVE_MESH_FILES,nglob, &
+ ANISOTROPY)
! create the different regions of the mesh
@@ -104,6 +105,8 @@
logical :: SAVE_MESH_FILES
integer :: nglob
+ logical :: ANISOTROPY
+
! local parameters
!-----------------------
@@ -135,7 +138,7 @@
etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore,jacobianstore
! for model density
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rhostore,kappastore,mustore,vpstore,vsstore
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rhostore,kappastore,mustore !,vpstore,vsstore
! mass matrix
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass
@@ -191,6 +194,15 @@
character(len=150) prname
character(len=150) prname_file
+! anisotropy
+ integer :: NSPEC_ANISO
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
+ c11store,c12store,c13store,c14store,c15store,c16store,&
+ c22store,c23store,c24store,c25store,c26store,c33store,&
+ c34store,c35store,c36store,c44store,c45store,c46store,&
+ c55store,c56store,c66store
+
+
! mask to sort ibool
! integer, dimension(:), allocatable :: mask_ibool
! integer, dimension(:,:,:,:), allocatable :: copy_ibool_ori
@@ -297,9 +309,9 @@
! array with model density
allocate(rhostore(NGLLX,NGLLY,NGLLZ,nspec), &
kappastore(NGLLX,NGLLY,NGLLZ,nspec), &
- mustore(NGLLX,NGLLY,NGLLZ,nspec), &
- vpstore(NGLLX,NGLLY,NGLLZ,nspec), &
- vsstore(NGLLX,NGLLY,NGLLZ,nspec),stat=ier) !pll
+ mustore(NGLLX,NGLLY,NGLLZ,nspec),stat=ier) !pll
+ !vpstore(NGLLX,NGLLY,NGLLZ,nspec), &
+ !vsstore(NGLLX,NGLLY,NGLLZ,nspec),
if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
! arrays with mesh parameters
@@ -356,6 +368,34 @@
absorbing_boundary_normal(NDIM,NGLLSQUARE,num_absorbing_boundary_faces),stat=ier)
if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
+! array with anisotropy
+ if( ANISOTROPY ) then
+ NSPEC_ANISO = nspec
+ else
+ NSPEC_ANISO = 1
+ endif
+ allocate(c11store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c12store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c13store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c14store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c15store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c16store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c22store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c23store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c24store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c25store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c26store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c33store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c34store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c35store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c36store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c44store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c45store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c46store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c55store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c56store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO), &
+ c66store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO),stat=ier)
+ if(ier /= 0) call exit_MPI(myrank,'not enough memory to allocate arrays')
! fills location and weights for Gauss-Lobatto-Legendre points, shape and derivations,
! returns jacobianstore,xixstore,...gammazstore
@@ -378,15 +418,21 @@
! sets material velocities
call sync_all()
if( myrank == 0) then
- write(IMAIN,*) ' ...determining kappa and mu parameters'
+ write(IMAIN,*) ' ...determining velocity model'
endif
- call create_regions_mesh_ext_mesh_determine_kappamu(nspec,mat_ext_mesh,nelmnts_ext_mesh,&
- materials_ext_mesh,nmat_ext_mesh,&
- undef_mat_prop,nundefMat_ext_mesh,&
- rhostore,kappastore,mustore,vpstore,vsstore,&
- iflag_attenuation_store,rho_vp,rho_vs)
-
+ call create_regions_mesh_ext_mesh_determine_velocity(nspec,mat_ext_mesh,nelmnts_ext_mesh, &
+ materials_ext_mesh,nmat_ext_mesh, &
+ undef_mat_prop,nundefMat_ext_mesh, &
+ rhostore,kappastore,mustore, &
+ iflag_attenuation_store,rho_vp,rho_vs, &
+ ANISOTROPY,NSPEC_ANISO, &
+ c11store,c12store,c13store,c14store,c15store,c16store, &
+ c22store,c23store,c24store,c25store,c26store,c33store, &
+ c34store,c35store,c36store,c44store,c45store,c46store, &
+ c55store,c56store,c66store)
+ !,vpstore,vsstore,
+
! creates ibool index array for projection from local to global points
call sync_all()
if( myrank == 0) then
@@ -479,7 +525,12 @@
num_absorbing_boundary_faces, &
num_interfaces_ext_mesh,my_neighbours_ext_mesh,nibool_interfaces_ext_mesh, &
max_interface_size_ext_mesh,ibool_interfaces_ext_mesh, &
- prname,SAVE_MESH_FILES)
+ prname,SAVE_MESH_FILES, &
+ ANISOTROPY,NSPEC_ANISO, &
+ c11store,c12store,c13store,c14store,c15store,c16store, &
+ c22store,c23store,c24store,c25store,c26store,c33store, &
+ c34store,c35store,c36store,c44store,c45store,c46store, &
+ c55store,c56store,c66store)
! computes the approximate amount of static memory needed to run the solver
call memory_eval(nspec,nglob,maxval(nibool_interfaces_ext_mesh),num_interfaces_ext_mesh,static_memory_size)
@@ -865,11 +916,17 @@
!----
!
-subroutine create_regions_mesh_ext_mesh_determine_kappamu(nspec,mat_ext_mesh,nelmnts_ext_mesh,&
- materials_ext_mesh,nmat_ext_mesh,&
- undef_mat_prop,nundefMat_ext_mesh,&
- rhostore,kappastore,mustore,vpstore,vsstore,&
- iflag_attenuation_store,rho_vp,rho_vs)
+subroutine create_regions_mesh_ext_mesh_determine_velocity(nspec,mat_ext_mesh,nelmnts_ext_mesh, &
+ materials_ext_mesh,nmat_ext_mesh, &
+ undef_mat_prop,nundefMat_ext_mesh, &
+ rhostore,kappastore,mustore, &
+ iflag_attenuation_store,rho_vp,rho_vs, &
+ ANISOTROPY,NSPEC_ANISO, &
+ c11store,c12store,c13store,c14store,c15store,c16store, &
+ c22store,c23store,c24store,c25store,c26store,c33store, &
+ c34store,c35store,c36store,c44store,c45store,c46store, &
+ c55store,c56store,c66store)
+ ! vpstore,vsstore,
implicit none
@@ -887,17 +944,32 @@
character (len=30), dimension(5,nundefMat_ext_mesh):: undef_mat_prop
! for model density
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: rhostore, &
- kappastore,mustore,vpstore,vsstore
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: &
+ rhostore,kappastore,mustore !,vpstore,vsstore
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: rho_vp,rho_vs
! attenuation
integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: iflag_attenuation_store
+! anisotropy
+ logical :: ANISOTROPY
+ integer :: NSPEC_ANISO
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO) :: &
+ c11store,c12store,c13store,c14store,c15store,c16store, &
+ c22store,c23store,c24store,c25store,c26store,c33store, &
+ c34store,c35store,c36store,c44store,c45store,c46store, &
+ c55store,c56store,c66store
+
+
! local parameters
- integer :: ispec,i,j,k,iundef
- integer :: iflag,flag_below,flag_above
-
+ real(kind=CUSTOM_REAL) :: vp,vs,rho
+ real(kind=CUSTOM_REAL) :: c11,c12,c13,c14,c15,c16,c22,c23,c24,c25, &
+ c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
+
+ integer :: ispec,i,j,k,iundef,iflag_atten
+ integer :: iflag,flag_below,flag_above
+ integer :: iflag_aniso
+
! ! Piero, read bedrock file
! allocate(ibedrock(NX_TOPO_ANT,NY_TOPO_ANT))
! if(myrank == 0) then
@@ -909,25 +981,38 @@
! ! broadcast the information read on the master to the nodes
! ! call MPI_BCAST(ibedrock,NX_TOPO_ANT*NY_TOPO_ANT,MPI_REAL,0,MPI_COMM_WORLD,ier)
! call bcast_all_cr(ibedrock,NX_TOPO_ANT*NY_TOPO_ANT)
-
-! kappastore and mustore
+
+! material properties on all GLL points: taken from material values defined for
+! each spectral element in input mesh
do ispec = 1, nspec
do k = 1, NGLLZ
do j = 1, NGLLY
do i = 1, NGLLX
! check if the material is known or unknown
if (mat_ext_mesh(1,ispec) > 0) then
+ ! density
! materials_ext_mesh format: #index1 = rho #index2 = vp #index3 = vs #index4 = Q_flag #index5 = 0
- rhostore(i,j,k,ispec) = materials_ext_mesh(1,mat_ext_mesh(1,ispec))
- vpstore(i,j,k,ispec) = materials_ext_mesh(2,mat_ext_mesh(1,ispec))
- vsstore(i,j,k,ispec) = materials_ext_mesh(3,mat_ext_mesh(1,ispec))
- iflag_attenuation_store(i,j,k,ispec) = materials_ext_mesh(4,mat_ext_mesh(1,ispec))
+ !rhostore(i,j,k,ispec) = materials_ext_mesh(1,mat_ext_mesh(1,ispec))
+ rho = materials_ext_mesh(1,mat_ext_mesh(1,ispec))
+
+ ! isotropic values: vp, vs
+ !vpstore(i,j,k,ispec) = materials_ext_mesh(2,mat_ext_mesh(1,ispec))
+ !vsstore(i,j,k,ispec) = materials_ext_mesh(3,mat_ext_mesh(1,ispec))
+ vp = materials_ext_mesh(2,mat_ext_mesh(1,ispec))
+ vs = materials_ext_mesh(3,mat_ext_mesh(1,ispec))
+
+ ! attenuation
+ iflag_atten = materials_ext_mesh(4,mat_ext_mesh(1,ispec))
!change for piero :
!if(mat_ext_mesh(1,ispec) == 1) then
! iflag_attenuation_store(i,j,k,ispec) = 1
!else
! iflag_attenuation_store(i,j,k,ispec) = 2
!endif
+
+ ! anisotropy
+ iflag_aniso = materials_ext_mesh(5,mat_ext_mesh(1,ispec))
+
else if (mat_ext_mesh(2,ispec) == 1) then
stop 'material: interface not implemented yet'
@@ -939,30 +1024,71 @@
enddo
!call interface(iflag,flag_below,flag_above,ispec,nspec,i,j,k,xstore,ystore,zstore,ibedrock)
iflag = 1
- rhostore(i,j,k,ispec) = materials_ext_mesh(1,iflag)
- vpstore(i,j,k,ispec) = materials_ext_mesh(2,iflag)
- vsstore(i,j,k,ispec) = materials_ext_mesh(3,iflag)
- iflag_attenuation_store(i,j,k,ispec) = materials_ext_mesh(4,iflag)
+ !rhostore(i,j,k,ispec) = materials_ext_mesh(1,iflag)
+ !vpstore(i,j,k,ispec) = materials_ext_mesh(2,iflag)
+ !vsstore(i,j,k,ispec) = materials_ext_mesh(3,iflag)
+ rho = materials_ext_mesh(1,iflag)
+ vp = materials_ext_mesh(2,iflag)
+ vs = materials_ext_mesh(3,iflag)
+ iflag_atten = materials_ext_mesh(4,iflag)
!change for piero :
! if(iflag == 1) then
! iflag_attenuation_store(i,j,k,ispec) = 1
! else
! iflag_attenuation_store(i,j,k,ispec) = 2
! endif
+ iflag_aniso = 0
else
stop 'material: tomography not implemented yet'
! call tomography()
end if
- kappastore(i,j,k,ispec) = rhostore(i,j,k,ispec)* &
- ( vpstore(i,j,k,ispec)*vpstore(i,j,k,ispec) &
- - FOUR_THIRDS*vsstore(i,j,k,ispec)*vsstore(i,j,k,ispec) )
-
- mustore(i,j,k,ispec) = rhostore(i,j,k,ispec)*vsstore(i,j,k,ispec)*vsstore(i,j,k,ispec)
-
+! adds anisotropic perturbation to vp, vs
+ if( ANISOTROPY ) then
+ call aniso_model(iflag_aniso,rho,vp,vs,c11,c12,c13,c14,c15,c16, &
+ c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66)
+
+ c11store(i,j,k,ispec) = c11
+ c12store(i,j,k,ispec) = c12
+ c13store(i,j,k,ispec) = c13
+ c14store(i,j,k,ispec) = c14
+ c15store(i,j,k,ispec) = c15
+ c16store(i,j,k,ispec) = c16
+ c22store(i,j,k,ispec) = c22
+ c23store(i,j,k,ispec) = c23
+ c24store(i,j,k,ispec) = c24
+ c25store(i,j,k,ispec) = c25
+ c26store(i,j,k,ispec) = c26
+ c33store(i,j,k,ispec) = c33
+ c34store(i,j,k,ispec) = c34
+ c35store(i,j,k,ispec) = c35
+ c36store(i,j,k,ispec) = c36
+ c44store(i,j,k,ispec) = c44
+ c45store(i,j,k,ispec) = c45
+ c46store(i,j,k,ispec) = c46
+ c55store(i,j,k,ispec) = c55
+ c56store(i,j,k,ispec) = c56
+ c66store(i,j,k,ispec) = c66
+
+ endif
+! density
+ rhostore(i,j,k,ispec) = rho
+
+! kappa, mu
+ !kappastore(i,j,k,ispec) = rhostore(i,j,k,ispec)* &
+ ! ( vpstore(i,j,k,ispec)*vpstore(i,j,k,ispec) &
+ ! - FOUR_THIRDS*vsstore(i,j,k,ispec)*vsstore(i,j,k,ispec) )
+ !mustore(i,j,k,ispec) = rhostore(i,j,k,ispec)*vsstore(i,j,k,ispec)*vsstore(i,j,k,ispec)
+ kappastore(i,j,k,ispec) = rho*( vp*vp - FOUR_THIRDS*vs*vs )
+ mustore(i,j,k,ispec) = rho*vs*vs
+
+! attenuation
+ iflag_attenuation_store(i,j,k,ispec) = iflag_atten
! Stacey, a completer par la suite
- rho_vp(i,j,k,ispec) = rhostore(i,j,k,ispec)*vpstore(i,j,k,ispec)
- rho_vs(i,j,k,ispec) = rhostore(i,j,k,ispec)*vsstore(i,j,k,ispec)
+ !rho_vp(i,j,k,ispec) = rhostore(i,j,k,ispec)*vpstore(i,j,k,ispec)
+ !rho_vs(i,j,k,ispec) = rhostore(i,j,k,ispec)*vsstore(i,j,k,ispec)
+ rho_vp(i,j,k,ispec) = rho*vp
+ rho_vs(i,j,k,ispec) = rho*vs
!end pll
enddo
@@ -1139,7 +1265,7 @@
! enddo
! enddo
-end subroutine create_regions_mesh_ext_mesh_determine_kappamu
+end subroutine create_regions_mesh_ext_mesh_determine_velocity
!
!----
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/OUTPUT_FILES/cubit2specfem3d.py
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/OUTPUT_FILES/cubit2specfem3d.py 2009-10-25 23:06:34 UTC (rev 15876)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/decompose_mesh_SCOTCH/OUTPUT_FILES/cubit2specfem3d.py 2009-10-25 23:54:42 UTC (rev 15877)
@@ -34,11 +34,13 @@
#or
# manually following the convention:
# - each material should have a block defined by name,flag of the material (integer),p velocity
-# (or the full description: name, flag, vp, vs, rho, Q ... if not present these last 3 parameters will be interpolated by module mat_parameter)
+# (or the full description: name, flag, vp, vs, rho, Q ... if not present these last 3 parameters will be
+# interpolated by module mat_parameter)
# - each mesh should have the block definition for the face on the free_surface (topography),
# the name of this block must be 'face_topo' or you can change the default name in mesh.topo defined in profile.
# - each mesh should have the block definition for the faces on the absorbing boundaries,
-# one block for each surface with x=Xmin,x=Xmax,y=Ymin,y=Ymax and z=bottom. The names of the blocks should contain the strings "xmin,xmax,ymin,ymax,bottom"
+# one block for each surface with x=Xmin,x=Xmax,y=Ymin,y=Ymax and z=bottom. The names of
+# the blocks should contain the strings "xmin,xmax,ymin,ymax,bottom"
#
#############################################################################
#RUN
@@ -79,14 +81,20 @@
# .....
# flag 'tomography' file_name #for interpolation with tomography
# .....
-# flag 'interface' file_name flag_for_the_gll_below_the_interface flag_for_the_gll_above_the_interface #for interpolation with interface
+# flag 'interface' file_name flag_for_the_gll_below_the_interface
+# flag_for_the_gll_above_the_interface #for interpolation with interface
#__________________________________________________________________________________________
##absname='absorbing_surface_file' -> this file contains all the face in all the absorbing boundaries
-##absname_local='absorbing_surface_file'+'_xmin' -> this file contains all the face in the absorbing boundary defined by x=Xmin
-##absname_local='absorbing_surface_file'+'_xmax' -> this file contains all the face in the absorbing boundary defined by x=Xmax
-##absname_local='absorbing_surface_file'+'_ymin' -> this file contains all the face in the absorbing boundary defined by y=Ymin
-##absname_local='absorbing_surface_file'+'_ymax' -> this file contains all the face in the absorbing boundary defined by y=Ymax
-##absname_local='absorbing_surface_file'+'_bottom' -> this file contains all the face in the absorbing boundary defined by z=bottom
+##absname_local='absorbing_surface_file'+'_xmin' -> this file contains all the face in the
+# absorbing boundary defined by x=Xmin
+##absname_local='absorbing_surface_file'+'_xmax' -> this file contains all the face in the
+# absorbing boundary defined by x=Xmax
+##absname_local='absorbing_surface_file'+'_ymin' -> this file contains all the face in the
+# absorbing boundary defined by y=Ymin
+##absname_local='absorbing_surface_file'+'_ymax' -> this file contains all the face in the
+# absorbing boundary defined by y=Ymax
+##absname_local='absorbing_surface_file'+'_bottom' -> this file contains all the face in the
+# absorbing boundary defined by z=bottom
# format:
# number of faces
# id_(element containg the face) id_node1_face id_node2_face id_node3_face id_node4_face
@@ -99,7 +107,8 @@
# id_(element containg the face) id_node1_face id_node2_face id_node3_face id_node4_face
#
#__________________________________________________________________________________________
-# it is possible save only one (or more) file singularly: for example if you want only the nodecoord_file call the module mesh.nodescoord_write(full path name)
+# it is possible save only one (or more) file singularly: for example if you want only the nodecoord_file
+# call the module mesh.nodescoord_write(full path name)
#
#############################################################################
@@ -120,7 +129,8 @@
def __repr__(self):
txt='Meshing for frequency up to '+str(self.frequency)+'Hz\n'
for surf,vp in zip(self.list_surf,self.list_vp):
- txt=txt+'surface '+str(surf)+', vp ='+str(vp)+' -> size '+str(self.freq2meshsize(vp)[0])+' -> dt '+str(self.freq2meshsize(vp)[0])+'\n'
+ txt=txt+'surface '+str(surf)+', vp ='+str(vp)+' -> size '+str(self.freq2meshsize(vp)[0])\
+ +' -> dt '+str(self.freq2meshsize(vp)[0])+'\n'
return txt
def freq2meshsize(self,vp):
velocity=vp*.5
@@ -182,8 +192,10 @@
"""Tools for the mesh
#########
dt,edge_dt,freq,edge_freq=seismic_resolution(edges,velocity,bins_d=None,bins_u=None,sidelist=None,ngll=5,np=8)
- Given the velocity of a list of edges, seismic_resolution provides the minimum Dt required for the stability condition (and the corrisponding edge).
- Furthermore, given the number of gll point in the element (ngll) and the number of GLL point for wavelength, it provide the maximum resolved frequency.
+ Given the velocity of a list of edges, seismic_resolution provides the minimum Dt
+ required for the stability condition (and the corrisponding edge).
+ Furthermore, given the number of gll point in the element (ngll) and the number
+ of GLL point for wavelength, it provide the maximum resolved frequency.
#########
length=edge_length(edge)
return the length of a edge
@@ -197,8 +209,10 @@
def seismic_resolution(self,edges,velocity,bins_d=None,bins_u=None,sidelist=None):
"""
dt,edge_dt,freq,edge_freq=seismic_resolution(edges,velocity,bins_d=None,bins_u=None,sidelist=None,ngll=5,np=8)
- Given the velocity of a list of edges, seismic_resolution provides the minimum Dt required for the stability condition (and the corrisponding edge).
- Furthermore, given the number of gll point in the element (ngll) and the number of GLL point for wavelength, it provide the maximum resolved frequency.
+ Given the velocity of a list of edges, seismic_resolution provides the minimum Dt
+ required for the stability condition (and the corrisponding edge).
+ Furthermore, given the number of gll point in the element (ngll) and the number
+ of GLL point for wavelength, it provide the maximum resolved frequency.
"""
ratiostore=1e10
dtstore=1e10
@@ -307,7 +321,8 @@
es=cubit.get_sub_elements("face", face, 1)
edges=edges+list(es)
edges=Set(edges)
- dtstore,edgedtstore,ratiostore,edgeratiostore=self.seismic_resolution(edges,velocity,bins_d,bins_u,sidelist)
+ dtstore,edgedtstore,ratiostore,edgeratiostore=self.seismic_resolution(edges,\
+ velocity,bins_d,bins_u,sidelist)
dt.append(dtstore)
ed_dt.append(edgedtstore)
r.append(ratiostore)
@@ -364,24 +379,28 @@
vs=None
rho=None
q=None
+ ani=None
if nattrib != 0:
flag=int(cubit.get_block_attribute_value(block,0))
if flag > 0 and nattrib >= 2:
- vel=cubit.get_block_attribute_value(block,1)
- if nattrib >= 3:
- vs=cubit.get_block_attribute_value(block,2)
- if nattrib >= 4:
- #density
- rho=cubit.get_block_attribute_value(block,3)
- if nattrib == 5:
- #Q_flag
- q=cubit.get_block_attribute_value(block,4)
- # for q to be valid: it must be an integer flag between 1 and 13
- # (see constants.h for IATTENUATION_SEDIMENT_40, etc. )
- if q < 0 or q > 13:
- print 'error, q flag invalid:', q
- print ' check with constants.h for IATTENUATION flags'
- break
+ vel=cubit.get_block_attribute_value(block,1)
+ if nattrib >= 3:
+ vs=cubit.get_block_attribute_value(block,2)
+ if nattrib >= 4:
+ #density
+ rho=cubit.get_block_attribute_value(block,3)
+ if nattrib >= 5:
+ #Q_flag
+ q=cubit.get_block_attribute_value(block,4)
+ # for q to be valid: it must be an integer flag between 1 and 13
+ # (see constants.h for IATTENUATION_SEDIMENT_40, etc. )
+ if q < 0 or q > 13:
+ print 'error, q flag invalid:', q
+ print ' check with constants.h for IATTENUATION flags'
+ break
+ if nattrib == 6:
+ #anisotropy_flag
+ ani=cubit.get_block_attribute_value(block,5)
elif flag < 0:
vel=name
attrib=cubit.get_block_attribute_value(block,1)
@@ -393,11 +412,11 @@
kind='tomography'
else:
flag=block
- vel,vs,rho,q=(name,0,0,0)
+ vel,vs,rho,q,ani=(name,0,0,0,0)
block_flag.append(int(flag))
block_mat.append(block)
if flag > 0:
- par=tuple([flag,vel,vs,rho,q])
+ par=tuple([flag,vel,vs,rho,q,ani])
elif flag < 0:
if kind=='interface':
par=tuple([flag,kind,name,flag_down,flag_up])
@@ -441,7 +460,8 @@
print bc
print topography
def mat_parameter(self,properties):
- #TODO: attenuation q .... where?
+ #TODO: material property acoustic/elastic/poroelastic ? .... where?
+ print "#material properties:"
print properties
flag=properties[0]
if flag > 0:
@@ -455,15 +475,17 @@
rho=(1.6612*vp-0.472*vp**2+0.0671*vp**3-0.0043*vp**4+0.000106*vp**4)*m2km
txt='%3i %20f %20f %20f %1i %1i\n' % (properties[0],rho,vel,vel/(3**.5),0,0)
elif type(vel) != str:
- #format nummaterials file: #material_id #rho #vp #vs #Q_flag #0
- txt='%3i %20f %20f %20f %2i %1i\n' % (properties[0],properties[3],properties[1],properties[2],properties[4],0)
+ #format nummaterials file: #material_id #rho #vp #vs #Q_flag #anisotropy_flag
+ txt='%3i %20f %20f %20f %2i %2i\n' % (properties[0],properties[3], \
+ properties[1],properties[2],properties[4],properties[5])
else:
txt='%3i %s \n' % (properties[0],properties[1])
elif flag < 0:
if properties[1] == 'tomography':
txt='%3i %s %s\n' % (properties[0],properties[1],properties[2])
elif properties[1] == 'interface':
- txt='%3i %s %s %1i %1i\n' % (properties[0],properties[1],properties[2],properties[3],properties[4])
+ txt='%3i %s %s %1i %1i\n' % (properties[0],properties[1],properties[2],\
+ properties[3],properties[4])
return txt
def nummaterial_write(self,nummaterial_name):
print 'Writing '+nummaterial_name+'.....'
@@ -542,7 +564,8 @@
#print f
nodes=cubit.get_connectivity('Face',f)
nodes_ok=self.normal_check(nodes,normal)
- txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],nodes_ok[1],nodes_ok[2],nodes_ok[3])
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],\
+ nodes_ok[1],nodes_ok[2],nodes_ok[3])
freehex.write(txt)
freehex.close()
print 'Ok'
@@ -606,9 +629,11 @@
nodes=cubit.get_connectivity('Face',f)
if not absflag:
nodes_ok=self.normal_check(nodes,normal)
- txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],nodes_ok[1],nodes_ok[2],nodes_ok[3])
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes_ok[0],\
+ nodes_ok[1],nodes_ok[2],nodes_ok[3])
else:
- txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],nodes[1],nodes[2],nodes[3])
+ txt='%10i %10i %10i %10i %10i\n' % (h,nodes[0],\
+ nodes[1],nodes[2],nodes[3])
abshex_local.write(txt)
abshex_local.close()
print 'Ok'
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/detect_mesh_surfaces.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/detect_mesh_surfaces.f90 2009-10-25 23:06:34 UTC (rev 15876)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/detect_mesh_surfaces.f90 2009-10-25 23:54:42 UTC (rev 15877)
@@ -39,7 +39,9 @@
if (.not. RECVS_CAN_BE_BURIED_EXT_MESH .or. EXTERNAL_MESH_MOVIE_SURFACE .or. EXTERNAL_MESH_CREATE_SHAKEMAP) then
- ! returns surace points/elements
+ ! returns surface points/elements
+ ! in ispec_is_surface_external_mesh / iglob_is_surface_external_mesh and
+ ! number of faces in nfaces_surface_external_mesh
call detect_surface(NPROC,NGLOB_AB,NSPEC_AB,ibool,&
ispec_is_surface_external_mesh, &
iglob_is_surface_external_mesh, &
@@ -49,75 +51,6 @@
nibool_interfaces_ext_mesh, &
my_neighbours_ext_mesh, &
ibool_interfaces_ext_mesh)
-
-!
-! valence_external_mesh(:) = 0
-! ispec_is_surface_external_mesh(:) = .false.
-! iglob_is_surface_external_mesh(:) = .false.
-! do ispec = 1, NSPEC_AB
-! do k = 1, NGLLZ
-! do j = 1, NGLLY
-! do i = 1, NGLLX
-! iglob = ibool(i,j,k,ispec)
-! valence_external_mesh(iglob) = valence_external_mesh(iglob) + 1
-! enddo
-! enddo
-! enddo
-! enddo
-!
-! allocate(buffer_send_scalar_i_ext_mesh(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
-! allocate(buffer_recv_scalar_i_ext_mesh(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
-!
-! ! adds contributions from different partitions to valence_external_mesh
-! call assemble_MPI_scalar_i_ext_mesh(NPROC,NGLOB_AB,valence_external_mesh, &
-! buffer_send_scalar_i_ext_mesh,buffer_recv_scalar_i_ext_mesh, &
-! num_interfaces_ext_mesh,max_nibool_interfaces_ext_mesh, &
-! nibool_interfaces_ext_mesh,ibool_interfaces_ext_mesh,my_neighbours_ext_mesh, &
-! request_send_scalar_ext_mesh,request_recv_scalar_ext_mesh)
-!
-! do ispec = 1, NSPEC_AB
-! do k = 1, NGLLZ
-! do j = 1, NGLLY
-! do i = 1, NGLLX
-! if ( &
-! (k == 1 .or. k == NGLLZ) .and. (j /= 1 .and. j /= NGLLY) .and. (i /= 1 .and. i /= NGLLX) .or. &
-! (j == 1 .or. j == NGLLY) .and. (k /= 1 .and. k /= NGLLZ) .and. (i /= 1 .and. i /= NGLLX) .or. &
-! (i == 1 .or. i == NGLLX) .and. (k /= 1 .and. k /= NGLLZ) .and. (j /= 1 .and. j /= NGLLY) &
-! ) then
-! iglob = ibool(i,j,k,ispec)
-! if (valence_external_mesh(iglob) == 1) then
-! ispec_is_surface_external_mesh(ispec) = .true.
-!
-! if (k == 1 .or. k == NGLLZ) then
-! do jj = 1, NGLLY
-! do ii = 1, NGLLX
-! iglob_is_surface_external_mesh(ibool(ii,jj,k,ispec)) = .true.
-! enddo
-! enddo
-! endif
-! if (j == 1 .or. j == NGLLY) then
-! do kk = 1, NGLLZ
-! do ii = 1, NGLLX
-! iglob_is_surface_external_mesh(ibool(ii,j,kk,ispec)) = .true.
-! enddo
-! enddo
-! endif
-! if (i == 1 .or. i == NGLLX) then
-! do kk = 1, NGLLZ
-! do jj = 1, NGLLY
-! iglob_is_surface_external_mesh(ibool(i,jj,kk,ispec)) = .true.
-! enddo
-! enddo
-! endif
-! endif
-!
-! endif
-! enddo
-! enddo
-! enddo
-!
-! enddo ! nspec
-
endif
! handles movies and shakemaps
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/generate_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/generate_databases.f90 2009-10-25 23:06:34 UTC (rev 15876)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/generate_databases.f90 2009-10-25 23:54:42 UTC (rev 15877)
@@ -306,7 +306,7 @@
integer :: nfaces_surface_external_mesh,nfaces_surface_glob_ext_mesh
integer :: i
- end module
+ end module generate_databases_par
!
!-------------------------------------------------------------------------------------------------
@@ -483,6 +483,13 @@
endif
write(IMAIN,*)
+ if(ANISOTROPY) then
+ write(IMAIN,*) 'incorporating anisotropy'
+ else
+ write(IMAIN,*) 'no anisotropy'
+ endif
+
+ write(IMAIN,*)
if(OCEANS) then
write(IMAIN,*) 'incorporating the oceans using equivalent load'
else
@@ -581,7 +588,7 @@
allocate(materials_ext_mesh(5,nmat_ext_mesh))
allocate(undef_mat_prop(5,nundefMat_ext_mesh))
do imat = 1, nmat_ext_mesh
- ! format: # rho # vp # vs # Q_flag # 0
+ ! format: #(1) rho #(2) vp #(3) vs #(4) Q_flag #(5) anisotropy_flag
read(IIN,*) materials_ext_mesh(1,imat), materials_ext_mesh(2,imat), materials_ext_mesh(3,imat), &
materials_ext_mesh(4,imat), materials_ext_mesh(5,imat)
end do
@@ -757,7 +764,8 @@
ibelm_xmin, ibelm_xmax, ibelm_ymin, ibelm_ymax, ibelm_bottom, ibelm_top, &
nodes_ibelm_xmin,nodes_ibelm_xmax,nodes_ibelm_ymin,nodes_ibelm_ymax, &
nodes_ibelm_bottom,nodes_ibelm_top, &
- SAVE_MESH_FILES,nglob)
+ SAVE_MESH_FILES,nglob, &
+ ANISOTROPY)
call sync_all()
@@ -835,7 +843,8 @@
do i = 1, num_interfaces_ext_mesh
ibool_interfaces_ext_mesh_dummy(:,:) = ibool_interfaces_ext_mesh(1:max_nibool_interfaces_ext_mesh,:)
enddo
-
+ call sync_all()
+
call detect_surface(NPROC,NGLOB_AB,NSPEC_AB,ibool, &
ispec_is_surface_external_mesh, &
iglob_is_surface_external_mesh, &
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/initialize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/initialize_simulation.f90 2009-10-25 23:06:34 UTC (rev 15876)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/initialize_simulation.f90 2009-10-25 23:54:42 UTC (rev 15877)
@@ -90,7 +90,7 @@
if( NGLLX /= NGLLY .and. NGLLY /= NGLLZ ) &
stop 'must have NGLLX = NGLLY = NGLLZ'
endif
-
+
! chris: DT_ext_mesh & NSTE_ext_mesh were in constants.h, I suppressed it, now it is Par_file & read in
! read_parameters_file.f90
! DT = DT_ext_mesh
@@ -102,6 +102,7 @@
read(27) NGLOB_AB
close(27)
+! attenuation arrays size
if( ATTENUATION ) then
!pll
NSPEC_ATTENUATION_AB = NSPEC_AB
@@ -110,6 +111,14 @@
NSPEC_ATTENUATION_AB = 1
endif
+! anisotropy arrays size
+ if( ANISOTROPY ) then
+ NSPEC_ANISO = NSPEC_AB
+ else
+ ! if off, set dummy size
+ NSPEC_ANISO = 1
+ endif
+
! open main output file, only written to by process 0
if(myrank == 0 .and. IMAIN /= ISTANDARD_OUTPUT) &
open(unit=IMAIN,file=trim(OUTPUT_FILES)//'/output_solver.txt',status='unknown')
@@ -154,11 +163,9 @@
endif
-
! check that we have at least one source
if(NSOURCES < 1) call exit_MPI(myrank,'need at least one source')
-
! allocate arrays for storing the databases
allocate(ibool(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
allocate(xix(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
@@ -171,22 +178,49 @@
allocate(gammay(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
allocate(gammaz(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
allocate(jacobian(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+! mesh node locations
allocate(xstore(NGLOB_AB))
allocate(ystore(NGLOB_AB))
allocate(zstore(NGLOB_AB))
+! material properties
allocate(kappastore(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
allocate(mustore(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
! allocate(not_fully_in_bedrock(NSPEC_AB))
! allocate(flag_sediments(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
allocate(rho_vp(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
allocate(rho_vs(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
+ allocate(c11store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
+ allocate(c12store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
+ allocate(c13store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
+ allocate(c14store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
+ allocate(c15store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
+ allocate(c16store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
+ allocate(c22store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
+ allocate(c23store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
+ allocate(c24store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
+ allocate(c25store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
+ allocate(c26store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
+ allocate(c33store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
+ allocate(c34store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
+ allocate(c35store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
+ allocate(c36store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
+ allocate(c44store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
+ allocate(c45store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
+ allocate(c46store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
+ allocate(c55store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
+ allocate(c56store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
+ allocate(c66store(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO))
+
! allocate(idoubling(NSPEC_AB))
+!mass matrix
allocate(rmass(NGLOB_AB))
- allocate(rmass_ocean_load(NGLOB_AB))
+ allocate(rmass_ocean_load(NGLOB_AB))
allocate(updated_dof_ocean_load(NGLOB_AB))
+! displacement,velocity,acceleration
allocate(displ(NDIM,NGLOB_AB))
allocate(veloc(NDIM,NGLOB_AB))
allocate(accel(NDIM,NGLOB_AB))
+! attenuation
allocate(iflag_attenuation_store(NGLLX,NGLLY,NGLLZ,NSPEC_AB))
end subroutine
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/iterate_time.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/iterate_time.f90 2009-10-25 23:06:34 UTC (rev 15876)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/iterate_time.f90 2009-10-25 23:54:42 UTC (rev 15877)
@@ -619,7 +619,7 @@
use specfem_par_movie
implicit none
-
+
! initializes arrays
if (it == 1) then
@@ -629,23 +629,30 @@
do ispec = 1,nfaces_surface_external_mesh
if (USE_HIGHRES_FOR_MOVIES) then
do ipoin = 1, NGLLX*NGLLY
+ ! x,y,z coordinates
store_val_x_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin) = xstore(faces_surface_external_mesh(ipoin,ispec))
store_val_y_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin) = ystore(faces_surface_external_mesh(ipoin,ispec))
store_val_z_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin) = zstore(faces_surface_external_mesh(ipoin,ispec))
enddo
else
- store_val_x_external_mesh(NGNOD2D*(ispec-1)+1) = xstore(faces_surface_external_mesh(1,ispec))
- store_val_x_external_mesh(NGNOD2D*(ispec-1)+2) = xstore(faces_surface_external_mesh(2,ispec))
- store_val_x_external_mesh(NGNOD2D*(ispec-1)+3) = xstore(faces_surface_external_mesh(3,ispec))
- store_val_x_external_mesh(NGNOD2D*(ispec-1)+4) = xstore(faces_surface_external_mesh(4,ispec))
- store_val_y_external_mesh(NGNOD2D*(ispec-1)+1) = ystore(faces_surface_external_mesh(1,ispec))
- store_val_y_external_mesh(NGNOD2D*(ispec-1)+2) = ystore(faces_surface_external_mesh(2,ispec))
- store_val_y_external_mesh(NGNOD2D*(ispec-1)+3) = ystore(faces_surface_external_mesh(3,ispec))
- store_val_y_external_mesh(NGNOD2D*(ispec-1)+4) = ystore(faces_surface_external_mesh(4,ispec))
- store_val_z_external_mesh(NGNOD2D*(ispec-1)+1) = zstore(faces_surface_external_mesh(1,ispec))
- store_val_z_external_mesh(NGNOD2D*(ispec-1)+2) = zstore(faces_surface_external_mesh(2,ispec))
- store_val_z_external_mesh(NGNOD2D*(ispec-1)+3) = zstore(faces_surface_external_mesh(3,ispec))
- store_val_z_external_mesh(NGNOD2D*(ispec-1)+4) = zstore(faces_surface_external_mesh(4,ispec))
+ do ipoin = 1, 4
+ ! x,y,z coordinates
+ store_val_x_external_mesh(NGNOD2D*(ispec-1)+ipoin) = xstore(faces_surface_external_mesh(ipoin,ispec))
+ store_val_y_external_mesh(NGNOD2D*(ispec-1)+ipoin) = ystore(faces_surface_external_mesh(ipoin,ispec))
+ store_val_z_external_mesh(NGNOD2D*(ispec-1)+ipoin) = zstore(faces_surface_external_mesh(ipoin,ispec))
+ enddo
+! store_val_x_external_mesh(NGNOD2D*(ispec-1)+1) = xstore(faces_surface_external_mesh(1,ispec))
+! store_val_x_external_mesh(NGNOD2D*(ispec-1)+2) = xstore(faces_surface_external_mesh(2,ispec))
+! store_val_x_external_mesh(NGNOD2D*(ispec-1)+3) = xstore(faces_surface_external_mesh(3,ispec))
+! store_val_x_external_mesh(NGNOD2D*(ispec-1)+4) = xstore(faces_surface_external_mesh(4,ispec))
+! store_val_y_external_mesh(NGNOD2D*(ispec-1)+1) = ystore(faces_surface_external_mesh(1,ispec))
+! store_val_y_external_mesh(NGNOD2D*(ispec-1)+2) = ystore(faces_surface_external_mesh(2,ispec))
+! store_val_y_external_mesh(NGNOD2D*(ispec-1)+3) = ystore(faces_surface_external_mesh(3,ispec))
+! store_val_y_external_mesh(NGNOD2D*(ispec-1)+4) = ystore(faces_surface_external_mesh(4,ispec))
+! store_val_z_external_mesh(NGNOD2D*(ispec-1)+1) = zstore(faces_surface_external_mesh(1,ispec))
+! store_val_z_external_mesh(NGNOD2D*(ispec-1)+2) = zstore(faces_surface_external_mesh(2,ispec))
+! store_val_z_external_mesh(NGNOD2D*(ispec-1)+3) = zstore(faces_surface_external_mesh(3,ispec))
+! store_val_z_external_mesh(NGNOD2D*(ispec-1)+4) = zstore(faces_surface_external_mesh(4,ispec))
endif
enddo
endif
@@ -654,16 +661,19 @@
do ispec = 1,nfaces_surface_external_mesh
if (USE_HIGHRES_FOR_MOVIES) then
do ipoin = 1, NGLLX*NGLLY
+ ! norm of displacement
store_val_ux_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin) = &
max(store_val_ux_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin), &
sqrt(displ(1,faces_surface_external_mesh(ipoin,ispec))**2 + &
displ(2,faces_surface_external_mesh(ipoin,ispec))**2 + &
displ(3,faces_surface_external_mesh(ipoin,ispec))**2))
+ ! norm of velocity
store_val_uy_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin) = &
max(store_val_uy_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin), &
sqrt(veloc(1,faces_surface_external_mesh(ipoin,ispec))**2 + &
veloc(2,faces_surface_external_mesh(ipoin,ispec))**2 + &
veloc(3,faces_surface_external_mesh(ipoin,ispec))**2))
+ ! norm of acceleration
store_val_uz_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin) = &
max(store_val_uz_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin), &
sqrt(accel(1,faces_surface_external_mesh(ipoin,ispec))**2 + &
@@ -672,70 +682,91 @@
enddo
else
- store_val_ux_external_mesh(NGNOD2D*(ispec-1)+1) = &
- max(store_val_ux_external_mesh(NGNOD2D*(ispec-1)+1), &
- sqrt(displ(1,faces_surface_external_mesh(1,ispec))**2 + &
- displ(2,faces_surface_external_mesh(1,ispec))**2 + &
- displ(3,faces_surface_external_mesh(1,ispec))**2))
- store_val_ux_external_mesh(NGNOD2D*(ispec-1)+2) = &
- max(store_val_ux_external_mesh(NGNOD2D*(ispec-1)+2), &
- sqrt(displ(1,faces_surface_external_mesh(2,ispec))**2 + &
- displ(2,faces_surface_external_mesh(2,ispec))**2 + &
- displ(3,faces_surface_external_mesh(2,ispec))**2))
- store_val_ux_external_mesh(NGNOD2D*(ispec-1)+3) = &
- max(store_val_ux_external_mesh(NGNOD2D*(ispec-1)+3), &
- sqrt(displ(1,faces_surface_external_mesh(3,ispec))**2 + &
- displ(2,faces_surface_external_mesh(3,ispec))**2 + &
- displ(3,faces_surface_external_mesh(3,ispec))**2))
- store_val_ux_external_mesh(NGNOD2D*(ispec-1)+4) = &
- max(store_val_ux_external_mesh(NGNOD2D*(ispec-1)+4), &
- sqrt(displ(1,faces_surface_external_mesh(4,ispec))**2 + &
- displ(2,faces_surface_external_mesh(4,ispec))**2 + &
- displ(3,faces_surface_external_mesh(4,ispec))**2))
- store_val_uy_external_mesh(NGNOD2D*(ispec-1)+1) = &
- max(store_val_uy_external_mesh(NGNOD2D*(ispec-1)+1), &
- sqrt(veloc(1,faces_surface_external_mesh(1,ispec))**2 + &
- veloc(2,faces_surface_external_mesh(1,ispec))**2 + &
- veloc(3,faces_surface_external_mesh(1,ispec))**2))
- store_val_uy_external_mesh(NGNOD2D*(ispec-1)+2) = &
- max(store_val_uy_external_mesh(NGNOD2D*(ispec-1)+2), &
- sqrt(veloc(1,faces_surface_external_mesh(2,ispec))**2 + &
- veloc(2,faces_surface_external_mesh(2,ispec))**2 + &
- veloc(3,faces_surface_external_mesh(2,ispec))**2))
- store_val_uy_external_mesh(NGNOD2D*(ispec-1)+3) = &
- max(store_val_uy_external_mesh(NGNOD2D*(ispec-1)+3), &
- sqrt(veloc(1,faces_surface_external_mesh(3,ispec))**2 + &
- veloc(2,faces_surface_external_mesh(3,ispec))**2 + &
- veloc(3,faces_surface_external_mesh(3,ispec))**2))
- store_val_uy_external_mesh(NGNOD2D*(ispec-1)+4) = &
- max(store_val_uy_external_mesh(NGNOD2D*(ispec-1)+4), &
- sqrt(veloc(1,faces_surface_external_mesh(4,ispec))**2 + &
- veloc(2,faces_surface_external_mesh(4,ispec))**2 + &
- veloc(3,faces_surface_external_mesh(4,ispec))**2))
- store_val_uz_external_mesh(NGNOD2D*(ispec-1)+1) = &
- max(store_val_uz_external_mesh(NGNOD2D*(ispec-1)+1), &
- sqrt(accel(1,faces_surface_external_mesh(1,ispec))**2 + &
- accel(2,faces_surface_external_mesh(1,ispec))**2 + &
- accel(3,faces_surface_external_mesh(1,ispec))**2))
- store_val_uz_external_mesh(NGNOD2D*(ispec-1)+2) = &
- max(store_val_uz_external_mesh(NGNOD2D*(ispec-1)+2), &
- sqrt(accel(1,faces_surface_external_mesh(2,ispec))**2 + &
- accel(2,faces_surface_external_mesh(2,ispec))**2 + &
- accel(3,faces_surface_external_mesh(2,ispec))**2))
- store_val_uz_external_mesh(NGNOD2D*(ispec-1)+3) = &
- max(store_val_uz_external_mesh(NGNOD2D*(ispec-1)+3), &
- sqrt(accel(1,faces_surface_external_mesh(3,ispec))**2 + &
- accel(2,faces_surface_external_mesh(3,ispec))**2 + &
- accel(3,faces_surface_external_mesh(3,ispec))**2))
- store_val_uz_external_mesh(NGNOD2D*(ispec-1)+4) = &
- max(store_val_uz_external_mesh(NGNOD2D*(ispec-1)+4), &
- sqrt(accel(1,faces_surface_external_mesh(4,ispec))**2 + &
- accel(2,faces_surface_external_mesh(4,ispec))**2 + &
- accel(3,faces_surface_external_mesh(4,ispec))**2))
+ do ipoin = 1, 4
+ ! norm of displacement
+ store_val_ux_external_mesh(NGNOD2D*(ispec-1)+ipoin) = &
+ max(store_val_ux_external_mesh(NGNOD2D*(ispec-1)+ipoin), &
+ sqrt(displ(1,faces_surface_external_mesh(ipoin,ispec))**2 + &
+ displ(2,faces_surface_external_mesh(ipoin,ispec))**2 + &
+ displ(3,faces_surface_external_mesh(ipoin,ispec))**2))
+ ! norm of velocity
+ store_val_uy_external_mesh(NGNOD2D*(ispec-1)+ipoin) = &
+ max(store_val_uy_external_mesh(NGNOD2D*(ispec-1)+ipoin), &
+ sqrt(veloc(1,faces_surface_external_mesh(ipoin,ispec))**2 + &
+ veloc(2,faces_surface_external_mesh(ipoin,ispec))**2 + &
+ veloc(3,faces_surface_external_mesh(ipoin,ispec))**2))
+ ! norm of acceleration
+ store_val_uz_external_mesh(NGNOD2D*(ispec-1)+ipoin) = &
+ max(store_val_uz_external_mesh(NGNOD2D*(ispec-1)+ipoin), &
+ sqrt(accel(1,faces_surface_external_mesh(ipoin,ispec))**2 + &
+ accel(2,faces_surface_external_mesh(ipoin,ispec))**2 + &
+ accel(3,faces_surface_external_mesh(ipoin,ispec))**2))
+ enddo
+
+! store_val_ux_external_mesh(NGNOD2D*(ispec-1)+1) = &
+! max(store_val_ux_external_mesh(NGNOD2D*(ispec-1)+1), &
+! sqrt(displ(1,faces_surface_external_mesh(1,ispec))**2 + &
+! displ(2,faces_surface_external_mesh(1,ispec))**2 + &
+! displ(3,faces_surface_external_mesh(1,ispec))**2))
+! store_val_ux_external_mesh(NGNOD2D*(ispec-1)+2) = &
+! max(store_val_ux_external_mesh(NGNOD2D*(ispec-1)+2), &
+! sqrt(displ(1,faces_surface_external_mesh(2,ispec))**2 + &
+! displ(2,faces_surface_external_mesh(2,ispec))**2 + &
+! displ(3,faces_surface_external_mesh(2,ispec))**2))
+! store_val_ux_external_mesh(NGNOD2D*(ispec-1)+3) = &
+! max(store_val_ux_external_mesh(NGNOD2D*(ispec-1)+3), &
+! sqrt(displ(1,faces_surface_external_mesh(3,ispec))**2 + &
+! displ(2,faces_surface_external_mesh(3,ispec))**2 + &
+! displ(3,faces_surface_external_mesh(3,ispec))**2))
+! store_val_ux_external_mesh(NGNOD2D*(ispec-1)+4) = &
+! max(store_val_ux_external_mesh(NGNOD2D*(ispec-1)+4), &
+! sqrt(displ(1,faces_surface_external_mesh(4,ispec))**2 + &
+! displ(2,faces_surface_external_mesh(4,ispec))**2 + &
+! displ(3,faces_surface_external_mesh(4,ispec))**2))
+! store_val_uy_external_mesh(NGNOD2D*(ispec-1)+1) = &
+! max(store_val_uy_external_mesh(NGNOD2D*(ispec-1)+1), &
+! sqrt(veloc(1,faces_surface_external_mesh(1,ispec))**2 + &
+! veloc(2,faces_surface_external_mesh(1,ispec))**2 + &
+! veloc(3,faces_surface_external_mesh(1,ispec))**2))
+! store_val_uy_external_mesh(NGNOD2D*(ispec-1)+2) = &
+! max(store_val_uy_external_mesh(NGNOD2D*(ispec-1)+2), &
+! sqrt(veloc(1,faces_surface_external_mesh(2,ispec))**2 + &
+! veloc(2,faces_surface_external_mesh(2,ispec))**2 + &
+! veloc(3,faces_surface_external_mesh(2,ispec))**2))
+! store_val_uy_external_mesh(NGNOD2D*(ispec-1)+3) = &
+! max(store_val_uy_external_mesh(NGNOD2D*(ispec-1)+3), &
+! sqrt(veloc(1,faces_surface_external_mesh(3,ispec))**2 + &
+! veloc(2,faces_surface_external_mesh(3,ispec))**2 + &
+! veloc(3,faces_surface_external_mesh(3,ispec))**2))
+! store_val_uy_external_mesh(NGNOD2D*(ispec-1)+4) = &
+! max(store_val_uy_external_mesh(NGNOD2D*(ispec-1)+4), &
+! sqrt(veloc(1,faces_surface_external_mesh(4,ispec))**2 + &
+! veloc(2,faces_surface_external_mesh(4,ispec))**2 + &
+! veloc(3,faces_surface_external_mesh(4,ispec))**2))
+! store_val_uz_external_mesh(NGNOD2D*(ispec-1)+1) = &
+! max(store_val_uz_external_mesh(NGNOD2D*(ispec-1)+1), &
+! sqrt(accel(1,faces_surface_external_mesh(1,ispec))**2 + &
+! accel(2,faces_surface_external_mesh(1,ispec))**2 + &
+! accel(3,faces_surface_external_mesh(1,ispec))**2))
+! store_val_uz_external_mesh(NGNOD2D*(ispec-1)+2) = &
+! max(store_val_uz_external_mesh(NGNOD2D*(ispec-1)+2), &
+! sqrt(accel(1,faces_surface_external_mesh(2,ispec))**2 + &
+! accel(2,faces_surface_external_mesh(2,ispec))**2 + &
+! accel(3,faces_surface_external_mesh(2,ispec))**2))
+! store_val_uz_external_mesh(NGNOD2D*(ispec-1)+3) = &
+! max(store_val_uz_external_mesh(NGNOD2D*(ispec-1)+3), &
+! sqrt(accel(1,faces_surface_external_mesh(3,ispec))**2 + &
+! accel(2,faces_surface_external_mesh(3,ispec))**2 + &
+! accel(3,faces_surface_external_mesh(3,ispec))**2))
+! store_val_uz_external_mesh(NGNOD2D*(ispec-1)+4) = &
+! max(store_val_uz_external_mesh(NGNOD2D*(ispec-1)+4), &
+! sqrt(accel(1,faces_surface_external_mesh(4,ispec))**2 + &
+! accel(2,faces_surface_external_mesh(4,ispec))**2 + &
+! accel(3,faces_surface_external_mesh(4,ispec))**2))
endif
enddo
-! finalizes shakemap
+! finalizes shakemap: master process collects all info
if (it == NSTEP) then
if (USE_HIGHRES_FOR_MOVIES) then
call gatherv_all_cr(store_val_x_external_mesh,nfaces_surface_external_mesh*NGLLX*NGLLY,&
@@ -780,12 +811,12 @@
! creates shakemap file
if(myrank == 0) then
open(unit=IOUT,file=trim(OUTPUT_FILES)//'/shakingdata',status='unknown',form='unformatted')
- write(IOUT) store_val_x_all_external_mesh
- write(IOUT) store_val_y_all_external_mesh
- write(IOUT) store_val_z_all_external_mesh
- write(IOUT) store_val_ux_all_external_mesh
- write(IOUT) store_val_uy_all_external_mesh
- write(IOUT) store_val_uz_all_external_mesh
+ write(IOUT) store_val_x_all_external_mesh ! x coordinates
+ write(IOUT) store_val_y_all_external_mesh ! y coordinates
+ write(IOUT) store_val_z_all_external_mesh ! z coordinates
+ write(IOUT) store_val_ux_all_external_mesh ! norm of displacement vector
+ write(IOUT) store_val_uy_all_external_mesh ! norm of velocity vector
+ write(IOUT) store_val_uz_all_external_mesh ! norm of acceleration vector
close(IOUT)
endif
endif
@@ -809,41 +840,55 @@
do ispec = 1,nfaces_surface_external_mesh
if (USE_HIGHRES_FOR_MOVIES) then
do ipoin = 1, NGLLX*NGLLY
+ ! x,y,z coordinates
store_val_x_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin) = xstore(faces_surface_external_mesh(ipoin,ispec))
store_val_y_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin) = ystore(faces_surface_external_mesh(ipoin,ispec))
store_val_z_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin) = zstore(faces_surface_external_mesh(ipoin,ispec))
+ ! velocity x,y,z-components
store_val_ux_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin) = veloc(1,faces_surface_external_mesh(ipoin,ispec))
store_val_uy_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin) = veloc(2,faces_surface_external_mesh(ipoin,ispec))
store_val_uz_external_mesh(NGLLX*NGLLY*(ispec-1)+ipoin) = veloc(3,faces_surface_external_mesh(ipoin,ispec))
enddo
else
- store_val_x_external_mesh(NGNOD2D*(ispec-1)+1) = xstore(faces_surface_external_mesh(1,ispec))
- store_val_x_external_mesh(NGNOD2D*(ispec-1)+2) = xstore(faces_surface_external_mesh(2,ispec))
- store_val_x_external_mesh(NGNOD2D*(ispec-1)+3) = xstore(faces_surface_external_mesh(3,ispec))
- store_val_x_external_mesh(NGNOD2D*(ispec-1)+4) = xstore(faces_surface_external_mesh(4,ispec))
- store_val_y_external_mesh(NGNOD2D*(ispec-1)+1) = ystore(faces_surface_external_mesh(1,ispec))
- store_val_y_external_mesh(NGNOD2D*(ispec-1)+2) = ystore(faces_surface_external_mesh(2,ispec))
- store_val_y_external_mesh(NGNOD2D*(ispec-1)+3) = ystore(faces_surface_external_mesh(3,ispec))
- store_val_y_external_mesh(NGNOD2D*(ispec-1)+4) = ystore(faces_surface_external_mesh(4,ispec))
- store_val_z_external_mesh(NGNOD2D*(ispec-1)+1) = zstore(faces_surface_external_mesh(1,ispec))
- store_val_z_external_mesh(NGNOD2D*(ispec-1)+2) = zstore(faces_surface_external_mesh(2,ispec))
- store_val_z_external_mesh(NGNOD2D*(ispec-1)+3) = zstore(faces_surface_external_mesh(3,ispec))
- store_val_z_external_mesh(NGNOD2D*(ispec-1)+4) = zstore(faces_surface_external_mesh(4,ispec))
- store_val_ux_external_mesh(NGNOD2D*(ispec-1)+1) = veloc(1,faces_surface_external_mesh(1,ispec))
- store_val_ux_external_mesh(NGNOD2D*(ispec-1)+2) = veloc(1,faces_surface_external_mesh(2,ispec))
- store_val_ux_external_mesh(NGNOD2D*(ispec-1)+3) = veloc(1,faces_surface_external_mesh(3,ispec))
- store_val_ux_external_mesh(NGNOD2D*(ispec-1)+4) = veloc(1,faces_surface_external_mesh(4,ispec))
- store_val_uy_external_mesh(NGNOD2D*(ispec-1)+1) = veloc(2,faces_surface_external_mesh(1,ispec))
- store_val_uy_external_mesh(NGNOD2D*(ispec-1)+2) = veloc(2,faces_surface_external_mesh(2,ispec))
- store_val_uy_external_mesh(NGNOD2D*(ispec-1)+3) = veloc(2,faces_surface_external_mesh(3,ispec))
- store_val_uy_external_mesh(NGNOD2D*(ispec-1)+4) = veloc(2,faces_surface_external_mesh(4,ispec))
- store_val_uz_external_mesh(NGNOD2D*(ispec-1)+1) = veloc(3,faces_surface_external_mesh(1,ispec))
- store_val_uz_external_mesh(NGNOD2D*(ispec-1)+2) = veloc(3,faces_surface_external_mesh(2,ispec))
- store_val_uz_external_mesh(NGNOD2D*(ispec-1)+3) = veloc(3,faces_surface_external_mesh(3,ispec))
- store_val_uz_external_mesh(NGNOD2D*(ispec-1)+4) = veloc(3,faces_surface_external_mesh(4,ispec))
+ do ipoin = 1, 4
+ ! x,y,z coordinates
+ store_val_x_external_mesh(NGNOD2D*(ispec-1)+ipoin) = xstore(faces_surface_external_mesh(ipoin,ispec))
+ store_val_y_external_mesh(NGNOD2D*(ispec-1)+ipoin) = ystore(faces_surface_external_mesh(ipoin,ispec))
+ store_val_z_external_mesh(NGNOD2D*(ispec-1)+ipoin) = zstore(faces_surface_external_mesh(ipoin,ispec))
+ ! velocity x,y,z-components
+ store_val_ux_external_mesh(NGNOD2D*(ispec-1)+ipoin) = veloc(1,faces_surface_external_mesh(ipoin,ispec))
+ store_val_uy_external_mesh(NGNOD2D*(ispec-1)+ipoin) = veloc(2,faces_surface_external_mesh(ipoin,ispec))
+ store_val_uz_external_mesh(NGNOD2D*(ispec-1)+ipoin) = veloc(3,faces_surface_external_mesh(ipoin,ispec))
+
+ enddo
+! store_val_x_external_mesh(NGNOD2D*(ispec-1)+1) = xstore(faces_surface_external_mesh(1,ispec))
+! store_val_x_external_mesh(NGNOD2D*(ispec-1)+2) = xstore(faces_surface_external_mesh(2,ispec))
+! store_val_x_external_mesh(NGNOD2D*(ispec-1)+3) = xstore(faces_surface_external_mesh(3,ispec))
+! store_val_x_external_mesh(NGNOD2D*(ispec-1)+4) = xstore(faces_surface_external_mesh(4,ispec))
+! store_val_y_external_mesh(NGNOD2D*(ispec-1)+1) = ystore(faces_surface_external_mesh(1,ispec))
+! store_val_y_external_mesh(NGNOD2D*(ispec-1)+2) = ystore(faces_surface_external_mesh(2,ispec))
+! store_val_y_external_mesh(NGNOD2D*(ispec-1)+3) = ystore(faces_surface_external_mesh(3,ispec))
+! store_val_y_external_mesh(NGNOD2D*(ispec-1)+4) = ystore(faces_surface_external_mesh(4,ispec))
+! store_val_z_external_mesh(NGNOD2D*(ispec-1)+1) = zstore(faces_surface_external_mesh(1,ispec))
+! store_val_z_external_mesh(NGNOD2D*(ispec-1)+2) = zstore(faces_surface_external_mesh(2,ispec))
+! store_val_z_external_mesh(NGNOD2D*(ispec-1)+3) = zstore(faces_surface_external_mesh(3,ispec))
+! store_val_z_external_mesh(NGNOD2D*(ispec-1)+4) = zstore(faces_surface_external_mesh(4,ispec))
+! store_val_ux_external_mesh(NGNOD2D*(ispec-1)+1) = veloc(1,faces_surface_external_mesh(1,ispec))
+! store_val_ux_external_mesh(NGNOD2D*(ispec-1)+2) = veloc(1,faces_surface_external_mesh(2,ispec))
+! store_val_ux_external_mesh(NGNOD2D*(ispec-1)+3) = veloc(1,faces_surface_external_mesh(3,ispec))
+! store_val_ux_external_mesh(NGNOD2D*(ispec-1)+4) = veloc(1,faces_surface_external_mesh(4,ispec))
+! store_val_uy_external_mesh(NGNOD2D*(ispec-1)+1) = veloc(2,faces_surface_external_mesh(1,ispec))
+! store_val_uy_external_mesh(NGNOD2D*(ispec-1)+2) = veloc(2,faces_surface_external_mesh(2,ispec))
+! store_val_uy_external_mesh(NGNOD2D*(ispec-1)+3) = veloc(2,faces_surface_external_mesh(3,ispec))
+! store_val_uy_external_mesh(NGNOD2D*(ispec-1)+4) = veloc(2,faces_surface_external_mesh(4,ispec))
+! store_val_uz_external_mesh(NGNOD2D*(ispec-1)+1) = veloc(3,faces_surface_external_mesh(1,ispec))
+! store_val_uz_external_mesh(NGNOD2D*(ispec-1)+2) = veloc(3,faces_surface_external_mesh(2,ispec))
+! store_val_uz_external_mesh(NGNOD2D*(ispec-1)+3) = veloc(3,faces_surface_external_mesh(3,ispec))
+! store_val_uz_external_mesh(NGNOD2D*(ispec-1)+4) = veloc(3,faces_surface_external_mesh(4,ispec))
endif
enddo
+! master process collects all info
if (USE_HIGHRES_FOR_MOVIES) then
call gatherv_all_cr(store_val_x_external_mesh,nfaces_surface_external_mesh*NGLLX*NGLLY,&
store_val_x_all_external_mesh,nfaces_perproc_surface_ext_mesh*NGLLX*NGLLY,faces_surface_offset_ext_mesh,&
@@ -884,15 +929,16 @@
nfaces_surface_glob_ext_mesh*NGNOD2D,NPROC)
endif
+! file output
if(myrank == 0) then
write(outputname,"('/moviedata',i6.6)") it
open(unit=IOUT,file=trim(OUTPUT_FILES)//outputname,status='unknown',form='unformatted')
- write(IOUT) store_val_x_all_external_mesh
- write(IOUT) store_val_y_all_external_mesh
- write(IOUT) store_val_z_all_external_mesh
- write(IOUT) store_val_ux_all_external_mesh
- write(IOUT) store_val_uy_all_external_mesh
- write(IOUT) store_val_uz_all_external_mesh
+ write(IOUT) store_val_x_all_external_mesh ! x coordinate
+ write(IOUT) store_val_y_all_external_mesh ! y coordinate
+ write(IOUT) store_val_z_all_external_mesh ! z coordinate
+ write(IOUT) store_val_ux_all_external_mesh ! velocity x-component
+ write(IOUT) store_val_uy_all_external_mesh ! velocity y-component
+ write(IOUT) store_val_uz_all_external_mesh ! velocity z-component
close(IOUT)
endif
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/prepare_timerun.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/prepare_timerun.f90 2009-10-25 23:06:34 UTC (rev 15876)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/prepare_timerun.f90 2009-10-25 23:54:42 UTC (rev 15877)
@@ -56,6 +56,13 @@
endif
write(IMAIN,*)
+ if(ANISOTROPY) then
+ write(IMAIN,*) 'incorporating anisotropy'
+ else
+ write(IMAIN,*) 'no anisotropy'
+ endif
+
+ write(IMAIN,*)
if(OCEANS) then
write(IMAIN,*) 'incorporating the oceans using equivalent load'
else
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/read_mesh_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/read_mesh_databases.f90 2009-10-25 23:06:34 UTC (rev 15876)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/read_mesh_databases.f90 2009-10-25 23:54:42 UTC (rev 15877)
@@ -179,6 +179,33 @@
read(27) nibool_interfaces_ext_mesh
read(27) ibool_interfaces_ext_mesh
+ if( ANISOTROPY ) then
+ read(27) c11store
+ read(27) c12store
+ read(27) c13store
+ read(27) c14store
+ read(27) c15store
+ read(27) c16store
+ read(27) c22store
+ read(27) c23store
+ read(27) c24store
+ read(27) c25store
+ read(27) c26store
+ read(27) c33store
+ read(27) c34store
+ read(27) c35store
+ read(27) c36store
+ read(27) c44store
+ read(27) c45store
+ read(27) c46store
+ read(27) c55store
+ read(27) c56store
+ read(27) c66store
+ endif
+
+ close(27)
+
+! MPI communications
allocate(buffer_send_vector_ext_mesh(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
allocate(buffer_recv_vector_ext_mesh(NDIM,max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
allocate(buffer_send_scalar_ext_mesh(max_nibool_interfaces_ext_mesh,num_interfaces_ext_mesh))
@@ -187,7 +214,6 @@
allocate(request_recv_vector_ext_mesh(num_interfaces_ext_mesh))
allocate(request_send_scalar_ext_mesh(num_interfaces_ext_mesh))
allocate(request_recv_scalar_ext_mesh(num_interfaces_ext_mesh))
- close(27)
! locate inner and outer elements
allocate(ispec_is_inner_ext_mesh(NSPEC_AB))
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/save_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/save_arrays_solver.f90 2009-10-25 23:06:34 UTC (rev 15876)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/save_arrays_solver.f90 2009-10-25 23:54:42 UTC (rev 15877)
@@ -27,16 +27,22 @@
! for external mesh
subroutine save_arrays_solver_ext_mesh(nspec,nglob, &
- xixstore,xiystore,xizstore,etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore, &
- jacobianstore, rho_vp,rho_vs,iflag_attenuation_store, &
- kappastore,mustore,rmass,ibool,xstore_dummy,ystore_dummy,zstore_dummy, &
- NSPEC2D_TOP,ibelm_top,normal_top,jacobian2D_top, &
- absorbing_boundary_normal,absorbing_boundary_jacobian2D, &
- absorbing_boundary_ijk,absorbing_boundary_ispec, &
- num_absorbing_boundary_faces, &
- num_interfaces_ext_mesh,my_neighbours_ext_mesh,nibool_interfaces_ext_mesh, &
- max_interface_size_ext_mesh,ibool_interfaces_ext_mesh, &
- prname,SAVE_MESH_FILES)
+ xixstore,xiystore,xizstore,etaxstore,etaystore,etazstore, &
+ gammaxstore,gammaystore,gammazstore, &
+ jacobianstore, rho_vp,rho_vs,iflag_attenuation_store, &
+ kappastore,mustore,rmass,ibool,xstore_dummy,ystore_dummy,zstore_dummy, &
+ NSPEC2D_TOP,ibelm_top,normal_top,jacobian2D_top, &
+ absorbing_boundary_normal,absorbing_boundary_jacobian2D, &
+ absorbing_boundary_ijk,absorbing_boundary_ispec, &
+ num_absorbing_boundary_faces, &
+ num_interfaces_ext_mesh,my_neighbours_ext_mesh,nibool_interfaces_ext_mesh, &
+ max_interface_size_ext_mesh,ibool_interfaces_ext_mesh, &
+ prname,SAVE_MESH_FILES, &
+ ANISOTROPY,NSPEC_ANISO, &
+ c11store,c12store,c13store,c14store,c15store,c16store, &
+ c22store,c23store,c24store,c25store,c26store,c33store, &
+ c34store,c35store,c36store,c44store,c45store,c46store, &
+ c55store,c56store,c66store)
implicit none
@@ -107,6 +113,15 @@
! file name
character(len=150) prname
logical :: SAVE_MESH_FILES
+
+! anisotropy
+ logical :: ANISOTROPY
+ integer :: NSPEC_ANISO
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_ANISO) :: &
+ c11store,c12store,c13store,c14store,c15store,c16store, &
+ c22store,c23store,c24store,c25store,c26store,c33store, &
+ c34store,c35store,c36store,c44store,c45store,c46store, &
+ c55store,c56store,c66store
! local parameters
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: v_tmp
@@ -219,9 +234,35 @@
enddo
write(IOUT) ibool_interfaces_ext_mesh_dummy
+ deallocate(ibool_interfaces_ext_mesh_dummy,stat=ier); if( ier /= 0 ) stop 'error deallocating array'
+
+! anisotropy
+ if( ANISOTROPY ) then
+ write(IOUT) c11store
+ write(IOUT) c12store
+ write(IOUT) c13store
+ write(IOUT) c14store
+ write(IOUT) c15store
+ write(IOUT) c16store
+ write(IOUT) c22store
+ write(IOUT) c23store
+ write(IOUT) c24store
+ write(IOUT) c25store
+ write(IOUT) c26store
+ write(IOUT) c33store
+ write(IOUT) c34store
+ write(IOUT) c35store
+ write(IOUT) c36store
+ write(IOUT) c44store
+ write(IOUT) c45store
+ write(IOUT) c46store
+ write(IOUT) c55store
+ write(IOUT) c56store
+ write(IOUT) c66store
+ endif
+
close(IOUT)
- deallocate(ibool_interfaces_ext_mesh_dummy,stat=ier); if( ier /= 0 ) stop 'error deallocating array'
! mesh arrays used for example in combine_vol_data.f90
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/save_header_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/save_header_file.f90 2009-10-25 23:06:34 UTC (rev 15876)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/save_header_file.f90 2009-10-25 23:54:42 UTC (rev 15877)
@@ -139,12 +139,14 @@
! anisotropy
if(ANISOTROPY) then
- stop 'ANISOTROPY not supported yet in the CUBIT + SCOTCH version because of arrays of constant size defined'
- ! write(IOUT,*) 'integer, parameter :: NSPEC_ANISO = ',NSPEC_AB
- write(IOUT,*) 'logical, parameter :: ANISOTROPY_VAL = .true.'
+ !stop 'ANISOTROPY not supported yet in the CUBIT + SCOTCH version because of arrays of constant size defined'
+ !write(IOUT,*) 'integer, parameter :: NSPEC_ANISO = ',NSPEC_AB
+ !write(IOUT,*) 'logical, parameter :: ANISOTROPY_VAL = .true.'
+ write(IOUT,*) '! with anisotropy'
else
- write(IOUT,*) 'integer, parameter :: NSPEC_ANISO = ', 1
- write(IOUT,*) 'logical, parameter :: ANISOTROPY_VAL = .false.'
+ !write(IOUT,*) 'integer, parameter :: NSPEC_ANISO = ', 1
+ !write(IOUT,*) 'logical, parameter :: ANISOTROPY_VAL = .false.'
+ write(IOUT,*) '! no anisotropy'
endif
write(IOUT,*)
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/setup_movie_meshes.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/setup_movie_meshes.f90 2009-10-25 23:06:34 UTC (rev 15876)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/setup_movie_meshes.f90 2009-10-25 23:54:42 UTC (rev 15877)
@@ -35,34 +35,6 @@
implicit none
! initializes mesh arrays for movies and shakemaps
-! nfaces_surface_external_mesh = 0
-! do ispec = 1, NSPEC_AB
-! iglob = ibool(2,2,1,ispec)
-! if (iglob_is_surface_external_mesh(iglob)) then
-! nfaces_surface_external_mesh = nfaces_surface_external_mesh + 1
-! endif
-! iglob = ibool(2,2,NGLLZ,ispec)
-! if (iglob_is_surface_external_mesh(iglob)) then
-! nfaces_surface_external_mesh = nfaces_surface_external_mesh + 1
-! endif
-! iglob = ibool(2,1,2,ispec)
-! if (iglob_is_surface_external_mesh(iglob)) then
-! nfaces_surface_external_mesh = nfaces_surface_external_mesh + 1
-! endif
-! iglob = ibool(2,NGLLY,2,ispec)
-! if (iglob_is_surface_external_mesh(iglob)) then
-! nfaces_surface_external_mesh = nfaces_surface_external_mesh + 1
-! endif
-! iglob = ibool(1,2,2,ispec)
-! if (iglob_is_surface_external_mesh(iglob)) then
-! nfaces_surface_external_mesh = nfaces_surface_external_mesh + 1
-! endif
-! iglob = ibool(NGLLX,2,2,ispec)
-! if (iglob_is_surface_external_mesh(iglob)) then
-! nfaces_surface_external_mesh = nfaces_surface_external_mesh + 1
-! endif
-! enddo ! NSPEC_AB
-
allocate(nfaces_perproc_surface_ext_mesh(NPROC))
allocate(faces_surface_offset_ext_mesh(NPROC))
if (nfaces_surface_external_mesh == 0) then
@@ -135,6 +107,7 @@
faces_surface_offset_ext_mesh(:) = faces_surface_offset_ext_mesh(:)*NGNOD2D
endif
+! stores global indices of GLL points on the surface to array faces_surface_external_mesh
nfaces_surface_external_mesh = 0
do ispec = 1, NSPEC_AB
if (ispec_is_surface_external_mesh(ispec)) then
Modified: seismo/3D/SPECFEM3D_SESAME/trunk/specfem3D_par.f90
===================================================================
--- seismo/3D/SPECFEM3D_SESAME/trunk/specfem3D_par.f90 2009-10-25 23:06:34 UTC (rev 15876)
+++ seismo/3D/SPECFEM3D_SESAME/trunk/specfem3D_par.f90 2009-10-25 23:54:42 UTC (rev 15877)
@@ -138,9 +138,12 @@
xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz,jacobian
real(kind=CUSTOM_REAL), dimension(:), allocatable :: xstore,ystore,zstore
+! material properties
+ ! isotropic
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
kappastore,mustore
+
! flag for sediments
! logical, dimension(:), allocatable :: not_fully_in_bedrock
! logical, dimension(:,:,:,:), allocatable :: flag_sediments
@@ -368,6 +371,14 @@
! Stacey
real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rho_vp,rho_vs
+ ! anisotropic
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
+ c11store,c12store,c13store,c14store,c15store,c16store,&
+ c22store,c23store,c24store,c25store,c26store,c33store,&
+ c34store,c35store,c36store,c44store,c45store,c46store,&
+ c55store,c56store,c66store
+ integer :: NSPEC_ANISO
+
! maximum of the norm of the displacement
real(kind=CUSTOM_REAL) Usolidnorm,Usolidnorm_all
integer:: Usolidnorm_index(1)
More information about the CIG-COMMITS
mailing list