[cig-commits] [commit] devel, master: removed src/auxiliaries/combine_vol_data_shared.f90, which contained routines duplicated from other files of the code (that is dangerous, and very uneady to maintain the same routines twice, in two different files of the code). Thus temporarily had to comment out xcombine_vol_data* from Makefile.in and from src/auxiliaries/rules.mk because I do not know how to fix the Makefile (probably not difficult, but I am not familiar with that). Also removed unused variables introduced by the CEM modification a few days ago. Added comments for ellipticity. Removed useless white spaces using a cleaning script. (b60514e)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Nov 6 08:30:06 PST 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branches: devel,master
Link : https://github.com/geodynamics/specfem3d_globe/compare/bc58e579b3b0838a0968725a076f5904845437ca...be63f20cbb6f462104e949894dbe205d2398cd7f
>---------------------------------------------------------------
commit b60514e5506d711b3f8464a34b7b9df5539a582c
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date: Sat Sep 6 14:10:50 2014 +0200
removed src/auxiliaries/combine_vol_data_shared.f90, which contained routines duplicated from other files of the code (that is dangerous, and very uneady to maintain the same routines twice, in two different files of the code).
Thus temporarily had to comment out xcombine_vol_data* from Makefile.in and from src/auxiliaries/rules.mk because I do not know how to fix the Makefile (probably not difficult, but I am not familiar with that).
Also removed unused variables introduced by the CEM modification a few days ago.
Added comments for ellipticity.
Removed useless white spaces using a cleaning script.
>---------------------------------------------------------------
b60514e5506d711b3f8464a34b7b9df5539a582c
Makefile.in | 25 +-
src/auxiliaries/combine_vol_data.F90 | 8 +-
src/auxiliaries/combine_vol_data_shared.f90 | 569 ---------------------
src/auxiliaries/rules.mk | 17 +-
src/meshfem3D/compute_coordinates_grid.f90 | 2 +-
src/meshfem3D/create_regions_mesh.F90 | 2 +-
src/meshfem3D/{get_model.f90 => get_model.F90} | 7 +-
src/meshfem3D/meshfem3D_models.F90 | 18 +-
src/meshfem3D/model_cem.f90 | 16 +-
src/meshfem3D/model_ppm.f90 | 2 +-
src/meshfem3D/write_AVS_DX_global_chunks_data.f90 | 3 +
.../write_AVS_DX_global_chunks_data_adios.f90 | 8 +-
src/meshfem3D/write_AVS_DX_global_faces_data.f90 | 3 +
.../write_AVS_DX_global_faces_data_adios.f90 | 3 +
src/meshfem3D/write_AVS_DX_surface_data.f90 | 3 +
src/meshfem3D/write_AVS_DX_surface_data_adios.f90 | 3 +
src/shared/asdf_helpers_definitions.f90 | 2 +-
src/shared/make_ellipticity.f90 | 45 ++
src/specfem3D/comp_source_time_function.f90 | 8 +-
src/specfem3D/compute_add_sources.f90 | 2 +
src/specfem3D/get_cmt.f90 | 2 +-
src/specfem3D/locate_receivers.f90 | 3 +
src/specfem3D/locate_sources.f90 | 3 +
src/specfem3D/setup_sources_receivers.f90 | 2 +-
src/specfem3D/write_output_SAC.f90 | 2 +-
25 files changed, 134 insertions(+), 624 deletions(-)
diff --git a/Makefile.in b/Makefile.in
index a471fa4..cd60d2c 100644
--- a/Makefile.in
+++ b/Makefile.in
@@ -225,21 +225,25 @@ DEFAULT = \
xcreate_header_file \
xmeshfem3D \
xspecfem3D \
- xcombine_vol_data \
- xcombine_vol_data_vtk \
xcombine_surf_data \
xcombine_AVS_DX \
xconvolve_source_timefunction \
xcreate_movie_AVS_DX \
xcreate_movie_GMT_global \
$(EMPTY_MACRO)
-
-ifeq ($(ADIOS), yes)
-DEFAULT += \
- xcombine_vol_data_adios \
- xcombine_vol_data_vtk_adios \
- $(EMPTY_MACRO)
-endif
+## DK DK comments out all the xcombine_vol_data* for now because they depend on $O/combine_vol_data_shared.aux.o, which contained duplicated code and which I have thus removed
+## DK DK put this back in the above list of default targets if you want to compile them again,
+## DK DK when adding comments I need to move them outside of the list to avoid getting an error from "make"
+# xcombine_vol_data \
+# xcombine_vol_data_vtk \
+
+## DK DK comments out all the xcombine_vol_data* for now because they depend on $O/combine_vol_data_shared.aux.o, which contained duplicated code and which I have thus removed
+#ifeq ($(ADIOS), yes)
+#DEFAULT += \
+# xcombine_vol_data_adios \
+# xcombine_vol_data_vtk_adios \
+# $(EMPTY_MACRO)
+#endif
all: default
@@ -255,6 +259,9 @@ clean:
-rm -f $(foreach dir, $(SUBDIRS), $($(dir)_OBJECTS) $($(dir)_MODULES) $($(dir)_TARGETS)) $O/*
endif
+realclean: clean
+mrproper: clean
+
help:
@echo "usage: make [executable]"
@echo ""
diff --git a/src/auxiliaries/combine_vol_data.F90 b/src/auxiliaries/combine_vol_data.F90
index 3e610b6..6974830 100644
--- a/src/auxiliaries/combine_vol_data.F90
+++ b/src/auxiliaries/combine_vol_data.F90
@@ -77,7 +77,7 @@ program combine_vol_data
!
! puts point locations back into a perfectly spherical shape by removing the ellipticity factor;
! useful for plotting spherical cuts at certain depths
- logical,parameter:: CORRECT_ELLIPTICITY = .false.
+ logical, parameter :: CORRECT_ELLIPTICITY = .false.
integer :: nspl
double precision :: rspl(NR),espl(NR),espl2(NR)
@@ -250,7 +250,7 @@ program combine_vol_data
endif
! sets up ellipticity splines in order to remove ellipticity from point coordinates
- if (CORRECT_ELLIPTICITY ) call make_ellipticity(nspl,rspl,espl,espl2,ONE_CRUST)
+ if (CORRECT_ELLIPTICITY) call make_ellipticity(nspl,rspl,espl,espl2,ONE_CRUST)
#ifdef ADIOS_INPUT
call init_adios(value_file_name, mesh_file_name, value_handle, mesh_handle)
@@ -481,7 +481,7 @@ program combine_vol_data
z = zstore(iglob)
! remove ellipticity
- if (CORRECT_ELLIPTICITY ) call reverse_ellipticity(x,y,z,nspl,rspl,espl,espl2)
+ if (CORRECT_ELLIPTICITY) call revert_ellipticity(x,y,z,nspl,rspl,espl,espl2)
!dat = data(i,j,k,ispec)
dat = ibool_dat(iglob)
@@ -508,7 +508,7 @@ program combine_vol_data
z = zstore(iglob)
! remove ellipticity
- if (CORRECT_ELLIPTICITY ) call reverse_ellipticity(x,y,z,nspl,rspl,espl,espl2)
+ if (CORRECT_ELLIPTICITY) call revert_ellipticity(x,y,z,nspl,rspl,espl,espl2)
dat = data(i,j,k,ispec)
diff --git a/src/auxiliaries/combine_vol_data_shared.f90 b/src/auxiliaries/combine_vol_data_shared.f90
deleted file mode 100644
index 3bc648a..0000000
--- a/src/auxiliaries/combine_vol_data_shared.f90
+++ /dev/null
@@ -1,569 +0,0 @@
-!=====================================================================
-!
-! S p e c f e m 3 D G l o b e V e r s i o n 6 . 0
-! --------------------------------------------------
-!
-! Main historical authors: Dimitri Komatitsch and Jeroen Tromp
-! Princeton University, USA
-! and CNRS / University of Marseille, France
-! (there are currently many more authors!)
-! (c) Princeton University and CNRS / University of Marseille, April 2014
-!
-! 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.
-!
-!=====================================================================
-
- subroutine reverse_ellipticity(x,y,z,nspl,rspl,espl,espl2)
-
- use constants
-
- implicit none
-
- real(kind=CUSTOM_REAL) :: x,y,z
- integer nspl
- double precision rspl(NR),espl(NR),espl2(NR)
- double precision x1,y1,z1
-
- double precision ell
- double precision r,theta,phi,factor
- double precision cost,p20
-
- ! gets spherical coordinates
- x1 = x
- y1 = y
- z1 = z
- call xyz_2_rthetaphi_dble(x1,y1,z1,r,theta,phi)
-
- cost=dcos(theta)
-! this is the Legendre polynomial of degree two, P2(cos(theta)), see the discussion above eq (14.4) in Dahlen and Tromp (1998)
- p20=0.5d0*(3.0d0*cost*cost-1.0d0)
-
- ! get ellipticity using spline evaluation
- call spline_evaluation(rspl,espl,espl2,nspl,r,ell)
-
-! this is eq (14.4) in Dahlen and Tromp (1998)
- factor=ONE-(TWO/3.0d0)*ell*p20
-
- ! removes ellipticity factor
- x = x / factor
- y = y / factor
- z = z / factor
-
- end subroutine reverse_ellipticity
-
-!
-! ------------------------------------------------------------------------------------------------
-!
-
-! copy from make_ellipticity.f90 to avoid compiling issues
-
- subroutine make_ellipticity(nspl,rspl,espl,espl2,ONE_CRUST)
-
-! creates a spline for the ellipticity profile in PREM
-! radius and density are non-dimensional
-
- use constants
-
- implicit none
-
- integer nspl
-
- logical ONE_CRUST
-
-! radius of the Earth for gravity calculation
- double precision, parameter :: R_EARTH_ELLIPTICITY = 6371000.d0
-! radius of the ocean floor for gravity calculation
- double precision, parameter :: ROCEAN_ELLIPTICITY = 6368000.d0
-
- double precision rspl(NR),espl(NR),espl2(NR)
-
- integer i
- double precision ROCEAN,RMIDDLE_CRUST,RMOHO,R80,R220,R400,R600,R670, &
- R771,RTOPDDOUBLEPRIME,RCMB,RICB
- double precision r_icb,r_cmb,r_topddoubleprime,r_771,r_670,r_600
- double precision r_400,r_220,r_80,r_moho,r_middle_crust,r_ocean,r_0
- double precision r(NR),rho(NR),epsilonval(NR),eta(NR)
- double precision radau(NR),z,k(NR),g_a,bom,exponentval,i_rho,i_radau
- double precision s1(NR),s2(NR),s3(NR)
- double precision yp1,ypn
-
-! PREM
- ROCEAN = 6368000.d0
- RMIDDLE_CRUST = 6356000.d0
- RMOHO = 6346600.d0
- R80 = 6291000.d0
- R220 = 6151000.d0
- R400 = 5971000.d0
- R600 = 5771000.d0
- R670 = 5701000.d0
- R771 = 5600000.d0
- RTOPDDOUBLEPRIME = 3630000.d0
- RCMB = 3480000.d0
- RICB = 1221000.d0
-
-! non-dimensionalize
- r_icb = RICB/R_EARTH_ELLIPTICITY
- r_cmb = RCMB/R_EARTH_ELLIPTICITY
- r_topddoubleprime = RTOPDDOUBLEPRIME/R_EARTH_ELLIPTICITY
- r_771 = R771/R_EARTH_ELLIPTICITY
- r_670 = R670/R_EARTH_ELLIPTICITY
- r_600 = R600/R_EARTH_ELLIPTICITY
- r_400 = R400/R_EARTH_ELLIPTICITY
- r_220 = R220/R_EARTH_ELLIPTICITY
- r_80 = R80/R_EARTH_ELLIPTICITY
- r_moho = RMOHO/R_EARTH_ELLIPTICITY
- r_middle_crust = RMIDDLE_CRUST/R_EARTH_ELLIPTICITY
- r_ocean = ROCEAN_ELLIPTICITY/R_EARTH_ELLIPTICITY
- r_0 = 1.d0
-
- do i=1,163
- r(i) = r_icb*dble(i-1)/dble(162)
- enddo
- do i=164,323
- r(i) = r_icb+(r_cmb-r_icb)*dble(i-164)/dble(159)
- enddo
- do i=324,336
- r(i) = r_cmb+(r_topddoubleprime-r_cmb)*dble(i-324)/dble(12)
- enddo
- do i=337,517
- r(i) = r_topddoubleprime+(r_771-r_topddoubleprime)*dble(i-337)/dble(180)
- enddo
- do i=518,530
- r(i) = r_771+(r_670-r_771)*dble(i-518)/dble(12)
- enddo
- do i=531,540
- r(i) = r_670+(r_600-r_670)*dble(i-531)/dble(9)
- enddo
- do i=541,565
- r(i) = r_600+(r_400-r_600)*dble(i-541)/dble(24)
- enddo
- do i=566,590
- r(i) = r_400+(r_220-r_400)*dble(i-566)/dble(24)
- enddo
- do i=591,609
- r(i) = r_220+(r_80-r_220)*dble(i-591)/dble(18)
- enddo
- do i=610,619
- r(i) = r_80+(r_moho-r_80)*dble(i-610)/dble(9)
- enddo
- do i=620,626
- r(i) = r_moho+(r_middle_crust-r_moho)*dble(i-620)/dble(6)
- enddo
- do i=627,633
- r(i) = r_middle_crust+(r_ocean-r_middle_crust)*dble(i-627)/dble(6)
- enddo
- do i=634,NR
- r(i) = r_ocean+(r_0-r_ocean)*dble(i-634)/dble(6)
- enddo
-
-! use PREM to get the density profile for ellipticity (fine for other 1D reference models)
- do i = 1,NR
- call prem_density(r(i),rho(i),ONE_CRUST,RICB,RCMB,RTOPDDOUBLEPRIME, &
- R600,R670,R220,R771,R400,R80,RMOHO,RMIDDLE_CRUST,ROCEAN)
- radau(i)=rho(i)*r(i)*r(i)
- enddo
-
- eta(1)=0.0d0
-
- k(1)=0.0d0
-
- do i=2,NR
- call intgrl(i_rho,r,1,i,rho,s1,s2,s3)
-! Radau approximation of Clairaut's equation for first-order terms of ellipticity, see e.g. Jeffreys H.,
-! The figures of rotating planets, Mon. Not. R. astr. Soc., vol. 113, p. 97-105 (1953).
-! The Radau approximation is mentioned on page 97.
-! For more details see Section 14.1.2 in Dahlen and Tromp (1998)
-! (see also in file ellipticity_equations_from_Dahlen_Tromp_1998.pdf in the "doc" directory of the code).
- call intgrl(i_radau,r,1,i,radau,s1,s2,s3)
- z=(2.0d0/3.0d0)*i_radau/(i_rho*r(i)*r(i))
-! this comes from equation (14.19) in Dahlen and Tromp (1998)
- eta(i)=(25.0d0/4.0d0)*((1.0d0-(3.0d0/2.0d0)*z)**2.0d0)-1.0d0
- k(i)=eta(i)/(r(i)**3.0d0)
- enddo
-
- bom=TWO_PI/(HOURS_PER_DAY*SECONDS_PER_HOUR)
-
-! non-dimensionalized version
- bom=bom/sqrt(PI*GRAV*RHOAV)
-
- g_a = 4.0d0*i_rho
-! this is the equation right above (14.21) in Dahlen and Tromp (1998)
- epsilonval(NR) = (5.0d0/2.d0)*(bom**2.0d0)*R_UNIT_SPHERE / (g_a * (eta(NR)+2.0d0))
-
- do i = 1,NR-1
- call intgrl(exponentval,r,i,NR,k,s1,s2,s3)
- epsilonval(i)=epsilonval(NR)*exp(-exponentval)
- enddo
-
-! get ready to spline epsilonval
- nspl = 1
- rspl(1)=r(1)
- espl(1)=epsilonval(1)
- do i=2,NR
- if (r(i) /= r(i-1)) then
- nspl=nspl+1
- rspl(nspl)=r(i)
- espl(nspl)=epsilonval(i)
- endif
- enddo
-
-! spline epsilonval
- yp1=0.0d0
- ypn=(5.0d0/2.0d0)*(bom**2)/g_a-2.0d0*epsilonval(NR)
- call spline_construction(rspl,espl,nspl,yp1,ypn,espl2)
-
- end subroutine make_ellipticity
-
-!
-! ------------------------------------------------------------------------------------------------
-!
-
-! copy from model_prem.f90 to avoid compiling issues
-
- subroutine prem_density(x,rho,ONE_CRUST,RICB,RCMB,RTOPDDOUBLEPRIME, &
- R600,R670,R220,R771,R400,R80,RMOHO,RMIDDLE_CRUST,ROCEAN)
-
- use constants
-
- implicit none
-
- double precision x,rho,RICB,RCMB,RTOPDDOUBLEPRIME, &
- R600,R670,R220,R771,R400,R80,RMOHO,RMIDDLE_CRUST,ROCEAN
-
- logical ONE_CRUST
-
- double precision r
-
- ! compute real physical radius in meters
- r = x * R_EARTH
-
- ! calculates density according to radius
- if (r <= RICB) then
- rho=13.0885d0-8.8381d0*x*x
- else if (r > RICB .and. r <= RCMB) then
- rho=12.5815d0-1.2638d0*x-3.6426d0*x*x-5.5281d0*x*x*x
- else if (r > RCMB .and. r <= RTOPDDOUBLEPRIME) then
- rho=7.9565d0-6.4761d0*x+5.5283d0*x*x-3.0807d0*x*x*x
- else if (r > RTOPDDOUBLEPRIME .and. r <= R771) then
- rho=7.9565d0-6.4761d0*x+5.5283d0*x*x-3.0807d0*x*x*x
- else if (r > R771 .and. r <= R670) then
- rho=7.9565d0-6.4761d0*x+5.5283d0*x*x-3.0807d0*x*x*x
- else if (r > R670 .and. r <= R600) then
- rho=5.3197d0-1.4836d0*x
- else if (r > R600 .and. r <= R400) then
- rho=11.2494d0-8.0298d0*x
- else if (r > R400 .and. r <= R220) then
- rho=7.1089d0-3.8045d0*x
- else if (r > R220 .and. r <= R80) then
- rho=2.6910d0+0.6924d0*x
- else
- if (r > R80 .and. r <= RMOHO) then
- rho=2.6910d0+0.6924d0*x
- else if (r > RMOHO .and. r <= RMIDDLE_CRUST) then
- if (ONE_CRUST) then
- rho=2.6d0
- else
- rho=2.9d0
- endif
- else if (r > RMIDDLE_CRUST .and. r <= ROCEAN) then
- rho=2.6d0
- else if (r > ROCEAN) then
- rho=2.6d0
- endif
- endif
-
- rho=rho*1000.0d0/RHOAV
-
- end subroutine prem_density
-
-!
-! ------------------------------------------------------------------------------------------------
-!
-
-! copy from intgrl.f90 to avoid compiling issues
-
-
- subroutine intgrl(sumval,r,nir,ner,f,s1,s2,s3)
-
-! Computes the integral of f[i]*r[i]*r[i] from i=nir to i=ner for
-! radii values as in model PREM_an640
-
- implicit none
-
-! Argument variables
- integer ner,nir
- double precision f(640),r(640),s1(640),s2(640)
- double precision s3(640),sumval
-
-! Local variables
- double precision, parameter :: third = 1.0d0/3.0d0
- double precision, parameter :: fifth = 1.0d0/5.0d0
- double precision, parameter :: sixth = 1.0d0/6.0d0
-
- double precision rji,yprime(640)
- double precision s1l,s2l,s3l
-
- integer i,j,n,kdis(28)
- integer ndis,nir1
-
- data kdis/163,323,336,517,530,540,565,590,609,619,626,633,16*0/
-
- ndis = 12
- n = 640
-
- call deriv(f,yprime,n,r,ndis,kdis,s1,s2,s3)
- nir1 = nir + 1
- sumval = 0.0d0
- do i=nir1,ner
- j = i-1
- rji = r(i) - r(j)
- s1l = s1(j)
- s2l = s2(j)
- s3l = s3(j)
- sumval = sumval + r(j)*r(j)*rji*(f(j) &
- + rji*(0.5d0*s1l + rji*(third*s2l + rji*0.25d0*s3l))) &
- + 2.0d0*r(j)*rji*rji*(0.5d0*f(j) + rji*(third*s1l + rji*(0.25d0*s2l + rji*fifth*s3l))) &
- + rji*rji*rji*(third*f(j) + rji*(0.25d0*s1l + rji*(fifth*s2l + rji*sixth*s3l)))
- enddo
-
- end subroutine intgrl
-
-! -------------------------------
-
- subroutine deriv(y,yprime,n,r,ndis,kdis,s1,s2,s3)
-
- implicit none
-
-! Argument variables
- integer kdis(28),n,ndis
- double precision r(n),s1(n),s2(n),s3(n)
- double precision y(n),yprime(n)
-
-! Local variables
- integer i,j,j1,j2
- integer k,nd,ndp
- double precision a0,b0,b1
- double precision f(3,1000),h,h2,h2a
- double precision h2b,h3a,ha,s13
- double precision s21,s32,yy(3)
-
- yy(1) = 0.d0
- yy(2) = 0.d0
- yy(3) = 0.d0
-
- ndp=ndis+1
- do 3 nd = 1,ndp
- if (nd == 1) goto 4
- if (nd == ndp) goto 5
- j1=kdis(nd-1)+1
- j2=kdis(nd)-2
- goto 6
- 4 j1=1
- j2=kdis(1)-2
- goto 6
- 5 j1=kdis(ndis)+1
- j2=n-2
- 6 if ((j2+1-j1)>0) goto 11
- j2=j2+2
- yy(1)=(y(j2)-y(j1))/(r(j2)-r(j1))
- s1(j1)=yy(1)
- s1(j2)=yy(1)
- s2(j1)=yy(2)
- s2(j2)=yy(2)
- s3(j1)=yy(3)
- s3(j2)=yy(3)
- goto 3
- 11 a0=0.0d0
- if (j1 == 1) goto 7
- h=r(j1+1)-r(j1)
- h2=r(j1+2)-r(j1)
- yy(1)=h*h2*(h2-h)
- h=h*h
- h2=h2*h2
- b0=(y(j1)*(h-h2)+y(j1+1)*h2-y(j1+2)*h)/yy(1)
- goto 8
- 7 b0=0.0d0
- 8 b1=b0
-
- if (j2 > 1000) stop 'Error in subroutine deriv for j2'
-
- do i=j1,j2
- h=r(i+1)-r(i)
- yy(1)=y(i+1)-y(i)
- h2=h*h
- ha=h-a0
- h2a=h-2.0d0*a0
- h3a=2.0d0*h-3.0d0*a0
- h2b=h2*b0
- s1(i)=h2/ha
- s2(i)=-ha/(h2a*h2)
- s3(i)=-h*h2a/h3a
- f(1,i)=(yy(1)-h*b0)/(h*ha)
- f(2,i)=(h2b-yy(1)*(2.0d0*h-a0))/(h*h2*h2a)
- f(3,i)=-(h2b-3.0d0*yy(1)*ha)/(h*h3a)
- a0=s3(i)
- b0=f(3,i)
- enddo
-
- i=j2+1
- h=r(i+1)-r(i)
- yy(1)=y(i+1)-y(i)
- h2=h*h
- ha=h-a0
- h2a=h*ha
- h2b=h2*b0-yy(1)*(2.d0*h-a0)
- s1(i)=h2/ha
- f(1,i)=(yy(1)-h*b0)/h2a
- ha=r(j2)-r(i+1)
- yy(1)=-h*ha*(ha+h)
- ha=ha*ha
- yy(1)=(y(i+1)*(h2-ha)+y(i)*ha-y(j2)*h2)/yy(1)
- s3(i)=(yy(1)*h2a+h2b)/(h*h2*(h-2.0d0*a0))
- s13=s1(i)*s3(i)
- s2(i)=f(1,i)-s13
-
- do j=j1,j2
- k=i-1
- s32=s3(k)*s2(i)
- s1(i)=f(3,k)-s32
- s21=s2(k)*s1(i)
- s3(k)=f(2,k)-s21
- s13=s1(k)*s3(k)
- s2(k)=f(1,k)-s13
- i=k
- enddo
-
- s1(i)=b1
- j2=j2+2
- s1(j2)=yy(1)
- s2(j2)=yy(2)
- s3(j2)=yy(3)
- 3 continue
-
- do i = 1,n
- yprime(i)=s1(i)
- enddo
-
- end subroutine deriv
-
-!
-! ------------------------------------------------------------------------------------------------
-!
-
-! copy from spline_routines.f90 to avoid compiling issues
-
-! compute spline coefficients
-
- subroutine spline_construction(xpoint,ypoint,npoint,tangent_first_point,tangent_last_point,spline_coefficients)
-
- implicit none
-
-! tangent to the spline imposed at the first and last points
- double precision, intent(in) :: tangent_first_point,tangent_last_point
-
-! number of input points and coordinates of the input points
- integer, intent(in) :: npoint
- double precision, dimension(npoint), intent(in) :: xpoint,ypoint
-
-! spline coefficients output by the routine
- double precision, dimension(npoint), intent(out) :: spline_coefficients
-
- integer :: i
-
- double precision, dimension(:), allocatable :: temporary_array
-
- allocate(temporary_array(npoint))
-
- spline_coefficients(1) = - 1.d0 / 2.d0
-
- temporary_array(1) = (3.d0/(xpoint(2)-xpoint(1)))*((ypoint(2)-ypoint(1))/(xpoint(2)-xpoint(1))-tangent_first_point)
-
- do i = 2,npoint-1
-
- spline_coefficients(i) = ((xpoint(i)-xpoint(i-1))/(xpoint(i+1)-xpoint(i-1))-1.d0) &
- / ((xpoint(i)-xpoint(i-1))/(xpoint(i+1)-xpoint(i-1))*spline_coefficients(i-1)+2.d0)
-
- temporary_array(i) = (6.d0*((ypoint(i+1)-ypoint(i))/(xpoint(i+1)-xpoint(i)) &
- - (ypoint(i)-ypoint(i-1))/(xpoint(i)-xpoint(i-1)))/(xpoint(i+1)-xpoint(i-1)) &
- - (xpoint(i)-xpoint(i-1))/(xpoint(i+1)-xpoint(i-1))*temporary_array(i-1)) &
- / ((xpoint(i)-xpoint(i-1))/(xpoint(i+1)-xpoint(i-1))*spline_coefficients(i-1)+2.d0)
-
- enddo
-
- spline_coefficients(npoint) = ((3.d0/(xpoint(npoint)-xpoint(npoint-1))) &
- * (tangent_last_point-(ypoint(npoint)-ypoint(npoint-1))/(xpoint(npoint)-xpoint(npoint-1))) &
- - 1.d0/2.d0*temporary_array(npoint-1))/(1.d0/2.d0*spline_coefficients(npoint-1)+1.d0)
-
- do i = npoint-1,1,-1
- spline_coefficients(i) = spline_coefficients(i)*spline_coefficients(i+1) + temporary_array(i)
- enddo
-
- deallocate(temporary_array)
-
- end subroutine spline_construction
-
-! --------------
-
-! evaluate a spline
-
- subroutine spline_evaluation(xpoint,ypoint,spline_coefficients,npoint,x_evaluate_spline,y_spline_obtained)
-
- implicit none
-
-! number of input points and coordinates of the input points
- integer, intent(in) :: npoint
- double precision, dimension(npoint), intent(in) :: xpoint,ypoint
-
-! spline coefficients to use
- double precision, dimension(npoint), intent(in) :: spline_coefficients
-
-! abscissa at which we need to evaluate the value of the spline
- double precision, intent(in):: x_evaluate_spline
-
-! ordinate evaluated by the routine for the spline at this abscissa
- double precision, intent(out):: y_spline_obtained
-
- integer :: index_loop,index_lower,index_higher
-
- double precision :: coef1,coef2
-
-! initialize to the whole interval
- index_lower = 1
- index_higher = npoint
-
-! determine the right interval to use, by dichotomy
- do while (index_higher - index_lower > 1)
-! compute the middle of the interval
- index_loop = (index_higher + index_lower) / 2
- if (xpoint(index_loop) > x_evaluate_spline) then
- index_higher = index_loop
- else
- index_lower = index_loop
- endif
- enddo
-
-! test that the interval obtained does not have a size of zero
-! (this could happen for instance in the case of duplicates in the input list of points)
- if (xpoint(index_higher) == xpoint(index_lower)) stop 'incorrect interval found in spline evaluation'
-
- coef1 = (xpoint(index_higher) - x_evaluate_spline) / (xpoint(index_higher) - xpoint(index_lower))
- coef2 = (x_evaluate_spline - xpoint(index_lower)) / (xpoint(index_higher) - xpoint(index_lower))
-
- y_spline_obtained = coef1*ypoint(index_lower) + coef2*ypoint(index_higher) + &
- ((coef1**3 - coef1)*spline_coefficients(index_lower) + &
- (coef2**3 - coef2)*spline_coefficients(index_higher))*((xpoint(index_higher) - xpoint(index_lower))**2)/6.d0
-
- end subroutine spline_evaluation
-
diff --git a/src/auxiliaries/rules.mk b/src/auxiliaries/rules.mk
index cbdde73..198be60 100644
--- a/src/auxiliaries/rules.mk
+++ b/src/auxiliaries/rules.mk
@@ -101,17 +101,18 @@ ${E}/xcombine_AVS_DX: $(auxiliaries_SHARED_OBJECTS) $O/get_cmt.solver.o $O/combi
${E}/xcombine_paraview_strain_data: $(auxiliaries_SHARED_OBJECTS) $O/combine_paraview_strain_data.auxsolver.o $O/binary_c_io.cc.o
${FCCOMPILE_CHECK} -o $@ $+
-${E}/xcombine_vol_data: $(auxiliaries_SHARED_OBJECTS) $O/combine_vol_data.auxsolver.o $O/binary_c_io.cc.o $O/combine_vol_data_shared.aux.o
- ${FCCOMPILE_CHECK} -o $@ $+
+## DK DK comments out all the ${E}/xcombine_vol_data* for now because they depend on $O/combine_vol_data_shared.aux.o, which contained duplicated code and which I have thus removed
+#${E}/xcombine_vol_data: $(auxiliaries_SHARED_OBJECTS) $O/combine_vol_data.auxsolver.o $O/binary_c_io.cc.o $O/combine_vol_data_shared.aux.o
+# ${FCCOMPILE_CHECK} -o $@ $+
-${E}/xcombine_vol_data_adios: $(auxiliaries_SHARED_OBJECTS) $O/combine_vol_data_adios_impl.auxmpi.o $O/combine_vol_data.auxadios.o $O/binary_c_io.cc.o $O/combine_vol_data_shared.aux.o $O/parallel.sharedmpi.o
- ${MPIFCCOMPILE_CHECK} -o $@ $+ $(MPILIBS)
+#${E}/xcombine_vol_data_adios: $(auxiliaries_SHARED_OBJECTS) $O/combine_vol_data_adios_impl.auxmpi.o $O/combine_vol_data.auxadios.o $O/binary_c_io.cc.o $O/combine_vol_data_shared.aux.o $O/parallel.sharedmpi.o
+# ${MPIFCCOMPILE_CHECK} -o $@ $+ $(MPILIBS)
-${E}/xcombine_vol_data_vtk_adios: $(auxiliaries_SHARED_OBJECTS) $O/combine_vol_data_adios_impl.auxmpi.o $O/combine_vol_data.auxadios_vtk.o $O/binary_c_io.cc.o $O/combine_vol_data_shared.aux.o $O/parallel.sharedmpi.o
- ${MPIFCCOMPILE_CHECK} -o $@ $+ $(MPILIBS)
+#${E}/xcombine_vol_data_vtk_adios: $(auxiliaries_SHARED_OBJECTS) $O/combine_vol_data_adios_impl.auxmpi.o $O/combine_vol_data.auxadios_vtk.o $O/binary_c_io.cc.o $O/combine_vol_data_shared.aux.o $O/parallel.sharedmpi.o
+# ${MPIFCCOMPILE_CHECK} -o $@ $+ $(MPILIBS)
-${E}/xcombine_vol_data_vtk: $(auxiliaries_SHARED_OBJECTS) $O/combine_vol_data.auxsolver_vtk.o $O/binary_c_io.cc.o $O/combine_vol_data_shared.aux.o
- ${FCCOMPILE_CHECK} -o $@ $+
+#${E}/xcombine_vol_data_vtk: $(auxiliaries_SHARED_OBJECTS) $O/combine_vol_data.auxsolver_vtk.o $O/binary_c_io.cc.o $O/combine_vol_data_shared.aux.o
+# ${FCCOMPILE_CHECK} -o $@ $+
${E}/xcombine_surf_data: $(auxiliaries_SHARED_OBJECTS) $O/combine_surf_data.auxsolver.o $O/binary_c_io.cc.o
${FCCOMPILE_CHECK} -o $@ $+
diff --git a/src/meshfem3D/compute_coordinates_grid.f90 b/src/meshfem3D/compute_coordinates_grid.f90
index fed71d4..36388de 100644
--- a/src/meshfem3D/compute_coordinates_grid.f90
+++ b/src/meshfem3D/compute_coordinates_grid.f90
@@ -59,7 +59,7 @@
double precision :: ratio_xi, ratio_eta, fact_xi, fact_eta, fact_xi_,fact_eta_
-! this to avoid compilation warnings
+! to avoid compiler warnings
x_ = 0
y_ = 0
diff --git a/src/meshfem3D/create_regions_mesh.F90 b/src/meshfem3D/create_regions_mesh.F90
index d6f7f26..31e9a06 100644
--- a/src/meshfem3D/create_regions_mesh.F90
+++ b/src/meshfem3D/create_regions_mesh.F90
@@ -200,7 +200,7 @@
iboolright_xi,iboolleft_eta,iboolright_eta,nimin,nimax, &
njmin, njmax,nkmin_xi,nkmin_eta,iboolfaces,iboolcorner)
- end if
+ endif
#endif
diff --git a/src/meshfem3D/get_model.f90 b/src/meshfem3D/get_model.F90
similarity index 99%
rename from src/meshfem3D/get_model.f90
rename to src/meshfem3D/get_model.F90
index c98465f..583769a 100644
--- a/src/meshfem3D/get_model.f90
+++ b/src/meshfem3D/get_model.F90
@@ -155,8 +155,11 @@
RCMB,R670,RMOHO, &
xmesh,ymesh,zmesh,r, &
c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,&
- c33,c34,c35,c36,c44,c45,c46,c55,c56,c66, &
- ispec,i,j,k)
+ c33,c34,c35,c36,c44,c45,c46,c55,c56,c66 &
+#if defined (CEM)
+ ,ispec,i,j,k
+#endif
+ )
! gets the 3-D crustal model
! M.A. don't overwrite crust if using CEM.
diff --git a/src/meshfem3D/meshfem3D_models.F90 b/src/meshfem3D/meshfem3D_models.F90
index 14fa012..d4d28d5 100644
--- a/src/meshfem3D/meshfem3D_models.F90
+++ b/src/meshfem3D/meshfem3D_models.F90
@@ -344,21 +344,26 @@
!-------------------------------------------------------------------------------------------------
!
-
subroutine meshfem3D_models_get3Dmntl_val(iregion_code,r_prem,rho,dvp,&
vpv,vph,vsv,vsh,eta_aniso, &
RCMB,R670,RMOHO, &
xmesh,ymesh,zmesh,r, &
c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,&
- c33,c34,c35,c36,c44,c45,c46,c55,c56,c66, &
- ispec,i,j,k)
+ c33,c34,c35,c36,c44,c45,c46,c55,c56,c66 &
+#if defined (CEM)
+ ,ispec,i,j,k
+#endif
+ )
use meshfem3D_models_par
implicit none
- intent (in) :: ispec, i, j, k
+#if defined (CEM)
+ ! CEM needs these to determine iglob
+ integer, intent (in) :: ispec, i, j, k
+#endif
integer iregion_code
double precision r_prem
double precision rho,dvp
@@ -377,9 +382,6 @@
real(kind=4) :: xcolat,xlon,xrad,dvpv,dvph,dvsv,dvsh
logical :: found_crust,suppress_mantle_extension
- ! CEM needs these to determine iglob
- integer :: ispec, i, j, k
-
! initializes perturbation values
dvs = ZERO
dvp = ZERO
@@ -611,7 +613,7 @@
#if defined (CEM)
if (CEM_ACCEPT) then
call request_cem (vsh, vsv, vph, vpv, rho, iregion_code, ispec, i, j, k)
- end if
+ endif
#endif
!> Hejun
diff --git a/src/meshfem3D/model_cem.f90 b/src/meshfem3D/model_cem.f90
index 7d6c363..7724e51 100644
--- a/src/meshfem3D/model_cem.f90
+++ b/src/meshfem3D/model_cem.f90
@@ -67,7 +67,7 @@ subroutine model_cem_broadcast ( myrank )
call synchronize_all ()
- end if
+ endif
end subroutine model_cem_broadcast
@@ -112,7 +112,7 @@ subroutine request_cem ( vsh, vsv, vph, vpv, rho, iregion_code, ispec, i, j, k )
vsv = reg3Bc%vsv(iglob) * 1000.0d0 / (R_EARTH * scaleval)
vsh = reg3Bc%vsh(iglob) * 1000.0d0 / (R_EARTH * scaleval)
- end if
+ endif
vph = vpv
@@ -316,15 +316,15 @@ subroutine build_global_coordinates ( nspec, nglob, iregion_code )
region = 2
else if ( rad >= R_020_KM ) then
region = 1
- end if
+ endif
iglob = ibool(i,j,k,ispec)
regCode(iglob) = region
- end do
- end do
- end do
- end do
+ enddo
+ enddo
+ enddo
+ enddo
else if (iregion_code == 2) then
@@ -334,6 +334,6 @@ subroutine build_global_coordinates ( nspec, nglob, iregion_code )
regCode(:) = 9
- end if
+ endif
end subroutine build_global_coordinates
diff --git a/src/meshfem3D/model_ppm.f90 b/src/meshfem3D/model_ppm.f90
index d843d64..cfb6935 100644
--- a/src/meshfem3D/model_ppm.f90
+++ b/src/meshfem3D/model_ppm.f90
@@ -1080,7 +1080,7 @@
! >>>>>
! uniform sigma
- ! just to avoid compiler warning
+ ! to avoid compiler warnings
ii = ispec2
!exp_val(:,:,:) = exp( -((xx(:,:,:,ispec2)-x0)**2+(yy(:,:,:,ispec2)-y0)**2 &
! +(zz(:,:,:,ispec2)-z0)**2 )/(2*sigma2) )*factor(:,:,:)
diff --git a/src/meshfem3D/write_AVS_DX_global_chunks_data.f90 b/src/meshfem3D/write_AVS_DX_global_chunks_data.f90
index 649c78d..267c48d 100644
--- a/src/meshfem3D/write_AVS_DX_global_chunks_data.f90
+++ b/src/meshfem3D/write_AVS_DX_global_chunks_data.f90
@@ -553,8 +553,11 @@
if (ELLIPTICITY) then
call xyz_2_rthetaphi_dble(x,y,z,r,theta,phi_dummy)
cost=dcos(theta)
+! this is the Legendre polynomial of degree two, P2(cos(theta)), see the discussion above eq (14.4) in Dahlen and Tromp (1998)
p20=0.5d0*(3.0d0*cost*cost-1.0d0)
+! get ellipticity using spline evaluation
call spline_evaluation(rspl,espl,espl2,nspl,r,ell)
+! this is eq (14.4) in Dahlen and Tromp (1998)
factor=ONE-(TWO/3.0d0)*ell*p20
r=r/factor
endif
diff --git a/src/meshfem3D/write_AVS_DX_global_chunks_data_adios.f90 b/src/meshfem3D/write_AVS_DX_global_chunks_data_adios.f90
index f34be69..c8cc0f6 100644
--- a/src/meshfem3D/write_AVS_DX_global_chunks_data_adios.f90
+++ b/src/meshfem3D/write_AVS_DX_global_chunks_data_adios.f90
@@ -80,11 +80,6 @@ contains
integer, dimension(8) :: iglobval
integer npoin,nspecface !,ispecface,numpoin
- !real(kind=CUSTOM_REAL) vmin,vmax
- !double precision r,rho,vp,vs,Qkappa,Qmu
- !double precision vpv,vph,vsv,vsh !,eta_aniso
- !double precision x,y,z,theta,phi_dummy,p20 !,ell,factor,cost
- !real(kind=CUSTOM_REAL) dvp,dvs
integer :: ierr
! Dummy arrays for type inference inside adios helpers
@@ -875,8 +870,11 @@ contains
if (ELLIPTICITY) then
call xyz_2_rthetaphi_dble(x,y,z,r,theta,phi_dummy)
cost=dcos(theta)
+! this is the Legendre polynomial of degree two, P2(cos(theta)), see the discussion above eq (14.4) in Dahlen and Tromp (1998)
p20=0.5d0*(3.0d0*cost*cost-1.0d0)
+! get ellipticity using spline evaluation
call spline_evaluation(rspl,espl,espl2,nspl,r,ell)
+! this is eq (14.4) in Dahlen and Tromp (1998)
factor=ONE-(TWO/3.0d0)*ell*p20
r=r/factor
endif
diff --git a/src/meshfem3D/write_AVS_DX_global_faces_data.f90 b/src/meshfem3D/write_AVS_DX_global_faces_data.f90
index 1c4f31a..b8f2c43 100644
--- a/src/meshfem3D/write_AVS_DX_global_faces_data.f90
+++ b/src/meshfem3D/write_AVS_DX_global_faces_data.f90
@@ -355,8 +355,11 @@
if (ELLIPTICITY) then
call xyz_2_rthetaphi_dble(x,y,z,r,theta,phi_dummy)
cost=dcos(theta)
+! this is the Legendre polynomial of degree two, P2(cos(theta)), see the discussion above eq (14.4) in Dahlen and Tromp (1998)
p20=0.5d0*(3.0d0*cost*cost-1.0d0)
+! get ellipticity using spline evaluation
call spline_evaluation(rspl,espl,espl2,nspl,r,ell)
+! this is eq (14.4) in Dahlen and Tromp (1998)
factor=ONE-(TWO/3.0d0)*ell*p20
r=r/factor
endif
diff --git a/src/meshfem3D/write_AVS_DX_global_faces_data_adios.f90 b/src/meshfem3D/write_AVS_DX_global_faces_data_adios.f90
index 5665128..eea3c1b 100644
--- a/src/meshfem3D/write_AVS_DX_global_faces_data_adios.f90
+++ b/src/meshfem3D/write_AVS_DX_global_faces_data_adios.f90
@@ -535,8 +535,11 @@ subroutine prepare_AVS_DX_global_faces_data_adios(myrank, nspec, &
if (ELLIPTICITY) then
call xyz_2_rthetaphi_dble(x,y,z,r,theta,phi_dummy)
cost=dcos(theta)
+! this is the Legendre polynomial of degree two, P2(cos(theta)), see the discussion above eq (14.4) in Dahlen and Tromp (1998)
p20=0.5d0*(3.0d0*cost*cost-1.0d0)
+! get ellipticity using spline evaluation
call spline_evaluation(rspl,espl,espl2,nspl,r,ell)
+! this is eq (14.4) in Dahlen and Tromp (1998)
factor=ONE-(TWO/3.0d0)*ell*p20
r=r/factor
endif
diff --git a/src/meshfem3D/write_AVS_DX_surface_data.f90 b/src/meshfem3D/write_AVS_DX_surface_data.f90
index 180f47b..3485879 100644
--- a/src/meshfem3D/write_AVS_DX_surface_data.f90
+++ b/src/meshfem3D/write_AVS_DX_surface_data.f90
@@ -219,8 +219,11 @@
if (ELLIPTICITY) then
call xyz_2_rthetaphi_dble(x,y,z,r,theta,phi_dummy)
cost=dcos(theta)
+! this is the Legendre polynomial of degree two, P2(cos(theta)), see the discussion above eq (14.4) in Dahlen and Tromp (1998)
p20=0.5d0*(3.0d0*cost*cost-1.0d0)
+! get ellipticity using spline evaluation
call spline_evaluation(rspl,espl,espl2,nspl,r,ell)
+! this is eq (14.4) in Dahlen and Tromp (1998)
factor=ONE-(TWO/3.0d0)*ell*p20
r=r/factor
endif
diff --git a/src/meshfem3D/write_AVS_DX_surface_data_adios.f90 b/src/meshfem3D/write_AVS_DX_surface_data_adios.f90
index 73515ee..2bd2ede 100644
--- a/src/meshfem3D/write_AVS_DX_surface_data_adios.f90
+++ b/src/meshfem3D/write_AVS_DX_surface_data_adios.f90
@@ -346,8 +346,11 @@ end subroutine define_AVS_DX_surfaces_data_adios
if (ELLIPTICITY) then
call xyz_2_rthetaphi_dble(x,y,z,r,theta,phi_dummy)
cost=dcos(theta)
+! this is the Legendre polynomial of degree two, P2(cos(theta)), see the discussion above eq (14.4) in Dahlen and Tromp (1998)
p20=0.5d0*(3.0d0*cost*cost-1.0d0)
+! get ellipticity using spline evaluation
call spline_evaluation(rspl,espl,espl2,nspl,r,ell)
+! this is eq (14.4) in Dahlen and Tromp (1998)
factor=ONE-(TWO/3.0d0)*ell*p20
r=r/factor
endif
diff --git a/src/shared/asdf_helpers_definitions.f90 b/src/shared/asdf_helpers_definitions.f90
index 45c510e..f9116d1 100644
--- a/src/shared/asdf_helpers_definitions.f90
+++ b/src/shared/asdf_helpers_definitions.f90
@@ -552,7 +552,7 @@
!
! call define_adios_global_1d_double_generic(adios_group, group_size_inc,full_name, local_dim)
!
-! ! to avoid compilation warnings
+! ! to avoid compiler warnings
! idummy = size(var)
!
!end subroutine define_adios_global_1d_double_1d
diff --git a/src/shared/make_ellipticity.f90 b/src/shared/make_ellipticity.f90
index 1422c82..980ae37 100644
--- a/src/shared/make_ellipticity.f90
+++ b/src/shared/make_ellipticity.f90
@@ -188,3 +188,48 @@
end subroutine make_ellipticity
+!
+!-----------------------------------------------------------------
+!
+
+ subroutine revert_ellipticity(x,y,z,nspl,rspl,espl,espl2)
+
+! this routine to revert ellipticity and go back to a spherical Earth
+! is currently used by src/auxiliaries/combine_vol_data.F90 only
+
+ use constants
+
+ implicit none
+
+ real(kind=CUSTOM_REAL) :: x,y,z
+ integer nspl
+ double precision rspl(NR),espl(NR),espl2(NR)
+ double precision x1,y1,z1
+
+ double precision ell
+ double precision r,theta,phi,factor
+ double precision cost,p20
+
+ ! gets spherical coordinates
+ x1 = x
+ y1 = y
+ z1 = z
+ call xyz_2_rthetaphi_dble(x1,y1,z1,r,theta,phi)
+
+ cost=dcos(theta)
+! this is the Legendre polynomial of degree two, P2(cos(theta)), see the discussion above eq (14.4) in Dahlen and Tromp (1998)
+ p20=0.5d0*(3.0d0*cost*cost-1.0d0)
+
+ ! get ellipticity using spline evaluation
+ call spline_evaluation(rspl,espl,espl2,nspl,r,ell)
+
+! this is eq (14.4) in Dahlen and Tromp (1998)
+ factor=ONE-(TWO/3.0d0)*ell*p20
+
+ ! removes ellipticity factor
+ x = x / factor
+ y = y / factor
+ z = z / factor
+
+ end subroutine revert_ellipticity
+
diff --git a/src/specfem3D/comp_source_time_function.f90 b/src/specfem3D/comp_source_time_function.f90
index 9fff1e2..7a07d3a 100644
--- a/src/specfem3D/comp_source_time_function.f90
+++ b/src/specfem3D/comp_source_time_function.f90
@@ -41,7 +41,7 @@
else
! quasi Heaviside
comp_source_time_function = 0.5d0*(1.0d0 + netlib_specfun_erf(t/hdur))
- end if
+ endif
end function comp_source_time_function
@@ -80,7 +80,7 @@
! On the first iteration, go get the ASCII file.
if ( .not. allocated (stfArray_external) ) then
call get_external_source_time_function ()
- end if
+ endif
comp_source_time_function_external = stfArray_external (it)
@@ -112,14 +112,14 @@
if ( RetCode /= 0 ) then
print *, "ERROR IN SOURCE TIME FUNCTION."
stop
- end if
+ endif
! Ignore lines with a hash (comments)
if ( index (line, "#") /= 0 ) cycle read_loop
read (line, *) stfArray_external(iterator)
- end do read_loop
+ enddo read_loop
close (10)
diff --git a/src/specfem3D/compute_add_sources.f90 b/src/specfem3D/compute_add_sources.f90
index ca48c62..8139bd2 100644
--- a/src/specfem3D/compute_add_sources.f90
+++ b/src/specfem3D/compute_add_sources.f90
@@ -57,6 +57,7 @@
nint(gamma_source(isource)), &
ispec_selected_source(isource))
+ !! question from DK DK: not sure how the line below works, how can a duration be used as a frequency???
f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing CMTSOLUTION file format
! This is the expression of a Ricker; should be changed according maybe to the Par_file.
@@ -356,6 +357,7 @@
nint(gamma_source(isource)), &
ispec_selected_source(isource))
+ !! question from DK DK: not sure how the line below works, how can a duration be used as a frequency???
f0 = hdur(isource) !! using hdur as a FREQUENCY just to avoid changing CMTSOLUTION file format
!if (it == 1 .and. myrank == 0) then
diff --git a/src/specfem3D/get_cmt.f90 b/src/specfem3D/get_cmt.f90
index e6932b1..176a1d3 100644
--- a/src/specfem3D/get_cmt.f90
+++ b/src/specfem3D/get_cmt.f90
@@ -148,7 +148,7 @@
! If we're using external stf, don't worry about hdur.
if (EXTERNAL_SOURCE_TIME_FUNCTION) then
hdur(:) = 0.d0
- end if
+ endif
close(1)
diff --git a/src/specfem3D/locate_receivers.f90 b/src/specfem3D/locate_receivers.f90
index 6fee35e..81def24 100644
--- a/src/specfem3D/locate_receivers.f90
+++ b/src/specfem3D/locate_receivers.f90
@@ -393,8 +393,11 @@
! ellipticity
if (ELLIPTICITY_VAL) then
cost=cos(theta)
+! this is the Legendre polynomial of degree two, P2(cos(theta)), see the discussion above eq (14.4) in Dahlen and Tromp (1998)
p20=0.5d0*(3.0d0*cost*cost-1.0d0)
+! get ellipticity using spline evaluation
call spline_evaluation(rspl,espl,espl2,nspl,r0,ell)
+! this is eq (14.4) in Dahlen and Tromp (1998)
r0=r0*(1.0d0-(2.0d0/3.0d0)*ell*p20)
endif
diff --git a/src/specfem3D/locate_sources.f90 b/src/specfem3D/locate_sources.f90
index d34ec25..6e5eb57 100644
--- a/src/specfem3D/locate_sources.f90
+++ b/src/specfem3D/locate_sources.f90
@@ -318,9 +318,12 @@
endif
if (ELLIPTICITY) then
dcost = dcos(theta)
+! this is the Legendre polynomial of degree two, P2(cos(theta)), see the discussion above eq (14.4) in Dahlen and Tromp (1998)
p20 = 0.5d0*(3.0d0*dcost*dcost-1.0d0)
radius = r0 - depth(isource)*1000.0d0/R_EARTH
+! get ellipticity using spline evaluation
call spline_evaluation(rspl,espl,espl2,nspl,radius,ell)
+! this is eq (14.4) in Dahlen and Tromp (1998)
r0 = r0*(1.0d0-(2.0d0/3.0d0)*ell*p20)
endif
diff --git a/src/specfem3D/setup_sources_receivers.f90 b/src/specfem3D/setup_sources_receivers.f90
index 0a522f4..6188c0a 100644
--- a/src/specfem3D/setup_sources_receivers.f90
+++ b/src/specfem3D/setup_sources_receivers.f90
@@ -181,7 +181,7 @@
if ( EXTERNAL_SOURCE_TIME_FUNCTION ) then
hdur(:) = 0._CUSTOM_REAL
t0 = 0.d0
- end if
+ endif
! point force sources will start depending on the frequency given by hdur
if (USE_FORCE_POINT_SOURCE) then
diff --git a/src/specfem3D/write_output_SAC.f90 b/src/specfem3D/write_output_SAC.f90
index 75ec52c..6accdca 100644
--- a/src/specfem3D/write_output_SAC.f90
+++ b/src/specfem3D/write_output_SAC.f90
@@ -213,7 +213,7 @@
!USER2 = sngl(depth)
!USER3 = sngl(cmt_hdur) !half duration from CMT if not changed to t0 = 0.d0 (point source)
- ! just to avoid compiler warning
+ ! to avoid compiler warnings
value1 = elat
value1 = elon
value1 = depth
More information about the CIG-COMMITS
mailing list