[cig-commits] r8601 - seismo/2D/SPECFEM2D/trunk
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 7 16:00:07 PST 2007
Author: walter
Date: 2007-12-07 16:00:07 -0800 (Fri, 07 Dec 2007)
New Revision: 8601
Added:
seismo/2D/SPECFEM2D/trunk/spline_routines.f90
Removed:
seismo/2D/SPECFEM2D/trunk/numerical_recipes.f90
Modified:
seismo/2D/SPECFEM2D/trunk/Makefile
seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
Log:
renamed numerical_recipes.f90 to spline_routines.f90 to avoid copyright problem
Modified: seismo/2D/SPECFEM2D/trunk/Makefile
===================================================================
--- seismo/2D/SPECFEM2D/trunk/Makefile 2007-11-21 03:30:13 UTC (rev 8600)
+++ seismo/2D/SPECFEM2D/trunk/Makefile 2007-12-08 00:00:07 UTC (rev 8601)
@@ -33,7 +33,7 @@
#LIB = /opt/metis-4.0/gcc64/lib/libmetis.a /opt/scotch-4.0/gcc64/lib/libscotch.a /opt/scotch-4.0/gcc64/lib/libscotcherr.a
LIB =
-OBJS_MESHFEM2D = $O/part_unstruct.o $O/meshfem2D.o $O/read_value_parameters.o $O/numerical_recipes.o
+OBJS_MESHFEM2D = $O/part_unstruct.o $O/meshfem2D.o $O/read_value_parameters.o $O/spline_routines.o
OBJS_SPECFEM2D = $O/checkgrid.o $O/datim.o $O/enforce_acoustic_free_surface.o\
$O/compute_forces_acoustic.o $O/compute_forces_elastic.o\
@@ -154,8 +154,8 @@
$O/create_color_image.o: create_color_image.f90 constants.h
${F90} $(FLAGS_CHECK) -c -o $O/create_color_image.o create_color_image.f90
-$O/numerical_recipes.o: numerical_recipes.f90 constants.h
- ${F90} $(FLAGS_CHECK) -c -o $O/numerical_recipes.o numerical_recipes.f90
+$O/spline_routines.o: spline_routines.f90 constants.h
+ ${F90} $(FLAGS_CHECK) -c -o $O/spline_routines.o spline_routines.f90
$O/netlib_specfun_erf.o: netlib_specfun_erf.f90
${F90} $(FLAGS_CHECK) -c -o $O/netlib_specfun_erf.o netlib_specfun_erf.f90
Modified: seismo/2D/SPECFEM2D/trunk/meshfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2007-11-21 03:30:13 UTC (rev 8600)
+++ seismo/2D/SPECFEM2D/trunk/meshfem2D.F90 2007-12-08 00:00:07 UTC (rev 8601)
@@ -657,13 +657,13 @@
tang1 = (zinterface_bottom(2)-zinterface_bottom(1)) / (xinterface_bottom(2)-xinterface_bottom(1))
tangN = (zinterface_bottom(npoints_interface_bottom)-zinterface_bottom(npoints_interface_bottom-1)) / &
(xinterface_bottom(npoints_interface_bottom)-xinterface_bottom(npoints_interface_bottom-1))
- call spline(xinterface_bottom,zinterface_bottom,npoints_interface_bottom,tang1,tangN,coefs_interface_bottom)
+ call spline_construction(xinterface_bottom,zinterface_bottom,npoints_interface_bottom,tang1,tangN,coefs_interface_bottom)
! compute the spline for the top interface, impose the tangent on both edges
tang1 = (zinterface_top(2)-zinterface_top(1)) / (xinterface_top(2)-xinterface_top(1))
tangN = (zinterface_top(npoints_interface_top)-zinterface_top(npoints_interface_top-1)) / &
(xinterface_top(npoints_interface_top)-xinterface_top(npoints_interface_top-1))
- call spline(xinterface_top,zinterface_top,npoints_interface_top,tang1,tangN,coefs_interface_top)
+ call spline_construction(xinterface_top,zinterface_top,npoints_interface_top,tang1,tangN,coefs_interface_top)
! check if we are in the last layer, which contains topography,
! and modify the position of the source accordingly if it is located exactly at the surface
@@ -1320,7 +1320,7 @@
if(xp < xinterface(1)) xp = xinterface(1)
if(xp > xinterface(npoints_interface)) xp = xinterface(npoints_interface)
- call splint(xinterface,zinterface,coefs_interface,npoints_interface,xp,value_spline)
+ call spline_evaluation(xinterface,zinterface,coefs_interface,npoints_interface,xp,value_spline)
end function value_spline
Deleted: seismo/2D/SPECFEM2D/trunk/numerical_recipes.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/numerical_recipes.f90 2007-11-21 03:30:13 UTC (rev 8600)
+++ seismo/2D/SPECFEM2D/trunk/numerical_recipes.f90 2007-12-08 00:00:07 UTC (rev 8601)
@@ -1,86 +0,0 @@
-
-!=====================================================================
-!
-! Routines from "Numerical Recipes: the Art of Scientific Computing"
-! W. H. Press et al., Cambridge University Press
-!
-!=====================================================================
-
-! --------------------------------------
-
-! compute spline coefficients (Numerical Recipes)
-! modified to use dynamic allocation
-
- subroutine spline(x,y,n,yp1,ypn,y2)
-
- implicit none
-
- integer n
- double precision, dimension(n) :: x,y,y2
- double precision, dimension(:), allocatable :: u
- double precision yp1,ypn
-
- integer i,k
- double precision sig,p,qn,un
-
- allocate(u(n))
-
- y2(1)=-0.5d0
- u(1)=(3.d0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
-
- do i=2,n-1
- sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
- p=sig*y2(i-1)+2.d0
- y2(i)=(sig-1.d0)/p
- u(i)=(6.d0*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
- enddo
-
- qn=0.5d0
- un=(3.d0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
- y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.d0)
-
- do k=n-1,1,-1
- y2(k)=y2(k)*y2(k+1)+u(k)
- enddo
-
- deallocate(u)
-
- end subroutine spline
-
-! --------------
-
-! evaluate spline (adapted from Numerical Recipes)
-
- subroutine splint(xa,ya,y2a,n,x,y)
-
- implicit none
-
- integer n
- double precision, dimension(n) :: XA,YA,Y2A
- double precision x,y
-
- integer k,klo,khi
- double precision h,a,b
-
- KLO = 1
- KHI = N
-
- do while (KHI-KLO > 1)
- K=(KHI+KLO)/2
- if(XA(K) > X) then
- KHI=K
- else
- KLO=K
- endif
- enddo
-
- H = XA(KHI) - XA(KLO)
- IF (H == 0.d0) stop 'bad input in spline evaluation'
-
- A = (XA(KHI)-X) / H
- B = (X-XA(KLO)) / H
-
- Y = A*YA(KLO) + B*YA(KHI) + ((A**3-A)*Y2A(KLO) + (B**3-B)*Y2A(KHI))*(H**2)/6.d0
-
- end subroutine splint
-
Copied: seismo/2D/SPECFEM2D/trunk/spline_routines.f90 (from rev 8600, seismo/2D/SPECFEM2D/trunk/numerical_recipes.f90)
===================================================================
--- seismo/2D/SPECFEM2D/trunk/numerical_recipes.f90 2007-11-21 03:30:13 UTC (rev 8600)
+++ seismo/2D/SPECFEM2D/trunk/spline_routines.f90 2007-12-08 00:00:07 UTC (rev 8601)
@@ -0,0 +1,86 @@
+
+!=====================================================================
+!
+! Routines from "Numerical Recipes: the Art of Scientific Computing"
+! W. H. Press et al., Cambridge University Press
+!
+!=====================================================================
+
+! --------------------------------------
+
+! compute spline coefficients (Numerical Recipes)
+! modified to use dynamic allocation
+
+ subroutine spline_construction(x,y,n,yp1,ypn,y2)
+
+ implicit none
+
+ integer n
+ double precision, dimension(n) :: x,y,y2
+ double precision, dimension(:), allocatable :: u
+ double precision yp1,ypn
+
+ integer i,k
+ double precision sig,p,qn,un
+
+ allocate(u(n))
+
+ y2(1)=-0.5d0
+ u(1)=(3.d0/(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1)
+
+ do i=2,n-1
+ sig=(x(i)-x(i-1))/(x(i+1)-x(i-1))
+ p=sig*y2(i-1)+2.d0
+ y2(i)=(sig-1.d0)/p
+ u(i)=(6.d0*((y(i+1)-y(i))/(x(i+1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig*u(i-1))/p
+ enddo
+
+ qn=0.5d0
+ un=(3.d0/(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1)))
+ y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.d0)
+
+ do k=n-1,1,-1
+ y2(k)=y2(k)*y2(k+1)+u(k)
+ enddo
+
+ deallocate(u)
+
+ end subroutine spline_construction
+
+! --------------
+
+! evaluate spline (adapted from Numerical Recipes)
+
+ subroutine spline_evaluation(xa,ya,y2a,n,x,y)
+
+ implicit none
+
+ integer n
+ double precision, dimension(n) :: XA,YA,Y2A
+ double precision x,y
+
+ integer k,klo,khi
+ double precision h,a,b
+
+ KLO = 1
+ KHI = N
+
+ do while (KHI-KLO > 1)
+ K=(KHI+KLO)/2
+ if(XA(K) > X) then
+ KHI=K
+ else
+ KLO=K
+ endif
+ enddo
+
+ H = XA(KHI) - XA(KLO)
+ IF (H == 0.d0) stop 'bad input in spline evaluation'
+
+ A = (XA(KHI)-X) / H
+ B = (X-XA(KLO)) / H
+
+ Y = A*YA(KLO) + B*YA(KHI) + ((A**3-A)*Y2A(KLO) + (B**3-B)*Y2A(KHI))*(H**2)/6.d0
+
+ end subroutine spline_evaluation
+
More information about the cig-commits
mailing list