[cig-commits] r18116 - in seismo/3D/SPECFEM3D_GLOBE/trunk: setup src/auxiliaries src/specfem3D
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Tue Mar 15 19:28:50 PDT 2011
Author: danielpeter
Date: 2011-03-15 19:28:49 -0700 (Tue, 15 Mar 2011)
New Revision: 18116
Added:
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_adj_source_frechet.f90
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/setup/config.h.in
seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_GMT_global.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/Makefile.in
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_arrays_source.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.f90
Log:
bug fix in array initialization; puts compute_adj_source_frechet routine into separate source file src/specfem3D/compute_adj_source_frechet.f90 (compilation of this routine with -O3 may take a while)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/setup/config.h.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/setup/config.h.in 2011-03-15 16:48:14 UTC (rev 18115)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/setup/config.h.in 2011-03-16 02:28:49 UTC (rev 18116)
@@ -31,3 +31,6 @@
/* Define to the version of this package. */
#undef PACKAGE_VERSION
+
+/* Uncomment and define to select optimized file i/o for regional simulations */
+#undef USE_MAP_FUNCTION
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_GMT_global.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_GMT_global.f90 2011-03-15 16:48:14 UTC (rev 18115)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_GMT_global.f90 2011-03-16 02:28:49 UTC (rev 18116)
@@ -50,17 +50,16 @@
! to avoid flickering in movies, the displacement field will get normalized with an
! averaged maximum value over the past few, available snapshots
logical,parameter :: USE_AVERAGED_MAXIMUM = .true.
-
! minimum number of frames to average maxima
integer,parameter :: AVERAGE_MINIMUM = 5
+ ! normalizes output values
+ logical, parameter :: NORMALIZE_VALUES = .true.
! muting source region
logical, parameter :: MUTE_SOURCE = .true.
- real(kind=CUSTOM_REAL) :: RADIUS_TO_MUTE = 1.0 ! start radius in degrees
- real(kind=CUSTOM_REAL) :: STARTTIME_TO_MUTE = 2.0 ! factor times hdur_movie
+ real(kind=CUSTOM_REAL) :: RADIUS_TO_MUTE = 0.5 ! start radius in degrees
+ real(kind=CUSTOM_REAL) :: STARTTIME_TO_MUTE = 0.5 ! adds seconds to shift starttime
- ! normalizes output values
- logical, parameter :: NORMALIZE_VALUES = .true.
!---------------------
@@ -154,9 +153,10 @@
integer, dimension(NB_SQUARE_EDGES_ONEDIR,NB_CUT_CASE) :: DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA
real(kind=CUSTOM_REAL) :: LAT_SOURCE,LON_SOURCE,DEP_SOURCE
- real(kind=CUSTOM_REAL) :: dist_lon,dist_lat,mute_factor
+ real(kind=CUSTOM_REAL) :: dist_lon,dist_lat,distance,mute_factor,val
character(len=256) line
-
+ real(kind=CUSTOM_REAL) :: cmt_hdur,cmt_t_shift,t0,hdur
+
! ************** PROGRAM STARTS HERE **************
print *
@@ -191,11 +191,14 @@
NSPEC,NSPEC2D_XI,NSPEC2D_ETA, &
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
NSPEC1D_RADIAL,NGLOB1D_RADIAL,NGLOB2DMAX_XMIN_XMAX,NGLOB2DMAX_YMIN_YMAX,NGLOB, &
- ratio_sampling_array, ner, doubling_index,r_bottom,r_top,this_region_has_a_doubling,rmins,rmaxs,CASE_3D, &
+ ratio_sampling_array, ner, doubling_index,r_bottom,r_top, &
+ this_region_has_a_doubling,rmins,rmaxs,CASE_3D, &
OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
- ROTATE_SEISMOGRAMS_RT,ratio_divide_central_cube,HONOR_1D_SPHERICAL_MOHO,CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,&
+ ROTATE_SEISMOGRAMS_RT,ratio_divide_central_cube,HONOR_1D_SPHERICAL_MOHO, &
+ CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,&
DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA,&
- WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,.false.,NOISE_TOMOGRAPHY)
+ WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE, &
+ USE_BINARY_FOR_LARGE_FILE,.false.,NOISE_TOMOGRAPHY)
if(.not. MOVIE_SURFACE) stop 'movie frames were not saved by the solver'
@@ -204,7 +207,8 @@
print *
if(MOVIE_COARSE) then
! note:
- ! nex_per_proc_xi*nex_per_proc_eta = nex_xi*nex_eta/nproc = nspec2d_top(iregion_crust_mantle) used in specfem3D.f90
+ ! nex_per_proc_xi*nex_per_proc_eta = nex_xi*nex_eta/nproc = nspec2d_top(iregion_crust_mantle)
+ ! used in specfem3D.f90
! and ilocnum = nmovie_points = 2 * 2 * NEX_XI * NEX_ETA / NPROC
ilocnum = 2 * 2 * NEX_PER_PROC_XI*NEX_PER_PROC_ETA
NIT =NGLLX-1
@@ -213,7 +217,8 @@
NIT = 1
endif
print *
- print *,'Allocating arrays for reading data of size ',ilocnum*NPROCTOT,'=',6*ilocnum*NPROCTOT*CUSTOM_REAL/1000000,'MB'
+ print *,'Allocating arrays for reading data of size ',ilocnum*NPROCTOT,'=', &
+ 6*ilocnum*NPROCTOT*CUSTOM_REAL/1000000,'MB'
print *
! allocates movie arrays
@@ -293,7 +298,8 @@
print *
- print *,'Allocating 4 outputdata arrays of size 4*CUSTOM_REAL',npointot,'=',4*npointot*CUSTOM_REAL/1000000,' MB'
+ print *,'Allocating 4 outputdata arrays of size 4*CUSTOM_REAL',npointot,'=', &
+ 4*npointot*CUSTOM_REAL/1000000,' MB'
print *
allocate(xp(npointot),stat=ierror)
@@ -309,7 +315,7 @@
if(ierror /= 0) stop 'error while allocating field_display'
- ! initializes maxima history
+ ! initializes maxima history
if( USE_AVERAGED_MAXIMUM ) then
! determines length of history
nmax_history = AVERAGE_MINIMUM + int( HDUR_MOVIE / (DT*NTSTEP_BETWEEN_FRAMES) * 1.5 )
@@ -318,63 +324,83 @@
allocate(max_history(nmax_history))
max_history(:) = 0.0d0
+ print *,'averages wavefield maxima'
print *
print *,'Movie half-duration: ',HDUR_MOVIE,'(s)'
+ print *,'DT per time step : ',DT,'(s)'
print *,'Frame step size : ',DT*NTSTEP_BETWEEN_FRAMES,'(s)'
print *,'Normalization by averaged maxima over ',nmax_history,'snapshots'
print *
+ endif
+
+ if( MUTE_SOURCE ) then
+ ! initializes
+ LAT_SOURCE = -1000.0
+ LON_SOURCE = -1000.0
+ DEP_SOURCE = 0.0
+ cmt_t_shift = 0.0
+ cmt_hdur = 0.0
+
+ ! reads in source lat/lon
+ open(22,file="DATA/CMTSOLUTION",status='old',action='read',iostat=ierror )
+ if( ierror == 0 ) then
+ ! skip first line, event name,timeshift,half duration
+ read(22,*,iostat=ierror ) line ! PDE line
+ read(22,*,iostat=ierror ) line ! event name
+ ! timeshift
+ read(22,'(a256)',iostat=ierror ) line
+ if( ierror == 0 ) read(line(12:len_trim(line)),*) cmt_t_shift
+ ! halfduration
+ read(22,'(a256)',iostat=ierror ) line
+ if( ierror == 0 ) read(line(15:len_trim(line)),*) cmt_hdur
+ ! latitude
+ read(22,'(a256)',iostat=ierror ) line
+ if( ierror == 0 ) read(line(10:len_trim(line)),*) LAT_SOURCE
+ ! longitude
+ read(22,'(a256)',iostat=ierror ) line
+ if( ierror == 0 ) read(line(11:len_trim(line)),*) LON_SOURCE
+ ! depth
+ read(22,'(a256)',iostat=ierror ) line
+ if( ierror == 0 ) read(line(7:len_trim(line)),*) DEP_SOURCE
+ close(22)
+ endif
+ ! effective half duration in movie runs
+ hdur = sqrt( cmt_hdur**2 + HDUR_MOVIE**2)
+ ! start time of simulation
+ t0 = - 1.5d0*( cmt_t_shift - hdur )
- if( MUTE_SOURCE ) then
- ! initializes
- LAT_SOURCE = -1000.0
- LON_SOURCE = -1000.0
+ ! becomes time (s) from hypocenter to reach surface (using average 8 km/s s-wave speed)
+ ! note: especially for deep sources, this helps determine a better starttime to mute
+ DEP_SOURCE = DEP_SOURCE / 8.0
- ! reads in source lat/lon
- open(22,file="DATA/CMTSOLUTION",status='old',action='read',iostat=ierror )
- if( ierror == 0 ) then
- ! skip first line, event name,timeshift,half duration
- read(22,*,iostat=ierror ) line ! PDE line
- read(22,*,iostat=ierror ) line ! event name
- read(22,*,iostat=ierror ) line ! timeshift
- read(22,*,iostat=ierror ) line ! halfduration
- ! latitude
- read(22,'(a256)',iostat=ierror ) line
- if( ierror == 0 ) read(line(10:len_trim(line)),*) LAT_SOURCE
- ! longitude
- read(22,'(a256)',iostat=ierror ) line
- if( ierror == 0 ) read(line(11:len_trim(line)),*) LON_SOURCE
- ! depth
- read(22,'(a256)',iostat=ierror ) line
- if( ierror == 0 ) read(line(11:len_trim(line)),*) DEP_SOURCE
- close(22)
- endif
+ ! time when muting starts
+ ! note: this starttime is supposed to be the time when displacements at the surface
+ ! can be observed;
+ ! it helps to mute out numerical noise before the source effects actually start showing up
+ STARTTIME_TO_MUTE = STARTTIME_TO_MUTE + DEP_SOURCE
+ if( STARTTIME_TO_MUTE < 0.0 ) STARTTIME_TO_MUTE = 0.0
+
+ print *,'mutes source area'
+ print *
+ print *,'source lat/lon/dep: ',LAT_SOURCE,LON_SOURCE,DEP_SOURCE
+ print *,'muting radius: ',RADIUS_TO_MUTE,'(degrees)'
+ print *,'muting starttime: ',STARTTIME_TO_MUTE,'(s)'
+ print *,'simulation starttime: ',-t0,'(s)'
+ print *
- print *,'muting source lat/lon/dep: ',LAT_SOURCE,LON_SOURCE,DEP_SOURCE
+ ! converts values into radians
+ ! colatitude [0, PI]
+ LAT_SOURCE = (90. - LAT_SOURCE)*PI/180.0
- ! becomes time (s) from hypocenter to reach surface (using average 8 km/s p-wave speed)
- DEP_SOURCE = DEP_SOURCE / 8.0
+ ! longitude [-PI, PI]
+ if( LON_SOURCE < -180.0 ) LON_SOURCE = LON_SOURCE + 360.0
+ if( LON_SOURCE > 180.0 ) LON_SOURCE = LON_SOURCE - 360.0
+ LON_SOURCE = LON_SOURCE *PI/180.0
- ! time when muting starts
- STARTTIME_TO_MUTE = STARTTIME_TO_MUTE * HDUR_MOVIE + DEP_SOURCE
+ ! mute radius in rad
+ RADIUS_TO_MUTE = RADIUS_TO_MUTE*PI/180.0
+ endif
- print *,'muting radius: ',RADIUS_TO_MUTE
- print *,'muting starttime: ',STARTTIME_TO_MUTE,'(s)'
- print *
-
- ! colatitude [0, PI]
- LAT_SOURCE = (90. - LAT_SOURCE)*PI/180.0
-
- ! longitude [-PI, PI]
- if( LON_SOURCE < -180.0 ) LON_SOURCE = LON_SOURCE + 360.0
- if( LON_SOURCE > 180.0 ) LON_SOURCE = LON_SOURCE - 360.0
- LON_SOURCE = LON_SOURCE *PI/180.0
-
- ! mute radius in rad
- RADIUS_TO_MUTE = RADIUS_TO_MUTE*PI/180.0
- endif
-
-
- endif
print *,'--------'
!--- ****** read data saved by solver ******
@@ -390,46 +416,6 @@
iframe = iframe + 1
- ! mutes source region
- if( MUTE_SOURCE ) then
-
- ! muting radius grows/shrinks with time
- if( (it-1)*DT > STARTTIME_TO_MUTE ) then
-
- ! approximate wavefront travel distance in degrees (~3.5 km/s wave speed for surface waves)
- mute_factor = 3.5 * (it-1)*DT / 6371. * 180./PI
-
- ! approximate distance to source (in degrees)
- do while ( mute_factor > 360. )
- mute_factor = mute_factor - 360.
- enddo
- if( mute_factor > 180. ) mute_factor = 360. - mute_factor
-
- ! limit size around source (in degrees)
- !if( mute_factor < 10. ) then
- ! mute_factor = 0.0
- !endif
- if( mute_factor > 80. ) then
- mute_factor = 80.0
- endif
-
- print*,'muting radius: ',0.7 * mute_factor
-
- RADIUS_TO_MUTE = 0.7 * mute_factor * PI/180.
-
- else
- ! mute_factor used at the beginning for scaling displacement values
- if( STARTTIME_TO_MUTE > TINYVAL ) then
- ! scales from 1 to 0
- mute_factor = ( STARTTIME_TO_MUTE - (it-1)*DT ) / STARTTIME_TO_MUTE
- if( mute_factor < TINYVAL ) mute_factor = TINYVAL
- else
- mute_factor = 1.0
- endif
- endif
-
- endif
-
! read all the elements from the same file
write(outputname,"('OUTPUT_FILES/moviedata',i6.6)") it
open(unit=IOUT,file=outputname,status='old',form='unformatted')
@@ -452,6 +438,56 @@
close(IOUT)
!print *, 'finished reading ',outputname
+ ! mutes source region
+ if( MUTE_SOURCE ) then
+ ! initialize factor
+ mute_factor = 1.0
+
+ print*,'simulation time: ',(it-1)*DT - t0,'(s)'
+
+ ! muting radius grows/shrinks with time
+ if( (it-1)*DT - t0 > STARTTIME_TO_MUTE ) then
+
+ ! approximate wavefront travel distance in degrees
+ ! (~3.5 km/s wave speed for surface waves)
+ distance = 3.5 * ((it-1)*DT-t0) / 6371.0 * 180./PI
+
+ ! approximate distance to source (in degrees)
+ ! (shrinks if waves travel back from antipode)
+ !do while ( distance > 360. )
+ ! distance = distance - 360.
+ !enddo
+ ! waves are back at origin, no source tapering anymore
+ if( distance > 360.0 ) distance = 0.0
+ ! shrinks when waves reached antipode
+ !if( distance > 180. ) distance = 360. - distance
+ ! shrinks when waves reached half-way to antipode
+ if( distance > 90.0 ) distance = 90.0 - distance
+
+ ! limit size around source (in degrees)
+ if( distance < 0.0 ) distance = 0.0
+ if( distance > 80.0 ) distance = 80.0
+
+ print*,'muting radius: ',0.7 * distance,'(degrees)'
+
+ ! new radius of mute area (in rad)
+ RADIUS_TO_MUTE = 0.7 * distance * PI/180.
+ else
+ ! mute_factor used at the beginning for scaling displacement values
+ if( STARTTIME_TO_MUTE > TINYVAL ) then
+ ! mute factor 1: no masking out
+ ! 0: masks out values (within mute radius)
+ ! linear scaling between [0,1]:
+ ! from 0 (simulation time equal to zero )
+ ! to 1 (simulation time equals starttime_to_mute)
+ mute_factor = 1.0 - ( STARTTIME_TO_MUTE - ((it-1)*DT-t0) ) / (STARTTIME_TO_MUTE+t0)
+ ! threshold value for mute_factor
+ if( mute_factor < TINYVAL ) mute_factor = TINYVAL
+ if( mute_factor > 1.0 ) mute_factor = 1.0
+ endif
+ endif
+ endif
+
! clear number of elements kept
ispec = 0
@@ -519,12 +555,13 @@
! mute values
if( MUTE_SOURCE ) then
- ! distance in colatitude
+ ! distance in colatitude (in rad)
! note: this mixes geocentric (point location) and geographic (source location) coordinates;
- ! since we only need approximate distances here, this should be fine for the muting region
+ ! since we only need approximate distances here,
+ ! this should be fine for the muting region
dist_lat = thetaval - LAT_SOURCE
- ! distance in longitude
+ ! distance in longitude (in rad)
! checks source longitude range
if( LON_SOURCE - RADIUS_TO_MUTE < -PI .or. LON_SOURCE + RADIUS_TO_MUTE > PI ) then
! source close to 180. longitudes, shifts range to [0, 2PI]
@@ -542,13 +579,21 @@
dist_lon = phival - LON_SOURCE
endif
-
+ ! distance of point to source (in rad)
+ distance = sqrt(dist_lat**2 + dist_lon**2)
+
! mutes source region values
- if ( ( dist_lat**2 + dist_lon**2 ) < RADIUS_TO_MUTE**2 ) then
+ if ( distance < RADIUS_TO_MUTE ) then
! muting takes account of the event time
- if( (it-1)*DT > STARTTIME_TO_MUTE ) then
- displn(i,j) = displn(i,j) * TINYVAL
+ if( (it-1)*DT-t0 > STARTTIME_TO_MUTE ) then
+ ! wavefield will be tapered to mask out noise in source area
+ ! factor from 0 to 1
+ mute_factor = ( 0.5*(1.0 - cos(distance/RADIUS_TO_MUTE*PI)) )**6
+ ! factor from 0.01 to 1
+ mute_factor = mute_factor * 0.99 + 0.01
+ displn(i,j) = displn(i,j) * mute_factor
else
+ ! wavefield will initially be scaled down to avoid noise being amplified at beginning
displn(i,j) = displn(i,j) * mute_factor
endif
endif
@@ -642,6 +687,7 @@
print *,'minimum amplitude in current snapshot = ',min_field_current
print *,'maximum amplitude in current snapshot = ',max_field_current
+
! takes average over last few snapshots available and uses it
! to normalize field values
if( USE_AVERAGED_MAXIMUM ) then
@@ -666,19 +712,61 @@
print *,'maximum amplitude over averaged last snapshots = ',max_average
+ ! thresholds positive & negative maximum values
+ where( field_display(:) > max_average ) field_display = max_average
+ where( field_display(:) < - max_average ) field_display = -max_average
+
+ ! updates current wavefield maxima
+ min_field_current = minval(field_display(:))
+ max_field_current = maxval(field_display(:))
+ max_absol = (abs(min_field_current)+abs(max_field_current))/2.0
+
! scales field values up to match average
if( abs(max_absol) > TINYVAL) &
field_display = field_display * max_average / max_absol
- ! thresholds positive & negative maximum values
- where( field_display(:) > max_average ) field_display = max_average
- where( field_display(:) < - max_average ) field_display = -max_average
-
! normalizes field values
if( NORMALIZE_VALUES ) then
- if( abs(max_average) > TINYVAL ) field_display = field_display / max_average
+ if( MUTE_SOURCE ) then
+ ! checks if source wavefield kicked in
+ if( (it-1)*DT - t0 > STARTTIME_TO_MUTE ) then
+ ! wavefield should be visible at surface now
+ ! normalizes wavefield
+ if( abs(max_average) > TINYVAL ) field_display = field_display / max_average
+ else
+ ! no wavefield yet assumed
+
+ ! we set two single field values (last in array)
+ ! to: +/- 100 * max_average
+ ! to avoid further amplifying when
+ ! a normalization routine is used for rendering images
+ ! (which for example is the case for shakemovies)
+ if( STARTTIME_TO_MUTE > TINYVAL ) then
+ ! with additional scale factor:
+ ! linear scaling between [0,1]:
+ ! from 0 (simulation time equal to -t0 )
+ ! to 1 (simulation time equals starttime_to_mute)
+ mute_factor = 1.0 - ( STARTTIME_TO_MUTE - ((it-1)*DT-t0) ) / (STARTTIME_TO_MUTE+t0)
+ ! takes complement and shifts scale to (1,100)
+ ! thus, mute factor is 100 at simulation start and 1.0 at starttime_to_mute
+ mute_factor = abs(1.0 - mute_factor) * 99.0 + 1.0
+ ! positive and negative maximum reach average when wavefield appears
+ val = mute_factor * max_average
+ else
+ ! uses a constant factor
+ val = 100.0 * max_average
+ endif
+ ! positive and negative maximum
+ field_display(ieoff) = + val
+ field_display(ieoff-1) = - val
+ if( abs(max_average) > TINYVAL ) field_display = field_display / val
+ endif
+ else
+ ! no source to mute
+ ! normalizes wavefield
+ if( abs(max_average) > TINYVAL ) field_display = field_display / max_average
+ endif
endif
-
endif
print *
@@ -696,8 +784,10 @@
elseif(USE_COMPONENT == 3) then
write(outputname,"('bin_movie_',i6.6,'.E')") it
endif
- open(unit=11,file='OUTPUT_FILES/'//trim(outputname),status='unknown',form='unformatted')
- if(iframe == 1) open(unit=12,file='OUTPUT_FILES/bin_movie.xy',status='unknown',form='unformatted')
+ open(unit=11,file='OUTPUT_FILES/'//trim(outputname),status='unknown', &
+ form='unformatted',action='write')
+ if(iframe == 1) open(unit=12,file='OUTPUT_FILES/bin_movie.xy', &
+ status='unknown',form='unformatted',action='write')
else
if(USE_COMPONENT == 1) then
write(outputname,"('ascii_movie_',i6.6,'.d')") it
@@ -706,8 +796,10 @@
elseif(USE_COMPONENT == 3) then
write(outputname,"('ascii_movie_',i6.6,'.E')") it
endif
- open(unit=11,file='OUTPUT_FILES/'//trim(outputname),status='unknown')
- if(iframe == 1) open(unit=12,file='OUTPUT_FILES/ascii_movie.xy',status='unknown')
+ open(unit=11,file='OUTPUT_FILES/'//trim(outputname),status='unknown', &
+ action='write')
+ if(iframe == 1) open(unit=12,file='OUTPUT_FILES/ascii_movie.xy', &
+ status='unknown',action='write')
endif
! clear number of elements kept
ispec = 0
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/Makefile.in 2011-03-15 16:48:14 UTC (rev 18115)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/Makefile.in 2011-03-16 02:28:49 UTC (rev 18116)
@@ -78,6 +78,7 @@
$O/calendar.o \
$O/comp_source_spectrum.o \
$O/comp_source_time_function.o \
+ $O/compute_adj_source_frechet.o \
$O/compute_arrays_source.o \
$O/convert_time.o \
$O/create_central_cube_buffers.o \
@@ -363,9 +364,11 @@
$O/comp_source_spectrum.o: ${SETUP}/constants.h $S/comp_source_spectrum.f90
${FCCOMPILE_CHECK} -c -o $O/comp_source_spectrum.o ${FCFLAGS_f90} $S/comp_source_spectrum.f90
-## DK DK added -O1 here for now, because compiling compute_arrays_source.f90 with -O3 currently takes forever...
+$O/compute_adj_source_frechet.o: ${SETUP}/constants.h $S/compute_adj_source_frechet.f90
+ ${FCCOMPILE_CHECK} -c -o $O/compute_adj_source_frechet.o ${FCFLAGS_f90} $S/compute_adj_source_frechet.f90
+
$O/compute_arrays_source.o: ${SETUP}/constants.h $S/compute_arrays_source.f90
- ${FCCOMPILE_CHECK} -O1 -c -o $O/compute_arrays_source.o ${FCFLAGS_f90} $S/compute_arrays_source.f90
+ ${FCCOMPILE_CHECK} -c -o $O/compute_arrays_source.o ${FCFLAGS_f90} $S/compute_arrays_source.f90
$O/compute_boundary_kernel.o: ${SETUP}/constants.h $S/compute_boundary_kernel.f90
${FCCOMPILE_CHECK} -c -o $O/compute_boundary_kernel.o ${FCFLAGS_f90} $S/compute_boundary_kernel.f90
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_adj_source_frechet.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_adj_source_frechet.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_adj_source_frechet.f90 2011-03-16 02:28:49 UTC (rev 18116)
@@ -0,0 +1,209 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 5 . 1
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Princeton University, USA
+! and University of Pau / CNRS / INRIA, France
+! (c) Princeton University / California Institute of Technology and University of Pau / CNRS / INRIA
+! February 2011
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+
+! compute the integrated derivatives of source parameters (M_jk and X_s)
+
+ subroutine compute_adj_source_frechet(displ_s,Mxx,Myy,Mzz,Mxy,Mxz,Myz, &
+ eps_s,eps_m_s,eps_m_l_s, &
+ hxir,hetar,hgammar,hpxir,hpetar,hpgammar, hprime_xx,hprime_yy,hprime_zz, &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
+
+ implicit none
+
+ include 'constants.h'
+
+ ! input
+ real(kind=CUSTOM_REAL) :: displ_s(NDIM,NGLLX,NGLLY,NGLLZ)
+ double precision :: Mxx, Myy, Mzz, Mxy, Mxz, Myz
+ ! output
+ real(kind=CUSTOM_REAL) :: eps_s(NDIM,NDIM), eps_m_s, eps_m_l_s(NDIM)
+
+ ! auxilliary
+ double precision :: hxir(NGLLX), hetar(NGLLY), hgammar(NGLLZ), &
+ hpxir(NGLLX),hpetar(NGLLY),hpgammar(NGLLZ)
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLY) :: hprime_yy
+ real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
+ xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+
+! local variables
+ real(kind=CUSTOM_REAL) :: tempx1l,tempx2l,tempx3l, tempy1l,tempy2l,tempy3l, &
+ tempz1l,tempz2l,tempz3l, hp1, hp2, hp3, &
+ xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl, &
+ duxdxl,duxdyl,duxdzl,duydxl,duydyl,duydzl,duzdxl,duzdyl,duzdzl, &
+ xix_s,xiy_s,xiz_s,etax_s,etay_s,etaz_s,gammax_s,gammay_s,gammaz_s, &
+ hlagrange_xi, hlagrange_eta, hlagrange_gamma, hlagrange
+
+ real(kind=CUSTOM_REAL) :: eps(NDIM,NDIM), eps_array(NDIM,NDIM,NGLLX,NGLLY,NGLLZ), &
+ eps_m_array(NGLLX,NGLLY,NGLLZ)
+
+ integer i,j,k,l
+
+
+! first compute the strain at all the GLL points of the source element
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+
+ tempx1l = 0._CUSTOM_REAL
+ tempx2l = 0._CUSTOM_REAL
+ tempx3l = 0._CUSTOM_REAL
+
+ tempy1l = 0._CUSTOM_REAL
+ tempy2l = 0._CUSTOM_REAL
+ tempy3l = 0._CUSTOM_REAL
+
+ tempz1l = 0._CUSTOM_REAL
+ tempz2l = 0._CUSTOM_REAL
+ tempz3l = 0._CUSTOM_REAL
+
+ do l=1,NGLLX
+ hp1 = hprime_xx(i,l)
+ tempx1l = tempx1l + displ_s(1,l,j,k)*hp1
+ tempy1l = tempy1l + displ_s(2,l,j,k)*hp1
+ tempz1l = tempz1l + displ_s(3,l,j,k)*hp1
+
+ hp2 = hprime_yy(j,l)
+ tempx2l = tempx2l + displ_s(1,i,l,k)*hp2
+ tempy2l = tempy2l + displ_s(2,i,l,k)*hp2
+ tempz2l = tempz2l + displ_s(3,i,l,k)*hp2
+
+ hp3 = hprime_zz(k,l)
+ tempx3l = tempx3l + displ_s(1,i,j,l)*hp3
+ tempy3l = tempy3l + displ_s(2,i,j,l)*hp3
+ tempz3l = tempz3l + displ_s(3,i,j,l)*hp3
+ enddo
+
+! dudx
+ xixl = xix(i,j,k)
+ xiyl = xiy(i,j,k)
+ xizl = xiz(i,j,k)
+ etaxl = etax(i,j,k)
+ etayl = etay(i,j,k)
+ etazl = etaz(i,j,k)
+ gammaxl = gammax(i,j,k)
+ gammayl = gammay(i,j,k)
+ gammazl = gammaz(i,j,k)
+
+ duxdxl = xixl*tempx1l + etaxl*tempx2l + gammaxl*tempx3l
+ duxdyl = xiyl*tempx1l + etayl*tempx2l + gammayl*tempx3l
+ duxdzl = xizl*tempx1l + etazl*tempx2l + gammazl*tempx3l
+
+ duydxl = xixl*tempy1l + etaxl*tempy2l + gammaxl*tempy3l
+ duydyl = xiyl*tempy1l + etayl*tempy2l + gammayl*tempy3l
+ duydzl = xizl*tempy1l + etazl*tempy2l + gammazl*tempy3l
+
+ duzdxl = xixl*tempz1l + etaxl*tempz2l + gammaxl*tempz3l
+ duzdyl = xiyl*tempz1l + etayl*tempz2l + gammayl*tempz3l
+ duzdzl = xizl*tempz1l + etazl*tempz2l + gammazl*tempz3l
+
+! strain eps_jk
+ eps(1,1) = duxdxl
+ eps(1,2) = (duxdyl + duydxl) / 2
+ eps(1,3) = (duxdzl + duzdxl) / 2
+ eps(2,2) = duydyl
+ eps(2,3) = (duydzl + duzdyl) / 2
+ eps(3,3) = duzdzl
+ eps(2,1) = eps(1,2)
+ eps(3,1) = eps(1,3)
+ eps(3,2) = eps(2,3)
+
+ eps_array(:,:,i,j,k) = eps(:,:)
+
+! Mjk eps_jk
+ eps_m_array(i,j,k) = Mxx * eps(1,1) + Myy * eps(2,2) + Mzz * eps(3,3) + &
+ 2 * (Mxy * eps(1,2) + Mxz * eps(1,3) + Myz * eps(2,3))
+
+ enddo
+ enddo
+ enddo
+
+ ! interpolate the strain eps_s(:,:) from eps_array(:,:,i,j,k)
+ eps_s = 0.; eps_m_s=0.;
+ xix_s = 0.; xiy_s = 0.; xiz_s = 0.
+ etax_s = 0.; etay_s = 0.; etaz_s = 0.
+ gammax_s = 0.; gammay_s = 0.; gammaz_s = 0.
+
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
+ hlagrange = hxir(i)*hetar(j)*hgammar(k)
+
+ eps_s(1,1) = eps_s(1,1) + eps_array(1,1,i,j,k)*hlagrange
+ eps_s(1,2) = eps_s(1,2) + eps_array(1,2,i,j,k)*hlagrange
+ eps_s(1,3) = eps_s(1,3) + eps_array(1,3,i,j,k)*hlagrange
+ eps_s(2,2) = eps_s(2,2) + eps_array(2,2,i,j,k)*hlagrange
+ eps_s(2,3) = eps_s(2,3) + eps_array(2,3,i,j,k)*hlagrange
+ eps_s(3,3) = eps_s(3,3) + eps_array(3,3,i,j,k)*hlagrange
+
+ xix_s = xix_s + xix(i,j,k)*hlagrange
+ xiy_s = xiy_s + xiy(i,j,k)*hlagrange
+ xiz_s = xiz_s + xiz(i,j,k)*hlagrange
+ etax_s = etax_s + etax(i,j,k)*hlagrange
+ etay_s = etay_s + etay(i,j,k)*hlagrange
+ etaz_s = etaz_s + etaz(i,j,k)*hlagrange
+ gammax_s = gammax_s + gammax(i,j,k)*hlagrange
+ gammay_s = gammay_s + gammay(i,j,k)*hlagrange
+ gammaz_s = gammaz_s + gammaz(i,j,k)*hlagrange
+
+ eps_m_s = eps_m_s + eps_m_array(i,j,k)*hlagrange
+ enddo
+ enddo
+ enddo
+
+! for completion purpose, not used in specfem3D.f90
+ eps_s(2,1) = eps_s(1,2)
+ eps_s(3,1) = eps_s(1,3)
+ eps_s(3,2) = eps_s(2,3)
+
+! compute the gradient of M_jk * eps_jk, and then interpolate it
+
+ eps_m_l_s = 0.
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
+ hlagrange_xi = hpxir(i)*hetar(j)*hgammar(k)
+ hlagrange_eta = hxir(i)*hpetar(j)*hgammar(k)
+ hlagrange_gamma = hxir(i)*hetar(j)*hpgammar(k)
+
+ eps_m_l_s(1) = eps_m_l_s(1) + eps_m_array(i,j,k) * (hlagrange_xi * xix_s &
+ + hlagrange_eta * etax_s + hlagrange_gamma * gammax_s)
+ eps_m_l_s(2) = eps_m_l_s(2) + eps_m_array(i,j,k) * (hlagrange_xi * xiy_s &
+ + hlagrange_eta * etay_s + hlagrange_gamma * gammay_s)
+ eps_m_l_s(3) = eps_m_l_s(3) + eps_m_array(i,j,k) * (hlagrange_xi * xiz_s &
+ + hlagrange_eta * etaz_s + hlagrange_gamma * gammaz_s)
+
+ enddo
+ enddo
+ enddo
+
+ end subroutine compute_adj_source_frechet
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_arrays_source.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_arrays_source.f90 2011-03-15 16:48:14 UTC (rev 18115)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_arrays_source.f90 2011-03-16 02:28:49 UTC (rev 18116)
@@ -166,7 +166,7 @@
!================================================================
-subroutine compute_arrays_source_adjoint(myrank, adj_source_file, &
+ subroutine compute_arrays_source_adjoint(myrank, adj_source_file, &
xi_receiver,eta_receiver,gamma_receiver, nu,adj_sourcearray, &
xigll,yigll,zigll,NSTEP_BLOCK,iadjsrc,it_sub_adj,NSTEP_SUB_ADJ, &
NTSTEP_BETWEEN_READ_ADJSRC,DT)
@@ -203,12 +203,14 @@
double precision :: hxir(NGLLX), hpxir(NGLLX), hetar(NGLLY), hpetar(NGLLY), &
hgammar(NGLLZ), hpgammar(NGLLZ)
- real(kind=CUSTOM_REAL) :: adj_src(NDIM,NSTEP_BLOCK),adj_src_u(NDIM,NSTEP_BLOCK)
+ double precision, dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: sourcearrayd
+
+ real(kind=CUSTOM_REAL) :: adj_src(NDIM,NSTEP_BLOCK)
+ double precision, dimension(NDIM,NSTEP_BLOCK) :: adj_src_u
integer icomp, itime, i, j, k, ios
integer it_start,it_end,index_i
real(kind=CUSTOM_REAL) :: junk
- !character(len=3),dimension(NDIM) :: comp = (/ "LHN", "LHE", "LHZ" /)
character(len=3),dimension(NDIM) :: comp
character(len=150) :: filename
character(len=2) :: bic
@@ -309,204 +311,64 @@
call lagrange_any(gamma_receiver,NGLLZ,zigll,hgammar,hpgammar)
! adds interpolated source contribution to all GLL points within this element
- do k = 1, NGLLZ
- do j = 1, NGLLY
- do i = 1, NGLLX
- do itime = 1, NSTEP_BLOCK
- adj_sourcearray(:,i,j,k,itime) = hxir(i) * hetar(j) * hgammar(k) * adj_src_u(:,itime)
- enddo
- enddo
- enddo
+ do itime = 1, NSTEP_BLOCK
+
+ ! multiply with interpolators
+ call multiply_arrays_adjoint(sourcearrayd,hxir,hetar,hgammar,adj_src_u(:,itime))
+
+ ! distinguish between single and double precision for reals
+ if(CUSTOM_REAL == SIZE_REAL) then
+ adj_sourcearray(:,:,:,:,itime) = sngl(sourcearrayd(:,:,:,:))
+ else
+ adj_sourcearray(:,:,:,:,itime) = sourcearrayd(:,:,:,:)
+ endif
+
enddo
+! do k = 1, NGLLZ
+! do j = 1, NGLLY
+! do i = 1, NGLLX
+! do itime = 1, NSTEP_BLOCK
+! adj_sourcearray(:,i,j,k,itime) = hxir(i) * hetar(j) * hgammar(k) * adj_src_u(:,itime)
+! enddo
+! enddo
+! enddo
+! enddo
-end subroutine compute_arrays_source_adjoint
+ end subroutine compute_arrays_source_adjoint
! =======================================================================
-! compute the integrated derivatives of source parameters (M_jk and X_s)
+! we put these multiplications in a separate routine because otherwise
+! some compilers try to unroll the four loops above and take forever to compile
+ subroutine multiply_arrays_adjoint(sourcearrayd,hxir,hetar,hgammar,adj_src_ud)
-subroutine compute_adj_source_frechet(displ_s,Mxx,Myy,Mzz,Mxy,Mxz,Myz, &
- eps_s,eps_m_s,eps_m_l_s, &
- hxir,hetar,hgammar,hpxir,hpetar,hpgammar, hprime_xx,hprime_yy,hprime_zz, &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
-
implicit none
- include 'constants.h'
+ include "constants.h"
- ! input
- real(kind=CUSTOM_REAL) :: displ_s(NDIM,NGLLX,NGLLY,NGLLZ)
- double precision :: Mxx, Myy, Mzz, Mxy, Mxz, Myz
- ! output
- real(kind=CUSTOM_REAL) :: eps_s(NDIM,NDIM), eps_m_s, eps_m_l_s(NDIM)
-
- ! auxilliary
- double precision :: hxir(NGLLX), hetar(NGLLY), hgammar(NGLLZ), &
- hpxir(NGLLX),hpetar(NGLLY),hpgammar(NGLLZ)
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx
- real(kind=CUSTOM_REAL), dimension(NGLLY,NGLLY) :: hprime_yy
- real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz
-
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: &
- xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
-
-! local variables
- real(kind=CUSTOM_REAL) :: tempx1l,tempx2l,tempx3l, tempy1l,tempy2l,tempy3l, &
- tempz1l,tempz2l,tempz3l, hp1, hp2, hp3, &
- xixl,xiyl,xizl,etaxl,etayl,etazl,gammaxl,gammayl,gammazl, &
- duxdxl,duxdyl,duxdzl,duydxl,duydyl,duydzl,duzdxl,duzdyl,duzdzl, &
- xix_s,xiy_s,xiz_s,etax_s,etay_s,etaz_s,gammax_s,gammay_s,gammaz_s, &
- hlagrange_xi, hlagrange_eta, hlagrange_gamma, hlagrange
-
- real(kind=CUSTOM_REAL) :: eps(NDIM,NDIM), eps_array(NDIM,NDIM,NGLLX,NGLLY,NGLLZ), &
- eps_m_array(NGLLX,NGLLY,NGLLZ)
-
- integer i,j,k,l
-
-
-! first compute the strain at all the GLL points of the source element
+ double precision, dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: sourcearrayd
+ double precision, dimension(NGLLX) :: hxir
+ double precision, dimension(NGLLY) :: hetar
+ double precision, dimension(NGLLZ) :: hgammar
+ double precision, dimension(NDIM) :: adj_src_ud
+
+ integer :: i,j,k
+
+ ! adds interpolated source contribution to all GLL points within this element
do k = 1, NGLLZ
do j = 1, NGLLY
do i = 1, NGLLX
-
- tempx1l = 0._CUSTOM_REAL
- tempx2l = 0._CUSTOM_REAL
- tempx3l = 0._CUSTOM_REAL
-
- tempy1l = 0._CUSTOM_REAL
- tempy2l = 0._CUSTOM_REAL
- tempy3l = 0._CUSTOM_REAL
-
- tempz1l = 0._CUSTOM_REAL
- tempz2l = 0._CUSTOM_REAL
- tempz3l = 0._CUSTOM_REAL
-
- do l=1,NGLLX
- hp1 = hprime_xx(i,l)
- tempx1l = tempx1l + displ_s(1,l,j,k)*hp1
- tempy1l = tempy1l + displ_s(2,l,j,k)*hp1
- tempz1l = tempz1l + displ_s(3,l,j,k)*hp1
-
- hp2 = hprime_yy(j,l)
- tempx2l = tempx2l + displ_s(1,i,l,k)*hp2
- tempy2l = tempy2l + displ_s(2,i,l,k)*hp2
- tempz2l = tempz2l + displ_s(3,i,l,k)*hp2
-
- hp3 = hprime_zz(k,l)
- tempx3l = tempx3l + displ_s(1,i,j,l)*hp3
- tempy3l = tempy3l + displ_s(2,i,j,l)*hp3
- tempz3l = tempz3l + displ_s(3,i,j,l)*hp3
- enddo
-
-! dudx
- xixl = xix(i,j,k)
- xiyl = xiy(i,j,k)
- xizl = xiz(i,j,k)
- etaxl = etax(i,j,k)
- etayl = etay(i,j,k)
- etazl = etaz(i,j,k)
- gammaxl = gammax(i,j,k)
- gammayl = gammay(i,j,k)
- gammazl = gammaz(i,j,k)
-
- duxdxl = xixl*tempx1l + etaxl*tempx2l + gammaxl*tempx3l
- duxdyl = xiyl*tempx1l + etayl*tempx2l + gammayl*tempx3l
- duxdzl = xizl*tempx1l + etazl*tempx2l + gammazl*tempx3l
-
- duydxl = xixl*tempy1l + etaxl*tempy2l + gammaxl*tempy3l
- duydyl = xiyl*tempy1l + etayl*tempy2l + gammayl*tempy3l
- duydzl = xizl*tempy1l + etazl*tempy2l + gammazl*tempy3l
-
- duzdxl = xixl*tempz1l + etaxl*tempz2l + gammaxl*tempz3l
- duzdyl = xiyl*tempz1l + etayl*tempz2l + gammayl*tempz3l
- duzdzl = xizl*tempz1l + etazl*tempz2l + gammazl*tempz3l
-
-! strain eps_jk
- eps(1,1) = duxdxl
- eps(1,2) = (duxdyl + duydxl) / 2
- eps(1,3) = (duxdzl + duzdxl) / 2
- eps(2,2) = duydyl
- eps(2,3) = (duydzl + duzdyl) / 2
- eps(3,3) = duzdzl
- eps(2,1) = eps(1,2)
- eps(3,1) = eps(1,3)
- eps(3,2) = eps(2,3)
-
- eps_array(:,:,i,j,k) = eps(:,:)
-
-! Mjk eps_jk
- eps_m_array(i,j,k) = Mxx * eps(1,1) + Myy * eps(2,2) + Mzz * eps(3,3) + &
- 2 * (Mxy * eps(1,2) + Mxz * eps(1,3) + Myz * eps(2,3))
-
+ sourcearrayd(:,i,j,k) = hxir(i) * hetar(j) * hgammar(k) * adj_src_ud(:)
enddo
enddo
enddo
- ! interpolate the strain eps_s(:,:) from eps_array(:,:,i,j,k)
- eps_s = 0.; eps_m_s=0.;
- xix_s = 0.; xiy_s = 0.; xiz_s = 0.
- etax_s = 0.; etay_s = 0.; etaz_s = 0.
- gammax_s = 0.; gammay_s = 0.; gammaz_s = 0.
+ end subroutine multiply_arrays_adjoint
- do k = 1,NGLLZ
- do j = 1,NGLLY
- do i = 1,NGLLX
- hlagrange = hxir(i)*hetar(j)*hgammar(k)
- eps_s(1,1) = eps_s(1,1) + eps_array(1,1,i,j,k)*hlagrange
- eps_s(1,2) = eps_s(1,2) + eps_array(1,2,i,j,k)*hlagrange
- eps_s(1,3) = eps_s(1,3) + eps_array(1,3,i,j,k)*hlagrange
- eps_s(2,2) = eps_s(2,2) + eps_array(2,2,i,j,k)*hlagrange
- eps_s(2,3) = eps_s(2,3) + eps_array(2,3,i,j,k)*hlagrange
- eps_s(3,3) = eps_s(3,3) + eps_array(3,3,i,j,k)*hlagrange
-
- xix_s = xix_s + xix(i,j,k)*hlagrange
- xiy_s = xiy_s + xiy(i,j,k)*hlagrange
- xiz_s = xiz_s + xiz(i,j,k)*hlagrange
- etax_s = etax_s + etax(i,j,k)*hlagrange
- etay_s = etay_s + etay(i,j,k)*hlagrange
- etaz_s = etaz_s + etaz(i,j,k)*hlagrange
- gammax_s = gammax_s + gammax(i,j,k)*hlagrange
- gammay_s = gammay_s + gammay(i,j,k)*hlagrange
- gammaz_s = gammaz_s + gammaz(i,j,k)*hlagrange
-
- eps_m_s = eps_m_s + eps_m_array(i,j,k)*hlagrange
- enddo
- enddo
- enddo
-
-! for completion purpose, not used in specfem3D.f90
- eps_s(2,1) = eps_s(1,2)
- eps_s(3,1) = eps_s(1,3)
- eps_s(3,2) = eps_s(2,3)
-
-! compute the gradient of M_jk * eps_jk, and then interpolate it
-
- eps_m_l_s = 0.
- do k = 1,NGLLZ
- do j = 1,NGLLY
- do i = 1,NGLLX
-
- hlagrange_xi = hpxir(i)*hetar(j)*hgammar(k)
- hlagrange_eta = hxir(i)*hpetar(j)*hgammar(k)
- hlagrange_gamma = hxir(i)*hetar(j)*hpgammar(k)
-
- eps_m_l_s(1) = eps_m_l_s(1) + eps_m_array(i,j,k) * (hlagrange_xi * xix_s &
- + hlagrange_eta * etax_s + hlagrange_gamma * gammax_s)
- eps_m_l_s(2) = eps_m_l_s(2) + eps_m_array(i,j,k) * (hlagrange_xi * xiy_s &
- + hlagrange_eta * etay_s + hlagrange_gamma * gammay_s)
- eps_m_l_s(3) = eps_m_l_s(3) + eps_m_array(i,j,k) * (hlagrange_xi * xiz_s &
- + hlagrange_eta * etaz_s + hlagrange_gamma * gammaz_s)
-
- enddo
- enddo
- enddo
-
-end subroutine compute_adj_source_frechet
-
-!================================================================
+! =======================================================================
!
! deprecated...
!
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.f90 2011-03-15 16:48:14 UTC (rev 18115)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.f90 2011-03-16 02:28:49 UTC (rev 18116)
@@ -1697,7 +1697,6 @@
rho_kl_outer_core(:,:,:,:) = 0._CUSTOM_REAL
alpha_kl_outer_core(:,:,:,:) = 0._CUSTOM_REAL
- beta_kl_outer_core(:,:,:,:) = 0._CUSTOM_REAL
rho_kl_inner_core(:,:,:,:) = 0._CUSTOM_REAL
beta_kl_inner_core(:,:,:,:) = 0._CUSTOM_REAL
@@ -1713,7 +1712,7 @@
nspec_beta_kl_outer_core = 1
endif
allocate(beta_kl_outer_core(NGLLX,NGLLY,NGLLZ,nspec_beta_kl_outer_core))
- beta_kl_outer_core = 0._CUSTOM_REAL
+ beta_kl_outer_core(:,:,:,:) = 0._CUSTOM_REAL
endif
! initialize to be on the save side for adjoint runs SIMULATION_TYPE==2
More information about the CIG-COMMITS
mailing list