[cig-commits] r13450 - in seismo/2D/SPECFEM2D/branches/BIOT: . DATA
cmorency at geodynamics.org
cmorency at geodynamics.org
Wed Dec 3 11:47:26 PST 2008
Author: cmorency
Date: 2008-12-03 11:47:26 -0800 (Wed, 03 Dec 2008)
New Revision: 13450
Modified:
seismo/2D/SPECFEM2D/branches/BIOT/DATA/Par_file
seismo/2D/SPECFEM2D/branches/BIOT/constants.h
seismo/2D/SPECFEM2D/branches/BIOT/meshfem2D.F90
Log:
creation of CMTSOLUTION for BIOT2D to handle multiple sources
Modified: seismo/2D/SPECFEM2D/branches/BIOT/DATA/Par_file
===================================================================
--- seismo/2D/SPECFEM2D/branches/BIOT/DATA/Par_file 2008-12-03 19:38:24 UTC (rev 13449)
+++ seismo/2D/SPECFEM2D/branches/BIOT/DATA/Par_file 2008-12-03 19:47:26 UTC (rev 13450)
@@ -43,20 +43,6 @@
#source parameters
NSOURCE = 1 # number of sources
-#source 1
-source_surf = .false. # source inside the medium or at the surface
-xs = 1600. # source location x in meters
-zs = 2900. # source location z in meters
-source_type = 1 # elastic force or acoustic pressure = 1 or moment tensor = 2
-time_function_type = 1 # Ricker = 1, first derivative = 2, Gaussian = 3, Dirac = 4, Heaviside = 5
-f0 = 15.0 # dominant source frequency (Hz) if not Dirac or Heaviside
-t0 = 0.0 # time shift when multi sources (if one source, must be zero)
-angleforce = 0. # angle of the source (for a force only)
-Mxx = 1. # Mxx component (for a moment tensor source only)
-Mzz = 1. # Mzz component (for a moment tensor source only)
-Mxz = 0. # Mxz component (for a moment tensor source only)
-factor = 1.d10 # amplification factor
-
# constants for attenuation
N_SLS = 2 # number of standard linear solids for attenuation
Qp_attenuation = 136.4376068115 # quality factor P for attenuation
Modified: seismo/2D/SPECFEM2D/branches/BIOT/constants.h
===================================================================
--- seismo/2D/SPECFEM2D/branches/BIOT/constants.h 2008-12-03 19:38:24 UTC (rev 13449)
+++ seismo/2D/SPECFEM2D/branches/BIOT/constants.h 2008-12-03 19:47:26 UTC (rev 13450)
@@ -39,6 +39,9 @@
! uncomment this to write to file instead
! integer, parameter :: IOUT = 41
+! number of lines per source in CMTSOLUTION file
+ integer, parameter :: NLINES_PER_CMTSOLUTION_SOURCE = 13
+
! flags for absorbing boundaries
integer, parameter :: IBOTTOM = 1
integer, parameter :: IRIGHT = 2
Modified: seismo/2D/SPECFEM2D/branches/BIOT/meshfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/branches/BIOT/meshfem2D.F90 2008-12-03 19:38:24 UTC (rev 13449)
+++ seismo/2D/SPECFEM2D/branches/BIOT/meshfem2D.F90 2008-12-03 19:47:26 UTC (rev 13450)
@@ -102,9 +102,12 @@
integer, dimension(:), allocatable :: source_type,time_function_type !yang
integer nrec_total,irec_global_number
double precision, dimension(:),allocatable :: xs,zs,f0,t0,angleforce,Mxx,Mzz,Mxz,factor !yang
- integer NSOURCE, i_source !yang
+ integer NSOURCE, NSOURCES, i_source, icounter, ios !yang
logical, dimension(:),allocatable :: source_surf
double precision xrec,zrec
+! file number for source file
+ integer, parameter :: IIN_SOURCE = 22
+ character(len=150) dummystring
character(len=50) interfacesfile,title
@@ -444,19 +447,36 @@
allocate(Mzz(NSOURCE))
allocate(factor(NSOURCE))
+!chris
+ open(unit=IIN_SOURCE,file='DATA/CMTSOLUTION',iostat=ios,status='old',action='read')
+ if(ios /= 0) stop 'error opening CMTSOLUTION file'
+ icounter = 0
+ do while(ios == 0)
+ read(IIN_SOURCE,"(a)",iostat=ios) dummystring
+ if(ios == 0) icounter = icounter + 1
+ enddo
+ close(IIN_SOURCE)
+ if(mod(icounter,NLINES_PER_CMTSOLUTION_SOURCE) /= 0) &
+ stop 'total number of lines in CMTSOLUTION file should be a multiple of NLINES_PER_CMTSOLUTION_SOURCE'
+ NSOURCES = icounter / NLINES_PER_CMTSOLUTION_SOURCE
+ if(NSOURCES < 1) stop 'need at least one source in CMTSOLUTION file'
+ if(NSOURCES /= NSOURCE) &
+ stop 'total number of sources read is different than declared in Par_file'
+
+ open(unit=IIN_SOURCE,file='DATA/CMTSOLUTION',status='old',action='read')
do i_source=1,NSOURCE
- call read_value_logical(IIN,IGNORE_JUNK,source_surf(i_source))
- call read_value_double_precision(IIN,IGNORE_JUNK,xs(i_source))
- call read_value_double_precision(IIN,IGNORE_JUNK,zs(i_source))
- call read_value_integer(IIN,IGNORE_JUNK,source_type(i_source))
- call read_value_integer(IIN,IGNORE_JUNK,time_function_type(i_source))
- call read_value_double_precision(IIN,IGNORE_JUNK,f0(i_source))
- call read_value_double_precision(IIN,IGNORE_JUNK,t0(i_source))
- call read_value_double_precision(IIN,IGNORE_JUNK,angleforce(i_source))
- call read_value_double_precision(IIN,IGNORE_JUNK,Mxx(i_source))
- call read_value_double_precision(IIN,IGNORE_JUNK,Mzz(i_source))
- call read_value_double_precision(IIN,IGNORE_JUNK,Mxz(i_source))
- call read_value_double_precision(IIN,IGNORE_JUNK,factor(i_source))
+ call read_value_logical(IIN_SOURCE,IGNORE_JUNK,source_surf(i_source))
+ call read_value_double_precision(IIN_SOURCE,IGNORE_JUNK,xs(i_source))
+ call read_value_double_precision(IIN_SOURCE,IGNORE_JUNK,zs(i_source))
+ call read_value_integer(IIN_SOURCE,IGNORE_JUNK,source_type(i_source))
+ call read_value_integer(IIN_SOURCE,IGNORE_JUNK,time_function_type(i_source))
+ call read_value_double_precision(IIN_SOURCE,IGNORE_JUNK,f0(i_source))
+ call read_value_double_precision(IIN_SOURCE,IGNORE_JUNK,t0(i_source))
+ call read_value_double_precision(IIN_SOURCE,IGNORE_JUNK,angleforce(i_source))
+ call read_value_double_precision(IIN_SOURCE,IGNORE_JUNK,Mxx(i_source))
+ call read_value_double_precision(IIN_SOURCE,IGNORE_JUNK,Mzz(i_source))
+ call read_value_double_precision(IIN_SOURCE,IGNORE_JUNK,Mxz(i_source))
+ call read_value_double_precision(IIN_SOURCE,IGNORE_JUNK,factor(i_source))
! if Dirac source time function, use a very thin Gaussian instead
! if Heaviside source time function, use a very thin error function instead
if(time_function_type(i_source) == 4 .or. time_function_type(i_source) == 5) f0(i_source) = 1.d0 / (10.d0 * deltat)
@@ -479,7 +499,8 @@
print *,'Mzz of the source if moment tensor = ',Mzz(i_source)
print *,'Mxz of the source if moment tensor = ',Mxz(i_source)
print *,'Multiplying factor = ',factor(i_source)
- enddo
+ enddo ! do i_source=1,NSOURCE
+ close(IIN_SOURCE)
! read constants for attenuation
call read_value_integer(IIN,IGNORE_JUNK,N_SLS)
More information about the CIG-COMMITS
mailing list