[cig-commits] r16060 - in seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate: . PLOTTING matlab src
carltape at geodynamics.org
carltape at geodynamics.org
Wed Dec 2 10:44:59 PST 2009
Author: carltape
Date: 2009-12-02 10:44:58 -0800 (Wed, 02 Dec 2009)
New Revision: 16060
Added:
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex00.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex01.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex02.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex03.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex04.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex05.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex06.f90
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex07.f90
Modified:
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/joint_inversion.pl
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/README
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_poly.m
seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_subspace.m
Log:
Tested all seven test simulations; updated README.
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/joint_inversion.pl
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/joint_inversion.pl 2009-11-30 23:43:20 UTC (rev 16059)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/PLOTTING/joint_inversion.pl 2009-12-02 18:44:58 UTC (rev 16060)
@@ -638,7 +638,8 @@
# third column is the target source
$shift = "-Y$dY1mB";
-$title = "($labs[$i]) $tlabs[$i]";
+$t = $t+1;
+$title = "($labs[$t]) $tlabs[$i]";
$source_error_file = $file1dat_src;
print CSH "awk '{print \$1,\$2,\$6,${ref_rad} + sqrt(\$7*\$7 + \$8*\$8)*$fac/1000}' $source_error_file > $dots_file\n";
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/README
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/README 2009-11-30 23:43:20 UTC (rev 16059)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/README 2009-12-02 18:44:58 UTC (rev 16060)
@@ -1,15 +1,20 @@
-Carl Tape, 12-Nov-2009
+Carl Tape, 01-Dec-2009
carltape at fas.harvard.edu
/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/README
-This is a developmental code for testing synthetic sesimic inversions using a 2D spectral-element method (SEM) code as the forward solver. The official 2D SEM solver, SPECFEM2D, is mantained by CIG. This version emphasizes the details in embedding the forward solver within the inverse problem, which is not contained in SPECFEM2D.
+This is a developmental code for testing synthetic sesimic inversions using a 2D spectral-element method (SEM) code as the forward solver. The official 2D SEM solver, SPECFEM2D, is mantained by CIG here:
+/cig/seismo/2D/SPECFEM2D/
+
+SEM2D_iterate emphasizes the details in embedding the forward solver within the inverse problem, which is not contained in SPECFEM2D.
+
The code can run in P-SV mode as a depth cross-section (Tromp, Tape, Liu, GJI, 2005).
Or it can run in SH mode as a membrane surface wave analogue (Tape, Liu, Tromp, GJI, 2007).
There are several perl-based GMT plotting scripts.
-The user needs access to gfortran, matlab, perl, and GMT, but not SAC.
+To compile and run, the user needs access to gfortran.
+For full plotting capabilities, the user needs perl, matlab, and GMT.
The best reference for the inverse method implemented here is:
Tarantola (2005), Inverse Problem Theory and Methods for Model Parameter Estimation
@@ -106,7 +111,8 @@
Run wave2d_cg_poly.m and then wave2d_cg_figs.m; this generates
data files that will be used in joint_inversion.pl.
cd ../PLOTTING/FIGURES
-Execute the example line in joint_inversion.pl --> joint_inversion_0100_h0.ps
+Open ../joint_inversion.pl
+Execute the example line --> joint_inversion_0100_h0.ps
------------------
EXAMPLE 2 (run_0200)
@@ -177,6 +183,9 @@
make clean ; make wave2d ; wave2d
Create a directory for the subspace update
mkdir OUTPUT/run_0700/HESSIAN
+ mkdir OUTPUT/run_0700/HESSIAN/model_m01
+ mkdir OUTPUT/run_0700/HESSIAN/model_m02
+ (make as many directories as you might need)
Open the matlab program matlab/wave2d_subspace.m
Enter the following values at the prompts:
-------------------------
@@ -192,20 +201,21 @@
Enter 1 if this balance looks okay : 1
Enter 0 for GCV, 1 for OCV, 2 for L-curve : 1
--------------------------------
-If everything has worked, you will see several sensible figures.
-You now have a model update computed using the subspace method.
-Now you want to compute the misfit function and event kernels using this model.
+If everything has worked, you will see several sensible figures. You now have a model update dm00 (and a new model m01) computed using the subspace method. Now you want to compute the misfit function and event kernels using m01.
Make the following modification to src/wave2d_constants.f90
HESSIAN = .false. --> HESSIAN = .true.
Then compile and execute:
make clean ; make wave2d ; wave2d
-After this is done, go back to wave2d_subspace.m and repeat as before.
-Only make this change
+After this is done, go back to wave2d_subspace.m and repeat as before, but make this change:
Enter next model number (hmod) : 2
+You should see more figures and write out model m02.
Then back to wave2d.f90 and modify
hmod = 1 --> hmod = 2
-This tells the code where to look for the new model.
+This tells the code where to look for model m02.
+Then compile and execute:
+ make clean ; make wave2d ; wave2d
+From PLOTTING/FIGURES/, execute ../joint_inversion.pl with $irun0_cg = 600
-As you can see, it is not ideal to go back-and-forth between Matlab and Fortran. However, many aspects are still in development, and Matlab provides quick graphics that are crucial for debugging, checking, and testing.
+Clearly for the subspace approach, it is not ideal to go back-and-forth between Matlab and Fortran. However, many aspects are still in development, and Matlab provides quick graphics that are crucial for debugging, checking, and testing.
==================================================================
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_poly.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_poly.m 2009-11-30 23:43:20 UTC (rev 16059)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_cg_poly.m 2009-12-02 18:44:58 UTC (rev 16060)
@@ -24,17 +24,17 @@
close all
clear
-% USER: change this
-dir0 = '/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_work/';
-if ~exist(dir0), error([dir0 ' does not exist']); end
-
colors;
npts = 100;
stfm = '%4.4i';
%---------------------------------------------------------
-% KEY COMMANDS
+% USER INPUT
+% base directory
+dir0 = '/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_work/';
+if ~exist(dir0), error([dir0 ' does not exist']); end
+
parms = [100 0 3]; % structure inversion: 1 source, Nfac=3
%parms = [200 0 3]; % source inversion (xs,ys,ts): single source
%parms = [300 0 3]; % joint inversion: 1 source, Nfac=3
Modified: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_subspace.m
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_subspace.m 2009-11-30 23:43:20 UTC (rev 16059)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/matlab/wave2d_subspace.m 2009-12-02 18:44:58 UTC (rev 16060)
@@ -19,9 +19,11 @@
ax1 = [-121 -114 31 37];
stfm = '%4.4i';
-% USER CHANGE
-dir_run = '/home/carltape/wave2d/SEM2D_iterate/OUTPUT/';
+%----------------------------------------------
+% USER INPUT
+dir_run = '/data1/cig/seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate_OUTPUT/';
+
inew = 0; % subspace AND parameter class inversions
%----------------------------------------------
@@ -840,7 +842,7 @@
% plot structure model update and new model
if and(INV_STRUCT == 1, ifigure == 1)
- tlabs = {{dir2lab,'data'},'m0',['dm (' stp ')'],['m1 (' stp ')']};
+ tlabs = {{dir2lab,'data'},'m00',['dm (' stp ')'],['m01 (' stp ')']};
cmax = 0.1;
figure; nr=3; nc=2;
for kk = 1:4
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex00.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex00.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex00.f90 2009-12-02 18:44:58 UTC (rev 16060)
@@ -0,0 +1,336 @@
+module wave2d_constants
+
+! This file is copied into the output directory when wave2d.f90 is run
+! There are two basic options: membrane (surface) waves and body waves.
+! Several lines need to be commented/uncommented to go from one to the other.
+
+! index of the reference directory for the simulation output
+ integer, parameter :: IRUNZ = 0
+
+!========================
+! GRID, TIME-STEP, AND SOURCE PARAMETERS
+
+! NFRAME : number of frames to save: membrane, 10; body, 20
+! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
+! NSTEP : number of timesteps
+ integer, parameter :: NFRAME = 10 ! membrane
+ !integer, parameter :: NFRAME = 11 ! body
+ integer, parameter :: NSAVE = 400 ! 200,400
+ integer, parameter :: NSTEP = NFRAME*NSAVE
+
+! time step in seconds
+ !double precision, parameter :: DT = 2.0e-2 ! body waves
+ double precision, parameter :: DT = 6.0e-2 ! membrane surface waves
+
+! temporal properties of source (source time function)
+ integer, parameter :: ISRC_TIME = 1 ! type (1)
+ double precision, parameter :: hdur = 10.0 ! HALF-duration (s), membrane
+ !double precision, parameter :: hdur = 2.0 ! HALF-duration (s), body
+ double precision, parameter :: tshift = 2.0*DT*dble(NSAVE) ! time shift (s)
+ !double precision, parameter :: tshift = 8.0*hdur
+ logical, parameter :: SRC_TAPER = .true. ! taper the endpoints of the time series
+
+! normalization factor of point source force
+ double precision, parameter :: FNORM = 1.0d10
+
+! forward wavefield source-time function
+! (For membrane waves, FOR_Y = 1 is all that matters.)
+ integer, parameter :: FOR_X = 1 ! boolean
+ integer, parameter :: FOR_Y = 1 ! boolean
+ integer, parameter :: FOR_Z = 1 ! boolean
+
+! adjoint wavefield source-time function
+! (For membrane waves, REV_Y = 1 is all that matters.)
+! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
+ integer, parameter :: REV_X = 1 ! boolean
+ integer, parameter :: REV_Y = 1 ! boolean
+ integer, parameter :: REV_Z = 1 ! boolean
+
+! spatial properties of sources
+! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
+! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
+! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
+ integer, parameter :: ISRC_SPACE = 6 ! see wave2d.f90
+
+! spatial properties of receivers
+! IREC_SPACE
+! (1) individual station(s)
+! (2) SoCal (used for GJI paper)
+! (3) regular mesh on land
+! (4) regular mesh
+! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
+! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
+! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
+ integer, parameter :: IREC_SPACE = 2 ! see wave2d.f90
+ integer, parameter :: NMESH_REC = 17
+ double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m ! 4km for membrane surface waves
+ double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m ! 15km for membrane surface waves
+ double precision, parameter :: STATION_COAST_BUFFER = 0.0e3 ! m
+
+! lower right corner for membrane surface waves plotting grid
+ double precision, parameter :: LAT_MIN = 32.0
+ double precision, parameter :: LON_MIN = -120.0
+ integer, parameter :: UTM_PROJECTION_ZONE = 11 ! southern California
+
+! mesh specifications: membrane surface waves
+ double precision, parameter :: LENGTH = 480.0e3 ! m
+ double precision, parameter :: HEIGHT = 480.0e3 ! m
+ double precision, parameter :: AREA = LENGTH*HEIGHT
+ integer, parameter :: NEX = 40
+ integer, parameter :: NEZ = 40
+!!$
+!!$! mesh specifications: body waves
+!!$ double precision, parameter :: LENGTH = 200.0e3 ! m ! 400 for 1D body waves
+!!$ double precision, parameter :: HEIGHT = 80.0e3 ! m
+!!$ integer, parameter :: NEX = 80 ! 160
+!!$ integer, parameter :: NEZ = 32 ! 32
+
+!========================
+! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
+
+! model perturbations for HOMOGENEOUS model (or perturbation)
+! scaling from beta to alpha
+! value is from Master et al. (2000), "The relative behavior of shear velocity..."
+ double precision, parameter :: R_BETA_OVER_ALPHA = 1.3
+ double precision, parameter :: PBETA = 10.0
+ double precision, parameter :: PALPHA = PBETA / R_BETA_OVER_ALPHA
+ double precision, parameter :: PRHO = 0.0
+
+ ! reference model and target model choice
+ integer, parameter :: IMODEL_SYN = 0
+ integer, parameter :: IMODEL_DAT = 2
+ !-----------------------------------------------------------------------------------------
+ ! 0 1 2 3
+ !-----------------------------------------------------------------------------------------
+ ! ISURFACE=1, IMODEL_SYN : homo checker het
+ ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het pert
+ !-----------------------------------------------------------------------------------------
+ ! ISURFACE=0, IMODEL_SYN : homo 1D model checker het
+ ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert het pert
+ !-----------------------------------------------------------------------------------------
+
+ ! whether you want to smooth the structure model for computing the synthetics
+ integer, parameter :: ISMOOTH_INITIAL_MODEL = 0
+ integer, parameter :: ISMOOTH_MODEL_UPDATE = 0
+ integer, parameter :: ISMOOTH_EVENT_KERNEL = 1
+ integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
+
+ ! smoothing scalelength
+ double precision, parameter :: GAMMA_SMOOTH_MODEL = 30.0e3
+
+!========================
+
+! boolean parameters
+! IKER: (0) waveform
+! (1) traveltime, cross-correlation, misfit
+! (2) amplitude, cross-correlation, misfit
+! (3) traveltime, multitaper
+! (4) amplitude, multitaper
+! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
+! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
+ integer, parameter :: IKER = 1
+ integer, parameter :: IAMP_VEL = 0 ! measure amplitudes between velocity traces
+
+! KEY: USE ONE LINE OR THE OTHER
+ integer, parameter :: ISURFACE = 1, NCOMP = 1, NABSORB = 4 ! surface waves
+! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
+
+! parameters controlling what to write to file
+! NOTE: for the tomography simulations, ALL of these can be .false.
+
+ logical, parameter :: WRITE_STF_F = .false.
+ logical, parameter :: WRITE_SEISMO_F = .false.
+! logical, parameter :: WRITE_SPECTRA_F = .false.
+! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
+
+ logical, parameter :: WRITE_SEISMO_RECONSTRUCT = .false. ! multitaper
+! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
+
+ logical, parameter :: WRITE_STF_A = .false.
+ logical, parameter :: WRITE_SEISMO_A = .false. ! source inversions
+! logical, parameter :: WRITE_SPECTRA_A = .false.
+! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
+
+ logical, parameter :: WRITE_KERNELS = .true. ! write all nine kernels
+ logical, parameter :: WRITE_KERNEL_SNAPSHOTS = .false. ! kernel snapshots
+ logical, parameter :: WRITE_WAVFIELD_SNAPSHOTS = .false. ! wavefield snapshots
+
+!--------------------------------------
+! INVERSION PARAMETERS
+
+ ! whether you want to compute kernels, or simply the misfit function
+ logical, parameter :: COMPUTE_KERNELS = .true.
+
+ ! whether to use the data subspace method, which has a Hessian
+ logical, parameter :: HESSIAN = .false.
+
+ ! stopping criteria
+ ! NITERATION : number of iterations
+ ! VAR_RED_MIN : minimum variance reduction (in percent)
+ ! SIGMA_FAC : stop if a model value exceeds SIGMA_FAC * sigma_m
+ ! CONV_STOP : stop when the misfit value is this fraction of the INITIAL misfit value
+ integer, parameter :: NITERATION = 1
+ double precision, parameter :: VAR_RED_MIN = 8.0
+
+ ! Gaussian errors containted in input file
+ ! see wave2d_sigmas.m and INPUT/sigma_0p1_pert.dat
+ double precision, parameter :: SIGMA_DT = 0.10
+ double precision, parameter :: SIGMA_DLNA = 1.0
+ double precision, parameter :: SIGMA_WAVEFORM = 1.0
+ logical, parameter :: ADD_DATA_ERRORS = .true.
+
+ !double precision, parameter :: SIGMA_FAC = 2.0
+ !double precision, parameter :: CONV_STOP = 1.0e-4
+
+ ! order of interpolating polynomial in conjugate gradient algorithm
+ ! using POLY_ORDER = 3 required computing the gradient of the test model
+ integer, parameter :: POLY_ORDER = 2 ! 2 (preferred) or 3
+
+ ! use inversion and scaling used for the GJI paper
+ !integer, parameter :: GJI_PAPER = 1
+
+ ! what to perturb, what to invert
+ integer, parameter :: PERT_STRUCT = 1
+ integer, parameter :: PERT_SOURCE_T = 0
+ integer, parameter :: PERT_SOURCE_X = 0
+
+ integer, parameter :: INV_STRUCT = 1
+ integer, parameter :: INV_SOURCE_T = 0
+ integer, parameter :: INV_SOURCE_X = 0
+
+ ! whether to include the model norm term in the misfit function, which acts like damping
+ logical, parameter :: INCLUDE_MODEL_NORM = .true.
+
+ ! log file showing source loactions
+ logical, parameter :: ISOURCE_LOG = .true.
+
+ ! DO NOT CHANGE THESE
+ integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ integer, parameter :: NVAR_SOURCE = 3 ! x position, z position, origin time
+ integer, parameter :: NVAR = NVAR_STRUCT + NVAR_SOURCE
+
+! parameterization for structure in the inversion
+! 1 : kappa-mu-rho
+! 2 : alpha-beta-rho
+! 3 : c-beta-rho
+ integer, parameter :: STRUCTURE_PARAMETER_TYPE = 2
+
+! scalelength of smoothing Gaussian
+! GJI paper : 30,60,90
+! body waves : 10
+ double precision, parameter :: GAMMA_SMOOTH = 60.0e3
+
+ ! parameters for smoothing
+ logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
+ logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
+
+! homogeneous background model (S.I. units)
+ double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
+ double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
+ double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
+
+!---------------------------------------------------------------
+! measurement windows
+
+ double precision, parameter :: HWIN1 = 2.5*hdur ! half-window width for pulse (2.5)
+ double precision, parameter :: HWIN2 = 5.0*hdur ! half-window with for MTM measurement (zero-padding)
+
+!---------------------------------------------------------------
+! CHT: do not change these
+
+! UTM zone for Southern California region
+! integer, parameter :: UTM_PROJECTION_ZONE = 11
+
+! to suppress UTM projection for SCEC benchmarks
+ logical, parameter :: SUPPRESS_UTM_PROJECTION = .false.
+
+! flag for projection from latitude/longitude to UTM, and back
+ integer, parameter :: ILONGLAT2UTM = 0, IUTM2LONGLAT = 1
+
+! flag for projection from latitude/longitude to mesh-UTM, and back
+ integer, parameter :: ILONLAT2MESH = 0, IMESH2LONLAT = 1
+
+! max number of fake receivers
+ integer, parameter :: MAX_SR_FAKE = 1000
+
+! max number of events, receivers, and phass
+ integer, parameter :: MAX_EVENT = 50
+ integer, parameter :: MAX_SR = 1400
+ integer, parameter :: MAX_PHASE = 1
+ integer, parameter :: MAX_COMP = NCOMP
+
+!========================
+! GRID AND GLL POINTS
+
+ integer, parameter :: NELE = MAX(NEX,NEZ)
+ integer, parameter :: NSPEC = NEX*NEZ
+
+! number of GLL points (polynomial degree plus one)
+ integer, parameter :: NGLLX = 5
+ integer, parameter :: NGLLZ = 5
+ integer, parameter :: NGLL = MAX(NGLLX,NGLLZ)
+
+! number of global points
+ integer, parameter :: NGLLSQUARE = NGLLX * NGLLZ ! element GLL points
+ integer, parameter :: NGLOB = ((NGLLX-1)*NEX + 1)*((NGLLZ-1)*NEZ +1) ! global GLL points
+ integer, parameter :: NLOCAL = NGLLX * NGLLZ * NSPEC ! local GLL points
+ integer, parameter :: NSPEC_CORNER = (NEX+1) * (NEZ+1) ! element corner points
+
+! number of nodes for 2D and 3D shape functions for hexahedra
+! we use 8-node mesh bricks, which are more stable than 27-node elements
+ integer, parameter :: NGNOD = 8, NGNOD2D = 4
+
+! number of iterations to solve the system for xi and eta
+ integer, parameter :: NUM_ITER = 1
+
+! very large and very small values
+ double precision, parameter :: HUGEVAL = 1.0e30, TINYVAL = 1.0e-9
+
+! for the Gauss-Lobatto-Legendre points and weights
+ double precision, parameter :: GAUSSALPHA = 0.0,GAUSSBETA = 0.0
+
+!
+! CONSTANTS
+!
+ double precision, parameter :: PI = 3.141592653589793
+ double precision, parameter :: FOUR_THIRDS = 4.0 / 3.0
+ double precision, parameter :: ONE_THIRD = 1.0 / 3.0
+ double precision, parameter :: ONEOVERTWO = 0.5
+ double precision, parameter :: DEG = 180.0/PI
+! double precision, parameter :: EPS = 1.0e-35
+
+!!$! parameter for FFTW
+!!$ integer, parameter :: NOUT = NSTEP/2 + 1
+!!$
+!!$! bounds for bandpass filter
+!!$ double precision, parameter :: hwid = 3.0 ! HALF-width of window (s)
+!!$ double precision, parameter :: tmin = 2.0*hdur-hwid
+!!$ double precision, parameter :: tmax = 2.0*hdur+hwid
+!!$ double precision, parameter :: fmin = 1.0/tmax, fmax = 1.0/tmin
+!!$ double precision, parameter :: trbdndw = 0.3, a = 30.0
+!!$ integer, parameter :: passes = 2, iord = 4
+!!$
+!!$!
+!!$! MULTI-TAPER PARAMETERS
+!!$!
+!!$! Ying Zhou: The fit between the recovered data and the data can be improved
+!!$! by either increasing the window width (HWIN above) or by decreasing NPI.
+!!$! In her experience, NPI = 2.5 is good for noisy data.
+!!$! For synthetic data, we can use a lower NPI.
+!!$! number of tapers should be fixed as twice NPI -- see Latex notes
+!!$ double precision, parameter :: WTR = 0.02
+!!$ double precision, parameter :: NPI = 2.5
+!!$ integer, parameter :: NTAPER = int(2.0*NPI)
+!!$
+!!$! FFT parameters
+!!$! FOR SOME REASON, I CANNOT SET NDIM = npt when lnpt <= 11 -- WHY NOT?
+!!$ integer, parameter :: lnpt = 13, npt = 2**lnpt, NDIM = 20000
+!!$ !integer, parameter :: lnpt = 14, npt = 2**lnpt, NDIM = npt
+!!$
+!!$ !double precision, parameter :: ZZIGN = -1.0 ! sign convention -- WATCH OUT
+!!$ double precision, parameter :: FORWARD_FFT = 1.0
+!!$ double precision, parameter :: REVERSE_FFT = -1.0
+
+ integer, parameter :: NDIM = 20000
+
+end module wave2d_constants
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex01.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex01.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex01.f90 2009-12-02 18:44:58 UTC (rev 16060)
@@ -0,0 +1,336 @@
+module wave2d_constants
+
+! This file is copied into the output directory when wave2d.f90 is run
+! There are two basic options: membrane (surface) waves and body waves.
+! Several lines need to be commented/uncommented to go from one to the other.
+
+! index of the reference directory for the simulation output
+ integer, parameter :: IRUNZ = 100
+
+!========================
+! GRID, TIME-STEP, AND SOURCE PARAMETERS
+
+! NFRAME : number of frames to save: membrane, 10; body, 20
+! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
+! NSTEP : number of timesteps
+ integer, parameter :: NFRAME = 10 ! membrane
+ !integer, parameter :: NFRAME = 11 ! body
+ integer, parameter :: NSAVE = 400 ! 200,400
+ integer, parameter :: NSTEP = NFRAME*NSAVE
+
+! time step in seconds
+ !double precision, parameter :: DT = 2.0e-2 ! body waves
+ double precision, parameter :: DT = 6.0e-2 ! membrane surface waves
+
+! temporal properties of source (source time function)
+ integer, parameter :: ISRC_TIME = 1 ! type (1)
+ double precision, parameter :: hdur = 10.0 ! HALF-duration (s), membrane
+ !double precision, parameter :: hdur = 2.0 ! HALF-duration (s), body
+ double precision, parameter :: tshift = 2.0*DT*dble(NSAVE) ! time shift (s)
+ !double precision, parameter :: tshift = 8.0*hdur
+ logical, parameter :: SRC_TAPER = .true. ! taper the endpoints of the time series
+
+! normalization factor of point source force
+ double precision, parameter :: FNORM = 1.0d10
+
+! forward wavefield source-time function
+! (For membrane waves, FOR_Y = 1 is all that matters.)
+ integer, parameter :: FOR_X = 1 ! boolean
+ integer, parameter :: FOR_Y = 1 ! boolean
+ integer, parameter :: FOR_Z = 1 ! boolean
+
+! adjoint wavefield source-time function
+! (For membrane waves, REV_Y = 1 is all that matters.)
+! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
+ integer, parameter :: REV_X = 1 ! boolean
+ integer, parameter :: REV_Y = 1 ! boolean
+ integer, parameter :: REV_Z = 1 ! boolean
+
+! spatial properties of sources
+! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
+! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
+! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
+ integer, parameter :: ISRC_SPACE = 6 ! see wave2d.f90
+
+! spatial properties of receivers
+! IREC_SPACE
+! (1) individual station(s)
+! (2) SoCal (used for GJI paper)
+! (3) regular mesh on land
+! (4) regular mesh
+! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
+! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
+! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
+ integer, parameter :: IREC_SPACE = 2 ! see wave2d.f90
+ integer, parameter :: NMESH_REC = 17
+ double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m ! 4km for membrane surface waves
+ double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m ! 15km for membrane surface waves
+ double precision, parameter :: STATION_COAST_BUFFER = 0.0e3 ! m
+
+! lower right corner for membrane surface waves plotting grid
+ double precision, parameter :: LAT_MIN = 32.0
+ double precision, parameter :: LON_MIN = -120.0
+ integer, parameter :: UTM_PROJECTION_ZONE = 11 ! southern California
+
+! mesh specifications: membrane surface waves
+ double precision, parameter :: LENGTH = 480.0e3 ! m
+ double precision, parameter :: HEIGHT = 480.0e3 ! m
+ double precision, parameter :: AREA = LENGTH*HEIGHT
+ integer, parameter :: NEX = 40
+ integer, parameter :: NEZ = 40
+!!$
+!!$! mesh specifications: body waves
+!!$ double precision, parameter :: LENGTH = 200.0e3 ! m ! 400 for 1D body waves
+!!$ double precision, parameter :: HEIGHT = 80.0e3 ! m
+!!$ integer, parameter :: NEX = 80 ! 160
+!!$ integer, parameter :: NEZ = 32 ! 32
+
+!========================
+! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
+
+! model perturbations for HOMOGENEOUS model (or perturbation)
+! scaling from beta to alpha
+! value is from Master et al. (2000), "The relative behavior of shear velocity..."
+ double precision, parameter :: R_BETA_OVER_ALPHA = 1.3
+ double precision, parameter :: PBETA = 10.0
+ double precision, parameter :: PALPHA = PBETA / R_BETA_OVER_ALPHA
+ double precision, parameter :: PRHO = 0.0
+
+ ! reference model and target model choice
+ integer, parameter :: IMODEL_SYN = 0
+ integer, parameter :: IMODEL_DAT = 2
+ !-----------------------------------------------------------------------------------------
+ ! 0 1 2 3
+ !-----------------------------------------------------------------------------------------
+ ! ISURFACE=1, IMODEL_SYN : homo checker het
+ ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het pert
+ !-----------------------------------------------------------------------------------------
+ ! ISURFACE=0, IMODEL_SYN : homo 1D model checker het
+ ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert het pert
+ !-----------------------------------------------------------------------------------------
+
+ ! whether you want to smooth the structure model for computing the synthetics
+ integer, parameter :: ISMOOTH_INITIAL_MODEL = 0
+ integer, parameter :: ISMOOTH_MODEL_UPDATE = 0
+ integer, parameter :: ISMOOTH_EVENT_KERNEL = 1
+ integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
+
+ ! smoothing scalelength
+ double precision, parameter :: GAMMA_SMOOTH_MODEL = 30.0e3
+
+!========================
+
+! boolean parameters
+! IKER: (0) waveform
+! (1) traveltime, cross-correlation, misfit
+! (2) amplitude, cross-correlation, misfit
+! (3) traveltime, multitaper
+! (4) amplitude, multitaper
+! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
+! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
+ integer, parameter :: IKER = 1
+ integer, parameter :: IAMP_VEL = 0 ! measure amplitudes between velocity traces
+
+! KEY: USE ONE LINE OR THE OTHER
+ integer, parameter :: ISURFACE = 1, NCOMP = 1, NABSORB = 4 ! surface waves
+! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
+
+! parameters controlling what to write to file
+! NOTE: for the tomography simulations, ALL of these can be .false.
+
+ logical, parameter :: WRITE_STF_F = .false.
+ logical, parameter :: WRITE_SEISMO_F = .false.
+! logical, parameter :: WRITE_SPECTRA_F = .false.
+! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
+
+ logical, parameter :: WRITE_SEISMO_RECONSTRUCT = .false. ! multitaper
+! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
+
+ logical, parameter :: WRITE_STF_A = .false.
+ logical, parameter :: WRITE_SEISMO_A = .false. ! source inversions
+! logical, parameter :: WRITE_SPECTRA_A = .false.
+! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
+
+ logical, parameter :: WRITE_KERNELS = .true. ! write all nine kernels
+ logical, parameter :: WRITE_KERNEL_SNAPSHOTS = .false. ! kernel snapshots
+ logical, parameter :: WRITE_WAVFIELD_SNAPSHOTS = .false. ! wavefield snapshots
+
+!--------------------------------------
+! INVERSION PARAMETERS
+
+ ! whether you want to compute kernels, or simply the misfit function
+ logical, parameter :: COMPUTE_KERNELS = .true.
+
+ ! whether to use the data subspace method, which has a Hessian
+ logical, parameter :: HESSIAN = .false.
+
+ ! stopping criteria
+ ! NITERATION : number of iterations
+ ! VAR_RED_MIN : minimum variance reduction (in percent)
+ ! SIGMA_FAC : stop if a model value exceeds SIGMA_FAC * sigma_m
+ ! CONV_STOP : stop when the misfit value is this fraction of the INITIAL misfit value
+ integer, parameter :: NITERATION = 16
+ double precision, parameter :: VAR_RED_MIN = 8.0
+
+ ! Gaussian errors containted in input file
+ ! see wave2d_sigmas.m and INPUT/sigma_0p1_pert.dat
+ double precision, parameter :: SIGMA_DT = 0.10
+ double precision, parameter :: SIGMA_DLNA = 1.0
+ double precision, parameter :: SIGMA_WAVEFORM = 1.0
+ logical, parameter :: ADD_DATA_ERRORS = .true.
+
+ !double precision, parameter :: SIGMA_FAC = 2.0
+ !double precision, parameter :: CONV_STOP = 1.0e-4
+
+ ! order of interpolating polynomial in conjugate gradient algorithm
+ ! using POLY_ORDER = 3 required computing the gradient of the test model
+ integer, parameter :: POLY_ORDER = 2 ! 2 (preferred) or 3
+
+ ! use inversion and scaling used for the GJI paper
+ !integer, parameter :: GJI_PAPER = 1
+
+ ! what to perturb, what to invert
+ integer, parameter :: PERT_STRUCT = 1
+ integer, parameter :: PERT_SOURCE_T = 0
+ integer, parameter :: PERT_SOURCE_X = 0
+
+ integer, parameter :: INV_STRUCT = 1
+ integer, parameter :: INV_SOURCE_T = 0
+ integer, parameter :: INV_SOURCE_X = 0
+
+ ! whether to include the model norm term in the misfit function, which acts like damping
+ logical, parameter :: INCLUDE_MODEL_NORM = .true.
+
+ ! log file showing source loactions
+ logical, parameter :: ISOURCE_LOG = .true.
+
+ ! DO NOT CHANGE THESE
+ integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ integer, parameter :: NVAR_SOURCE = 3 ! x position, z position, origin time
+ integer, parameter :: NVAR = NVAR_STRUCT + NVAR_SOURCE
+
+! parameterization for structure in the inversion
+! 1 : kappa-mu-rho
+! 2 : alpha-beta-rho
+! 3 : c-beta-rho
+ integer, parameter :: STRUCTURE_PARAMETER_TYPE = 2
+
+! scalelength of smoothing Gaussian
+! GJI paper : 30,60,90
+! body waves : 10
+ double precision, parameter :: GAMMA_SMOOTH = 60.0e3
+
+ ! parameters for smoothing
+ logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
+ logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
+
+! homogeneous background model (S.I. units)
+ double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
+ double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
+ double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
+
+!---------------------------------------------------------------
+! measurement windows
+
+ double precision, parameter :: HWIN1 = 2.5*hdur ! half-window width for pulse (2.5)
+ double precision, parameter :: HWIN2 = 5.0*hdur ! half-window with for MTM measurement (zero-padding)
+
+!---------------------------------------------------------------
+! CHT: do not change these
+
+! UTM zone for Southern California region
+! integer, parameter :: UTM_PROJECTION_ZONE = 11
+
+! to suppress UTM projection for SCEC benchmarks
+ logical, parameter :: SUPPRESS_UTM_PROJECTION = .false.
+
+! flag for projection from latitude/longitude to UTM, and back
+ integer, parameter :: ILONGLAT2UTM = 0, IUTM2LONGLAT = 1
+
+! flag for projection from latitude/longitude to mesh-UTM, and back
+ integer, parameter :: ILONLAT2MESH = 0, IMESH2LONLAT = 1
+
+! max number of fake receivers
+ integer, parameter :: MAX_SR_FAKE = 1000
+
+! max number of events, receivers, and phass
+ integer, parameter :: MAX_EVENT = 50
+ integer, parameter :: MAX_SR = 1400
+ integer, parameter :: MAX_PHASE = 1
+ integer, parameter :: MAX_COMP = NCOMP
+
+!========================
+! GRID AND GLL POINTS
+
+ integer, parameter :: NELE = MAX(NEX,NEZ)
+ integer, parameter :: NSPEC = NEX*NEZ
+
+! number of GLL points (polynomial degree plus one)
+ integer, parameter :: NGLLX = 5
+ integer, parameter :: NGLLZ = 5
+ integer, parameter :: NGLL = MAX(NGLLX,NGLLZ)
+
+! number of global points
+ integer, parameter :: NGLLSQUARE = NGLLX * NGLLZ ! element GLL points
+ integer, parameter :: NGLOB = ((NGLLX-1)*NEX + 1)*((NGLLZ-1)*NEZ +1) ! global GLL points
+ integer, parameter :: NLOCAL = NGLLX * NGLLZ * NSPEC ! local GLL points
+ integer, parameter :: NSPEC_CORNER = (NEX+1) * (NEZ+1) ! element corner points
+
+! number of nodes for 2D and 3D shape functions for hexahedra
+! we use 8-node mesh bricks, which are more stable than 27-node elements
+ integer, parameter :: NGNOD = 8, NGNOD2D = 4
+
+! number of iterations to solve the system for xi and eta
+ integer, parameter :: NUM_ITER = 1
+
+! very large and very small values
+ double precision, parameter :: HUGEVAL = 1.0e30, TINYVAL = 1.0e-9
+
+! for the Gauss-Lobatto-Legendre points and weights
+ double precision, parameter :: GAUSSALPHA = 0.0,GAUSSBETA = 0.0
+
+!
+! CONSTANTS
+!
+ double precision, parameter :: PI = 3.141592653589793
+ double precision, parameter :: FOUR_THIRDS = 4.0 / 3.0
+ double precision, parameter :: ONE_THIRD = 1.0 / 3.0
+ double precision, parameter :: ONEOVERTWO = 0.5
+ double precision, parameter :: DEG = 180.0/PI
+! double precision, parameter :: EPS = 1.0e-35
+
+!!$! parameter for FFTW
+!!$ integer, parameter :: NOUT = NSTEP/2 + 1
+!!$
+!!$! bounds for bandpass filter
+!!$ double precision, parameter :: hwid = 3.0 ! HALF-width of window (s)
+!!$ double precision, parameter :: tmin = 2.0*hdur-hwid
+!!$ double precision, parameter :: tmax = 2.0*hdur+hwid
+!!$ double precision, parameter :: fmin = 1.0/tmax, fmax = 1.0/tmin
+!!$ double precision, parameter :: trbdndw = 0.3, a = 30.0
+!!$ integer, parameter :: passes = 2, iord = 4
+!!$
+!!$!
+!!$! MULTI-TAPER PARAMETERS
+!!$!
+!!$! Ying Zhou: The fit between the recovered data and the data can be improved
+!!$! by either increasing the window width (HWIN above) or by decreasing NPI.
+!!$! In her experience, NPI = 2.5 is good for noisy data.
+!!$! For synthetic data, we can use a lower NPI.
+!!$! number of tapers should be fixed as twice NPI -- see Latex notes
+!!$ double precision, parameter :: WTR = 0.02
+!!$ double precision, parameter :: NPI = 2.5
+!!$ integer, parameter :: NTAPER = int(2.0*NPI)
+!!$
+!!$! FFT parameters
+!!$! FOR SOME REASON, I CANNOT SET NDIM = npt when lnpt <= 11 -- WHY NOT?
+!!$ integer, parameter :: lnpt = 13, npt = 2**lnpt, NDIM = 20000
+!!$ !integer, parameter :: lnpt = 14, npt = 2**lnpt, NDIM = npt
+!!$
+!!$ !double precision, parameter :: ZZIGN = -1.0 ! sign convention -- WATCH OUT
+!!$ double precision, parameter :: FORWARD_FFT = 1.0
+!!$ double precision, parameter :: REVERSE_FFT = -1.0
+
+ integer, parameter :: NDIM = 20000
+
+end module wave2d_constants
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex02.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex02.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex02.f90 2009-12-02 18:44:58 UTC (rev 16060)
@@ -0,0 +1,336 @@
+module wave2d_constants
+
+! This file is copied into the output directory when wave2d.f90 is run
+! There are two basic options: membrane (surface) waves and body waves.
+! Several lines need to be commented/uncommented to go from one to the other.
+
+! index of the reference directory for the simulation output
+ integer, parameter :: IRUNZ = 200
+
+!========================
+! GRID, TIME-STEP, AND SOURCE PARAMETERS
+
+! NFRAME : number of frames to save: membrane, 10; body, 20
+! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
+! NSTEP : number of timesteps
+ integer, parameter :: NFRAME = 10 ! membrane
+ !integer, parameter :: NFRAME = 11 ! body
+ integer, parameter :: NSAVE = 400 ! 200,400
+ integer, parameter :: NSTEP = NFRAME*NSAVE
+
+! time step in seconds
+ !double precision, parameter :: DT = 2.0e-2 ! body waves
+ double precision, parameter :: DT = 6.0e-2 ! membrane surface waves
+
+! temporal properties of source (source time function)
+ integer, parameter :: ISRC_TIME = 1 ! type (1)
+ double precision, parameter :: hdur = 10.0 ! HALF-duration (s), membrane
+ !double precision, parameter :: hdur = 2.0 ! HALF-duration (s), body
+ double precision, parameter :: tshift = 2.0*DT*dble(NSAVE) ! time shift (s)
+ !double precision, parameter :: tshift = 8.0*hdur
+ logical, parameter :: SRC_TAPER = .true. ! taper the endpoints of the time series
+
+! normalization factor of point source force
+ double precision, parameter :: FNORM = 1.0d10
+
+! forward wavefield source-time function
+! (For membrane waves, FOR_Y = 1 is all that matters.)
+ integer, parameter :: FOR_X = 1 ! boolean
+ integer, parameter :: FOR_Y = 1 ! boolean
+ integer, parameter :: FOR_Z = 1 ! boolean
+
+! adjoint wavefield source-time function
+! (For membrane waves, REV_Y = 1 is all that matters.)
+! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
+ integer, parameter :: REV_X = 1 ! boolean
+ integer, parameter :: REV_Y = 1 ! boolean
+ integer, parameter :: REV_Z = 1 ! boolean
+
+! spatial properties of sources
+! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
+! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
+! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
+ integer, parameter :: ISRC_SPACE = 6 ! see wave2d.f90
+
+! spatial properties of receivers
+! IREC_SPACE
+! (1) individual station(s)
+! (2) SoCal (used for GJI paper)
+! (3) regular mesh on land
+! (4) regular mesh
+! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
+! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
+! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
+ integer, parameter :: IREC_SPACE = 2 ! see wave2d.f90
+ integer, parameter :: NMESH_REC = 17
+ double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m ! 4km for membrane surface waves
+ double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m ! 15km for membrane surface waves
+ double precision, parameter :: STATION_COAST_BUFFER = 0.0e3 ! m
+
+! lower right corner for membrane surface waves plotting grid
+ double precision, parameter :: LAT_MIN = 32.0
+ double precision, parameter :: LON_MIN = -120.0
+ integer, parameter :: UTM_PROJECTION_ZONE = 11 ! southern California
+
+! mesh specifications: membrane surface waves
+ double precision, parameter :: LENGTH = 480.0e3 ! m
+ double precision, parameter :: HEIGHT = 480.0e3 ! m
+ double precision, parameter :: AREA = LENGTH*HEIGHT
+ integer, parameter :: NEX = 40
+ integer, parameter :: NEZ = 40
+!!$
+!!$! mesh specifications: body waves
+!!$ double precision, parameter :: LENGTH = 200.0e3 ! m ! 400 for 1D body waves
+!!$ double precision, parameter :: HEIGHT = 80.0e3 ! m
+!!$ integer, parameter :: NEX = 80 ! 160
+!!$ integer, parameter :: NEZ = 32 ! 32
+
+!========================
+! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
+
+! model perturbations for HOMOGENEOUS model (or perturbation)
+! scaling from beta to alpha
+! value is from Master et al. (2000), "The relative behavior of shear velocity..."
+ double precision, parameter :: R_BETA_OVER_ALPHA = 1.3
+ double precision, parameter :: PBETA = 10.0
+ double precision, parameter :: PALPHA = PBETA / R_BETA_OVER_ALPHA
+ double precision, parameter :: PRHO = 0.0
+
+ ! reference model and target model choice
+ integer, parameter :: IMODEL_SYN = 0
+ integer, parameter :: IMODEL_DAT = 2
+ !-----------------------------------------------------------------------------------------
+ ! 0 1 2 3
+ !-----------------------------------------------------------------------------------------
+ ! ISURFACE=1, IMODEL_SYN : homo checker het
+ ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het pert
+ !-----------------------------------------------------------------------------------------
+ ! ISURFACE=0, IMODEL_SYN : homo 1D model checker het
+ ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert het pert
+ !-----------------------------------------------------------------------------------------
+
+ ! whether you want to smooth the structure model for computing the synthetics
+ integer, parameter :: ISMOOTH_INITIAL_MODEL = 0
+ integer, parameter :: ISMOOTH_MODEL_UPDATE = 0
+ integer, parameter :: ISMOOTH_EVENT_KERNEL = 1
+ integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
+
+ ! smoothing scalelength
+ double precision, parameter :: GAMMA_SMOOTH_MODEL = 30.0e3
+
+!========================
+
+! boolean parameters
+! IKER: (0) waveform
+! (1) traveltime, cross-correlation, misfit
+! (2) amplitude, cross-correlation, misfit
+! (3) traveltime, multitaper
+! (4) amplitude, multitaper
+! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
+! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
+ integer, parameter :: IKER = 1
+ integer, parameter :: IAMP_VEL = 0 ! measure amplitudes between velocity traces
+
+! KEY: USE ONE LINE OR THE OTHER
+ integer, parameter :: ISURFACE = 1, NCOMP = 1, NABSORB = 4 ! surface waves
+! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
+
+! parameters controlling what to write to file
+! NOTE: for the tomography simulations, ALL of these can be .false.
+
+ logical, parameter :: WRITE_STF_F = .false.
+ logical, parameter :: WRITE_SEISMO_F = .false.
+! logical, parameter :: WRITE_SPECTRA_F = .false.
+! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
+
+ logical, parameter :: WRITE_SEISMO_RECONSTRUCT = .false. ! multitaper
+! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
+
+ logical, parameter :: WRITE_STF_A = .false.
+ logical, parameter :: WRITE_SEISMO_A = .false. ! source inversions
+! logical, parameter :: WRITE_SPECTRA_A = .false.
+! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
+
+ logical, parameter :: WRITE_KERNELS = .true. ! write all nine kernels
+ logical, parameter :: WRITE_KERNEL_SNAPSHOTS = .false. ! kernel snapshots
+ logical, parameter :: WRITE_WAVFIELD_SNAPSHOTS = .false. ! wavefield snapshots
+
+!--------------------------------------
+! INVERSION PARAMETERS
+
+ ! whether you want to compute kernels, or simply the misfit function
+ logical, parameter :: COMPUTE_KERNELS = .true.
+
+ ! whether to use the data subspace method, which has a Hessian
+ logical, parameter :: HESSIAN = .false.
+
+ ! stopping criteria
+ ! NITERATION : number of iterations
+ ! VAR_RED_MIN : minimum variance reduction (in percent)
+ ! SIGMA_FAC : stop if a model value exceeds SIGMA_FAC * sigma_m
+ ! CONV_STOP : stop when the misfit value is this fraction of the INITIAL misfit value
+ integer, parameter :: NITERATION = 16
+ double precision, parameter :: VAR_RED_MIN = 8.0
+
+ ! Gaussian errors containted in input file
+ ! see wave2d_sigmas.m and INPUT/sigma_0p1_pert.dat
+ double precision, parameter :: SIGMA_DT = 0.10
+ double precision, parameter :: SIGMA_DLNA = 1.0
+ double precision, parameter :: SIGMA_WAVEFORM = 1.0
+ logical, parameter :: ADD_DATA_ERRORS = .true.
+
+ !double precision, parameter :: SIGMA_FAC = 2.0
+ !double precision, parameter :: CONV_STOP = 1.0e-4
+
+ ! order of interpolating polynomial in conjugate gradient algorithm
+ ! using POLY_ORDER = 3 required computing the gradient of the test model
+ integer, parameter :: POLY_ORDER = 2 ! 2 (preferred) or 3
+
+ ! use inversion and scaling used for the GJI paper
+ !integer, parameter :: GJI_PAPER = 1
+
+ ! what to perturb, what to invert
+ integer, parameter :: PERT_STRUCT = 0
+ integer, parameter :: PERT_SOURCE_T = 1
+ integer, parameter :: PERT_SOURCE_X = 1
+
+ integer, parameter :: INV_STRUCT = 0
+ integer, parameter :: INV_SOURCE_T = 1
+ integer, parameter :: INV_SOURCE_X = 1
+
+ ! whether to include the model norm term in the misfit function, which acts like damping
+ logical, parameter :: INCLUDE_MODEL_NORM = .true.
+
+ ! log file showing source loactions
+ logical, parameter :: ISOURCE_LOG = .true.
+
+ ! DO NOT CHANGE THESE
+ integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ integer, parameter :: NVAR_SOURCE = 3 ! x position, z position, origin time
+ integer, parameter :: NVAR = NVAR_STRUCT + NVAR_SOURCE
+
+! parameterization for structure in the inversion
+! 1 : kappa-mu-rho
+! 2 : alpha-beta-rho
+! 3 : c-beta-rho
+ integer, parameter :: STRUCTURE_PARAMETER_TYPE = 2
+
+! scalelength of smoothing Gaussian
+! GJI paper : 30,60,90
+! body waves : 10
+ double precision, parameter :: GAMMA_SMOOTH = 60.0e3
+
+ ! parameters for smoothing
+ logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
+ logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
+
+! homogeneous background model (S.I. units)
+ double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
+ double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
+ double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
+
+!---------------------------------------------------------------
+! measurement windows
+
+ double precision, parameter :: HWIN1 = 2.5*hdur ! half-window width for pulse (2.5)
+ double precision, parameter :: HWIN2 = 5.0*hdur ! half-window with for MTM measurement (zero-padding)
+
+!---------------------------------------------------------------
+! CHT: do not change these
+
+! UTM zone for Southern California region
+! integer, parameter :: UTM_PROJECTION_ZONE = 11
+
+! to suppress UTM projection for SCEC benchmarks
+ logical, parameter :: SUPPRESS_UTM_PROJECTION = .false.
+
+! flag for projection from latitude/longitude to UTM, and back
+ integer, parameter :: ILONGLAT2UTM = 0, IUTM2LONGLAT = 1
+
+! flag for projection from latitude/longitude to mesh-UTM, and back
+ integer, parameter :: ILONLAT2MESH = 0, IMESH2LONLAT = 1
+
+! max number of fake receivers
+ integer, parameter :: MAX_SR_FAKE = 1000
+
+! max number of events, receivers, and phass
+ integer, parameter :: MAX_EVENT = 50
+ integer, parameter :: MAX_SR = 1400
+ integer, parameter :: MAX_PHASE = 1
+ integer, parameter :: MAX_COMP = NCOMP
+
+!========================
+! GRID AND GLL POINTS
+
+ integer, parameter :: NELE = MAX(NEX,NEZ)
+ integer, parameter :: NSPEC = NEX*NEZ
+
+! number of GLL points (polynomial degree plus one)
+ integer, parameter :: NGLLX = 5
+ integer, parameter :: NGLLZ = 5
+ integer, parameter :: NGLL = MAX(NGLLX,NGLLZ)
+
+! number of global points
+ integer, parameter :: NGLLSQUARE = NGLLX * NGLLZ ! element GLL points
+ integer, parameter :: NGLOB = ((NGLLX-1)*NEX + 1)*((NGLLZ-1)*NEZ +1) ! global GLL points
+ integer, parameter :: NLOCAL = NGLLX * NGLLZ * NSPEC ! local GLL points
+ integer, parameter :: NSPEC_CORNER = (NEX+1) * (NEZ+1) ! element corner points
+
+! number of nodes for 2D and 3D shape functions for hexahedra
+! we use 8-node mesh bricks, which are more stable than 27-node elements
+ integer, parameter :: NGNOD = 8, NGNOD2D = 4
+
+! number of iterations to solve the system for xi and eta
+ integer, parameter :: NUM_ITER = 1
+
+! very large and very small values
+ double precision, parameter :: HUGEVAL = 1.0e30, TINYVAL = 1.0e-9
+
+! for the Gauss-Lobatto-Legendre points and weights
+ double precision, parameter :: GAUSSALPHA = 0.0,GAUSSBETA = 0.0
+
+!
+! CONSTANTS
+!
+ double precision, parameter :: PI = 3.141592653589793
+ double precision, parameter :: FOUR_THIRDS = 4.0 / 3.0
+ double precision, parameter :: ONE_THIRD = 1.0 / 3.0
+ double precision, parameter :: ONEOVERTWO = 0.5
+ double precision, parameter :: DEG = 180.0/PI
+! double precision, parameter :: EPS = 1.0e-35
+
+!!$! parameter for FFTW
+!!$ integer, parameter :: NOUT = NSTEP/2 + 1
+!!$
+!!$! bounds for bandpass filter
+!!$ double precision, parameter :: hwid = 3.0 ! HALF-width of window (s)
+!!$ double precision, parameter :: tmin = 2.0*hdur-hwid
+!!$ double precision, parameter :: tmax = 2.0*hdur+hwid
+!!$ double precision, parameter :: fmin = 1.0/tmax, fmax = 1.0/tmin
+!!$ double precision, parameter :: trbdndw = 0.3, a = 30.0
+!!$ integer, parameter :: passes = 2, iord = 4
+!!$
+!!$!
+!!$! MULTI-TAPER PARAMETERS
+!!$!
+!!$! Ying Zhou: The fit between the recovered data and the data can be improved
+!!$! by either increasing the window width (HWIN above) or by decreasing NPI.
+!!$! In her experience, NPI = 2.5 is good for noisy data.
+!!$! For synthetic data, we can use a lower NPI.
+!!$! number of tapers should be fixed as twice NPI -- see Latex notes
+!!$ double precision, parameter :: WTR = 0.02
+!!$ double precision, parameter :: NPI = 2.5
+!!$ integer, parameter :: NTAPER = int(2.0*NPI)
+!!$
+!!$! FFT parameters
+!!$! FOR SOME REASON, I CANNOT SET NDIM = npt when lnpt <= 11 -- WHY NOT?
+!!$ integer, parameter :: lnpt = 13, npt = 2**lnpt, NDIM = 20000
+!!$ !integer, parameter :: lnpt = 14, npt = 2**lnpt, NDIM = npt
+!!$
+!!$ !double precision, parameter :: ZZIGN = -1.0 ! sign convention -- WATCH OUT
+!!$ double precision, parameter :: FORWARD_FFT = 1.0
+!!$ double precision, parameter :: REVERSE_FFT = -1.0
+
+ integer, parameter :: NDIM = 20000
+
+end module wave2d_constants
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex03.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex03.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex03.f90 2009-12-02 18:44:58 UTC (rev 16060)
@@ -0,0 +1,336 @@
+module wave2d_constants
+
+! This file is copied into the output directory when wave2d.f90 is run
+! There are two basic options: membrane (surface) waves and body waves.
+! Several lines need to be commented/uncommented to go from one to the other.
+
+! index of the reference directory for the simulation output
+ integer, parameter :: IRUNZ = 300
+
+!========================
+! GRID, TIME-STEP, AND SOURCE PARAMETERS
+
+! NFRAME : number of frames to save: membrane, 10; body, 20
+! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
+! NSTEP : number of timesteps
+ integer, parameter :: NFRAME = 10 ! membrane
+ !integer, parameter :: NFRAME = 11 ! body
+ integer, parameter :: NSAVE = 400 ! 200,400
+ integer, parameter :: NSTEP = NFRAME*NSAVE
+
+! time step in seconds
+ !double precision, parameter :: DT = 2.0e-2 ! body waves
+ double precision, parameter :: DT = 6.0e-2 ! membrane surface waves
+
+! temporal properties of source (source time function)
+ integer, parameter :: ISRC_TIME = 1 ! type (1)
+ double precision, parameter :: hdur = 10.0 ! HALF-duration (s), membrane
+ !double precision, parameter :: hdur = 2.0 ! HALF-duration (s), body
+ double precision, parameter :: tshift = 2.0*DT*dble(NSAVE) ! time shift (s)
+ !double precision, parameter :: tshift = 8.0*hdur
+ logical, parameter :: SRC_TAPER = .true. ! taper the endpoints of the time series
+
+! normalization factor of point source force
+ double precision, parameter :: FNORM = 1.0d10
+
+! forward wavefield source-time function
+! (For membrane waves, FOR_Y = 1 is all that matters.)
+ integer, parameter :: FOR_X = 1 ! boolean
+ integer, parameter :: FOR_Y = 1 ! boolean
+ integer, parameter :: FOR_Z = 1 ! boolean
+
+! adjoint wavefield source-time function
+! (For membrane waves, REV_Y = 1 is all that matters.)
+! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
+ integer, parameter :: REV_X = 1 ! boolean
+ integer, parameter :: REV_Y = 1 ! boolean
+ integer, parameter :: REV_Z = 1 ! boolean
+
+! spatial properties of sources
+! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
+! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
+! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
+ integer, parameter :: ISRC_SPACE = 6 ! see wave2d.f90
+
+! spatial properties of receivers
+! IREC_SPACE
+! (1) individual station(s)
+! (2) SoCal (used for GJI paper)
+! (3) regular mesh on land
+! (4) regular mesh
+! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
+! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
+! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
+ integer, parameter :: IREC_SPACE = 2 ! see wave2d.f90
+ integer, parameter :: NMESH_REC = 17
+ double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m ! 4km for membrane surface waves
+ double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m ! 15km for membrane surface waves
+ double precision, parameter :: STATION_COAST_BUFFER = 0.0e3 ! m
+
+! lower right corner for membrane surface waves plotting grid
+ double precision, parameter :: LAT_MIN = 32.0
+ double precision, parameter :: LON_MIN = -120.0
+ integer, parameter :: UTM_PROJECTION_ZONE = 11 ! southern California
+
+! mesh specifications: membrane surface waves
+ double precision, parameter :: LENGTH = 480.0e3 ! m
+ double precision, parameter :: HEIGHT = 480.0e3 ! m
+ double precision, parameter :: AREA = LENGTH*HEIGHT
+ integer, parameter :: NEX = 40
+ integer, parameter :: NEZ = 40
+!!$
+!!$! mesh specifications: body waves
+!!$ double precision, parameter :: LENGTH = 200.0e3 ! m ! 400 for 1D body waves
+!!$ double precision, parameter :: HEIGHT = 80.0e3 ! m
+!!$ integer, parameter :: NEX = 80 ! 160
+!!$ integer, parameter :: NEZ = 32 ! 32
+
+!========================
+! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
+
+! model perturbations for HOMOGENEOUS model (or perturbation)
+! scaling from beta to alpha
+! value is from Master et al. (2000), "The relative behavior of shear velocity..."
+ double precision, parameter :: R_BETA_OVER_ALPHA = 1.3
+ double precision, parameter :: PBETA = 10.0
+ double precision, parameter :: PALPHA = PBETA / R_BETA_OVER_ALPHA
+ double precision, parameter :: PRHO = 0.0
+
+ ! reference model and target model choice
+ integer, parameter :: IMODEL_SYN = 0
+ integer, parameter :: IMODEL_DAT = 2
+ !-----------------------------------------------------------------------------------------
+ ! 0 1 2 3
+ !-----------------------------------------------------------------------------------------
+ ! ISURFACE=1, IMODEL_SYN : homo checker het
+ ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het pert
+ !-----------------------------------------------------------------------------------------
+ ! ISURFACE=0, IMODEL_SYN : homo 1D model checker het
+ ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert het pert
+ !-----------------------------------------------------------------------------------------
+
+ ! whether you want to smooth the structure model for computing the synthetics
+ integer, parameter :: ISMOOTH_INITIAL_MODEL = 0
+ integer, parameter :: ISMOOTH_MODEL_UPDATE = 0
+ integer, parameter :: ISMOOTH_EVENT_KERNEL = 1
+ integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
+
+ ! smoothing scalelength
+ double precision, parameter :: GAMMA_SMOOTH_MODEL = 30.0e3
+
+!========================
+
+! boolean parameters
+! IKER: (0) waveform
+! (1) traveltime, cross-correlation, misfit
+! (2) amplitude, cross-correlation, misfit
+! (3) traveltime, multitaper
+! (4) amplitude, multitaper
+! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
+! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
+ integer, parameter :: IKER = 1
+ integer, parameter :: IAMP_VEL = 0 ! measure amplitudes between velocity traces
+
+! KEY: USE ONE LINE OR THE OTHER
+ integer, parameter :: ISURFACE = 1, NCOMP = 1, NABSORB = 4 ! surface waves
+! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
+
+! parameters controlling what to write to file
+! NOTE: for the tomography simulations, ALL of these can be .false.
+
+ logical, parameter :: WRITE_STF_F = .false.
+ logical, parameter :: WRITE_SEISMO_F = .false.
+! logical, parameter :: WRITE_SPECTRA_F = .false.
+! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
+
+ logical, parameter :: WRITE_SEISMO_RECONSTRUCT = .false. ! multitaper
+! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
+
+ logical, parameter :: WRITE_STF_A = .false.
+ logical, parameter :: WRITE_SEISMO_A = .false. ! source inversions
+! logical, parameter :: WRITE_SPECTRA_A = .false.
+! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
+
+ logical, parameter :: WRITE_KERNELS = .true. ! write all nine kernels
+ logical, parameter :: WRITE_KERNEL_SNAPSHOTS = .false. ! kernel snapshots
+ logical, parameter :: WRITE_WAVFIELD_SNAPSHOTS = .false. ! wavefield snapshots
+
+!--------------------------------------
+! INVERSION PARAMETERS
+
+ ! whether you want to compute kernels, or simply the misfit function
+ logical, parameter :: COMPUTE_KERNELS = .true.
+
+ ! whether to use the data subspace method, which has a Hessian
+ logical, parameter :: HESSIAN = .false.
+
+ ! stopping criteria
+ ! NITERATION : number of iterations
+ ! VAR_RED_MIN : minimum variance reduction (in percent)
+ ! SIGMA_FAC : stop if a model value exceeds SIGMA_FAC * sigma_m
+ ! CONV_STOP : stop when the misfit value is this fraction of the INITIAL misfit value
+ integer, parameter :: NITERATION = 16
+ double precision, parameter :: VAR_RED_MIN = 8.0
+
+ ! Gaussian errors containted in input file
+ ! see wave2d_sigmas.m and INPUT/sigma_0p1_pert.dat
+ double precision, parameter :: SIGMA_DT = 0.10
+ double precision, parameter :: SIGMA_DLNA = 1.0
+ double precision, parameter :: SIGMA_WAVEFORM = 1.0
+ logical, parameter :: ADD_DATA_ERRORS = .true.
+
+ !double precision, parameter :: SIGMA_FAC = 2.0
+ !double precision, parameter :: CONV_STOP = 1.0e-4
+
+ ! order of interpolating polynomial in conjugate gradient algorithm
+ ! using POLY_ORDER = 3 required computing the gradient of the test model
+ integer, parameter :: POLY_ORDER = 2 ! 2 (preferred) or 3
+
+ ! use inversion and scaling used for the GJI paper
+ !integer, parameter :: GJI_PAPER = 1
+
+ ! what to perturb, what to invert
+ integer, parameter :: PERT_STRUCT = 1
+ integer, parameter :: PERT_SOURCE_T = 1
+ integer, parameter :: PERT_SOURCE_X = 1
+
+ integer, parameter :: INV_STRUCT = 1
+ integer, parameter :: INV_SOURCE_T = 1
+ integer, parameter :: INV_SOURCE_X = 1
+
+ ! whether to include the model norm term in the misfit function, which acts like damping
+ logical, parameter :: INCLUDE_MODEL_NORM = .true.
+
+ ! log file showing source loactions
+ logical, parameter :: ISOURCE_LOG = .true.
+
+ ! DO NOT CHANGE THESE
+ integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ integer, parameter :: NVAR_SOURCE = 3 ! x position, z position, origin time
+ integer, parameter :: NVAR = NVAR_STRUCT + NVAR_SOURCE
+
+! parameterization for structure in the inversion
+! 1 : kappa-mu-rho
+! 2 : alpha-beta-rho
+! 3 : c-beta-rho
+ integer, parameter :: STRUCTURE_PARAMETER_TYPE = 2
+
+! scalelength of smoothing Gaussian
+! GJI paper : 30,60,90
+! body waves : 10
+ double precision, parameter :: GAMMA_SMOOTH = 60.0e3
+
+ ! parameters for smoothing
+ logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
+ logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
+
+! homogeneous background model (S.I. units)
+ double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
+ double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
+ double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
+
+!---------------------------------------------------------------
+! measurement windows
+
+ double precision, parameter :: HWIN1 = 2.5*hdur ! half-window width for pulse (2.5)
+ double precision, parameter :: HWIN2 = 5.0*hdur ! half-window with for MTM measurement (zero-padding)
+
+!---------------------------------------------------------------
+! CHT: do not change these
+
+! UTM zone for Southern California region
+! integer, parameter :: UTM_PROJECTION_ZONE = 11
+
+! to suppress UTM projection for SCEC benchmarks
+ logical, parameter :: SUPPRESS_UTM_PROJECTION = .false.
+
+! flag for projection from latitude/longitude to UTM, and back
+ integer, parameter :: ILONGLAT2UTM = 0, IUTM2LONGLAT = 1
+
+! flag for projection from latitude/longitude to mesh-UTM, and back
+ integer, parameter :: ILONLAT2MESH = 0, IMESH2LONLAT = 1
+
+! max number of fake receivers
+ integer, parameter :: MAX_SR_FAKE = 1000
+
+! max number of events, receivers, and phass
+ integer, parameter :: MAX_EVENT = 50
+ integer, parameter :: MAX_SR = 1400
+ integer, parameter :: MAX_PHASE = 1
+ integer, parameter :: MAX_COMP = NCOMP
+
+!========================
+! GRID AND GLL POINTS
+
+ integer, parameter :: NELE = MAX(NEX,NEZ)
+ integer, parameter :: NSPEC = NEX*NEZ
+
+! number of GLL points (polynomial degree plus one)
+ integer, parameter :: NGLLX = 5
+ integer, parameter :: NGLLZ = 5
+ integer, parameter :: NGLL = MAX(NGLLX,NGLLZ)
+
+! number of global points
+ integer, parameter :: NGLLSQUARE = NGLLX * NGLLZ ! element GLL points
+ integer, parameter :: NGLOB = ((NGLLX-1)*NEX + 1)*((NGLLZ-1)*NEZ +1) ! global GLL points
+ integer, parameter :: NLOCAL = NGLLX * NGLLZ * NSPEC ! local GLL points
+ integer, parameter :: NSPEC_CORNER = (NEX+1) * (NEZ+1) ! element corner points
+
+! number of nodes for 2D and 3D shape functions for hexahedra
+! we use 8-node mesh bricks, which are more stable than 27-node elements
+ integer, parameter :: NGNOD = 8, NGNOD2D = 4
+
+! number of iterations to solve the system for xi and eta
+ integer, parameter :: NUM_ITER = 1
+
+! very large and very small values
+ double precision, parameter :: HUGEVAL = 1.0e30, TINYVAL = 1.0e-9
+
+! for the Gauss-Lobatto-Legendre points and weights
+ double precision, parameter :: GAUSSALPHA = 0.0,GAUSSBETA = 0.0
+
+!
+! CONSTANTS
+!
+ double precision, parameter :: PI = 3.141592653589793
+ double precision, parameter :: FOUR_THIRDS = 4.0 / 3.0
+ double precision, parameter :: ONE_THIRD = 1.0 / 3.0
+ double precision, parameter :: ONEOVERTWO = 0.5
+ double precision, parameter :: DEG = 180.0/PI
+! double precision, parameter :: EPS = 1.0e-35
+
+!!$! parameter for FFTW
+!!$ integer, parameter :: NOUT = NSTEP/2 + 1
+!!$
+!!$! bounds for bandpass filter
+!!$ double precision, parameter :: hwid = 3.0 ! HALF-width of window (s)
+!!$ double precision, parameter :: tmin = 2.0*hdur-hwid
+!!$ double precision, parameter :: tmax = 2.0*hdur+hwid
+!!$ double precision, parameter :: fmin = 1.0/tmax, fmax = 1.0/tmin
+!!$ double precision, parameter :: trbdndw = 0.3, a = 30.0
+!!$ integer, parameter :: passes = 2, iord = 4
+!!$
+!!$!
+!!$! MULTI-TAPER PARAMETERS
+!!$!
+!!$! Ying Zhou: The fit between the recovered data and the data can be improved
+!!$! by either increasing the window width (HWIN above) or by decreasing NPI.
+!!$! In her experience, NPI = 2.5 is good for noisy data.
+!!$! For synthetic data, we can use a lower NPI.
+!!$! number of tapers should be fixed as twice NPI -- see Latex notes
+!!$ double precision, parameter :: WTR = 0.02
+!!$ double precision, parameter :: NPI = 2.5
+!!$ integer, parameter :: NTAPER = int(2.0*NPI)
+!!$
+!!$! FFT parameters
+!!$! FOR SOME REASON, I CANNOT SET NDIM = npt when lnpt <= 11 -- WHY NOT?
+!!$ integer, parameter :: lnpt = 13, npt = 2**lnpt, NDIM = 20000
+!!$ !integer, parameter :: lnpt = 14, npt = 2**lnpt, NDIM = npt
+!!$
+!!$ !double precision, parameter :: ZZIGN = -1.0 ! sign convention -- WATCH OUT
+!!$ double precision, parameter :: FORWARD_FFT = 1.0
+!!$ double precision, parameter :: REVERSE_FFT = -1.0
+
+ integer, parameter :: NDIM = 20000
+
+end module wave2d_constants
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex04.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex04.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex04.f90 2009-12-02 18:44:58 UTC (rev 16060)
@@ -0,0 +1,336 @@
+module wave2d_constants
+
+! This file is copied into the output directory when wave2d.f90 is run
+! There are two basic options: membrane (surface) waves and body waves.
+! Several lines need to be commented/uncommented to go from one to the other.
+
+! index of the reference directory for the simulation output
+ integer, parameter :: IRUNZ = 400
+
+!========================
+! GRID, TIME-STEP, AND SOURCE PARAMETERS
+
+! NFRAME : number of frames to save: membrane, 10; body, 20
+! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
+! NSTEP : number of timesteps
+ integer, parameter :: NFRAME = 10 ! membrane
+ !integer, parameter :: NFRAME = 11 ! body
+ integer, parameter :: NSAVE = 400 ! 200,400
+ integer, parameter :: NSTEP = NFRAME*NSAVE
+
+! time step in seconds
+ !double precision, parameter :: DT = 2.0e-2 ! body waves
+ double precision, parameter :: DT = 6.0e-2 ! membrane surface waves
+
+! temporal properties of source (source time function)
+ integer, parameter :: ISRC_TIME = 1 ! type (1)
+ double precision, parameter :: hdur = 10.0 ! HALF-duration (s), membrane
+ !double precision, parameter :: hdur = 2.0 ! HALF-duration (s), body
+ double precision, parameter :: tshift = 2.0*DT*dble(NSAVE) ! time shift (s)
+ !double precision, parameter :: tshift = 8.0*hdur
+ logical, parameter :: SRC_TAPER = .true. ! taper the endpoints of the time series
+
+! normalization factor of point source force
+ double precision, parameter :: FNORM = 1.0d10
+
+! forward wavefield source-time function
+! (For membrane waves, FOR_Y = 1 is all that matters.)
+ integer, parameter :: FOR_X = 1 ! boolean
+ integer, parameter :: FOR_Y = 1 ! boolean
+ integer, parameter :: FOR_Z = 1 ! boolean
+
+! adjoint wavefield source-time function
+! (For membrane waves, REV_Y = 1 is all that matters.)
+! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
+ integer, parameter :: REV_X = 1 ! boolean
+ integer, parameter :: REV_Y = 1 ! boolean
+ integer, parameter :: REV_Z = 1 ! boolean
+
+! spatial properties of sources
+! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
+! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
+! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
+ integer, parameter :: ISRC_SPACE = 6 ! see wave2d.f90
+
+! spatial properties of receivers
+! IREC_SPACE
+! (1) individual station(s)
+! (2) SoCal (used for GJI paper)
+! (3) regular mesh on land
+! (4) regular mesh
+! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
+! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
+! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
+ integer, parameter :: IREC_SPACE = 2 ! see wave2d.f90
+ integer, parameter :: NMESH_REC = 17
+ double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m ! 4km for membrane surface waves
+ double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m ! 15km for membrane surface waves
+ double precision, parameter :: STATION_COAST_BUFFER = 0.0e3 ! m
+
+! lower right corner for membrane surface waves plotting grid
+ double precision, parameter :: LAT_MIN = 32.0
+ double precision, parameter :: LON_MIN = -120.0
+ integer, parameter :: UTM_PROJECTION_ZONE = 11 ! southern California
+
+! mesh specifications: membrane surface waves
+ double precision, parameter :: LENGTH = 480.0e3 ! m
+ double precision, parameter :: HEIGHT = 480.0e3 ! m
+ double precision, parameter :: AREA = LENGTH*HEIGHT
+ integer, parameter :: NEX = 40
+ integer, parameter :: NEZ = 40
+!!$
+!!$! mesh specifications: body waves
+!!$ double precision, parameter :: LENGTH = 200.0e3 ! m ! 400 for 1D body waves
+!!$ double precision, parameter :: HEIGHT = 80.0e3 ! m
+!!$ integer, parameter :: NEX = 80 ! 160
+!!$ integer, parameter :: NEZ = 32 ! 32
+
+!========================
+! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
+
+! model perturbations for HOMOGENEOUS model (or perturbation)
+! scaling from beta to alpha
+! value is from Master et al. (2000), "The relative behavior of shear velocity..."
+ double precision, parameter :: R_BETA_OVER_ALPHA = 1.3
+ double precision, parameter :: PBETA = 10.0
+ double precision, parameter :: PALPHA = PBETA / R_BETA_OVER_ALPHA
+ double precision, parameter :: PRHO = 0.0
+
+ ! reference model and target model choice
+ integer, parameter :: IMODEL_SYN = 0
+ integer, parameter :: IMODEL_DAT = 2
+ !-----------------------------------------------------------------------------------------
+ ! 0 1 2 3
+ !-----------------------------------------------------------------------------------------
+ ! ISURFACE=1, IMODEL_SYN : homo checker het
+ ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het pert
+ !-----------------------------------------------------------------------------------------
+ ! ISURFACE=0, IMODEL_SYN : homo 1D model checker het
+ ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert het pert
+ !-----------------------------------------------------------------------------------------
+
+ ! whether you want to smooth the structure model for computing the synthetics
+ integer, parameter :: ISMOOTH_INITIAL_MODEL = 0
+ integer, parameter :: ISMOOTH_MODEL_UPDATE = 0
+ integer, parameter :: ISMOOTH_EVENT_KERNEL = 1
+ integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
+
+ ! smoothing scalelength
+ double precision, parameter :: GAMMA_SMOOTH_MODEL = 30.0e3
+
+!========================
+
+! boolean parameters
+! IKER: (0) waveform
+! (1) traveltime, cross-correlation, misfit
+! (2) amplitude, cross-correlation, misfit
+! (3) traveltime, multitaper
+! (4) amplitude, multitaper
+! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
+! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
+ integer, parameter :: IKER = 1
+ integer, parameter :: IAMP_VEL = 0 ! measure amplitudes between velocity traces
+
+! KEY: USE ONE LINE OR THE OTHER
+ integer, parameter :: ISURFACE = 1, NCOMP = 1, NABSORB = 4 ! surface waves
+! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
+
+! parameters controlling what to write to file
+! NOTE: for the tomography simulations, ALL of these can be .false.
+
+ logical, parameter :: WRITE_STF_F = .false.
+ logical, parameter :: WRITE_SEISMO_F = .false.
+! logical, parameter :: WRITE_SPECTRA_F = .false.
+! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
+
+ logical, parameter :: WRITE_SEISMO_RECONSTRUCT = .false. ! multitaper
+! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
+
+ logical, parameter :: WRITE_STF_A = .false.
+ logical, parameter :: WRITE_SEISMO_A = .false. ! source inversions
+! logical, parameter :: WRITE_SPECTRA_A = .false.
+! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
+
+ logical, parameter :: WRITE_KERNELS = .true. ! write all nine kernels
+ logical, parameter :: WRITE_KERNEL_SNAPSHOTS = .false. ! kernel snapshots
+ logical, parameter :: WRITE_WAVFIELD_SNAPSHOTS = .false. ! wavefield snapshots
+
+!--------------------------------------
+! INVERSION PARAMETERS
+
+ ! whether you want to compute kernels, or simply the misfit function
+ logical, parameter :: COMPUTE_KERNELS = .true.
+
+ ! whether to use the data subspace method, which has a Hessian
+ logical, parameter :: HESSIAN = .false.
+
+ ! stopping criteria
+ ! NITERATION : number of iterations
+ ! VAR_RED_MIN : minimum variance reduction (in percent)
+ ! SIGMA_FAC : stop if a model value exceeds SIGMA_FAC * sigma_m
+ ! CONV_STOP : stop when the misfit value is this fraction of the INITIAL misfit value
+ integer, parameter :: NITERATION = 16
+ double precision, parameter :: VAR_RED_MIN = 8.0
+
+ ! Gaussian errors containted in input file
+ ! see wave2d_sigmas.m and INPUT/sigma_0p1_pert.dat
+ double precision, parameter :: SIGMA_DT = 0.10
+ double precision, parameter :: SIGMA_DLNA = 1.0
+ double precision, parameter :: SIGMA_WAVEFORM = 1.0
+ logical, parameter :: ADD_DATA_ERRORS = .true.
+
+ !double precision, parameter :: SIGMA_FAC = 2.0
+ !double precision, parameter :: CONV_STOP = 1.0e-4
+
+ ! order of interpolating polynomial in conjugate gradient algorithm
+ ! using POLY_ORDER = 3 required computing the gradient of the test model
+ integer, parameter :: POLY_ORDER = 2 ! 2 (preferred) or 3
+
+ ! use inversion and scaling used for the GJI paper
+ !integer, parameter :: GJI_PAPER = 1
+
+ ! what to perturb, what to invert
+ integer, parameter :: PERT_STRUCT = 1
+ integer, parameter :: PERT_SOURCE_T = 1
+ integer, parameter :: PERT_SOURCE_X = 1
+
+ integer, parameter :: INV_STRUCT = 1
+ integer, parameter :: INV_SOURCE_T = 1
+ integer, parameter :: INV_SOURCE_X = 1
+
+ ! whether to include the model norm term in the misfit function, which acts like damping
+ logical, parameter :: INCLUDE_MODEL_NORM = .true.
+
+ ! log file showing source loactions
+ logical, parameter :: ISOURCE_LOG = .true.
+
+ ! DO NOT CHANGE THESE
+ integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ integer, parameter :: NVAR_SOURCE = 3 ! x position, z position, origin time
+ integer, parameter :: NVAR = NVAR_STRUCT + NVAR_SOURCE
+
+! parameterization for structure in the inversion
+! 1 : kappa-mu-rho
+! 2 : alpha-beta-rho
+! 3 : c-beta-rho
+ integer, parameter :: STRUCTURE_PARAMETER_TYPE = 2
+
+! scalelength of smoothing Gaussian
+! GJI paper : 30,60,90
+! body waves : 10
+ double precision, parameter :: GAMMA_SMOOTH = 60.0e3
+
+ ! parameters for smoothing
+ logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
+ logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
+
+! homogeneous background model (S.I. units)
+ double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
+ double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
+ double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
+
+!---------------------------------------------------------------
+! measurement windows
+
+ double precision, parameter :: HWIN1 = 2.5*hdur ! half-window width for pulse (2.5)
+ double precision, parameter :: HWIN2 = 5.0*hdur ! half-window with for MTM measurement (zero-padding)
+
+!---------------------------------------------------------------
+! CHT: do not change these
+
+! UTM zone for Southern California region
+! integer, parameter :: UTM_PROJECTION_ZONE = 11
+
+! to suppress UTM projection for SCEC benchmarks
+ logical, parameter :: SUPPRESS_UTM_PROJECTION = .false.
+
+! flag for projection from latitude/longitude to UTM, and back
+ integer, parameter :: ILONGLAT2UTM = 0, IUTM2LONGLAT = 1
+
+! flag for projection from latitude/longitude to mesh-UTM, and back
+ integer, parameter :: ILONLAT2MESH = 0, IMESH2LONLAT = 1
+
+! max number of fake receivers
+ integer, parameter :: MAX_SR_FAKE = 1000
+
+! max number of events, receivers, and phass
+ integer, parameter :: MAX_EVENT = 50
+ integer, parameter :: MAX_SR = 1400
+ integer, parameter :: MAX_PHASE = 1
+ integer, parameter :: MAX_COMP = NCOMP
+
+!========================
+! GRID AND GLL POINTS
+
+ integer, parameter :: NELE = MAX(NEX,NEZ)
+ integer, parameter :: NSPEC = NEX*NEZ
+
+! number of GLL points (polynomial degree plus one)
+ integer, parameter :: NGLLX = 5
+ integer, parameter :: NGLLZ = 5
+ integer, parameter :: NGLL = MAX(NGLLX,NGLLZ)
+
+! number of global points
+ integer, parameter :: NGLLSQUARE = NGLLX * NGLLZ ! element GLL points
+ integer, parameter :: NGLOB = ((NGLLX-1)*NEX + 1)*((NGLLZ-1)*NEZ +1) ! global GLL points
+ integer, parameter :: NLOCAL = NGLLX * NGLLZ * NSPEC ! local GLL points
+ integer, parameter :: NSPEC_CORNER = (NEX+1) * (NEZ+1) ! element corner points
+
+! number of nodes for 2D and 3D shape functions for hexahedra
+! we use 8-node mesh bricks, which are more stable than 27-node elements
+ integer, parameter :: NGNOD = 8, NGNOD2D = 4
+
+! number of iterations to solve the system for xi and eta
+ integer, parameter :: NUM_ITER = 1
+
+! very large and very small values
+ double precision, parameter :: HUGEVAL = 1.0e30, TINYVAL = 1.0e-9
+
+! for the Gauss-Lobatto-Legendre points and weights
+ double precision, parameter :: GAUSSALPHA = 0.0,GAUSSBETA = 0.0
+
+!
+! CONSTANTS
+!
+ double precision, parameter :: PI = 3.141592653589793
+ double precision, parameter :: FOUR_THIRDS = 4.0 / 3.0
+ double precision, parameter :: ONE_THIRD = 1.0 / 3.0
+ double precision, parameter :: ONEOVERTWO = 0.5
+ double precision, parameter :: DEG = 180.0/PI
+! double precision, parameter :: EPS = 1.0e-35
+
+!!$! parameter for FFTW
+!!$ integer, parameter :: NOUT = NSTEP/2 + 1
+!!$
+!!$! bounds for bandpass filter
+!!$ double precision, parameter :: hwid = 3.0 ! HALF-width of window (s)
+!!$ double precision, parameter :: tmin = 2.0*hdur-hwid
+!!$ double precision, parameter :: tmax = 2.0*hdur+hwid
+!!$ double precision, parameter :: fmin = 1.0/tmax, fmax = 1.0/tmin
+!!$ double precision, parameter :: trbdndw = 0.3, a = 30.0
+!!$ integer, parameter :: passes = 2, iord = 4
+!!$
+!!$!
+!!$! MULTI-TAPER PARAMETERS
+!!$!
+!!$! Ying Zhou: The fit between the recovered data and the data can be improved
+!!$! by either increasing the window width (HWIN above) or by decreasing NPI.
+!!$! In her experience, NPI = 2.5 is good for noisy data.
+!!$! For synthetic data, we can use a lower NPI.
+!!$! number of tapers should be fixed as twice NPI -- see Latex notes
+!!$ double precision, parameter :: WTR = 0.02
+!!$ double precision, parameter :: NPI = 2.5
+!!$ integer, parameter :: NTAPER = int(2.0*NPI)
+!!$
+!!$! FFT parameters
+!!$! FOR SOME REASON, I CANNOT SET NDIM = npt when lnpt <= 11 -- WHY NOT?
+!!$ integer, parameter :: lnpt = 13, npt = 2**lnpt, NDIM = 20000
+!!$ !integer, parameter :: lnpt = 14, npt = 2**lnpt, NDIM = npt
+!!$
+!!$ !double precision, parameter :: ZZIGN = -1.0 ! sign convention -- WATCH OUT
+!!$ double precision, parameter :: FORWARD_FFT = 1.0
+!!$ double precision, parameter :: REVERSE_FFT = -1.0
+
+ integer, parameter :: NDIM = 20000
+
+end module wave2d_constants
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex05.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex05.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex05.f90 2009-12-02 18:44:58 UTC (rev 16060)
@@ -0,0 +1,336 @@
+module wave2d_constants
+
+! This file is copied into the output directory when wave2d.f90 is run
+! There are two basic options: membrane (surface) waves and body waves.
+! Several lines need to be commented/uncommented to go from one to the other.
+
+! index of the reference directory for the simulation output
+ integer, parameter :: IRUNZ = 500
+
+!========================
+! GRID, TIME-STEP, AND SOURCE PARAMETERS
+
+! NFRAME : number of frames to save: membrane, 10; body, 20
+! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
+! NSTEP : number of timesteps
+ integer, parameter :: NFRAME = 10 ! membrane
+ !integer, parameter :: NFRAME = 11 ! body
+ integer, parameter :: NSAVE = 400 ! 200,400
+ integer, parameter :: NSTEP = NFRAME*NSAVE
+
+! time step in seconds
+ !double precision, parameter :: DT = 2.0e-2 ! body waves
+ double precision, parameter :: DT = 6.0e-2 ! membrane surface waves
+
+! temporal properties of source (source time function)
+ integer, parameter :: ISRC_TIME = 1 ! type (1)
+ double precision, parameter :: hdur = 10.0 ! HALF-duration (s), membrane
+ !double precision, parameter :: hdur = 2.0 ! HALF-duration (s), body
+ double precision, parameter :: tshift = 2.0*DT*dble(NSAVE) ! time shift (s)
+ !double precision, parameter :: tshift = 8.0*hdur
+ logical, parameter :: SRC_TAPER = .true. ! taper the endpoints of the time series
+
+! normalization factor of point source force
+ double precision, parameter :: FNORM = 1.0d10
+
+! forward wavefield source-time function
+! (For membrane waves, FOR_Y = 1 is all that matters.)
+ integer, parameter :: FOR_X = 1 ! boolean
+ integer, parameter :: FOR_Y = 1 ! boolean
+ integer, parameter :: FOR_Z = 1 ! boolean
+
+! adjoint wavefield source-time function
+! (For membrane waves, REV_Y = 1 is all that matters.)
+! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
+ integer, parameter :: REV_X = 1 ! boolean
+ integer, parameter :: REV_Y = 1 ! boolean
+ integer, parameter :: REV_Z = 1 ! boolean
+
+! spatial properties of sources
+! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
+! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
+! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
+ integer, parameter :: ISRC_SPACE = 6 ! see wave2d.f90
+
+! spatial properties of receivers
+! IREC_SPACE
+! (1) individual station(s)
+! (2) SoCal (used for GJI paper)
+! (3) regular mesh on land
+! (4) regular mesh
+! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
+! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
+! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
+ integer, parameter :: IREC_SPACE = 2 ! see wave2d.f90
+ integer, parameter :: NMESH_REC = 17
+ double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m ! 4km for membrane surface waves
+ double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m ! 15km for membrane surface waves
+ double precision, parameter :: STATION_COAST_BUFFER = 0.0e3 ! m
+
+! lower right corner for membrane surface waves plotting grid
+ double precision, parameter :: LAT_MIN = 32.0
+ double precision, parameter :: LON_MIN = -120.0
+ integer, parameter :: UTM_PROJECTION_ZONE = 11 ! southern California
+
+! mesh specifications: membrane surface waves
+ double precision, parameter :: LENGTH = 480.0e3 ! m
+ double precision, parameter :: HEIGHT = 480.0e3 ! m
+ double precision, parameter :: AREA = LENGTH*HEIGHT
+ integer, parameter :: NEX = 40
+ integer, parameter :: NEZ = 40
+!!$
+!!$! mesh specifications: body waves
+!!$ double precision, parameter :: LENGTH = 200.0e3 ! m ! 400 for 1D body waves
+!!$ double precision, parameter :: HEIGHT = 80.0e3 ! m
+!!$ integer, parameter :: NEX = 80 ! 160
+!!$ integer, parameter :: NEZ = 32 ! 32
+
+!========================
+! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
+
+! model perturbations for HOMOGENEOUS model (or perturbation)
+! scaling from beta to alpha
+! value is from Master et al. (2000), "The relative behavior of shear velocity..."
+ double precision, parameter :: R_BETA_OVER_ALPHA = 1.3
+ double precision, parameter :: PBETA = 10.0
+ double precision, parameter :: PALPHA = PBETA / R_BETA_OVER_ALPHA
+ double precision, parameter :: PRHO = 0.0
+
+ ! reference model and target model choice
+ integer, parameter :: IMODEL_SYN = 0
+ integer, parameter :: IMODEL_DAT = 2
+ !-----------------------------------------------------------------------------------------
+ ! 0 1 2 3
+ !-----------------------------------------------------------------------------------------
+ ! ISURFACE=1, IMODEL_SYN : homo checker het
+ ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het pert
+ !-----------------------------------------------------------------------------------------
+ ! ISURFACE=0, IMODEL_SYN : homo 1D model checker het
+ ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert het pert
+ !-----------------------------------------------------------------------------------------
+
+ ! whether you want to smooth the structure model for computing the synthetics
+ integer, parameter :: ISMOOTH_INITIAL_MODEL = 0
+ integer, parameter :: ISMOOTH_MODEL_UPDATE = 0
+ integer, parameter :: ISMOOTH_EVENT_KERNEL = 1
+ integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
+
+ ! smoothing scalelength
+ double precision, parameter :: GAMMA_SMOOTH_MODEL = 30.0e3
+
+!========================
+
+! boolean parameters
+! IKER: (0) waveform
+! (1) traveltime, cross-correlation, misfit
+! (2) amplitude, cross-correlation, misfit
+! (3) traveltime, multitaper
+! (4) amplitude, multitaper
+! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
+! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
+ integer, parameter :: IKER = 1
+ integer, parameter :: IAMP_VEL = 0 ! measure amplitudes between velocity traces
+
+! KEY: USE ONE LINE OR THE OTHER
+ integer, parameter :: ISURFACE = 1, NCOMP = 1, NABSORB = 4 ! surface waves
+! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
+
+! parameters controlling what to write to file
+! NOTE: for the tomography simulations, ALL of these can be .false.
+
+ logical, parameter :: WRITE_STF_F = .false.
+ logical, parameter :: WRITE_SEISMO_F = .false.
+! logical, parameter :: WRITE_SPECTRA_F = .false.
+! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
+
+ logical, parameter :: WRITE_SEISMO_RECONSTRUCT = .false. ! multitaper
+! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
+
+ logical, parameter :: WRITE_STF_A = .false.
+ logical, parameter :: WRITE_SEISMO_A = .false. ! source inversions
+! logical, parameter :: WRITE_SPECTRA_A = .false.
+! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
+
+ logical, parameter :: WRITE_KERNELS = .true. ! write all nine kernels
+ logical, parameter :: WRITE_KERNEL_SNAPSHOTS = .false. ! kernel snapshots
+ logical, parameter :: WRITE_WAVFIELD_SNAPSHOTS = .false. ! wavefield snapshots
+
+!--------------------------------------
+! INVERSION PARAMETERS
+
+ ! whether you want to compute kernels, or simply the misfit function
+ logical, parameter :: COMPUTE_KERNELS = .true.
+
+ ! whether to use the data subspace method, which has a Hessian
+ logical, parameter :: HESSIAN = .false.
+
+ ! stopping criteria
+ ! NITERATION : number of iterations
+ ! VAR_RED_MIN : minimum variance reduction (in percent)
+ ! SIGMA_FAC : stop if a model value exceeds SIGMA_FAC * sigma_m
+ ! CONV_STOP : stop when the misfit value is this fraction of the INITIAL misfit value
+ integer, parameter :: NITERATION = 16
+ double precision, parameter :: VAR_RED_MIN = 8.0
+
+ ! Gaussian errors containted in input file
+ ! see wave2d_sigmas.m and INPUT/sigma_0p1_pert.dat
+ double precision, parameter :: SIGMA_DT = 0.10
+ double precision, parameter :: SIGMA_DLNA = 1.0
+ double precision, parameter :: SIGMA_WAVEFORM = 1.0
+ logical, parameter :: ADD_DATA_ERRORS = .true.
+
+ !double precision, parameter :: SIGMA_FAC = 2.0
+ !double precision, parameter :: CONV_STOP = 1.0e-4
+
+ ! order of interpolating polynomial in conjugate gradient algorithm
+ ! using POLY_ORDER = 3 required computing the gradient of the test model
+ integer, parameter :: POLY_ORDER = 2 ! 2 (preferred) or 3
+
+ ! use inversion and scaling used for the GJI paper
+ !integer, parameter :: GJI_PAPER = 1
+
+ ! what to perturb, what to invert
+ integer, parameter :: PERT_STRUCT = 1
+ integer, parameter :: PERT_SOURCE_T = 1
+ integer, parameter :: PERT_SOURCE_X = 1
+
+ integer, parameter :: INV_STRUCT = 1
+ integer, parameter :: INV_SOURCE_T = 1
+ integer, parameter :: INV_SOURCE_X = 1
+
+ ! whether to include the model norm term in the misfit function, which acts like damping
+ logical, parameter :: INCLUDE_MODEL_NORM = .true.
+
+ ! log file showing source loactions
+ logical, parameter :: ISOURCE_LOG = .true.
+
+ ! DO NOT CHANGE THESE
+ integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ integer, parameter :: NVAR_SOURCE = 3 ! x position, z position, origin time
+ integer, parameter :: NVAR = NVAR_STRUCT + NVAR_SOURCE
+
+! parameterization for structure in the inversion
+! 1 : kappa-mu-rho
+! 2 : alpha-beta-rho
+! 3 : c-beta-rho
+ integer, parameter :: STRUCTURE_PARAMETER_TYPE = 2
+
+! scalelength of smoothing Gaussian
+! GJI paper : 30,60,90
+! body waves : 10
+ double precision, parameter :: GAMMA_SMOOTH = 60.0e3
+
+ ! parameters for smoothing
+ logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
+ logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
+
+! homogeneous background model (S.I. units)
+ double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
+ double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
+ double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
+
+!---------------------------------------------------------------
+! measurement windows
+
+ double precision, parameter :: HWIN1 = 2.5*hdur ! half-window width for pulse (2.5)
+ double precision, parameter :: HWIN2 = 5.0*hdur ! half-window with for MTM measurement (zero-padding)
+
+!---------------------------------------------------------------
+! CHT: do not change these
+
+! UTM zone for Southern California region
+! integer, parameter :: UTM_PROJECTION_ZONE = 11
+
+! to suppress UTM projection for SCEC benchmarks
+ logical, parameter :: SUPPRESS_UTM_PROJECTION = .false.
+
+! flag for projection from latitude/longitude to UTM, and back
+ integer, parameter :: ILONGLAT2UTM = 0, IUTM2LONGLAT = 1
+
+! flag for projection from latitude/longitude to mesh-UTM, and back
+ integer, parameter :: ILONLAT2MESH = 0, IMESH2LONLAT = 1
+
+! max number of fake receivers
+ integer, parameter :: MAX_SR_FAKE = 1000
+
+! max number of events, receivers, and phass
+ integer, parameter :: MAX_EVENT = 50
+ integer, parameter :: MAX_SR = 1400
+ integer, parameter :: MAX_PHASE = 1
+ integer, parameter :: MAX_COMP = NCOMP
+
+!========================
+! GRID AND GLL POINTS
+
+ integer, parameter :: NELE = MAX(NEX,NEZ)
+ integer, parameter :: NSPEC = NEX*NEZ
+
+! number of GLL points (polynomial degree plus one)
+ integer, parameter :: NGLLX = 5
+ integer, parameter :: NGLLZ = 5
+ integer, parameter :: NGLL = MAX(NGLLX,NGLLZ)
+
+! number of global points
+ integer, parameter :: NGLLSQUARE = NGLLX * NGLLZ ! element GLL points
+ integer, parameter :: NGLOB = ((NGLLX-1)*NEX + 1)*((NGLLZ-1)*NEZ +1) ! global GLL points
+ integer, parameter :: NLOCAL = NGLLX * NGLLZ * NSPEC ! local GLL points
+ integer, parameter :: NSPEC_CORNER = (NEX+1) * (NEZ+1) ! element corner points
+
+! number of nodes for 2D and 3D shape functions for hexahedra
+! we use 8-node mesh bricks, which are more stable than 27-node elements
+ integer, parameter :: NGNOD = 8, NGNOD2D = 4
+
+! number of iterations to solve the system for xi and eta
+ integer, parameter :: NUM_ITER = 1
+
+! very large and very small values
+ double precision, parameter :: HUGEVAL = 1.0e30, TINYVAL = 1.0e-9
+
+! for the Gauss-Lobatto-Legendre points and weights
+ double precision, parameter :: GAUSSALPHA = 0.0,GAUSSBETA = 0.0
+
+!
+! CONSTANTS
+!
+ double precision, parameter :: PI = 3.141592653589793
+ double precision, parameter :: FOUR_THIRDS = 4.0 / 3.0
+ double precision, parameter :: ONE_THIRD = 1.0 / 3.0
+ double precision, parameter :: ONEOVERTWO = 0.5
+ double precision, parameter :: DEG = 180.0/PI
+! double precision, parameter :: EPS = 1.0e-35
+
+!!$! parameter for FFTW
+!!$ integer, parameter :: NOUT = NSTEP/2 + 1
+!!$
+!!$! bounds for bandpass filter
+!!$ double precision, parameter :: hwid = 3.0 ! HALF-width of window (s)
+!!$ double precision, parameter :: tmin = 2.0*hdur-hwid
+!!$ double precision, parameter :: tmax = 2.0*hdur+hwid
+!!$ double precision, parameter :: fmin = 1.0/tmax, fmax = 1.0/tmin
+!!$ double precision, parameter :: trbdndw = 0.3, a = 30.0
+!!$ integer, parameter :: passes = 2, iord = 4
+!!$
+!!$!
+!!$! MULTI-TAPER PARAMETERS
+!!$!
+!!$! Ying Zhou: The fit between the recovered data and the data can be improved
+!!$! by either increasing the window width (HWIN above) or by decreasing NPI.
+!!$! In her experience, NPI = 2.5 is good for noisy data.
+!!$! For synthetic data, we can use a lower NPI.
+!!$! number of tapers should be fixed as twice NPI -- see Latex notes
+!!$ double precision, parameter :: WTR = 0.02
+!!$ double precision, parameter :: NPI = 2.5
+!!$ integer, parameter :: NTAPER = int(2.0*NPI)
+!!$
+!!$! FFT parameters
+!!$! FOR SOME REASON, I CANNOT SET NDIM = npt when lnpt <= 11 -- WHY NOT?
+!!$ integer, parameter :: lnpt = 13, npt = 2**lnpt, NDIM = 20000
+!!$ !integer, parameter :: lnpt = 14, npt = 2**lnpt, NDIM = npt
+!!$
+!!$ !double precision, parameter :: ZZIGN = -1.0 ! sign convention -- WATCH OUT
+!!$ double precision, parameter :: FORWARD_FFT = 1.0
+!!$ double precision, parameter :: REVERSE_FFT = -1.0
+
+ integer, parameter :: NDIM = 20000
+
+end module wave2d_constants
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex06.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex06.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex06.f90 2009-12-02 18:44:58 UTC (rev 16060)
@@ -0,0 +1,336 @@
+module wave2d_constants
+
+! This file is copied into the output directory when wave2d.f90 is run
+! There are two basic options: membrane (surface) waves and body waves.
+! Several lines need to be commented/uncommented to go from one to the other.
+
+! index of the reference directory for the simulation output
+ integer, parameter :: IRUNZ = 600
+
+!========================
+! GRID, TIME-STEP, AND SOURCE PARAMETERS
+
+! NFRAME : number of frames to save: membrane, 10; body, 20
+! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
+! NSTEP : number of timesteps
+ integer, parameter :: NFRAME = 10 ! membrane
+ !integer, parameter :: NFRAME = 11 ! body
+ integer, parameter :: NSAVE = 400 ! 200,400
+ integer, parameter :: NSTEP = NFRAME*NSAVE
+
+! time step in seconds
+ !double precision, parameter :: DT = 2.0e-2 ! body waves
+ double precision, parameter :: DT = 6.0e-2 ! membrane surface waves
+
+! temporal properties of source (source time function)
+ integer, parameter :: ISRC_TIME = 1 ! type (1)
+ double precision, parameter :: hdur = 10.0 ! HALF-duration (s), membrane
+ !double precision, parameter :: hdur = 2.0 ! HALF-duration (s), body
+ double precision, parameter :: tshift = 2.0*DT*dble(NSAVE) ! time shift (s)
+ !double precision, parameter :: tshift = 8.0*hdur
+ logical, parameter :: SRC_TAPER = .true. ! taper the endpoints of the time series
+
+! normalization factor of point source force
+ double precision, parameter :: FNORM = 1.0d10
+
+! forward wavefield source-time function
+! (For membrane waves, FOR_Y = 1 is all that matters.)
+ integer, parameter :: FOR_X = 1 ! boolean
+ integer, parameter :: FOR_Y = 1 ! boolean
+ integer, parameter :: FOR_Z = 1 ! boolean
+
+! adjoint wavefield source-time function
+! (For membrane waves, REV_Y = 1 is all that matters.)
+! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
+ integer, parameter :: REV_X = 1 ! boolean
+ integer, parameter :: REV_Y = 1 ! boolean
+ integer, parameter :: REV_Z = 1 ! boolean
+
+! spatial properties of sources
+! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
+! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
+! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
+ integer, parameter :: ISRC_SPACE = 6 ! see wave2d.f90
+
+! spatial properties of receivers
+! IREC_SPACE
+! (1) individual station(s)
+! (2) SoCal (used for GJI paper)
+! (3) regular mesh on land
+! (4) regular mesh
+! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
+! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
+! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
+ integer, parameter :: IREC_SPACE = 2 ! see wave2d.f90
+ integer, parameter :: NMESH_REC = 17
+ double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m ! 4km for membrane surface waves
+ double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m ! 15km for membrane surface waves
+ double precision, parameter :: STATION_COAST_BUFFER = 0.0e3 ! m
+
+! lower right corner for membrane surface waves plotting grid
+ double precision, parameter :: LAT_MIN = 32.0
+ double precision, parameter :: LON_MIN = -120.0
+ integer, parameter :: UTM_PROJECTION_ZONE = 11 ! southern California
+
+! mesh specifications: membrane surface waves
+ double precision, parameter :: LENGTH = 480.0e3 ! m
+ double precision, parameter :: HEIGHT = 480.0e3 ! m
+ double precision, parameter :: AREA = LENGTH*HEIGHT
+ integer, parameter :: NEX = 40
+ integer, parameter :: NEZ = 40
+!!$
+!!$! mesh specifications: body waves
+!!$ double precision, parameter :: LENGTH = 200.0e3 ! m ! 400 for 1D body waves
+!!$ double precision, parameter :: HEIGHT = 80.0e3 ! m
+!!$ integer, parameter :: NEX = 80 ! 160
+!!$ integer, parameter :: NEZ = 32 ! 32
+
+!========================
+! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
+
+! model perturbations for HOMOGENEOUS model (or perturbation)
+! scaling from beta to alpha
+! value is from Master et al. (2000), "The relative behavior of shear velocity..."
+ double precision, parameter :: R_BETA_OVER_ALPHA = 1.3
+ double precision, parameter :: PBETA = 10.0
+ double precision, parameter :: PALPHA = PBETA / R_BETA_OVER_ALPHA
+ double precision, parameter :: PRHO = 0.0
+
+ ! reference model and target model choice
+ integer, parameter :: IMODEL_SYN = 0
+ integer, parameter :: IMODEL_DAT = 2
+ !-----------------------------------------------------------------------------------------
+ ! 0 1 2 3
+ !-----------------------------------------------------------------------------------------
+ ! ISURFACE=1, IMODEL_SYN : homo checker het
+ ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het pert
+ !-----------------------------------------------------------------------------------------
+ ! ISURFACE=0, IMODEL_SYN : homo 1D model checker het
+ ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert het pert
+ !-----------------------------------------------------------------------------------------
+
+ ! whether you want to smooth the structure model for computing the synthetics
+ integer, parameter :: ISMOOTH_INITIAL_MODEL = 0
+ integer, parameter :: ISMOOTH_MODEL_UPDATE = 0
+ integer, parameter :: ISMOOTH_EVENT_KERNEL = 1
+ integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
+
+ ! smoothing scalelength
+ double precision, parameter :: GAMMA_SMOOTH_MODEL = 30.0e3
+
+!========================
+
+! boolean parameters
+! IKER: (0) waveform
+! (1) traveltime, cross-correlation, misfit
+! (2) amplitude, cross-correlation, misfit
+! (3) traveltime, multitaper
+! (4) amplitude, multitaper
+! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
+! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
+ integer, parameter :: IKER = 1
+ integer, parameter :: IAMP_VEL = 0 ! measure amplitudes between velocity traces
+
+! KEY: USE ONE LINE OR THE OTHER
+ integer, parameter :: ISURFACE = 1, NCOMP = 1, NABSORB = 4 ! surface waves
+! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
+
+! parameters controlling what to write to file
+! NOTE: for the tomography simulations, ALL of these can be .false.
+
+ logical, parameter :: WRITE_STF_F = .false.
+ logical, parameter :: WRITE_SEISMO_F = .false.
+! logical, parameter :: WRITE_SPECTRA_F = .false.
+! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
+
+ logical, parameter :: WRITE_SEISMO_RECONSTRUCT = .false. ! multitaper
+! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
+
+ logical, parameter :: WRITE_STF_A = .false.
+ logical, parameter :: WRITE_SEISMO_A = .false. ! source inversions
+! logical, parameter :: WRITE_SPECTRA_A = .false.
+! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
+
+ logical, parameter :: WRITE_KERNELS = .true. ! write all nine kernels
+ logical, parameter :: WRITE_KERNEL_SNAPSHOTS = .false. ! kernel snapshots
+ logical, parameter :: WRITE_WAVFIELD_SNAPSHOTS = .false. ! wavefield snapshots
+
+!--------------------------------------
+! INVERSION PARAMETERS
+
+ ! whether you want to compute kernels, or simply the misfit function
+ logical, parameter :: COMPUTE_KERNELS = .true.
+
+ ! whether to use the data subspace method, which has a Hessian
+ logical, parameter :: HESSIAN = .false.
+
+ ! stopping criteria
+ ! NITERATION : number of iterations
+ ! VAR_RED_MIN : minimum variance reduction (in percent)
+ ! SIGMA_FAC : stop if a model value exceeds SIGMA_FAC * sigma_m
+ ! CONV_STOP : stop when the misfit value is this fraction of the INITIAL misfit value
+ integer, parameter :: NITERATION = 16
+ double precision, parameter :: VAR_RED_MIN = 8.0
+
+ ! Gaussian errors containted in input file
+ ! see wave2d_sigmas.m and INPUT/sigma_0p1_pert.dat
+ double precision, parameter :: SIGMA_DT = 0.10
+ double precision, parameter :: SIGMA_DLNA = 1.0
+ double precision, parameter :: SIGMA_WAVEFORM = 1.0
+ logical, parameter :: ADD_DATA_ERRORS = .true.
+
+ !double precision, parameter :: SIGMA_FAC = 2.0
+ !double precision, parameter :: CONV_STOP = 1.0e-4
+
+ ! order of interpolating polynomial in conjugate gradient algorithm
+ ! using POLY_ORDER = 3 required computing the gradient of the test model
+ integer, parameter :: POLY_ORDER = 2 ! 2 (preferred) or 3
+
+ ! use inversion and scaling used for the GJI paper
+ !integer, parameter :: GJI_PAPER = 1
+
+ ! what to perturb, what to invert
+ integer, parameter :: PERT_STRUCT = 1
+ integer, parameter :: PERT_SOURCE_T = 0
+ integer, parameter :: PERT_SOURCE_X = 0
+
+ integer, parameter :: INV_STRUCT = 1
+ integer, parameter :: INV_SOURCE_T = 0
+ integer, parameter :: INV_SOURCE_X = 0
+
+ ! whether to include the model norm term in the misfit function, which acts like damping
+ logical, parameter :: INCLUDE_MODEL_NORM = .true.
+
+ ! log file showing source loactions
+ logical, parameter :: ISOURCE_LOG = .true.
+
+ ! DO NOT CHANGE THESE
+ integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ integer, parameter :: NVAR_SOURCE = 3 ! x position, z position, origin time
+ integer, parameter :: NVAR = NVAR_STRUCT + NVAR_SOURCE
+
+! parameterization for structure in the inversion
+! 1 : kappa-mu-rho
+! 2 : alpha-beta-rho
+! 3 : c-beta-rho
+ integer, parameter :: STRUCTURE_PARAMETER_TYPE = 2
+
+! scalelength of smoothing Gaussian
+! GJI paper : 30,60,90
+! body waves : 10
+ double precision, parameter :: GAMMA_SMOOTH = 60.0e3
+
+ ! parameters for smoothing
+ logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
+ logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
+
+! homogeneous background model (S.I. units)
+ double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
+ double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
+ double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
+
+!---------------------------------------------------------------
+! measurement windows
+
+ double precision, parameter :: HWIN1 = 2.5*hdur ! half-window width for pulse (2.5)
+ double precision, parameter :: HWIN2 = 5.0*hdur ! half-window with for MTM measurement (zero-padding)
+
+!---------------------------------------------------------------
+! CHT: do not change these
+
+! UTM zone for Southern California region
+! integer, parameter :: UTM_PROJECTION_ZONE = 11
+
+! to suppress UTM projection for SCEC benchmarks
+ logical, parameter :: SUPPRESS_UTM_PROJECTION = .false.
+
+! flag for projection from latitude/longitude to UTM, and back
+ integer, parameter :: ILONGLAT2UTM = 0, IUTM2LONGLAT = 1
+
+! flag for projection from latitude/longitude to mesh-UTM, and back
+ integer, parameter :: ILONLAT2MESH = 0, IMESH2LONLAT = 1
+
+! max number of fake receivers
+ integer, parameter :: MAX_SR_FAKE = 1000
+
+! max number of events, receivers, and phass
+ integer, parameter :: MAX_EVENT = 50
+ integer, parameter :: MAX_SR = 1400
+ integer, parameter :: MAX_PHASE = 1
+ integer, parameter :: MAX_COMP = NCOMP
+
+!========================
+! GRID AND GLL POINTS
+
+ integer, parameter :: NELE = MAX(NEX,NEZ)
+ integer, parameter :: NSPEC = NEX*NEZ
+
+! number of GLL points (polynomial degree plus one)
+ integer, parameter :: NGLLX = 5
+ integer, parameter :: NGLLZ = 5
+ integer, parameter :: NGLL = MAX(NGLLX,NGLLZ)
+
+! number of global points
+ integer, parameter :: NGLLSQUARE = NGLLX * NGLLZ ! element GLL points
+ integer, parameter :: NGLOB = ((NGLLX-1)*NEX + 1)*((NGLLZ-1)*NEZ +1) ! global GLL points
+ integer, parameter :: NLOCAL = NGLLX * NGLLZ * NSPEC ! local GLL points
+ integer, parameter :: NSPEC_CORNER = (NEX+1) * (NEZ+1) ! element corner points
+
+! number of nodes for 2D and 3D shape functions for hexahedra
+! we use 8-node mesh bricks, which are more stable than 27-node elements
+ integer, parameter :: NGNOD = 8, NGNOD2D = 4
+
+! number of iterations to solve the system for xi and eta
+ integer, parameter :: NUM_ITER = 1
+
+! very large and very small values
+ double precision, parameter :: HUGEVAL = 1.0e30, TINYVAL = 1.0e-9
+
+! for the Gauss-Lobatto-Legendre points and weights
+ double precision, parameter :: GAUSSALPHA = 0.0,GAUSSBETA = 0.0
+
+!
+! CONSTANTS
+!
+ double precision, parameter :: PI = 3.141592653589793
+ double precision, parameter :: FOUR_THIRDS = 4.0 / 3.0
+ double precision, parameter :: ONE_THIRD = 1.0 / 3.0
+ double precision, parameter :: ONEOVERTWO = 0.5
+ double precision, parameter :: DEG = 180.0/PI
+! double precision, parameter :: EPS = 1.0e-35
+
+!!$! parameter for FFTW
+!!$ integer, parameter :: NOUT = NSTEP/2 + 1
+!!$
+!!$! bounds for bandpass filter
+!!$ double precision, parameter :: hwid = 3.0 ! HALF-width of window (s)
+!!$ double precision, parameter :: tmin = 2.0*hdur-hwid
+!!$ double precision, parameter :: tmax = 2.0*hdur+hwid
+!!$ double precision, parameter :: fmin = 1.0/tmax, fmax = 1.0/tmin
+!!$ double precision, parameter :: trbdndw = 0.3, a = 30.0
+!!$ integer, parameter :: passes = 2, iord = 4
+!!$
+!!$!
+!!$! MULTI-TAPER PARAMETERS
+!!$!
+!!$! Ying Zhou: The fit between the recovered data and the data can be improved
+!!$! by either increasing the window width (HWIN above) or by decreasing NPI.
+!!$! In her experience, NPI = 2.5 is good for noisy data.
+!!$! For synthetic data, we can use a lower NPI.
+!!$! number of tapers should be fixed as twice NPI -- see Latex notes
+!!$ double precision, parameter :: WTR = 0.02
+!!$ double precision, parameter :: NPI = 2.5
+!!$ integer, parameter :: NTAPER = int(2.0*NPI)
+!!$
+!!$! FFT parameters
+!!$! FOR SOME REASON, I CANNOT SET NDIM = npt when lnpt <= 11 -- WHY NOT?
+!!$ integer, parameter :: lnpt = 13, npt = 2**lnpt, NDIM = 20000
+!!$ !integer, parameter :: lnpt = 14, npt = 2**lnpt, NDIM = npt
+!!$
+!!$ !double precision, parameter :: ZZIGN = -1.0 ! sign convention -- WATCH OUT
+!!$ double precision, parameter :: FORWARD_FFT = 1.0
+!!$ double precision, parameter :: REVERSE_FFT = -1.0
+
+ integer, parameter :: NDIM = 20000
+
+end module wave2d_constants
Added: seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex07.f90
===================================================================
--- seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex07.f90 (rev 0)
+++ seismo/3D/ADJOINT_TOMO/iterate_adj/SEM2D_iterate/src/wave2d_constants_ex07.f90 2009-12-02 18:44:58 UTC (rev 16060)
@@ -0,0 +1,336 @@
+module wave2d_constants
+
+! This file is copied into the output directory when wave2d.f90 is run
+! There are two basic options: membrane (surface) waves and body waves.
+! Several lines need to be commented/uncommented to go from one to the other.
+
+! index of the reference directory for the simulation output
+ integer, parameter :: IRUNZ = 700
+
+!========================
+! GRID, TIME-STEP, AND SOURCE PARAMETERS
+
+! NFRAME : number of frames to save: membrane, 10; body, 20
+! NSAVE : timestep increment to save the wavefield: membrane, 400; body, 400
+! NSTEP : number of timesteps
+ integer, parameter :: NFRAME = 10 ! membrane
+ !integer, parameter :: NFRAME = 11 ! body
+ integer, parameter :: NSAVE = 400 ! 200,400
+ integer, parameter :: NSTEP = NFRAME*NSAVE
+
+! time step in seconds
+ !double precision, parameter :: DT = 2.0e-2 ! body waves
+ double precision, parameter :: DT = 6.0e-2 ! membrane surface waves
+
+! temporal properties of source (source time function)
+ integer, parameter :: ISRC_TIME = 1 ! type (1)
+ double precision, parameter :: hdur = 10.0 ! HALF-duration (s), membrane
+ !double precision, parameter :: hdur = 2.0 ! HALF-duration (s), body
+ double precision, parameter :: tshift = 2.0*DT*dble(NSAVE) ! time shift (s)
+ !double precision, parameter :: tshift = 8.0*hdur
+ logical, parameter :: SRC_TAPER = .true. ! taper the endpoints of the time series
+
+! normalization factor of point source force
+ double precision, parameter :: FNORM = 1.0d10
+
+! forward wavefield source-time function
+! (For membrane waves, FOR_Y = 1 is all that matters.)
+ integer, parameter :: FOR_X = 1 ! boolean
+ integer, parameter :: FOR_Y = 1 ! boolean
+ integer, parameter :: FOR_Z = 1 ! boolean
+
+! adjoint wavefield source-time function
+! (For membrane waves, REV_Y = 1 is all that matters.)
+! integer, parameter :: IPSV = 0 ! boolean: plot PSV or SH kernels
+ integer, parameter :: REV_X = 1 ! boolean
+ integer, parameter :: REV_Y = 1 ! boolean
+ integer, parameter :: REV_Z = 1 ! boolean
+
+! spatial properties of sources
+! BODY WAVE OPTIONS : (6) a point source EVENT (GJI paper)
+! SURFACE WAVE OPTIONS : (1) point source, (2) finite segment, (3) CA shelf boundary
+! (4) CA coast, (5) finite circle, (6) a point source EVENT (GJI paper)
+ integer, parameter :: ISRC_SPACE = 6 ! see wave2d.f90
+
+! spatial properties of receivers
+! IREC_SPACE
+! (1) individual station(s)
+! (2) SoCal (used for GJI paper)
+! (3) regular mesh on land
+! (4) regular mesh
+! NMESH_REC : determines the number of receivers in a regular mesh (IREC_SPACE > 3)
+! STATION_GRID_BUFFER : exclude stations within this distance from edge of grid
+! STATION_COAST_BUFFER : exclude stations within this distance from edge of coast
+ integer, parameter :: IREC_SPACE = 2 ! see wave2d.f90
+ integer, parameter :: NMESH_REC = 17
+ double precision, parameter :: SOURCE_GRID_BUFFER = 4.0e3 ! m ! 4km for membrane surface waves
+ double precision, parameter :: STATION_GRID_BUFFER = 15.0e3 ! m ! 15km for membrane surface waves
+ double precision, parameter :: STATION_COAST_BUFFER = 0.0e3 ! m
+
+! lower right corner for membrane surface waves plotting grid
+ double precision, parameter :: LAT_MIN = 32.0
+ double precision, parameter :: LON_MIN = -120.0
+ integer, parameter :: UTM_PROJECTION_ZONE = 11 ! southern California
+
+! mesh specifications: membrane surface waves
+ double precision, parameter :: LENGTH = 480.0e3 ! m
+ double precision, parameter :: HEIGHT = 480.0e3 ! m
+ double precision, parameter :: AREA = LENGTH*HEIGHT
+ integer, parameter :: NEX = 40
+ integer, parameter :: NEZ = 40
+!!$
+!!$! mesh specifications: body waves
+!!$ double precision, parameter :: LENGTH = 200.0e3 ! m ! 400 for 1D body waves
+!!$ double precision, parameter :: HEIGHT = 80.0e3 ! m
+!!$ integer, parameter :: NEX = 80 ! 160
+!!$ integer, parameter :: NEZ = 32 ! 32
+
+!========================
+! MODEL SPECIFICATIONS (REFERENCE MODEL AND TARGET MODEL)
+
+! model perturbations for HOMOGENEOUS model (or perturbation)
+! scaling from beta to alpha
+! value is from Master et al. (2000), "The relative behavior of shear velocity..."
+ double precision, parameter :: R_BETA_OVER_ALPHA = 1.3
+ double precision, parameter :: PBETA = 10.0
+ double precision, parameter :: PALPHA = PBETA / R_BETA_OVER_ALPHA
+ double precision, parameter :: PRHO = 0.0
+
+ ! reference model and target model choice
+ integer, parameter :: IMODEL_SYN = 0
+ integer, parameter :: IMODEL_DAT = 2
+ !-----------------------------------------------------------------------------------------
+ ! 0 1 2 3
+ !-----------------------------------------------------------------------------------------
+ ! ISURFACE=1, IMODEL_SYN : homo checker het
+ ! ISURFACE=1, IMODEL_DAT : homo pert checker pert het pert
+ !-----------------------------------------------------------------------------------------
+ ! ISURFACE=0, IMODEL_SYN : homo 1D model checker het
+ ! ISURFACE=0, IMODEL_DAT : homo pert 1D model checker pert het pert
+ !-----------------------------------------------------------------------------------------
+
+ ! whether you want to smooth the structure model for computing the synthetics
+ integer, parameter :: ISMOOTH_INITIAL_MODEL = 0
+ integer, parameter :: ISMOOTH_MODEL_UPDATE = 0
+ integer, parameter :: ISMOOTH_EVENT_KERNEL = 1
+ integer, parameter :: ISMOOTH_MISFIT_KERNEL = 1
+
+ ! smoothing scalelength
+ double precision, parameter :: GAMMA_SMOOTH_MODEL = 30.0e3
+
+!========================
+
+! boolean parameters
+! IKER: (0) waveform
+! (1) traveltime, cross-correlation, misfit
+! (2) amplitude, cross-correlation, misfit
+! (3) traveltime, multitaper
+! (4) amplitude, multitaper
+! (5) traveltime, cross-correlation, sampling -- NOT AN OPTION
+! (6) amplitude, cross-correlation, sampling -- NOT AN OPTION
+ integer, parameter :: IKER = 1
+ integer, parameter :: IAMP_VEL = 0 ! measure amplitudes between velocity traces
+
+! KEY: USE ONE LINE OR THE OTHER
+ integer, parameter :: ISURFACE = 1, NCOMP = 1, NABSORB = 4 ! surface waves
+! integer, parameter :: ISURFACE = 0, NCOMP = 3, NABSORB = 3 ! body waves
+
+! parameters controlling what to write to file
+! NOTE: for the tomography simulations, ALL of these can be .false.
+
+ logical, parameter :: WRITE_STF_F = .false.
+ logical, parameter :: WRITE_SEISMO_F = .false.
+! logical, parameter :: WRITE_SPECTRA_F = .false.
+! logical, parameter :: WRITE_SPECTRAL_MAP_F = .false.
+
+ logical, parameter :: WRITE_SEISMO_RECONSTRUCT = .false. ! multitaper
+! logical, parameter :: WRITE_MTM_FILES = .false. ! multitaper
+
+ logical, parameter :: WRITE_STF_A = .false.
+ logical, parameter :: WRITE_SEISMO_A = .false. ! source inversions
+! logical, parameter :: WRITE_SPECTRA_A = .false.
+! logical, parameter :: WRITE_SPECTRAL_MAP_A = .false.
+
+ logical, parameter :: WRITE_KERNELS = .true. ! write all nine kernels
+ logical, parameter :: WRITE_KERNEL_SNAPSHOTS = .false. ! kernel snapshots
+ logical, parameter :: WRITE_WAVFIELD_SNAPSHOTS = .false. ! wavefield snapshots
+
+!--------------------------------------
+! INVERSION PARAMETERS
+
+ ! whether you want to compute kernels, or simply the misfit function
+ logical, parameter :: COMPUTE_KERNELS = .true.
+
+ ! whether to use the data subspace method, which has a Hessian
+ logical, parameter :: HESSIAN = .false.
+
+ ! stopping criteria
+ ! NITERATION : number of iterations
+ ! VAR_RED_MIN : minimum variance reduction (in percent)
+ ! SIGMA_FAC : stop if a model value exceeds SIGMA_FAC * sigma_m
+ ! CONV_STOP : stop when the misfit value is this fraction of the INITIAL misfit value
+ integer, parameter :: NITERATION = 0
+ double precision, parameter :: VAR_RED_MIN = 8.0
+
+ ! Gaussian errors containted in input file
+ ! see wave2d_sigmas.m and INPUT/sigma_0p1_pert.dat
+ double precision, parameter :: SIGMA_DT = 0.10
+ double precision, parameter :: SIGMA_DLNA = 1.0
+ double precision, parameter :: SIGMA_WAVEFORM = 1.0
+ logical, parameter :: ADD_DATA_ERRORS = .true.
+
+ !double precision, parameter :: SIGMA_FAC = 2.0
+ !double precision, parameter :: CONV_STOP = 1.0e-4
+
+ ! order of interpolating polynomial in conjugate gradient algorithm
+ ! using POLY_ORDER = 3 required computing the gradient of the test model
+ integer, parameter :: POLY_ORDER = 2 ! 2 (preferred) or 3
+
+ ! use inversion and scaling used for the GJI paper
+ !integer, parameter :: GJI_PAPER = 1
+
+ ! what to perturb, what to invert
+ integer, parameter :: PERT_STRUCT = 1
+ integer, parameter :: PERT_SOURCE_T = 0
+ integer, parameter :: PERT_SOURCE_X = 0
+
+ integer, parameter :: INV_STRUCT = 1
+ integer, parameter :: INV_SOURCE_T = 0
+ integer, parameter :: INV_SOURCE_X = 0
+
+ ! whether to include the model norm term in the misfit function, which acts like damping
+ logical, parameter :: INCLUDE_MODEL_NORM = .true.
+
+ ! log file showing source loactions
+ logical, parameter :: ISOURCE_LOG = .true.
+
+ ! DO NOT CHANGE THESE
+ integer, parameter :: NVAR_STRUCT = 2 ! alpha, beta (seismic velocities only)
+ integer, parameter :: NVAR_SOURCE = 3 ! x position, z position, origin time
+ integer, parameter :: NVAR = NVAR_STRUCT + NVAR_SOURCE
+
+! parameterization for structure in the inversion
+! 1 : kappa-mu-rho
+! 2 : alpha-beta-rho
+! 3 : c-beta-rho
+ integer, parameter :: STRUCTURE_PARAMETER_TYPE = 2
+
+! scalelength of smoothing Gaussian
+! GJI paper : 30,60,90
+! body waves : 10
+ double precision, parameter :: GAMMA_SMOOTH = 60.0e3
+
+ ! parameters for smoothing
+ logical, parameter :: HIGH_RES_SMOOTHING = .true. ! smooth at high resolution
+ logical, parameter :: EXAMPLE_GAUSSIAN = .false. ! plot an example Gaussian
+
+! homogeneous background model (S.I. units)
+ double precision, parameter :: DENSITY = 2.60e3 ! kg/m^3
+ double precision, parameter :: INCOMPRESSIBILITY = 5.20e10 ! Pa
+ double precision, parameter :: RIGIDITY = 2.66e10 ! Pa
+
+!---------------------------------------------------------------
+! measurement windows
+
+ double precision, parameter :: HWIN1 = 2.5*hdur ! half-window width for pulse (2.5)
+ double precision, parameter :: HWIN2 = 5.0*hdur ! half-window with for MTM measurement (zero-padding)
+
+!---------------------------------------------------------------
+! CHT: do not change these
+
+! UTM zone for Southern California region
+! integer, parameter :: UTM_PROJECTION_ZONE = 11
+
+! to suppress UTM projection for SCEC benchmarks
+ logical, parameter :: SUPPRESS_UTM_PROJECTION = .false.
+
+! flag for projection from latitude/longitude to UTM, and back
+ integer, parameter :: ILONGLAT2UTM = 0, IUTM2LONGLAT = 1
+
+! flag for projection from latitude/longitude to mesh-UTM, and back
+ integer, parameter :: ILONLAT2MESH = 0, IMESH2LONLAT = 1
+
+! max number of fake receivers
+ integer, parameter :: MAX_SR_FAKE = 1000
+
+! max number of events, receivers, and phass
+ integer, parameter :: MAX_EVENT = 50
+ integer, parameter :: MAX_SR = 1400
+ integer, parameter :: MAX_PHASE = 1
+ integer, parameter :: MAX_COMP = NCOMP
+
+!========================
+! GRID AND GLL POINTS
+
+ integer, parameter :: NELE = MAX(NEX,NEZ)
+ integer, parameter :: NSPEC = NEX*NEZ
+
+! number of GLL points (polynomial degree plus one)
+ integer, parameter :: NGLLX = 5
+ integer, parameter :: NGLLZ = 5
+ integer, parameter :: NGLL = MAX(NGLLX,NGLLZ)
+
+! number of global points
+ integer, parameter :: NGLLSQUARE = NGLLX * NGLLZ ! element GLL points
+ integer, parameter :: NGLOB = ((NGLLX-1)*NEX + 1)*((NGLLZ-1)*NEZ +1) ! global GLL points
+ integer, parameter :: NLOCAL = NGLLX * NGLLZ * NSPEC ! local GLL points
+ integer, parameter :: NSPEC_CORNER = (NEX+1) * (NEZ+1) ! element corner points
+
+! number of nodes for 2D and 3D shape functions for hexahedra
+! we use 8-node mesh bricks, which are more stable than 27-node elements
+ integer, parameter :: NGNOD = 8, NGNOD2D = 4
+
+! number of iterations to solve the system for xi and eta
+ integer, parameter :: NUM_ITER = 1
+
+! very large and very small values
+ double precision, parameter :: HUGEVAL = 1.0e30, TINYVAL = 1.0e-9
+
+! for the Gauss-Lobatto-Legendre points and weights
+ double precision, parameter :: GAUSSALPHA = 0.0,GAUSSBETA = 0.0
+
+!
+! CONSTANTS
+!
+ double precision, parameter :: PI = 3.141592653589793
+ double precision, parameter :: FOUR_THIRDS = 4.0 / 3.0
+ double precision, parameter :: ONE_THIRD = 1.0 / 3.0
+ double precision, parameter :: ONEOVERTWO = 0.5
+ double precision, parameter :: DEG = 180.0/PI
+! double precision, parameter :: EPS = 1.0e-35
+
+!!$! parameter for FFTW
+!!$ integer, parameter :: NOUT = NSTEP/2 + 1
+!!$
+!!$! bounds for bandpass filter
+!!$ double precision, parameter :: hwid = 3.0 ! HALF-width of window (s)
+!!$ double precision, parameter :: tmin = 2.0*hdur-hwid
+!!$ double precision, parameter :: tmax = 2.0*hdur+hwid
+!!$ double precision, parameter :: fmin = 1.0/tmax, fmax = 1.0/tmin
+!!$ double precision, parameter :: trbdndw = 0.3, a = 30.0
+!!$ integer, parameter :: passes = 2, iord = 4
+!!$
+!!$!
+!!$! MULTI-TAPER PARAMETERS
+!!$!
+!!$! Ying Zhou: The fit between the recovered data and the data can be improved
+!!$! by either increasing the window width (HWIN above) or by decreasing NPI.
+!!$! In her experience, NPI = 2.5 is good for noisy data.
+!!$! For synthetic data, we can use a lower NPI.
+!!$! number of tapers should be fixed as twice NPI -- see Latex notes
+!!$ double precision, parameter :: WTR = 0.02
+!!$ double precision, parameter :: NPI = 2.5
+!!$ integer, parameter :: NTAPER = int(2.0*NPI)
+!!$
+!!$! FFT parameters
+!!$! FOR SOME REASON, I CANNOT SET NDIM = npt when lnpt <= 11 -- WHY NOT?
+!!$ integer, parameter :: lnpt = 13, npt = 2**lnpt, NDIM = 20000
+!!$ !integer, parameter :: lnpt = 14, npt = 2**lnpt, NDIM = npt
+!!$
+!!$ !double precision, parameter :: ZZIGN = -1.0 ! sign convention -- WATCH OUT
+!!$ double precision, parameter :: FORWARD_FFT = 1.0
+!!$ double precision, parameter :: REVERSE_FFT = -1.0
+
+ integer, parameter :: NDIM = 20000
+
+end module wave2d_constants
More information about the CIG-COMMITS
mailing list