[cig-commits] r22460 - in seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils: . Profiles Visualization Visualization/VTK_ParaView Visualization/VTK_ParaView/global_slice_util/mod Visualization/VTK_ParaView/global_slice_util/obj Visualization/VTK_ParaView/mesh2vtu adjoint_sources/amplitude adjoint_sources/traveltime attenuation extract_database funaro_routines_bernhard
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sun Jun 30 07:11:15 PDT 2013
Author: dkomati1
Date: 2013-06-30 07:11:14 -0700 (Sun, 30 Jun 2013)
New Revision: 22460
Added:
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/Visualization/VTK_ParaView/
Removed:
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/Visualization/Paraview/
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/Visualization/VTK/
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/chunk_notes_scanned/
Modified:
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/Profiles/write_profile.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/Visualization/VTK_ParaView/global_slice_util/mod/
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/Visualization/VTK_ParaView/global_slice_util/obj/
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/Visualization/VTK_ParaView/mesh2vtu/
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/adjoint_sources/amplitude/
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/adjoint_sources/traveltime/
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/attenuation/attenuation.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/attenuation/attenuation_output.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/attenuation/attenuation_simplex.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/extract_database/extract_database.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/funaro_routines_bernhard/routines_SEM_Bernhard_diff_lagrange.f90
Log:
merged "utils" from the trunk
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/Profiles/write_profile.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/Profiles/write_profile.f90 2013-06-30 14:00:04 UTC (rev 22459)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/Profiles/write_profile.f90 2013-06-30 14:11:14 UTC (rev 22460)
@@ -431,7 +431,7 @@
! make sure that the first point of the profile is at zero and the last is at the surface
if(islice == 1) then
r_prem = rmin
- elseif(islice == nit+1) then
+ else if(islice == nit+1) then
r_prem = rmax
endif
@@ -446,7 +446,7 @@
! print *,'rprem before: ',r_prem*R_EARTH
r_prem = r_prem*(ONE + gamma * (elevation/R_EARTH) /r_prem)
! print *,'r_prem after: ',r_prem*R_EARTH
-elseif((.not. CRUSTAL) .and. (ROCEAN < R_EARTH)) then
+else if((.not. CRUSTAL) .and. (ROCEAN < R_EARTH)) then
r_prem = ROCEAN/R_EARTH
endif
! end add_topography
@@ -481,7 +481,7 @@
if(OCEANS .and. elevation < -500.0) then
! iline_ocean = iline
nlayers_ocean = floor(-elevation/500.0d0)
- elseif((.not. CRUSTAL) .and. (ROCEAN < R_EARTH)) then
+ else if((.not. CRUSTAL) .and. (ROCEAN < R_EARTH)) then
nlayers_ocean = floor((R_EARTH - ROCEAN)/500.0d0)
endif
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/Visualization/VTK_ParaView
___________________________________________________________________
Added: svn:ignore
+ xglobal_slice_number
xnormal_plane
xmake_az_stations
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/adjoint_sources/amplitude
___________________________________________________________________
Added: svn:ignore
+ *.o
xcreate_adjsrc_amplitude
Property changes on: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/adjoint_sources/traveltime
___________________________________________________________________
Added: svn:ignore
+ *.o
xcreate_adjsrc_traveltime
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/attenuation/attenuation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/attenuation/attenuation.f90 2013-06-30 14:00:04 UTC (rev 22459)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/attenuation/attenuation.f90 2013-06-30 14:11:14 UTC (rev 22460)
@@ -9,7 +9,7 @@
!
! Brian Savage 19/01/05
! This code should not produce the exact values that attenuation_prem.c
-! It has been updated to use a more robust inversion for the
+! It has been updated to use a more robust inversion for the
! the stress and strain relaxation times
! A simplex inversion is used
!
@@ -36,14 +36,14 @@
write(*,*)
! write(*,*)'number of mechanisms: '
read(5,*, end=13)n
-42 continue
+42 continue
! write(*,*)'Q: '
read(5,*, end=13)Q
-
+
tau_e(:) = 0.0d0
tau_s(:) = 0.0d0
omega_not = 0.0d0
-
+
! call attenuation_liu(t1, t2, n, Q, omega_not, tau_s, tau_e)
call attenuation_simplex(t1, t2, n, Q, omega_not, tau_s, tau_e)
if(write_central_period == 0) then
@@ -53,7 +53,7 @@
write(*,*)'! tau sigma evenly spaced in log frequency, does not depend on value of Q'
do i = 1,n
write(*,'(A13,I1,A4,F30.20,A2)')' tau_sigma(',i,') = ', tau_s(i), 'd0'
- end do
+ enddo
write(*,*)
write(*,*)"! check in which region we are based upon doubling flag"
write(*,*)
@@ -67,41 +67,41 @@
write(*,*)
write(*,*)' case(IREGION_ATTENUATION_INNER_CORE)'
write(*,*)
- end if
+ endif
if(Q == 312.0d0) then
write(*,*)'!--- CMB -> d670 (no attenuation in fluid outer core), target Q_mu = 312.'
write(*,*)
write(*,*)' case(IREGION_ATTENUATION_CMB_670)'
write(*,*)
- end if
+ endif
if(Q == 143.0d0) then
write(*,*)'!--- d670 -> d220, target Q_mu: 143.'
write(*,*)
write(*,*)' case(IREGION_ATTENUATION_670_220)'
write(*,*)
- end if
+ endif
if(Q == 80.0d0) then
write(*,*)'!--- d220 -> depth of 80 km, target Q_mu: 80.'
write(*,*)
write(*,*)' case(IREGION_ATTENUATION_220_80)'
write(*,*)
- end if
+ endif
if(Q == 600.0d0) then
write(*,*)'!--- depth of 80 km -> surface, target Q_mu: 600.'
write(*,*)
write(*,*)' case(IREGION_ATTENUATION_80_SURFACE)'
write(*,*)
- end if
-
+ endif
+
do i = 1,n
write(*,'(A12,I1,A4,F30.20,A2)')' tau_mu(',i,') = ', tau_e(i), 'd0'
- end do
+ enddo
write(*,'(A17,F20.10,A2)')' Q_mu = ', Q, 'd0'
write(*,*)
goto 42
-
-
+
+
13 continue
write(*,*)'!--- do nothing for fluid outer core (no attenuation there)'
write(*,*)
@@ -124,12 +124,12 @@
real(8), dimension(N_SLS) :: tauinv
tauinv(:) = - 1.0 / tau_s(:)
-
+
alphaval(:) = 1 + deltat*tauinv(:) + deltat**2*tauinv(:)**2 / 2. + &
deltat**3*tauinv(:)**3 / 6. + deltat**4*tauinv(:)**4 / 24.
betaval(:) = deltat / 2. + deltat**2*tauinv(:) / 3. + deltat**3*tauinv(:)**2 / 8. + deltat**4*tauinv(:)**3 / 24.
gammaval(:) = deltat / 2. + deltat**2*tauinv(:) / 6. + deltat**3*tauinv(:)**2 / 24.0
-
+
end subroutine attenuation_memory_values
subroutine attenuation_scale_factor(myrank, T_c_source, tau_mu, tau_sigma, Q_mu, scale_factor)
@@ -160,7 +160,7 @@
scale_t = ONE/dsqrt(PI*GRAV*RHOAV)
T_c_source_nondim = T_c_source / scale_t
-
+
!--- compute central angular frequency of source (non dimensionalized)
f_c_source = ONE / T_c_source_nondim
w_c_source = TWO_PI * f_c_source
@@ -174,14 +174,14 @@
!--- compute a, b and Omega parameters, also compute one minus sum of betas
a_val = ONE
b_val = ZERO
-
+
do i = 1,N_SLS
a_val = a_val - w_c_source * w_c_source * tau_mu(i) * &
(tau_mu(i) - tau_sigma(i)) / (1.d0 + w_c_source * w_c_source * tau_mu(i) * tau_mu(i))
b_val = b_val + w_c_source * (tau_mu(i) - tau_sigma(i)) / &
(1.d0 + w_c_source * w_c_source * tau_mu(i) * tau_mu(i))
enddo
-
+
big_omega = a_val*(sqrt(1.d0 + b_val*b_val/(a_val*a_val))-1.d0);
!--- quantity by which to scale mu to get mu_relaxed
@@ -286,7 +286,7 @@
expo = exp1 - dexp
do i=1,100
expo = expo + dexp
- df = 10.0d0**(expo+dexp) - 10.0d0**(expo)
+ df = 10.0d0**(expo+dexp) - 10.0d0**(expo)
omega = PI2 * 10.0d0**(expo)
write(*,*) 'df,expo,omega: ',df,expo,omega
a = real(1.0d0 - n)
@@ -332,7 +332,7 @@
subroutine invert(x,b,A,n)
implicit none
-
+
integer n
real(8), dimension(n) :: x, b
real(8), dimension(n,n) :: A
@@ -368,10 +368,10 @@
enddo
end subroutine invert
-
-
+
+
FUNCTION pythag_dp(a,b)
IMPLICIT NONE
INTEGER, PARAMETER :: DP = KIND(1.0D0)
@@ -387,8 +387,8 @@
pythag_dp=0.0d0
else
pythag_dp=absb*sqrt(1.0d0+(absa/absb)**2)
- end if
- end if
+ endif
+ endif
END FUNCTION pythag_dp
SUBROUTINE svdcmp_dp(a,w,v,p)
@@ -432,8 +432,8 @@
a(i:m,l:n)=a(i:m,l:n)+spread(a(i:m,i),dim=2,ncopies=size(tempn(l:n))) * &
spread(tempn(l:n),dim=1,ncopies=size(a(i:m,i)))
a(i:m,i)=scale*a(i:m,i)
- end if
- end if
+ endif
+ endif
w(i)=scale*g
g=0.0d0
scale=0.0d0
@@ -452,9 +452,9 @@
a(l:m,l:n)=a(l:m,l:n)+spread(tempm(l:m),dim=2,ncopies=size(rv1(l:n))) * &
spread(rv1(l:n),dim=1,ncopies=size(tempm(l:m)))
a(i,l:n)=scale*a(i,l:n)
- end if
- end if
- end do
+ endif
+ endif
+ enddo
anorm=maxval(abs(w)+abs(rv1))
! write(*,*)W
do i=n,1,-1
@@ -465,14 +465,14 @@
! v(l:n,l:n)=v(l:n,l:n)+outerprod_d(v(l:n,i),n-1+1,tempn(l:n),n-l+1)
v(l:n,l:n)=v(l:n,l:n)+spread(v(l:n,i),dim=2,ncopies=size(tempn(l:n))) * &
spread(tempn(l:n), dim=1, ncopies=size(v(l:n,i)))
- end if
+ endif
v(i,l:n)=0.0d0
v(l:n,i)=0.0d0
- end if
+ endif
v(i,i)=1.0d0
g=rv1(i)
l=i
- end do
+ enddo
do i=min(m,n),1,-1
l=i+1
g=w(i)
@@ -486,9 +486,9 @@
a(i:m,i)=a(i:m,i)*g
else
a(i:m,i)=0.0d0
- end if
+ endif
a(i,i)=a(i,i)+1.0d0
- end do
+ enddo
do k=n,1,-1
do its=1,30
do l=k,1,-1
@@ -510,19 +510,19 @@
tempm(1:m)=a(1:m,nm)
a(1:m,nm)=a(1:m,nm)*c+a(1:m,i)*s
a(1:m,i)=-tempm(1:m)*s+a(1:m,i)*c
- end do
+ enddo
exit
- end if
- end do
+ endif
+ enddo
z=w(k)
if (l == k) then
if (z < 0.0d0) then
w(k)=-z
v(1:n,k)=-v(1:n,k)
- end if
+ endif
exit
- end if
- if (its == 30) then
+ endif
+ if (its == 30) then
write(*,*) 'svdcmp_dp: no convergence in svdcmp'
call exit(-1)
endif
@@ -559,17 +559,17 @@
z=1.0d0/z
c=f*z
s=h*z
- end if
+ endif
f= (c*g)+(s*y)
x=-(s*g)+(c*y)
tempm(1:m)=a(1:m,j)
a(1:m,j)=a(1:m,j)*c+a(1:m,i)*s
a(1:m,i)=-tempm(1:m)*s+a(1:m,i)*c
- end do
+ enddo
rv1(l)=0.0d0
rv1(k)=f
w(k)=x
- end do
- end do
+ enddo
+ enddo
END SUBROUTINE svdcmp_dp
-
+
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/attenuation/attenuation_output.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/attenuation/attenuation_output.f90 2013-06-30 14:00:04 UTC (rev 22459)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/attenuation/attenuation_output.f90 2013-06-30 14:11:14 UTC (rev 22460)
@@ -12,10 +12,10 @@
! However, the values here should match those which are output from attenuation_prem.c
! This code was used for testing the implementation of 3D attenuation
! Brian Savage, 22/03/04
- implicit none
+ implicit none
include "OUTPUT_FILES/values_from_mesher.h"
include "constants.h"
-
+
integer i,j,k,ispec
integer myrank, vnspec, process, iregion
character(len=150) prname, LOCAL_PATH
@@ -23,7 +23,7 @@
double precision, dimension(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_AC) :: factor_common
double precision, dimension(N_SLS) :: tau_s
double precision T_c_source
-
+
LOCAL_PATH = "/scratch/DATABASES_MPI_BRIAN/att"
process = 42
iregion = IREGION_CRUST_MANTLE
@@ -33,11 +33,11 @@
open(unit=27, file=prname(1:len_trim(prname))//'tau_s.bin',status='old',form='unformatted')
read(27) tau_s
close(27)
-
+
open(unit=27, file=prname(1:len_trim(prname))//'T_c_source.bin',status='old',form='unformatted')
read(27) T_c_source
close(27);
-
+
open(unit=27, file=prname(1:len_trim(prname))//'Q.bin',status='old',form='unformatted')
read(27) scale_factor
close(27)
@@ -50,7 +50,7 @@
write(*,*)' tau_sigma(1) = ', tau_s(1)
write(*,*)' tau_sigma(2) = ', tau_s(2)
write(*,*)' tau_sigma(3) = ', tau_s(3)
-
+
do ispec = 1, NSPEC_CRUST_MANTLE_AC
do k = 1, NGLLZ
do j = 1, NGLLY
@@ -65,5 +65,5 @@
enddo
enddo
end program attenuation_output
-
-
+
+
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/attenuation/attenuation_simplex.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/attenuation/attenuation_simplex.f90 2013-06-30 14:00:04 UTC (rev 22459)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/attenuation/attenuation_simplex.f90 2013-06-30 14:11:14 UTC (rev 22460)
@@ -2,7 +2,7 @@
module attenuation_simplex_variables
implicit none
-
+
! nf = Number of Frequencies
! nsls = Number of Standard Linear Solids
integer nf, nsls
@@ -14,7 +14,7 @@
! iQ = 1/Q
real(8) Q, iQ
- ! tau_s = Tau_sigma defined by the frequency range and
+ ! tau_s = Tau_sigma defined by the frequency range and
! number of standard linear solids
real(8), allocatable, dimension(:) :: tau_s
@@ -49,11 +49,11 @@
! Determine the min and max frequencies
f1 = 1.0d0 / t1
f2 = 1.0d0 / t2
-
+
! Determine the exponents of the frequencies
exp1 = log10(f1);
exp2 = log10(f2);
-
+
if(f2 < f1 .OR. Q_real < 0.0d0 .OR. n < 1) then
write(*,*)'bad parameters'
call exit(-1)
@@ -61,30 +61,30 @@
! Determine the Source frequency
omega_not = 1.0e+03 * 10.0d0**(0.5 * (log10(f1) + log10(f2)))
-
- ! Determine the Frequencies at which to compare solutions
+
+ ! Determine the Frequencies at which to compare solutions
! The frequencies should be equally spaced in log10 frequency
do i = 1,nf
f(i) = exp1 + ((i-1)*1.0d0 * (exp2-exp1) / ((nf-1)*1.0d0))
- end do
+ enddo
! Set the Tau_sigma (tau_s) to be equally spaced in log10 frequency
dexp = (exp2-exp1) / ((n*1.0d0) - 1)
do i = 1,n
tau_s(i) = 1.0 / (PI * 2.0d0 * 10**(exp1 + (i - 1)* 1.0d0 *dexp))
- end do
+ enddo
- ! Shove the paramters into the module
+ ! Shove the paramters into the module
call attenuation_simplex_setup(nf,n,f,Q_real,tau_s)
-
+
! Set the Tau_epsilon (tau_e) to an initial value
! at omega*tau = 1; tan_delta = 1/Q = (tau_e - tau_s)/(2 * sqrt(tau e*tau_s))
! if we assume tau_e =~ tau_s
! we get the equation below
do i = 1,n
tau_e(i) = tau_s(i) + (tau_s(i) * 2.0d0/Q_real)
- end do
-
+ enddo
+
! Run a simplex search to determine the optimum values of tau_e
call fminsearch(attenuation_eval, tau_e, n, iterations, min_value, prnt, err)
if(err > 0) then
@@ -93,11 +93,11 @@
write(*,*)' Min Value: ', min_value
write(*,*)' Aborting program'
call exit(-1)
- end if
-
+ endif
+
deallocate(f)
call attenuation_simplex_finish()
-
+
end subroutine attenuation_simplex
subroutine attenuation_simplex_finish()
@@ -116,7 +116,7 @@
subroutine attenuation_simplex_setup(nf_in,nsls_in,f_in,Q_in,tau_s_in)
use attenuation_simplex_variables
implicit none
-
+
integer nf_in, nsls_in
real(8) Q_in
real(8), dimension(nf_in) :: f_in
@@ -136,14 +136,14 @@
!!!!!!!!
! subroutine attenuation_maxwell
-! - Computes the Moduli (Maxwell Solid) for a series of
+! - Computes the Moduli (Maxwell Solid) for a series of
! Standard Linear Solids
! - Computes M1 and M2 parameters after Dahlen and Tromp pp.203
! here called B and A after Liu et al. 1976
! - Another formulation uses Kelvin-Voigt Solids and computes
! Compliences J1 and J2 after Dahlen and Tromp pp.203
!
-! Input
+! Input
! nf = Number of Frequencies
! nsls = Number of Standard Linear Solids
! f = Frequencies (in log10 of frequencies)
@@ -188,7 +188,7 @@
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)
- end do
+ enddo
! write(*,*)A(i),B(i),10**f(i)
enddo
@@ -196,42 +196,42 @@
!!!!!!!!
! subroutine attenuation_eval
-! - Computes the misfit from a set of relaxation paramters
+! - Computes the misfit from a set of relaxation paramters
! given a set of frequencies and target attenuation
! - Evaluates only at the given frequencies
! - Evaluation is done with an L2 norm
!
! Input
! Xin = Tau_epsilon, Strain Relaxation Time
-! Note: Tau_sigma the Stress Relaxation Time is loaded
+! Note: Tau_sigma the Stress Relaxation Time is loaded
! with attenuation_simplex_setup and stored in
! attenuation_simplex_variables
-!
+!
! Xi = Sum_i^N sqrt [ (1/Qc_i - 1/Qt_i)^2 / 1/Qt_i^2 ]
-!
+!
! where Qc_i is the computed attenuation at a specific frequency
! Qt_i is the desired attenuaiton at that frequency
-!
+!
! Uses attenuation_simplex_variables to store constant values
!
! See atteunation_simplex_setup
-!
+!
real(8) function attenuation_eval(Xin)
use attenuation_simplex_variables
implicit none
! Input
real(8), dimension(nsls) :: Xin
real(8), dimension(nsls) :: tau_e
-
+
real(8), dimension(nf) :: A, B, tan_delta
-
+
integer i
real(8) xi, iQ2
tau_e = Xin
-
+
call attenuation_maxwell(nf,nsls,f,tau_s,tau_e,B,A)
-
+
tan_delta = B / A
attenuation_eval = 0.0d0
@@ -239,7 +239,7 @@
do i = 1,nf
xi = sqrt(( ( (tan_delta(i) - iQ) ** 2 ) / iQ2 ))
attenuation_eval = attenuation_eval + xi
- end do
+ enddo
end function attenuation_eval
@@ -259,20 +259,20 @@
! Output: Mimimized Value
! n = number of variables
! itercount = Input/Output
-! Input: maximum number of iterations
+! Input: maximum number of iterations
! if < 0 default is used (200 * n)
! Output: total number of iterations on output
! tolf = Input/Output
! Input: minimium tolerance of the function funk(x)
! Output: minimium value of funk(x)(i.e. "a" solution)
! prnt = Input
-! 3 => report every iteration
+! 3 => report every iteration
! 4 => report every iteration, total simplex
! err = Output
! 0 => Normal exeecution, converged within desired range
! 1 => Function Evaluation exceeded limit
! 2 => Iterations exceeded limit
-!
+!
! See Matlab fminsearch
!
subroutine fminsearch(funk, x, n, itercount, tolf, prnt, err)
@@ -296,7 +296,7 @@
integer, parameter :: contract_outside = 4
integer, parameter :: contract_inside = 5
integer, parameter :: shrink = 6
-
+
integer maxiter, maxfun
integer func_evals
real(8) tolx
@@ -318,86 +318,86 @@
if(itercount > 0) then
maxiter = itercount
- else
+ else
maxiter = 200 * n
- end if
+ endif
itercount = 0
maxfun = 200 * n
-
+
if(tolf > 0.0d0) then
tolx = 1.0e-4
else
tolx = 1.0e-4
tolf = 1.0e-4
- end if
+ endif
err = 0
xin = x
v(:,:) = 0.0d0
fv(:) = 0.0d0
-
+
v(:,1) = xin
x = xin
-
+
fv(1) = funk(xin)
-
+
usual_delta = 0.05
zero_term_delta = 0.00025
do j = 1,n
y = xin
- if(y(j) .NE. 0.0d0) then
+ if(y(j) /= 0.0d0) then
y(j) = (1.0d0 + usual_delta) * y(j)
else
y(j) = zero_term_delta
- end if
+ endif
v(:,j+1) = y
x(:) = y
fv(j+1) = funk(x)
- end do
-
- call qsort(fv,n+1,place)
+ enddo
+ call bubble_sort(fv,n+1,place)
+
do i = 1,n+1
vtmp(:,i) = v(:,place(i))
- end do
+ enddo
v = vtmp
-
+
how = initial
itercount = 1
func_evals = n+1
- if(prnt .EQ. 3) then
+ if(prnt == 3) then
write(*,*)'Iterations Funk Evals Value How'
write(*,*)itercount, func_evals, fv(1), how
endif
- if(prnt .EQ. 4) then
+ if(prnt == 4) then
write(*,*)'How: ',how
write(*,*)'V: ', v
write(*,*)'fv: ',fv
write(*,*)'evals: ',func_evals
- end if
+ endif
- do while (func_evals < maxfun .AND. itercount < maxiter)
+ do while (func_evals < maxfun .AND. itercount < maxiter)
-! if(max(max(abs(v(:,2:n+1) - v(:,1)))) .LE. tolx .AND. &
-! max(abs(fv(1) - fv(2:n+1))) .LE. tolf) then
+! if(max(max(abs(v(:,2:n+1) - v(:,1)))) <= tolx .AND. &
+! max(abs(fv(1) - fv(2:n+1))) <= tolf) then
- if(max_size_simplex(v,n) .LE. tolx .AND. &
- max_value(fv,n+1) .LE. tolf) then
+ if(max_size_simplex(v,n) <= tolx .AND. &
+ max_value(fv,n+1) <= tolf) then
goto 666
- end if
+ endif
how = none
-
+
! xbar = average of the n (NOT n+1) best points
! xbar = sum(v(:,1:n), 2)/n
xbar(:) = 0.0d0
do i = 1,n
do j = 1,n
xbar(i) = xbar(i) + v(i,j)
- end do
+ enddo
xbar(i) = xbar(i) / (n*1.0d0)
- end do
+ enddo
xr = (1 + rho)*xbar - rho*v(:,n+1)
x(:) = xr
fxr = funk(x)
@@ -413,36 +413,36 @@
fv(n+1) = fxe
how = expand
else
- v(:,n+1) = xr
+ v(:,n+1) = xr
fv(n+1) = fxr
how = reflect
- end if
+ endif
else ! fv(:,1) <= fxr
- if (fxr < fv(n)) then
- v(:,n+1) = xr
+ if (fxr < fv(n)) then
+ v(:,n+1) = xr
fv(n+1) = fxr
how = reflect
- else ! fxr >= fv(:,n)
+ else ! fxr >= fv(:,n)
! Perform contraction
if (fxr < fv(n+1)) then
! Perform an outside contraction
xc = (1 + psi*rho)*xbar - psi*rho*v(:,n+1)
- x(:) = xc
+ x(:) = xc
fxc = funk(x)
func_evals = func_evals+1
-
+
if (fxc <= fxr) then
- v(:,n+1) = xc
+ v(:,n+1) = xc
fv(n+1) = fxc
how = contract_outside
else
! perform a shrink
how = shrink
- end if
+ endif
else
! Perform an inside contraction
xcc = (1-psi)*xbar + psi*v(:,n+1)
- x(:) = xcc
+ x(:) = xcc
fxcc = funk(x)
func_evals = func_evals+1
@@ -453,45 +453,45 @@
else
! perform a shrink
how = shrink
- end if
- end if
- if (how .EQ. shrink) then
+ endif
+ endif
+ if (how == shrink) then
do j=2,n+1
v(:,j)=v(:,1)+sigma*(v(:,j) - v(:,1))
- x(:) = v(:,j)
+ x(:) = v(:,j)
fv(j) = funk(x)
- end do
+ enddo
func_evals = func_evals + n
- end if
- end if
- end if
+ endif
+ endif
+ endif
- call qsort(fv,n+1,place)
+ call bubble_sort(fv,n+1,place)
do i = 1,n+1
vtmp(:,i) = v(:,place(i))
- end do
+ enddo
v = vtmp
-
+
itercount = itercount + 1
if (prnt == 3) then
write(*,*)itercount, func_evals, fv(1), how
- elseif (prnt == 4) then
+ else if (prnt == 4) then
write(*,*)
write(*,*)'How: ',how
write(*,*)'v: ',v
write(*,*)'fv: ',fv
write(*,*)'evals: ',func_evals
- end if
- end do
+ endif
+ enddo
if(func_evals > maxfun) then
write(*,*)'function evaluations exceeded prescribed limit', maxfun
err = 1
- end if
+ endif
if(itercount > maxiter) then
write(*,*)'iterations exceeded prescribed limit', maxiter
err = 2
- end if
+ endif
666 continue
x = v(:,1)
@@ -510,7 +510,7 @@
! dimension(n)
! n = Input
! Length of fv
-!
+!
! Returns:
! Xi = max( || fv(1)- fv(i) || ); i=2:n
!
@@ -521,14 +521,14 @@
integer i
real(8) m, z
-
+
m = 0.0d0
do i = 2,n
z = abs(fv(1) - fv(i))
if(z > m) then
m = z
- end if
- end do
+ endif
+ enddo
max_value = m
@@ -542,7 +542,7 @@
! Simplex Verticies
! dimension(n, n+1)
! n = Pseudo Length of n
-!
+!
! Returns:
! Xi = max( max( || v(:,1) - v(:,i) || ) ) ; i=2:n+1
!
@@ -553,24 +553,23 @@
integer i,j
real(8) m, z
-
+
m = 0.0d0
do i = 1,n
do j = 2,n+1
z = abs(v(i,j) - v(i,1))
if(z > m) then
m = z
- end if
- end do
- end do
-
+ endif
+ enddo
+ enddo
+
max_size_simplex = m
end function max_size_simplex
!!!!!!!
-! subroutine qsort
! - Implementation of a Bubble Sort Routine
! Input
! X = Input/Output
@@ -588,19 +587,19 @@
! X = [ 1 2 3 4 ] on Output
! I = [ 3 4 2 1 ] on Output
!
-subroutine qsort(X,n,I)
+subroutine bubble_sort(X,n,I)
implicit none
integer n
real(8) X(n)
integer I(n)
-
+
integer j,k
real(8) rtmp
integer itmp
do j = 1,n
I(j) = j
- end do
+ enddo
do j = 1,n
do k = 1,n-j
@@ -608,15 +607,15 @@
rtmp = X(k)
X(k) = X(k+1)
X(k+1) = rtmp
-
+
itmp = I(k)
I(k) = I(k+1)
I(k+1) = itmp
- end if
+ endif
enddo
enddo
-end subroutine qsort
+end subroutine bubble_sort
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/extract_database/extract_database.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/extract_database/extract_database.f90 2013-06-30 14:00:04 UTC (rev 22459)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/extract_database/extract_database.f90 2013-06-30 14:11:14 UTC (rev 22460)
@@ -16,12 +16,12 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: junk
real(kind=CUSTOM_REAL) :: junk2
- call getarg(1,infile)
- call getarg(2,s_ireg)
- call getarg(3,s_num)
- call getarg(4,outfile)
- if (trim(infile) == '' .or. trim(s_ireg) == '' .or. trim(s_num) == '' &
- .or. trim(outfile) == '') then
+ call get_command_argument(1,infile)
+ call get_command_argument(2,s_ireg)
+ call get_command_argument(3,s_num)
+ call get_command_argument(4,outfile)
+ if (len_trim(infile) == 0 .or. len_trim(s_ireg) == 0 .or. len_trim(s_num) == 0 &
+ .or. len_trim(outfile) == 0) then
print *, 'Usage: extract_databases infile ireg num outfile'
print *, ' ireg = 1, 2, 3'
print *, ' num = 10 for rho, 11 for kappav, 12 for muv '
@@ -47,7 +47,7 @@
enddo
read(11) junk(:,:,:,1:nspec)
close(11)
-
+
open(12,file=outfile,status='unknown',form='unformatted')
write(12) junk(:,:,:,1:nspec)
close(12)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/funaro_routines_bernhard/routines_SEM_Bernhard_diff_lagrange.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/funaro_routines_bernhard/routines_SEM_Bernhard_diff_lagrange.f90 2013-06-30 14:00:04 UTC (rev 22459)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/utils/funaro_routines_bernhard/routines_SEM_Bernhard_diff_lagrange.f90 2013-06-30 14:11:14 UTC (rev 22460)
@@ -1,6 +1,6 @@
SUBROUTINE diff_lagrange(N,DL,rho,xi)
!
-!calculates the derivatives of the Lagrange Polynomials at the GLL quadrature
+!calculates the derivatives of the Lagrange Polynomials at the GLL quadrature
!points and stores them in Matrix DL (coloumn = d (L)/d (xi) = dL, row = dL at
!point xi
@@ -8,12 +8,12 @@
double precision xi(0:N),LL(0:N),VN(0:N),rho(0:N)
double precision DL(0:N,0:N)
-integer N,i,j,k
+integer N,i,j,k
!Calculate the Gauß-Lobatto-Legendre Quadrature points and the values of the
-!Legendre polynomials VN at the gll nodes
+!Legendre polynomials VN at the gll nodes
call main_pol(xi,N,VN)
@@ -49,7 +49,7 @@
ENDDO
DO j=0,N
- xi(j)=ET(j)
+ xi(j)=ET(j)
!WRITE(6,*)ET(j),VN(j)
ENDDO
@@ -61,66 +61,66 @@
!****************************************************************************
SUBROUTINE VALEPO(N,X,Y,DY,D2Y)
-
+
!* COMPUTES THE VALUE OF THE LEGENDRE POLYNOMIAL OF DEGREE N
-!* AND ITS FIRST AND SECOND DERIVATIVES AT A GIVEN POINT
+!* AND ITS FIRST AND SECOND DERIVATIVES AT A GIVEN POINT
!* N = DEGREE OF THE POLYNOMIAL
-!* X = POINT IN WHICH THE COMPUTATION IS PERFORMED
+!* X = POINT IN WHICH THE COMPUTATION IS PERFORMED
!* Y = VALUE OF THE POLYNOMIAL IN X
!* DY = VALUE OF THE FIRST DERIVATIVE IN X
-!* D2Y= VALUE OF THE SECOND DERIVATIVE IN X
+!* D2Y= VALUE OF THE SECOND DERIVATIVE IN X
-! IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+! IMPLICIT DOUBLE PRECISION (A-H,O-Z)
IMPLICIT NONE
DOUBLE PRECISION Y,DY,D2Y,X,YP,DYP,D2YP,C1,C2,C3,C4,YM,DYM,D2YM
- INTEGER N,I
+ INTEGER N,I
- Y = 1.
+ Y = 1.
DY = 0.
- D2Y = 0.
+ D2Y = 0.
- IF (N .EQ. 0) RETURN
+ IF (N == 0) RETURN
- Y = X
+ Y = X
DY = 1.
- D2Y = 0.
+ D2Y = 0.
- IF(N .EQ. 1) RETURN
+ IF(N == 1) RETURN
- YP = 1.
+ YP = 1.
DYP = 0.
- D2YP = 0.
+ D2YP = 0.
- DO I=2,N
+ DO I=2,N
C1 = DFLOAT(I)
- C2 = 2.*C1-1.
+ C2 = 2.*C1-1.
C4 = C1-1.
- YM = Y
+ YM = Y
Y = (C2*X*Y-C4*YP)/C1
- YP = YM
+ YP = YM
DYM = DY
- DY = (C2*X*DY-C4*DYP+C2*YP)/C1
+ DY = (C2*X*DY-C4*DYP+C2*YP)/C1
DYP = DYM
- D2YM = D2Y
+ D2YM = D2Y
D2Y = (C2*X*D2Y-C4*D2YP+2.D0*C2*DYP)/C1
- D2YP = D2YM
+ D2YP = D2YM
ENDDO
-
+
RETURN
-END SUBROUTINE VALEPO
+END SUBROUTINE VALEPO
!****************************************************************************
SUBROUTINE ZELEGL(N,ET,VN)
-
+
!* COMPUTES THE NODES RELATIVE TO THE LEGENDRE GAUSS-LOBATTO FORMULA
!* N = ORDER OF THE FORMULA
!* ET = VECTOR OF THE NODES, ET(I), I=0,N
!* VN = VALUES OF THE LEGENDRE POLYNOMIAL AT THE NODES, VN(I), I=0,N
-
+
! IMPLICIT DOUBLE PRECISION (A-H,O-Z)
IMPLICIT NONE
@@ -130,21 +130,21 @@
! DIMENSION ET(0:*), VN(0:*)
- IF (N .EQ. 0) RETURN
+ IF (N == 0) RETURN
- N2 = (N-1)/2
+ N2 = (N-1)/2
SN = DFLOAT(2*N-4*N2-3)
ET(0) = -1.
ET(N) = 1.
VN(0) = SN
VN(N) = 1.
- IF (N .EQ. 1) RETURN
+ IF (N == 1) RETURN
ET(N2+1) = 0.
X = 0.
CALL VALEPO(N,X,Y,DY,D2Y)
VN(N2+1) = Y
- IF(N .EQ. 2) RETURN
+ IF(N == 2) RETURN
PI = 3.14159265358979323846
C = PI/DFLOAT(N)
@@ -164,44 +164,44 @@
END SUBROUTINE ZELEGL
-!***********************************************************************
-
- SUBROUTINE DELEGL(N,ET,VN,QN,DQN)
+!***********************************************************************
-!************************************************************************
-!* COMPUTES THE DERIVATIVE OF A POLYNOMIAL AT THE LEGENDRE GAUSS-LOBATTO
-!* NODES FROM THE VALUES OF THE POLYNOMIAL ATTAINED AT THE SAME POINTS
-!* N = THE DEGREE OF THE POLYNOMIAL
-!* ET = VECTOR OF THE NODES, ET(I), I=0,N
-!* VN = VALUES OF THE LEGENDRE POLYNOMIAL AT THE NODES, VN(I), I=0,N
-!* QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N
-!* DQN = DERIVATIVES OF THE POLYNOMIAL AT THE NODES, DQZ(I), I=0,N
-!************************************************************************
- IMPLICIT DOUBLE PRECISION (A-H,O-Z)
-! DIMENSION ET(0:*), VN(0:*), QN(0:*), DQN(0:*)
+ SUBROUTINE DELEGL(N,ET,VN,QN,DQN)
+
+!************************************************************************
+!* COMPUTES THE DERIVATIVE OF A POLYNOMIAL AT THE LEGENDRE GAUSS-LOBATTO
+!* NODES FROM THE VALUES OF THE POLYNOMIAL ATTAINED AT THE SAME POINTS
+!* N = THE DEGREE OF THE POLYNOMIAL
+!* ET = VECTOR OF THE NODES, ET(I), I=0,N
+!* VN = VALUES OF THE LEGENDRE POLYNOMIAL AT THE NODES, VN(I), I=0,N
+!* QN = VALUES OF THE POLYNOMIAL AT THE NODES, QN(I), I=0,N
+!* DQN = DERIVATIVES OF THE POLYNOMIAL AT THE NODES, DQZ(I), I=0,N
+!************************************************************************
+ IMPLICIT DOUBLE PRECISION (A-H,O-Z)
+! DIMENSION ET(0:*), VN(0:*), QN(0:*), DQN(0:*)
double precision ET(0:*), VN(0:*), QN(0:*), DQN(0:*)
- DQN(0) = 0.D0
- IF (N .EQ. 0) RETURN
-
- DO 1 I=0,N
- SU = 0.D0
- VI = VN(I)
- EI = ET(I)
- DO 2 J=0,N
- IF (I .EQ. J) GOTO 2
- VJ = VN(J)
- EJ = ET(J)
- SU = SU+QN(J)/(VJ*(EI-EJ))
-2 CONTINUE
- DQN(I) = VI*SU
-1 CONTINUE
-
- DN = DFLOAT(N)
- C = .25D0*DN*(DN+1.D0)
- DQN(0) = DQN(0)-C*QN(0)
- DQN(N) = DQN(N)+C*QN(N)
-
- RETURN
- END SUBROUTINE DELEGL
-
+ DQN(0) = 0.D0
+ IF (N == 0) RETURN
+ DO 1 I=0,N
+ SU = 0.D0
+ VI = VN(I)
+ EI = ET(I)
+ DO 2 J=0,N
+ IF (I == J) GOTO 2
+ VJ = VN(J)
+ EJ = ET(J)
+ SU = SU+QN(J)/(VJ*(EI-EJ))
+2 CONTINUE
+ DQN(I) = VI*SU
+1 CONTINUE
+
+ DN = DFLOAT(N)
+ C = .25D0*DN*(DN+1.D0)
+ DQN(0) = DQN(0)-C*QN(0)
+ DQN(N) = DQN(N)+C*QN(N)
+
+ RETURN
+ END SUBROUTINE DELEGL
+
+
More information about the CIG-COMMITS
mailing list