[cig-commits] r22551 - in seismo/3D/SPECFEM3D_GLOBE/trunk/src: meshfem3D specfem3D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Mon Jul 8 12:07:27 PDT 2013
Author: dkomati1
Date: 2013-07-08 12:07:27 -0700 (Mon, 08 Jul 2013)
New Revision: 22551
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_attenuation.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_s362ani.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_backazimuth.f90
Log:
improved splcon() in model_s362ani.f90 and added a comment about memory copies that can potentially slow down the code;
also removed some constants that were already defined in setup/constants.h
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_attenuation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_attenuation.f90 2013-07-08 17:56:54 UTC (rev 22550)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_attenuation.f90 2013-07-08 19:07:27 UTC (rev 22551)
@@ -572,6 +572,8 @@
implicit none
+ include "constants.h"
+
integer n
double precision tau_s(n)
double precision min_period, max_period
@@ -579,7 +581,6 @@
double precision exp1, exp2
double precision dexp
integer i
- double precision, parameter :: PI = 3.14159265358979d0
f1 = 1.0d0 / max_period
f2 = 1.0d0 / min_period
@@ -604,6 +605,8 @@
implicit none
+ include "constants.h"
+
! attenuation_simplex_variables
type attenuation_simplex_variables
sequence
@@ -633,7 +636,6 @@
integer i, iterations, err,prnt
double precision f1, f2, exp1,exp2,dexp, min_value
double precision, allocatable, dimension(:) :: f
- double precision, parameter :: PI = 3.14159265358979d0
integer, parameter :: nf = 100
double precision, external :: attenuation_eval
@@ -805,6 +807,8 @@
implicit none
+ include "constants.h"
+
! Input
integer nf, nsls
double precision, dimension(nf) :: f
@@ -813,21 +817,17 @@
double precision, dimension(nf) :: A,B
integer i,j
- double precision w, pi, demon
+ double precision w, demon
- PI = 3.14159265358979d0
-
A(:) = 1.0d0 - nsls*1.0d0
B(:) = 0.0d0
do i = 1,nf
w = 2.0d0 * PI * 10**f(i)
do j = 1,nsls
-! write(*,*)j,tau_s(j),tau_e(j)
demon = 1.0d0 + w**2 * tau_s(j)**2
A(i) = A(i) + ((1.0d0 + (w**2 * tau_e(j) * tau_s(j)))/ demon)
B(i) = B(i) + ((w * (tau_e(j) - tau_s(j))) / demon)
enddo
-! write(*,*)A(i),B(i),10**f(i)
enddo
end subroutine attenuation_maxwell
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_s362ani.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_s362ani.f90 2013-07-08 17:56:54 UTC (rev 22550)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_s362ani.f90 2013-07-08 19:07:27 UTC (rev 22551)
@@ -1095,13 +1095,15 @@
!-------------------------------------------------------------------------------------------------
!
-
subroutine splcon(xlat,xlon,nver,verlat,verlon,verrad,ncon,icon,con)
implicit none
- integer :: ncon,nver
+ include "constants.h"
+ integer, intent(in) :: nver
+ integer, intent(out) :: ncon
+
! Daniel Peter: original define
!
! real(kind=4) verlat(1)
@@ -1112,39 +1114,39 @@
! real(kind=4) con(1)
! Daniel Peter: avoiding out-of-bounds errors
- real(kind=4) verlat(nver)
- real(kind=4) verlon(nver)
- real(kind=4) verrad(nver)
+ real(kind=4), intent(in) :: verlat(nver)
+ real(kind=4), intent(in) :: verlon(nver)
+ real(kind=4), intent(in) :: verrad(nver)
integer, parameter :: maxver=1000
- integer icon(maxver)
- real(kind=4) con(maxver)
+ integer, intent(out) :: icon(maxver)
+ real(kind=4), intent(out) :: con(maxver)
double precision dd
double precision rn
double precision dr
- double precision xrad
double precision ver8
double precision xla8
integer :: iver
real(kind=4) :: xlat,xlon
- xrad=3.14159265358979/180.d0
-
ncon=0
do iver=1,nver
if(xlat > verlat(iver)-2.*verrad(iver)) then
if(xlat < verlat(iver)+2.*verrad(iver)) then
- ver8=xrad*(verlat(iver))
- xla8=xrad*(xlat)
- dd=sin(ver8)*sin(xla8)
- dd=dd+cos(ver8)*cos(xla8)* cos(xrad*(xlon-verlon(iver)))
- dd=acos(dd)/xrad
+ ver8=DEGREES_TO_RADIANS*(verlat(iver))
+ xla8=DEGREES_TO_RADIANS*(xlat)
+ dd=sin(ver8)*sin(xla8) + cos(ver8)*cos(xla8)* cos(DEGREES_TO_RADIANS*(xlon-verlon(iver)))
+ dd=acos(dd) * RADIANS_TO_DEGREES
if(dd > (verrad(iver))*2.d0) then
else
ncon=ncon+1
+
+!! DK DK added this safety test
+ if(ncon > maxver) stop 'error: ncon > maxver in splcon()'
+
icon(ncon)=iver
rn=dd/(verrad(iver))
dr=rn-1.d0
@@ -1261,6 +1263,11 @@
! nconpt(ihpa),iconpt(1,ihpa),conpt(1,ihpa))
! making sure of array bounds
+!! note from DK DK, July 2013: the code below is correct but it seems to be exactly the same as the code
+!! note from DK DK, July 2013: commented out above because xlaspl(1:numcof,ihpa) is contiguous in memory, and thus
+!! note from DK DK, July 2013: using a pointer to it such as xlaspl(1,ihpa) in the call seems fine as well;
+!! note from DK DK, July 2013: and some compilers could create expensive and useless memory copies for each call
+!! note from DK DK, July 2013: with the new syntax below
call splcon(y,x,numcof,xlaspl(1:numcof,ihpa), &
xlospl(1:numcof,ihpa),radspl(1:numcof,ihpa), &
nconpt(ihpa),iconpt(1,ihpa),conpt(1,ihpa))
@@ -1423,6 +1430,11 @@
! nconpt(ihpa),iconpt(1,ihpa),conpt(1,ihpa))
! making sure of array bounds
+!! note from DK DK, July 2013: the code below is correct but it seems to be exactly the same as the code
+!! note from DK DK, July 2013: commented out above because xlaspl(1:numcof,ihpa) is contiguous in memory, and thus
+!! note from DK DK, July 2013: using a pointer to it such as xlaspl(1,ihpa) in the call seems fine as well;
+!! note from DK DK, July 2013: and some compilers could create expensive and useless memory copies for each call
+!! note from DK DK, July 2013: with the new syntax below
call splcon(y,x,numcof,xlaspl(1:numcof,ihpa), &
xlospl(1:numcof,ihpa),radspl(1:numcof,ihpa), &
nconpt(ihpa),iconpt(1:maxver,ihpa),conpt(1:maxver,ihpa))
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_backazimuth.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_backazimuth.f90 2013-07-08 17:56:54 UTC (rev 22550)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_backazimuth.f90 2013-07-08 19:07:27 UTC (rev 22551)
@@ -30,6 +30,8 @@
implicit none
+ include "constants.h"
+
double precision the, phe
double precision ths, phs
double precision az,baz,xdeg
@@ -40,24 +42,22 @@
double precision phsrad, sc, sd, ss
double precision temp, therad, thg, thsrad
- double precision, parameter :: rad = 6378.160
- double precision, parameter :: fl = 0.00335293
- double precision, parameter :: twopideg = 360.
- double precision, parameter :: c00 = 1.
- double precision, parameter :: c01 = 0.25
- double precision, parameter :: c02 = -4.6875e-02
- double precision, parameter :: c03 = 1.953125e-02
- double precision, parameter :: c21 = -0.125
- double precision, parameter :: c22 = 3.125e-02
- double precision, parameter :: c23 = -1.46484375e-02
- double precision, parameter :: c42 = -3.90625e-03
- double precision, parameter :: c43 = 2.9296875e-03
- double precision, parameter :: degtokm = 111.3199
- double precision, parameter :: pi = 3.141592654
- double precision, parameter :: TORAD = pi/180.
- double precision, parameter :: TODEG = 1./TORAD
+ double precision, parameter :: rad = 6378.160d0
+ double precision, parameter :: fl = 0.00335293d0
+ double precision, parameter :: twopideg = 360.d0
+ double precision, parameter :: c00 = 1.d0
+ double precision, parameter :: c01 = 0.25d0
+ double precision, parameter :: c02 = -4.6875d-02
+ double precision, parameter :: c03 = 1.953125d-02
+ double precision, parameter :: c21 = -0.125d0
+ double precision, parameter :: c22 = 3.125d-02
+ double precision, parameter :: c23 = -1.46484375d-02
+ double precision, parameter :: c42 = -3.90625d-03
+ double precision, parameter :: c43 = 2.9296875d-03
+ double precision, parameter :: degtokm = 111.3199d0
+ double precision, parameter :: TORAD = DEGREES_TO_RADIANS
+ double precision, parameter :: TODEG = RADIANS_TO_DEGREES
-
!=====================================================================
! PURPOSE: To compute the distance and azimuth between locations.
!=====================================================================
@@ -104,7 +104,7 @@
! (Equations are unstable for latidudes of exactly 0 degrees.)
temp = the
- if( temp == 0. ) temp = 1.0e-08
+ if( temp == 0. ) temp = 1.0d-08
therad = TORAD*temp
pherad = TORAD*phe
@@ -130,7 +130,7 @@
! -- Convert to radians.
temp = Ths
- if( temp == 0. ) temp = 1.0e-08
+ if( temp == 0. ) temp = 1.0d-08
thsrad = TORAD*temp
phsrad = TORAD*Phs
More information about the CIG-COMMITS
mailing list