[cig-commits] r21899 - in seismo/3D/SPECFEM3D_GLOBE/trunk: src/auxiliaries src/meshfem3D src/shared src/specfem3D utils/Profiles utils/attenuation utils/extract_database utils/funaro_routines_bernhard
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Thu Apr 18 12:23:18 PDT 2013
Author: dkomati1
Date: 2013-04-18 12:23:17 -0700 (Thu, 18 Apr 2013)
New Revision: 21899
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_paraview_strain_data.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_AVS_DX.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_GMT_global.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/add_topography.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/add_topography_410_650.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/add_topography_cmb.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/add_topography_icb.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/calc_jacobian.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_central_cube.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regions_mesh.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_ellipticity.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_jacobian_boundaries.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/meshfem3D_models.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_atten3D_QRFSI12.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_attenuation.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_crustmaps.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_epcrust.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_eucrust.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_gapp2.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_heterogen_mantle.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_jp1d.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_jp3d.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_s362ani.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_sea99_s.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/moho_stretching.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/write_AVS_DX_global_data.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/auto_ner.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/calendar.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/read_compute_parameters.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/read_value_parameters.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_Dev.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/convert_time.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_regular_points.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_sources.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/prepare_timerun.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_output_SAC.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_seismograms.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/utils/Profiles/write_profile.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/utils/attenuation/attenuation.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/utils/attenuation/attenuation_output.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/utils/attenuation/attenuation_simplex.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/utils/extract_database/extract_database.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/utils/funaro_routines_bernhard/routines_SEM_Bernhard_diff_lagrange.f90
Log:
switched to unified syntax for endif, enddo, else if;
also removed extra white spaces
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_paraview_strain_data.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_paraview_strain_data.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_paraview_strain_data.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -157,7 +157,7 @@
read(27) datstore(1:npoint(iproc))
endif
close(27)
- elseif(comp == 'SI1' .or. comp == 'SI2') then
+ else if(comp == 'SI1' .or. comp == 'SI2') then
write(local_data_file,'(a,i6.6,a)') 'movie3D_SEE',it,'.bin'
!print *, iproc,'reading from file:'//trim(prname)//trim(local_data_file)
!print *, 'reading from file:',local_data_file
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/combine_vol_data.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -368,7 +368,7 @@
if (HIGH_RESOLUTION_MESH ) then
if( ir==3 ) then
npoint(it) = numpoin
- elseif( numpoin /= npoint(it)) then
+ else if( numpoin /= npoint(it)) then
print*,'region:',ir
print*,'error number of points:',numpoin,npoint(it)
stop 'different number of points (high-res)'
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_AVS_DX.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_AVS_DX.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_AVS_DX.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -395,7 +395,7 @@
displn(i,j) = displx*normal_x + disply*normal_y + displz*normal_z
- elseif(USE_COMPONENT == 2) then
+ else if(USE_COMPONENT == 2) then
! compute unit tangent vector to the surface (N-S)
RRval = sqrt(xcoord**2 + ycoord**2 + zcoord**2)
@@ -413,7 +413,7 @@
displn(i,j) = - (displx*thetahat_x + disply*thetahat_y + displz*thetahat_z)
- elseif(USE_COMPONENT == 3) then
+ else if(USE_COMPONENT == 3) then
! compute unit tangent to the surface (E-W)
rhoval = sqrt(xcoord**2 + ycoord**2)
@@ -445,12 +445,12 @@
yp(ieoff+ilocnum) = dble(y(i,j))
zp(ieoff+ilocnum) = dble(z(i,j))
field_display(ieoff+ilocnum) = dble(displn(i,j))
- elseif(ilocnum == 2) then
+ else if(ilocnum == 2) then
xp(ieoff+ilocnum) = dble(x(i+1,j))
yp(ieoff+ilocnum) = dble(y(i+1,j))
zp(ieoff+ilocnum) = dble(z(i+1,j))
field_display(ieoff+ilocnum) = dble(displn(i+1,j))
- elseif(ilocnum == 3) then
+ else if(ilocnum == 3) then
xp(ieoff+ilocnum) = dble(x(i+1,j+1))
yp(ieoff+ilocnum) = dble(y(i+1,j+1))
zp(ieoff+ilocnum) = dble(z(i+1,j+1))
@@ -554,9 +554,9 @@
if(USE_OPENDX) then
if( USE_COMPONENT == 1) then
write(outputname,"('/DX_movie_',i6.6,'.Z.dx')") it
- elseif( USE_COMPONENT == 2) then
+ else if( USE_COMPONENT == 2) then
write(outputname,"('/DX_movie_',i6.6,'.N.dx')") it
- elseif( USE_COMPONENT == 3) then
+ else if( USE_COMPONENT == 3) then
write(outputname,"('/DX_movie_',i6.6,'.E.dx')") it
endif
open(unit=11,file=trim(OUTPUT_FILES)//trim(outputname),status='unknown')
@@ -565,9 +565,9 @@
if(UNIQUE_FILE .and. iframe == 1) then
if( USE_COMPONENT == 1) then
outputname = '/AVS_movie_all.Z.inp'
- elseif( USE_COMPONENT == 2) then
+ else if( USE_COMPONENT == 2) then
outputname = '/AVS_movie_all.N.inp'
- elseif( USE_COMPONENT == 3) then
+ else if( USE_COMPONENT == 3) then
outputname = '/AVS_movie_all.E.inp'
endif
open(unit=11,file=trim(OUTPUT_FILES)//trim(outputname),status='unknown')
@@ -578,9 +578,9 @@
else if(.not. UNIQUE_FILE) then
if( USE_COMPONENT == 1) then
write(outputname,"('/AVS_movie_',i6.6,'.Z.inp')") it
- elseif( USE_COMPONENT == 2) then
+ else if( USE_COMPONENT == 2) then
write(outputname,"('/AVS_movie_',i6.6,'.N.inp')") it
- elseif( USE_COMPONENT == 3) then
+ else if( USE_COMPONENT == 3) then
write(outputname,"('/AVS_movie_',i6.6,'.E.inp')") it
endif
open(unit=11,file=trim(OUTPUT_FILES)//trim(outputname),status='unknown')
@@ -589,9 +589,9 @@
else if(USE_GMT) then
if( USE_COMPONENT == 1) then
write(outputname,"('/gmt_movie_',i6.6,'.Z.xyz')") it
- elseif( USE_COMPONENT == 2) then
+ else if( USE_COMPONENT == 2) then
write(outputname,"('/gmt_movie_',i6.6,'.N.xyz')") it
- elseif( USE_COMPONENT == 3) then
+ else if( USE_COMPONENT == 3) then
write(outputname,"('/gmt_movie_',i6.6,'.E.xyz')") it
endif
open(unit=11,file=trim(OUTPUT_FILES)//trim(outputname),status='unknown')
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 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/auxiliaries/create_movie_GMT_global.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -536,7 +536,7 @@
displn(i,j) = displx*normal_x + disply*normal_y + displz*normal_z
- elseif(USE_COMPONENT == 2) then
+ else if(USE_COMPONENT == 2) then
! compute unit tangent vector to the surface (N-S)
RRval = sqrt(xcoord**2 + ycoord**2 + zcoord**2)
@@ -554,7 +554,7 @@
displn(i,j) = - (displx*thetahat_x + disply*thetahat_y + displz*thetahat_z)
- elseif(USE_COMPONENT == 3) then
+ else if(USE_COMPONENT == 3) then
! compute unit tangent to the surface (E-W)
rhoval = sqrt(xcoord**2 + ycoord**2)
@@ -850,9 +850,9 @@
if(OUTPUT_BINARY) then
if(USE_COMPONENT == 1) then
write(outputname,"('bin_movie_',i6.6,'.d')") it
- elseif(USE_COMPONENT == 2) then
+ else if(USE_COMPONENT == 2) then
write(outputname,"('bin_movie_',i6.6,'.N')") it
- elseif(USE_COMPONENT == 3) then
+ else if(USE_COMPONENT == 3) then
write(outputname,"('bin_movie_',i6.6,'.E')") it
endif
open(unit=11,file='OUTPUT_FILES/'//trim(outputname),status='unknown', &
@@ -862,9 +862,9 @@
else
if(USE_COMPONENT == 1) then
write(outputname,"('ascii_movie_',i6.6,'.d')") it
- elseif(USE_COMPONENT == 2) then
+ else if(USE_COMPONENT == 2) then
write(outputname,"('ascii_movie_',i6.6,'.N')") it
- elseif(USE_COMPONENT == 3) then
+ else if(USE_COMPONENT == 3) then
write(outputname,"('ascii_movie_',i6.6,'.E')") it
endif
open(unit=11,file='OUTPUT_FILES/'//trim(outputname),status='unknown', &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/add_topography.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/add_topography.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/add_topography.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -154,19 +154,19 @@
! also make sure factor makes sense
if(gamma < -0.02 .or. gamma > 1.02) then
call exit_MPI(myrank,'incorrect value of factor for topography gll points')
- end if
+ endif
!
! since not all GLL points are exactlly at R220, use a small
! tolerance for R220 detection
if (abs(gamma) < SMALLVAL) then
gamma = 0.0
- end if
+ endif
xstore(i,j,k,ispec) = xstore(i,j,k,ispec)*(ONE + gamma * elevation / r)
ystore(i,j,k,ispec) = ystore(i,j,k,ispec)*(ONE + gamma * elevation / r)
zstore(i,j,k,ispec) = zstore(i,j,k,ispec)*(ONE + gamma * elevation / r)
- end do
- end do
- end do
+ enddo
+ enddo
+ enddo
end subroutine add_topography_gll
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/add_topography_410_650.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/add_topography_410_650.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/add_topography_410_650.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -113,13 +113,13 @@
xelm(ia) = xelm(ia)*(ONE + gamma * topo410 / r)
yelm(ia) = yelm(ia)*(ONE + gamma * topo410 / r)
zelm(ia) = zelm(ia)*(ONE + gamma * topo410 / r)
- elseif(r>= R771/R_EARTH .and. r <= R670/R_EARTH) then
+ else if(r>= R771/R_EARTH .and. r <= R670/R_EARTH) then
! stretching between R771 and R670
gamma = (r - R771/R_EARTH) / (R670/R_EARTH - R771/R_EARTH)
xelm(ia) = xelm(ia)*(ONE + gamma * topo650 / r)
yelm(ia) = yelm(ia)*(ONE + gamma * topo650 / r)
zelm(ia) = zelm(ia)*(ONE + gamma * topo650 / r)
- elseif(r > R670/R_EARTH .and. r < R400/R_EARTH) then
+ else if(r > R670/R_EARTH .and. r < R400/R_EARTH) then
! stretching between R670 and R400
gamma = (R400/R_EARTH - r) / (R400/R_EARTH - R670/R_EARTH)
xelm(ia) = xelm(ia)*(ONE + (topo410 + gamma * (topo650 - topo410)) / r)
@@ -227,13 +227,13 @@
xstore(i,j,k,ispec) = xstore(i,j,k,ispec)*(ONE + gamma * topo410 / r)
ystore(i,j,k,ispec) = ystore(i,j,k,ispec)*(ONE + gamma * topo410 / r)
zstore(i,j,k,ispec) = zstore(i,j,k,ispec)*(ONE + gamma * topo410 / r)
- elseif(r>= R771/R_EARTH .and. r <= R670/R_EARTH) then
+ else if(r>= R771/R_EARTH .and. r <= R670/R_EARTH) then
! stretching between R771 and R670
gamma = (r - R771/R_EARTH) / (R670/R_EARTH - R771/R_EARTH)
xstore(i,j,k,ispec) = xstore(i,j,k,ispec)*(ONE + gamma * topo650 / r)
ystore(i,j,k,ispec) = ystore(i,j,k,ispec)*(ONE + gamma * topo650 / r)
zstore(i,j,k,ispec) = zstore(i,j,k,ispec)*(ONE + gamma * topo650 / r)
- elseif(r > R670/R_EARTH .and. r < R400/R_EARTH) then
+ else if(r > R670/R_EARTH .and. r < R400/R_EARTH) then
! stretching between R670 and R400
gamma = (R400/R_EARTH - r) / (R400/R_EARTH - R670/R_EARTH)
xstore(i,j,k,ispec) = xstore(i,j,k,ispec)*(ONE + (topo410 + gamma * (topo650 - topo410)) / r)
@@ -243,7 +243,7 @@
if(gamma < -0.0001 .or. gamma > 1.0001) call exit_MPI(myrank,'incorrect value of gamma for 410-650 topography')
enddo
- end do
- end do
+ enddo
+ enddo
end subroutine add_topography_410_650_gll
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/add_topography_cmb.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/add_topography_cmb.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/add_topography_cmb.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -68,7 +68,7 @@
if(r >= RCMB/R_EARTH .and. r <= RTOPDDOUBLEPRIME/R_EARTH) then
! stretching between RCMB and RTOPDDOUBLEPRIME
gamma = (RTOPDDOUBLEPRIME/R_EARTH - r) / (RTOPDDOUBLEPRIME/R_EARTH - RCMB/R_EARTH)
- elseif(r>= r_start .and. r <= RCMB/R_EARTH) then
+ else if(r>= r_start .and. r <= RCMB/R_EARTH) then
! stretching between r_start and RCMB
gamma = (r - r_start) / (RCMB/R_EARTH - r_start)
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/add_topography_icb.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/add_topography_icb.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/add_topography_icb.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -65,7 +65,7 @@
if(r > 0.0d0 .and. r <= RICB/R_EARTH) then
! stretching between center and RICB
gamma = r/(RICB/R_EARTH)
- elseif(r>= RICB/R_EARTH .and. r <= RCMB/R_EARTH) then
+ else if(r>= RICB/R_EARTH .and. r <= RCMB/R_EARTH) then
! stretching between RICB and RCMB
gamma = (r - RCMB/R_EARTH) / (RICB/R_EARTH - RCMB/R_EARTH)
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/calc_jacobian.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/calc_jacobian.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/calc_jacobian.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -148,28 +148,28 @@
sumdershapeeta = sumdershapeeta + hlagrange_eta
sumdershapegamma = sumdershapegamma + hlagrange_gamma
- end do
- end do
- end do
+ enddo
+ enddo
+ enddo
! Check the lagrange polynomial and its derivative
if (abs(xmesh - xstore(i,j,k,ispec)) > TINYVAL &
.or. abs(ymesh - ystore(i,j,k,ispec)) > TINYVAL &
.or. abs(zmesh - zstore(i,j,k,ispec)) > TINYVAL ) then
call exit_MPI(myrank,'new mesh are wrong in recalc_jacobian_gall3D.f90')
- end if
+ endif
if(abs(sumshape-one) > TINYVAL) then
call exit_MPI(myrank,'error shape functions in recalc_jacobian_gll3D.f90')
- end if
+ endif
if(abs(sumdershapexi) > TINYVAL) then
call exit_MPI(myrank,'error derivative xi in recalc_jacobian_gll3D.f90')
- end if
+ endif
if(abs(sumdershapeeta) > TINYVAL) then
call exit_MPI(myrank,'error derivative eta in recalc_jacobian_gll3D.f90')
- end if
+ endif
if(abs(sumdershapegamma) > TINYVAL) then
call exit_MPI(myrank,'error derivative gamma in recalc_jacobian_gll3D.f90')
- end if
+ endif
jacobian = xxi*(yeta*zgamma-ygamma*zeta) - &
@@ -181,7 +181,7 @@
call xyz_2_rthetaphi_dble(xmesh,ymesh,zmesh,r,theta,phi)
print*,'r/lat/lon:',r*R_EARTH_KM,90.0-theta*180./PI,phi*180./PI
call exit_MPI(myrank,'3D Jacobian undefined in recalc_jacobian_gll3D.f90')
- end if
+ endif
! invert the relation (Fletcher p. 50 vol. 2)
xix = (yeta*zgamma-ygamma*zeta) / jacobian
@@ -219,7 +219,7 @@
gammaystore(i,j,k,ispec) = gammay
gammazstore(i,j,k,ispec) = gammaz
endif
- end if
+ endif
enddo
enddo
enddo
@@ -308,8 +308,8 @@
sumshape = sumshape + hlagrange
sumdershapexi = sumdershapexi + hlagrange_xi
sumdershapeeta = sumdershapeeta + hlagrange_eta
- end do
- end do
+ enddo
+ enddo
! Check the lagrange polynomial
@@ -317,17 +317,17 @@
.or. abs(ymesh - yelm2D(i,j)) > TINYVAL &
.or. abs(zmesh - zelm2D(i,j)) > TINYVAL ) then
call exit_MPI(myrank,'new boundary mesh is wrong in recalc_jacobian_gll2D')
- end if
+ endif
if (abs(sumshape-one) > TINYVAL) then
call exit_MPI(myrank,'error shape functions in recalc_jacobian_gll2D')
- end if
+ endif
if (abs(sumdershapexi) > TINYVAL) then
call exit_MPI(myrank,'error derivative xi in recalc_jacobian_gll2D')
- end if
+ endif
if (abs(sumdershapeeta) > TINYVAL) then
call exit_MPI(myrank,'error derivative eta in recalc_jacobian_gll2D')
- end if
+ endif
unx = yxi*zeta - yeta*zxi
uny = zxi*xeta - zeta*xxi
@@ -346,8 +346,8 @@
normal(2,i,j,ispecb)=uny/jacobian
normal(3,i,j,ispecb)=unz/jacobian
endif
- end do
- end do
+ enddo
+ enddo
end subroutine recalc_jacobian_gll2D
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_central_cube.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_central_cube.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_central_cube.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -235,7 +235,7 @@
if (mod(NPROC_XI,2)/=0) then
if (ichunk == CHUNK_AB) then
nz_inf_limit = ((nz_central_cube*2)/NPROC_XI)*floor(NPROC_XI/2.d0)
- elseif (ichunk == CHUNK_AB_ANTIPODE) then
+ else if (ichunk == CHUNK_AB_ANTIPODE) then
nz_inf_limit = ((nz_central_cube*2)/NPROC_XI)*ceiling(NPROC_XI/2.d0)
endif
else
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regions_mesh.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regions_mesh.F90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/create_regions_mesh.F90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -280,7 +280,7 @@
nspec_att = nspec
else
nspec_att = 1
- end if
+ endif
allocate(Qmu_store(NGLLX,NGLLY,NGLLZ,nspec_att), &
tau_e_store(N_SLS,NGLLX,NGLLY,NGLLZ,nspec_att),stat=ier)
if(ier /= 0) stop 'error in allocate 1'
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_ellipticity.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_ellipticity.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_ellipticity.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -105,8 +105,8 @@
ystore(i,j,k,ispec)=ystore(i,j,k,ispec)*factor
zstore(i,j,k,ispec)=zstore(i,j,k,ispec)*factor
- end do
- end do
- end do
+ enddo
+ enddo
+ enddo
end subroutine get_ellipticity_gll
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_jacobian_boundaries.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_jacobian_boundaries.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/get_jacobian_boundaries.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -151,13 +151,13 @@
xelm2D(j,k) = xstore(1,j,k,ispec)
yelm2D(j,k) = ystore(1,j,k,ispec)
zelm2D(j,k) = zstore(1,j,k,ispec)
- end do
- end do
+ enddo
+ enddo
! recalculate jacobian according to 2D GLL points
call recalc_jacobian_gll2D(myrank,ispecb1,xelm2D,yelm2D,zelm2D, &
yigll,zigll,jacobian2D_xmin,normal_xmin,&
NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
- end if
+ endif
endif
! on boundary: xmax
@@ -207,13 +207,13 @@
xelm2D(j,k) = xstore(NGLLX,j,k,ispec)
yelm2D(j,k) = ystore(NGLLX,j,k,ispec)
zelm2D(j,k) = zstore(NGLLX,j,k,ispec)
- end do
- end do
+ enddo
+ enddo
! recalculate jacobian according to 2D GLL points
call recalc_jacobian_gll2D(myrank,ispecb2,xelm2D,yelm2D,zelm2D,&
yigll,zigll,jacobian2D_xmax,normal_xmax,&
NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
- end if
+ endif
endif
! on boundary: ymin
@@ -263,13 +263,13 @@
xelm2D(i,k) = xstore(i,1,k,ispec)
yelm2D(i,k) = ystore(i,1,k,ispec)
zelm2D(i,k) = zstore(i,1,k,ispec)
- end do
- end do
+ enddo
+ enddo
! recalcualte 2D jacobian according to GLL points
call recalc_jacobian_gll2D(myrank,ispecb3,xelm2D,yelm2D,zelm2D,&
xigll,zigll,jacobian2D_ymin,normal_ymin,&
NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
- end if
+ endif
endif
! on boundary: ymax
@@ -319,13 +319,13 @@
xelm2D(i,k) = xstore(i,NGLLY,k,ispec)
yelm2D(i,k) = ystore(i,NGLLY,k,ispec)
zelm2D(i,k) = zstore(i,NGLLY,k,ispec)
- end do
- end do
+ enddo
+ enddo
! recalculate jacobian for 2D GLL points
call recalc_jacobian_gll2D(myrank,ispecb4,xelm2D,yelm2D,zelm2D,&
xigll,zigll,jacobian2D_ymax,normal_ymax,&
NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
- end if
+ endif
endif
! on boundary: bottom
@@ -374,13 +374,13 @@
xelm2D(i,j) = xstore(i,j,1,ispec)
yelm2D(i,j) = ystore(i,j,1,ispec)
zelm2D(i,j) = zstore(i,j,1,ispec)
- end do
- end do
+ enddo
+ enddo
! recalcuate 2D jacobian according to GLL points
call recalc_jacobian_gll2D(myrank,ispecb5,xelm2D,yelm2D,zelm2D,&
xigll,yigll,jacobian2D_bottom,normal_bottom,&
NGLLX,NGLLY,NSPEC2D_BOTTOM)
- end if
+ endif
endif
@@ -429,14 +429,14 @@
xelm2D(i,j) = xstore(i,j,NGLLZ,ispec)
yelm2D(i,j) = ystore(i,j,NGLLZ,ispec)
zelm2D(i,j) = zstore(i,j,NGLLZ,ispec)
- end do
- end do
+ enddo
+ enddo
! recalcuate jacobian according to 2D gll points
call recalc_jacobian_gll2D(myrank,ispecb6,xelm2D,yelm2D,zelm2D,&
xigll,yigll,jacobian2D_top,normal_top,&
NGLLX,NGLLY,NSPEC2D_TOP)
- end if
+ endif
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/meshfem3D_models.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/meshfem3D_models.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/meshfem3D_models.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -1212,7 +1212,7 @@
! moho=0.0d0
! found_crust = .false.
call model_epcrust(lat,lon,r,vpc,vsc,rhoc,moho,found_crust,EPCRUST,elem_in_crust)
-! end if
+! endif
case default
stop 'crustal model type not defined'
@@ -1331,7 +1331,7 @@
end select
- end if
+ endif
! Get tau_e from tau_s and Qmu
call model_attenuation_getstored_tau(Qmu, T_c_source, tau_s, tau_e, AM_V, AM_S, AS_V)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_atten3D_QRFSI12.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_atten3D_QRFSI12.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_atten3D_QRFSI12.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -119,7 +119,7 @@
j=0
do l=0,MAXL_Q
do m=0,l
- if(m.eq.0)then
+ if(m==0)then
j=j+1
read(10,*)ll,mm,v1
QRFSI12_Q%dqmu(k,j)=v1
@@ -684,15 +684,15 @@
!!$ XCOSEC(I)=0.
!!$ XP(I)=0.
!!$ enddo
-!!$ IF(L.GT.1.AND.ABS(THETA).GT.1.E-5) GO TO 3
+!!$ IF(L>1.AND.ABS(THETA)>1.E-5) GO TO 3
!!$ X(1)=FCT
-!!$ IF(L.EQ.0) RETURN
+!!$ IF(L==0) RETURN
!!$ X(1)=CT*FCT
!!$ X(2)=-ST*FCT/DSFL3
!!$ XP(1)=-ST*FCT
!!$ XP(2)=-.5D0*CT*FCT*DSFL3
-!!$ IF(ABS(THETA).LT.1.E-5) XCOSEC(2)=XP(2)
-!!$ IF(ABS(THETA).GE.1.E-5) XCOSEC(2)=X(2)/ST
+!!$ IF(ABS(THETA)<1.E-5) XCOSEC(2)=XP(2)
+!!$ IF(ABS(THETA)>=1.E-5) XCOSEC(2)=X(2)/ST
!!$ RETURN
!!$ 3 X1=1.D0
!!$ X2=CT
@@ -714,7 +714,7 @@
!!$ XCOSEC(2)=X(2)*COSEC
!!$ XP(2)=-XP(2)/SFL3
!!$ SUM=SUM+2.D0*X(2)*X(2)
-!!$ IF(SUM-COMPAR.GT.SMALL) RETURN
+!!$ IF(SUM-COMPAR>SMALL) RETURN
!!$ X1=X3
!!$ X2=-X2/DSQRT(dble(FLOAT(L*(L+1))))
!!$ DO I=3,MP1
@@ -724,7 +724,7 @@
!!$ XM=K
!!$ X3=-(2.D0*COT*(XM-1.D0)*X2+F2*X1)/F1
!!$ SUM=SUM+2.D0*X3*X3
-!!$ IF(SUM-COMPAR.GT.SMALL.AND.I.NE.LP1) RETURN
+!!$ IF(SUM-COMPAR>SMALL.AND.I/=LP1) RETURN
!!$ X(I)=X3
!!$ XCOSEC(I)=X(I)*COSEC
!!$ X1=X2
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_attenuation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_attenuation.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_attenuation.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -353,12 +353,12 @@
else if(REFERENCE_1D_MODEL == REFERENCE_MODEL_SEA1D) then
AM_V%Qr(:) = SEA1DM_V%radius_sea1d(:)
AM_V%Qmu(:) = SEA1DM_V%Qmu_sea1d(:)
- end if
+ endif
do i = 1, AM_V%Qn
call model_attenuation_getstored_tau(AM_V%Qmu(i), AM_V%QT_c_source, AM_V%Qtau_s, tau_e, AM_V, AM_S,AS_V)
AM_V%Qtau_e(:,i) = tau_e(:)
- end do
+ enddo
end subroutine model_attenuation_setup
@@ -816,7 +816,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
@@ -1138,7 +1138,7 @@
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
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_crustmaps.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_crustmaps.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_crustmaps.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -626,7 +626,7 @@
velsl(i)=weightul*GC_V%velocsnp(i)+weightur*GC_V%velocsnp(i)+&
weightll*GC_V%velocs(1,ileftlng,i)+weightlr*GC_V%velocs(1,irightlng,i)
enddo
- elseif(iupcolat==180*CRUSTMAP_RESOLUTION) then
+ else if(iupcolat==180*CRUSTMAP_RESOLUTION) then
! south pole
do i=1,NLAYERS_CRUSTMAP
thickl(i)=weightul*GC_V%thickness(iupcolat,ileftlng,i)+weightur*GC_V%thickness(iupcolat,irightlng,i)+&
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_epcrust.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_epcrust.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_epcrust.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -111,8 +111,8 @@
EPCRUST%vp_ep(ilon,jlat,1:3) = tmp(7:9)
EPCRUST%vs_ep(ilon,jlat,1:3) = tmp(10:12)
EPCRUST%rho_ep(ilon,jlat,1:3) = tmp(13:15)
- end do
- end do
+ enddo
+ enddo
close(1001)
end subroutine read_epcrust_model
@@ -150,7 +150,7 @@
!if ( lat < EPCRUST_LAT_MIN .or. lat > EPCRUST_LAT_MAX &
! .or. lon < EPCRUST_LON_MIN .or. lon > EPCRUST_LON_MAX ) then
! stop 'incorrect enter EPCRUST model, check lat and lon'
- !end if
+ !endif
vp=0.0d0
vs=0.0d0
@@ -180,8 +180,8 @@
vpsmooth(:)=vpsmooth(:)+weightl*EPCRUST%vp_ep(ilon,jlat,:)
vssmooth(:)=vssmooth(:)+weightl*EPCRUST%vs_ep(ilon,jlat,:)
rhosmooth(:)=rhosmooth(:)+weightl*EPCRUST%rho_ep(ilon,jlat,:)
- end do
- end if
+ enddo
+ endif
!topo=(R_EARTH_KM+z0)/R_EARTH_KM
!basement=(R_EARTH_KM+z0-zsmooth(1))/R_EARTH_KM
@@ -211,14 +211,14 @@
rho=rhosmooth(3)
else
found_crust=.false.
- end if
+ endif
if (found_crust ) then
scaleval=dsqrt(PI*GRAV*RHOAV)
vp=vp*1000.d0/(R_EARTH*scaleval)
vs=vs*1000.d0/(R_EARTH*scaleval)
rho=rho*1000.d0/RHOAV
- end if
+ endif
!moho = -(z0-zsmooth(1)-zsmooth(2)-zsmooth(3))/R_EARTH_KM
moho = (zsmooth(1)+zsmooth(2)+zsmooth(3))/R_EARTH_KM
@@ -227,7 +227,7 @@
cut=7.0/R_EARTH_KM
if ( moho < cut ) then
moho = cut
- end if
+ endif
end subroutine model_epcrust
@@ -263,7 +263,7 @@
print*, 'error cap:', cap_degree_EP
print*, 'lat/lon:', x,y
stop 'error cap_degree too small'
- end if
+ endif
CAP=cap_degree_EP*PI/180.0d0
dtheta=0.5d0*CAP/dble(NTHETA_EP)
@@ -313,20 +313,20 @@
xx(j)=0.0d0
do k = 1,3
xx(j)=xx(j)+rotation_matrix(j,k)*xc(k)
- end do
- end do
+ enddo
+ enddo
call xyz_2_rthetaphi_dble(xx(1),xx(2),xx(3),r_rot,theta_rot,phi_rot)
call reduce(theta_rot,phi_rot)
x1(i)=phi_rot*RADIANS_TO_DEGREES
y1(i)=(PI_OVER_TWO-theta_rot)*RADIANS_TO_DEGREES
if (x1(i) > 180.d0) x1(i)=x1(i)-360.d0
- end do
- end do
+ enddo
+ enddo
if (abs(total-1.0d0) > 0.001d0) then
print*,'error cap:',total,cap_degree_EP
stop
- end if
+ endif
end subroutine epcrust_smooth_base
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_eucrust.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_eucrust.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_eucrust.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -131,7 +131,7 @@
do i=1,36058
read(11,'(a80)',iostat=ierror) line
- if(ierror .ne. 0 ) stop
+ if(ierror /= 0 ) stop
read(line,*)lon,lat,vp_uppercrust,vp_lowercrust,vp_avg,topo,basement,upper_lower_depth,moho_depth
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_gapp2.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_gapp2.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_gapp2.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -201,21 +201,21 @@
d=R_EARTH_-radius*R_EARTH_
call d2id(d,nnr,dep,id,icon)
- if(icon.ne.0) then
+ if(icon/=0) then
write(6,*)icon
write(6,*) radius,theta,phi,dvp,dvs,drho
endif
! latitude
- if(theta.ge.PI) then
+ if(theta>=PI) then
ia = na
else
ia = theta / dtheta + 1
endif
! longitude
- if(phi .lt. 0.0d0) phi = phi + 2.*PI
+ if(phi < 0.0d0) phi = phi + 2.*PI
io=phi / dphi + 1
- if(io.gt.no) io=io-no
+ if(io>no) io=io-no
! velocity and density perturbations
dvp = vp3(ia,io,id)/100.d0
@@ -246,20 +246,20 @@
icon=0
dmax=di(mr)
dmin=di(0)
- if(d.gt.dmax) then
+ if(d>dmax) then
icon=99
- else if(d.lt.dmin) then
+ else if(d<dmin) then
icon=-99
- else if(d.eq.dmax) then
+ else if(d==dmax) then
id=mr+1
else
do i = 0, mr
- if(d.lt.di(i)) then
+ if(d<di(i)) then
id=i
goto 900
endif
enddo
- end if
+ endif
900 continue
!..................................................................
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_heterogen_mantle.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_heterogen_mantle.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_heterogen_mantle.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -105,7 +105,7 @@
do i = 1,j
read(10,rec=i,fmt='(F20.15)') HMM%rho_in(i)
- end do
+ enddo
close(10)
@@ -215,6 +215,6 @@
drho = 0.
dvp = 0.
dvs = 0.
- end if
+ endif
end subroutine model_heterogen_mantle
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_jp1d.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_jp1d.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_jp1d.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -195,7 +195,7 @@
vs = 3.5d0
Qmu=600.0d0
Qkappa=57827.0d0
- end if
+ endif
! non-dimensionalize
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_jp3d.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_jp3d.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_jp3d.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -393,17 +393,17 @@
! when LAY = 2, the focus is in the lower crust;
! when LAY = 3, the focus is in the mantle wedge;
! when LAY = 4, the focus is beneath the plate boundary.
- IF(HE.LE.H1) THEN
+ IF(HE<=H1) THEN
LAY = 1
found_crust = .true.
- ELSE IF(HE.GT.H1.AND.HE.LE.H2) THEN
+ ELSE IF(HE>H1.AND.HE<=H2) THEN
LAY = 2
found_crust = .true.
- ELSE IF(HE.GT.H2.AND.HE.LE.H3) THEN
+ ELSE IF(HE>H2.AND.HE<=H3) THEN
LAY = 3
ELSE
LAY = 4
- END IF
+ endif
CALL VEL1D(HE,vp,LAY,1,JP3DM_V)
CALL VEL1D(HE,vs,LAY,2,JP3DM_V)
@@ -415,13 +415,13 @@
vs = vs*(1.0d0+dvs)
! determine rho
- if(LAY .eq. 1) then
+ if(LAY == 1) then
rho=2.6
endif
- if(LAY .eq. 2) then
+ if(LAY == 2) then
rho=2.9
endif
- if(LAY .GT. 2) then
+ if(LAY > 2) then
rho=3.3+(vs-4.4)*0.66667
endif
! non-dimensionalize
@@ -735,7 +735,7 @@
DO 10 I = 1,IPMAX
IP1 = IP+1
PNOW = (FLOAT(I)-PLX)/100.0
- IF(PNOW.GE.PNX(IP1)) IP = IP1
+ IF(PNOW>=PNX(IP1)) IP = IP1
IPLOCX(I)= IP
10 CONTINUE
RLX = 1.0-RNX(1)*100.0
@@ -744,7 +744,7 @@
DO 20 I = 1,IRMAX
IR1 = IR+1
RNOW = (FLOAT(I)-RLX)/100.0
- IF(RNOW.GE.RNX(IR1)) IR = IR1
+ IF(RNOW>=RNX(IR1)) IR = IR1
IRLOCX(I)= IR
20 CONTINUE
HLX = 1.0-HNX(1)
@@ -753,7 +753,7 @@
DO 30 I = 1,IHMAX
IH1 = IH+1
HNOW = FLOAT(I)-HLX
- IF(HNOW.GE.HNX(IH1)) IH = IH1
+ IF(HNOW>=HNX(IH1)) IH = IH1
IHLOCX(I)= IH
30 CONTINUE
RETURN
@@ -845,14 +845,14 @@
JP3DM_V%P = 90.0-PE/DEGREES_TO_RADIANS
JP3DM_V%R = RE/DEGREES_TO_RADIANS
JP3DM_V%H = HE
- IF(LAY.LE.3) THEN
+ IF(LAY<=3) THEN
CALL PRHF(JP3DM_V%IPLOCA,JP3DM_V%IRLOCA,JP3DM_V%IHLOCA,JP3DM_V%PLA,JP3DM_V%RLA,JP3DM_V%HLA, &
JP3DM_V%PNA,JP3DM_V%RNA,JP3DM_V%HNA,MPA,MRA,MHA,MKA,JP3DM_V)
- ELSE IF(LAY.EQ.4) THEN
+ ELSE IF(LAY==4) THEN
CALL PRHF(JP3DM_V%IPLOCB,JP3DM_V%IRLOCB,JP3DM_V%IHLOCB,JP3DM_V%PLB,JP3DM_V%RLB,JP3DM_V%HLB, &
JP3DM_V%PNB,JP3DM_V%RNB,JP3DM_V%HNB,MPB,MRB,MHB,MKB,JP3DM_V)
ELSE
- END IF
+ endif
JP3DM_V%WV(1) = JP3DM_V%PF1*JP3DM_V%RF1*JP3DM_V%HF1
JP3DM_V%WV(2) = JP3DM_V%PF*JP3DM_V%RF1*JP3DM_V%HF1
JP3DM_V%WV(3) = JP3DM_V%PF1*JP3DM_V%RF*JP3DM_V%HF1
@@ -862,12 +862,12 @@
JP3DM_V%WV(7) = JP3DM_V%PF1*JP3DM_V%RF*JP3DM_V%HF
JP3DM_V%WV(8) = JP3DM_V%PF*JP3DM_V%RF*JP3DM_V%HF
! calculate velocity
- IF(LAY.LE.3) THEN
+ IF(LAY<=3) THEN
CALL VABPS(MPA,MRA,MHA,JP3DM_V%VELAP,V,JP3DM_V)
- ELSE IF(LAY.EQ.4) THEN
+ ELSE IF(LAY==4) THEN
CALL VABPS(MPB,MRB,MHB,JP3DM_V%VELBP,V,JP3DM_V)
ELSE
- END IF
+ endif
RETURN
END SUBROUTINE VEL3
@@ -1162,12 +1162,12 @@
CALL LIMIT(JP3DM_V%RRN(1),JP3DM_V%RRN(63),R)
DO 1 I = 1,50
I1 = I+1
- IF(P.GE.JP3DM_V%PN(I).AND.P.LT.JP3DM_V%PN(I1)) GO TO 11
+ IF(P>=JP3DM_V%PN(I).AND.P<JP3DM_V%PN(I1)) GO TO 11
1 CONTINUE
11 CONTINUE
DO 2 J = 1,62
J1 = J+1
- IF(R.GE.JP3DM_V%RRN(J).AND.R.LT.JP3DM_V%RRN(J1)) GO TO 22
+ IF(R>=JP3DM_V%RRN(J).AND.R<JP3DM_V%RRN(J1)) GO TO 22
2 CONTINUE
22 CONTINUE
PF = (P-JP3DM_V%PN(I))/(JP3DM_V%PN(I1)-JP3DM_V%PN(I))
@@ -1178,17 +1178,17 @@
WV2 = PF*RF1
WV3 = PF1*RF
WV4 = PF*RF
- IF(IJK.EQ.1) THEN
+ IF(IJK==1) THEN
HE = WV1*JP3DM_V%DEPA(I,J) + WV2*JP3DM_V%DEPA(I1,J) &
+ WV3*JP3DM_V%DEPA(I,J1) + WV4*JP3DM_V%DEPA(I1,J1)
- ELSE IF(IJK.EQ.2) THEN
+ ELSE IF(IJK==2) THEN
HE = WV1*JP3DM_V%DEPB(I,J) + WV2*JP3DM_V%DEPB(I1,J) &
+ WV3*JP3DM_V%DEPB(I,J1) + WV4*JP3DM_V%DEPB(I1,J1)
- ELSE IF(IJK.EQ.3) THEN
+ ELSE IF(IJK==3) THEN
HE = WV1*JP3DM_V%DEPC(I,J) + WV2*JP3DM_V%DEPC(I1,J) &
+ WV3*JP3DM_V%DEPC(I,J1) + WV4*JP3DM_V%DEPC(I1,J1)
ELSE
- END IF
+ endif
RETURN
END SUBROUTINE HLAY
@@ -1196,8 +1196,8 @@
double precision :: A1,A2,C1,C2,C
A1 = dmin1(C1,C2)
A2 = dmax1(C1,C2)
- IF(C.LT.A1) C = A1
- IF(C.GT.A2) C = A2
+ IF(C<A1) C = A1
+ IF(C>A2) C = A2
END SUBROUTINE LIMIT
!
@@ -1280,22 +1280,22 @@
integer :: IPS,LAY
double precision :: HE,V,VM,HM
- IF(LAY.EQ.1) THEN
+ IF(LAY==1) THEN
V = 6.0
- IF(IPS.EQ.2) V = 3.5
- ELSE IF(LAY.EQ.2) THEN
+ IF(IPS==2) V = 3.5
+ ELSE IF(LAY==2) THEN
V = 6.7
- IF(IPS.EQ.2) V = 3.8
- ELSE IF(LAY.GE.3) THEN
+ IF(IPS==2) V = 3.8
+ ELSE IF(LAY>=3) THEN
HM = 40.0
- IF(HE.LT.HM) THEN
+ IF(HE<HM) THEN
CALL JPMODEL(IPS,HM,VM,JP3DM_V)
V = VM-(HM-HE)*0.003
ELSE
CALL JPMODEL(IPS,HE,V,JP3DM_V)
- END IF
+ endif
ELSE
- END IF
+ endif
RETURN
END
@@ -1479,15 +1479,15 @@
K1 = K+1
H1 = JP3DM_V%DEPJ(K)
H2 = JP3DM_V%DEPJ(K1)
- IF(H.GE.H1.AND.H.LT.H2) GO TO 3
+ IF(H>=H1.AND.H<H2) GO TO 3
2 CONTINUE
3 CONTINUE
H12 = (H-H1)/(H2-H1)
- IF(IPS.EQ.1) THEN
+ IF(IPS==1) THEN
V = (JP3DM_V%VP(K1)-JP3DM_V%VP(K))*H12+JP3DM_V%VP(K)
ELSE
V = (JP3DM_V%VS(K1)-JP3DM_V%VS(K))*H12+JP3DM_V%VS(K)
- END IF
+ endif
RETURN
END
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_s362ani.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_s362ani.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_s362ani.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -185,11 +185,11 @@
lu=1 ! --- log unit: input 3-D model
if(THREE_D_MODEL == THREE_D_MODEL_S362ANI) then
modeldef='DATA/s362ani/S362ANI'
- elseif(THREE_D_MODEL == THREE_D_MODEL_S362WMANI) then
+ else if(THREE_D_MODEL == THREE_D_MODEL_S362WMANI) then
modeldef='DATA/s362ani/S362WMANI'
- elseif(THREE_D_MODEL == THREE_D_MODEL_S362ANI_PREM) then
+ else if(THREE_D_MODEL == THREE_D_MODEL_S362ANI_PREM) then
modeldef='DATA/s362ani/S362ANI_PREM'
- elseif(THREE_D_MODEL == THREE_D_MODEL_S29EA) then
+ else if(THREE_D_MODEL == THREE_D_MODEL_S29EA) then
modeldef='DATA/s362ani/S2.9EA'
else
stop 'unknown 3D model in read_model_s362ani'
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_sea99_s.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_sea99_s.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/model_sea99_s.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -196,15 +196,15 @@
!----------------------- depth in the model ------------------
dep=R_EARTH_KM*(R_UNIT_SPHERE - radius)
- if (dep .le. SEA99M_V%sea99_depth(1)) then
+ if (dep <= SEA99M_V%sea99_depth(1)) then
id1 = 1
xd1 = 0
- else if (dep .ge. SEA99M_V%sea99_depth(SEA99M_V%sea99_ndep)) then
+ else if (dep >= SEA99M_V%sea99_depth(SEA99M_V%sea99_ndep)) then
id1 = SEA99M_V%sea99_ndep
xd1 = 0
else
do i = 2, SEA99M_V%sea99_ndep
- if (dep .le. SEA99M_V%sea99_depth(i)) then
+ if (dep <= SEA99M_V%sea99_depth(i)) then
id1 = i-1
xd1 = (dep-SEA99M_V%sea99_depth(i-1)) / (SEA99M_V%sea99_depth(i) - SEA99M_V%sea99_depth(i-1))
exit
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/moho_stretching.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/moho_stretching.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/moho_stretching.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -120,8 +120,8 @@
! tolerance for R220 detection, fix R220
if (abs(gamma) < SMALLVAL) then
gamma = 0.0d0
- end if
- end if
+ endif
+ endif
if(gamma < -0.0001d0 .or. gamma > 1.0001d0) &
call exit_MPI(myrank,'incorrect value of gamma for moho from crust 2.0')
@@ -147,15 +147,15 @@
! tolerance for R220 detection, fix R220
if (abs(gamma) < SMALLVAL) then
gamma = 0.0d0
- end if
- end if
+ endif
+ endif
if(gamma < -0.0001d0 .or. gamma > 1.0001d0) &
call exit_MPI(myrank,'incorrect value of gamma for moho from crust 2.0')
call move_point(ia,xelm,yelm,zelm,x,y,z,gamma,elevation,r)
- end if
+ endif
endif ! TOPOGRAPHY
@@ -170,16 +170,16 @@
count_mantle = count_mantle + 1
endif
- end do
+ enddo
! sets flag when all corners are above moho
if( count_crust == NGNOD) then
elem_in_crust = .true.
- end if
+ endif
! sets flag when all corners are below moho
if( count_mantle == NGNOD) then
elem_in_mantle = .true.
- end if
+ endif
! small stretch check: stretching should affect only points above R220
if( r*R_EARTH < R220 ) then
@@ -283,16 +283,16 @@
count_mantle = count_mantle + 1
endif
- end do
+ enddo
! sets flag when all corners are above moho
if( count_crust == NGNOD) then
elem_in_crust = .true.
- end if
+ endif
! sets flag when all corners are below moho
if( count_mantle == NGNOD) then
elem_in_mantle = .true.
- end if
+ endif
! small stretch check: stretching should affect only points above R220
if( r*R_EARTH < R220 ) then
@@ -363,10 +363,10 @@
gamma = (( r - R60)/( R35 - R60)) ! keeps r60 fixed
if (abs(gamma)<SMALLVAL) then
gamma=0.0d0
- end if
+ endif
else
gamma=0.0d0
- end if
+ endif
if(gamma < -0.0001d0 .or. gamma > 1.0001d0) &
stop 'incorrect value of gamma for moho from crust 2.0'
@@ -383,10 +383,10 @@
gamma=((r-R60)/(R35-R60)) ! keeps r60 fixed
if (abs(gamma)<SMALLVAL) then
gamma=0.0d0
- end if
+ endif
else
gamma=0.0d0
- end if
+ endif
if(gamma < -0.0001d0 .or. gamma > 1.0001d0) &
stop 'incorrect value of gamma for moho from crust 2.0'
@@ -404,13 +404,13 @@
gamma=(r-R220/R_EARTH)/(R60-R220/R_EARTH)
if (abs(gamma)<SMALLVAL) then
gamma=0.0d0
- end if
+ endif
else
gamma=0.0d0
- end if
+ endif
call move_point(ia,xelm,yelm,zelm,x,y,z,gamma,elevation,r)
- end if
+ endif
else if (moho > R25) then
! moho above r25
@@ -423,10 +423,10 @@
gamma=(r-R60)/(R35-R60) ! keeps r60 fixed
if (abs(gamma)<SMALLVAL) then
gamma=0.0d0
- end if
+ endif
else
gamma=0.0d0
- end if
+ endif
if(gamma < -0.0001d0 .or. gamma > 1.0001d0) &
stop 'incorrect value of gamma for moho from crust 2.0'
@@ -444,14 +444,14 @@
gamma=(r-R25)/(R15-R25) ! keeps mesh at r25 fixed
if (abs(gamma)<SMALLVAL) then
gamma=0.0d0
- end if
+ endif
else
gamma=0.0d0
- end if
+ endif
call move_point(ia,xelm,yelm,zelm,x,y,z,gamma,elevation,r)
- end if
- end if
+ endif
+ endif
end subroutine stretch_deep_moho
@@ -512,10 +512,10 @@
gamma = ((r-R220/R_EARTH)/(R35-R220/R_EARTH))
if (abs(gamma)<SMALLVAL) then
gamma=0.0d0
- end if
+ endif
else
gamma=0.0d0
- end if
+ endif
if(gamma < -0.0001d0 .or. gamma > 1.0001d0) &
stop 'incorrect value of gamma for moho from crust 2.0'
@@ -532,10 +532,10 @@
gamma=((r-R220/R_EARTH)/(R35-R220/R_EARTH))
if (abs(gamma)<SMALLVAL) then
gamma=0.0d0
- end if
+ endif
else
gamma=0.0d0
- end if
+ endif
if(gamma < -0.0001d0 .or. gamma > 1.0001d0) &
stop 'incorrect value of gamma for moho from crust 2.0'
@@ -552,10 +552,10 @@
gamma=(r-R220/R_EARTH)/(R35-R220/R_EARTH)
if (abs(gamma)<SMALLVAL) then
gamma=0.0d0
- end if
+ endif
else
gamma=0.0d0
- end if
+ endif
if(gamma < -0.0001d0 .or. gamma > 1.0001d0) &
stop 'incorrect value of gamma for moho from crust 2.0'
@@ -570,13 +570,13 @@
gamma=(r-R25)/(R15-R25)
if (abs(gamma)<SMALLVAL) then
gamma=0.0d0
- end if
+ endif
else
gamma=0.0d0
- end if
+ endif
call move_point(ia,xelm,yelm,zelm,x,y,z,gamma,elevation,r)
- end if
+ endif
endif
end subroutine stretch_moho
@@ -687,7 +687,7 @@
! if(r >= RMOHO/R_EARTH) then
!! stretching above the Moho
! gamma = (1.0d0 - r) / (1.0d0 - RMOHO/R_EARTH)
-! elseif(r>= R220/R_EARTH .and. r< RMOHO/R_EARTH) then
+! else if(r>= R220/R_EARTH .and. r< RMOHO/R_EARTH) then
!! stretching between R220 and RMOHO
! gamma = (r - R220/R_EARTH) / (RMOHO/R_EARTH - R220/R_EARTH)
! endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/write_AVS_DX_global_data.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/write_AVS_DX_global_data.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/write_AVS_DX_global_data.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -249,9 +249,9 @@
write(10,*) numpoin,sngl(xstore(i,j,k,ispec)),&
sngl(ystore(i,j,k,ispec)),sngl(zstore(i,j,k,ispec))
flag(i,j,k,ispec) = numpoin
- end do
- end do
- end do
+ enddo
+ enddo
+ enddo
enddo
close(10)
@@ -282,9 +282,9 @@
write(10,*) nelem,iglob1, &
iglob2,iglob3,iglob4,&
iglob5,iglob6,iglob7,iglob8
- end do
- end do
- end do
+ enddo
+ enddo
+ enddo
enddo
close(10)
@@ -344,7 +344,7 @@
Qmu(7)=dble(Qmustore(i+1,j+1,k+1,ispec))
Qmu(8)=dble(Qmustore(i,j+1,k+1,ispec))
Qmu_average=Qmu(1)
- end if
+ endif
!rho_average=sum(rho(1:4))/4.d0
!vp_average=sum(vp(1:4))/4.d0
!vs_average=sum(vs(1:4))/4.d0
@@ -356,11 +356,11 @@
write(1001,*) nelem,rho_average,vp_average,vs_average,Qmu_average
else
write(1001,*) nelem,rho_average,vp_average,vs_average
- end if
+ endif
- end do
- end do
- end do
+ enddo
+ enddo
+ enddo
enddo
close(1001)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/auto_ner.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/auto_ner.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/auto_ner.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -288,10 +288,10 @@
ner_test = ner_test + 1 ! Increment ner_test and
ratio = (dr / ner_test) / w ! look for a better
xi = dabs(ratio - 1.0d0) ! solution
- end do
+ enddo
rt(i) = dr / NER(i) / wt ! Find the Ratio of Top
rb(i) = dr / NER(i) / wb ! and Bottom for completeness
- end do
+ enddo
end subroutine auto_optimal_ner
@@ -344,7 +344,7 @@
max_edgemax = MAX(max_edgemax, edgemax)
min_edgemin = MIN(min_edgemin, edgemin)
max_aspect_ratio = MAX(max_aspect_ratio, aspect_ratio)
- end do
+ enddo
xi = (max_edgemax / min_edgemin)
! xi = abs(rcube_test - 981.0d0) / 45.0d0
! write(*,'(a,5(f14.4,2x))')'rcube, xi, ximin:-',rcube_test, xi, min_edgemin,max_edgemax,max_aspect_ratio
@@ -395,7 +395,7 @@
dist_cc_icb = dist_cc_icb * 2
endif
somme = somme + dist_cc_icb
- end do
+ enddo
dist_moy = somme / (nex_xi + 1)
ner = nint(dist_moy / ((PI * RICB_KM) / (2*nex_xi)))
@@ -444,7 +444,7 @@
sqrt( (pts(ix2,1) - pts(ix3,1))**2 + (pts(ix2,2) - pts(ix3,2))**2 )
edgemax = MAX(edgemax, edge)
edgemin = MIN(edgemin, edge)
- end do
+ enddo
end subroutine get_size_min_max
@@ -482,11 +482,11 @@
points(k,1) = x
points(k,2) = y
k = k + 1
- end do
+ enddo
nspec_chunks = nspec_chunks + 1
- end do
- end do
- end do
+ enddo
+ enddo
+ enddo
nspec_cube = 0
do ix = 0,(nex_xi-1)*2,2
@@ -496,10 +496,10 @@
points(k,1) = x
points(k,2) = y
k = k + 1
- end do
+ enddo
nspec_cube = nspec_cube + 1
- end do
- end do
+ enddo
+ enddo
end subroutine compute_IC_mesh
@@ -586,6 +586,6 @@
temp = x
x = -y
y = temp
- end if
+ endif
end subroutine compute_coordinate
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/calendar.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/calendar.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/calendar.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -455,10 +455,10 @@
!
! Explanation of all internal variables.
! jdref Julian Day on which 1 March begins in the reference year.
-! jmonth Month counter which equals month+1 if month .gt. 2
-! or month+13 if month .le. 2.
-! jyear Year index, jyear=iyear if iyear .gt. 0, jyear=iyear+1
-! if iyear .lt. 0. Thus, jyear does not skip year 0
+! jmonth Month counter which equals month+1 if month > 2
+! or month+13 if month <= 2.
+! jyear Year index, jyear=iyear if iyear > 0, jyear=iyear+1
+! if iyear < 0. Thus, jyear does not skip year 0
! like iyear does between BC and AD years.
! leap =1 if the year is a leap year, =0 if not.
! n1yr Number of complete individual years between iyear and
@@ -482,7 +482,7 @@
! this is 100*365 + 25 = 36525.
! nyrs Number of years from the beginning of yr400
! to the beginning of jyear. (Used for option +/-3).
-! yr400 The largest multiple of 400 years that is .le. jyear.
+! yr400 The largest multiple of 400 years that is <= jyear.
!
!
!----------------------------------------------------------------
@@ -502,7 +502,7 @@
if (abs(ioptn) <= 3) then
if (iyear > 0) then
jyear = iyear
- elseif (iyear == 0) then
+ else if (iyear == 0) then
write(*,*) 'For calndr(), you specified the nonexistent year 0'
stop
else
@@ -585,7 +585,7 @@
! OPTIONS -2 and +2:
! Given the day number of the year (idayct) and the year (iyear),
! compute the day of the month (iday) and the month (month).
- elseif (abs(ioptn) == 2) then
+ else if (abs(ioptn) == 2) then
!
if (idayct < 60+leap) then
month = (idayct-1)/31
@@ -604,7 +604,7 @@
! OPTIONS -3 and +3:
! Given a calendar date (iday,month,iyear), compute the Julian Day
! number (idayct) that starts at noon.
- elseif (abs(ioptn) == 3) then
+ else if (abs(ioptn) == 3) then
!
! Shift to a system where the year starts on 1 March, so January
! and February belong to the preceding year.
@@ -616,7 +616,7 @@
jmonth = month + 1
endif
!
-! Find the closest multiple of 400 years that is .le. jyear.
+! Find the closest multiple of 400 years that is <= jyear.
yr400 = (jyear/400)*400
! = multiple of 400 years at or less than jyear.
if (jyear < yr400) then
@@ -743,7 +743,7 @@
endif
endif
!
-! Adjust the year if it is .le. 0, and hence BC (Before Christ).
+! Adjust the year if it is <= 0, and hence BC (Before Christ).
if (iyear <= 0) then
iyear = iyear - 1
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/read_compute_parameters.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/read_compute_parameters.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/read_compute_parameters.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -263,7 +263,7 @@
if(mod(NEX_ETA/4,NPROC_ETA) /= 0) stop 'NEX_ETA must be a multiple of 4*NPROC_ETA'
if(mod(NEX_XI/8,NPROC_XI) /=0) CUT_SUPERBRICK_XI = .true.
if(mod(NEX_ETA/8,NPROC_ETA) /=0) CUT_SUPERBRICK_ETA = .true.
- elseif (SUPPRESS_CRUSTAL_MESH .or. .not. ADD_4TH_DOUBLING) then
+ else if (SUPPRESS_CRUSTAL_MESH .or. .not. ADD_4TH_DOUBLING) then
if(mod(NEX_XI,16) /= 0) stop 'NEX_XI must be a multiple of 16'
if(mod(NEX_ETA,16) /= 0) stop 'NEX_ETA must be a multiple of 16'
if(mod(NEX_XI/8,NPROC_XI) /= 0) stop 'NEX_XI must be a multiple of 8*NPROC_XI'
@@ -614,7 +614,7 @@
! the 670-discontinuity is moved up to 650 km depth.
if (REFERENCE_1D_MODEL == REFERENCE_MODEL_1DREF) then
NER_771_670 = NER_771_670 + 1
- end if
+ endif
!----
!---- change some values in the case of regular PREM with two crustal layers or of 3D models
@@ -1117,7 +1117,7 @@
rmaxs(14) = RICB / R_EARTH
rmins(14) = R_CENTRAL_CUBE / R_EARTH
- elseif (ONE_CRUST) then
+ else if (ONE_CRUST) then
! 1D models:
! in order to increase stability and therefore to allow cheaper
@@ -1521,7 +1521,7 @@
rmaxs(15) = RICB / R_EARTH
rmins(15) = R_CENTRAL_CUBE / R_EARTH
- elseif (ONE_CRUST) then
+ else if (ONE_CRUST) then
! 1D models:
! in order to increase stability and therefore to allow cheaper
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/read_value_parameters.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/read_value_parameters.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/read_value_parameters.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -37,7 +37,7 @@
integer ierr
call param_read(string_read, len(string_read), name, len(name), ierr)
- if (ierr .ne. 0) return
+ if (ierr /= 0) return
read(string_read,*,iostat=ierr) value_to_read
end subroutine read_value_integer
@@ -54,7 +54,7 @@
integer ierr
call param_read(string_read, len(string_read), name, len(name), ierr)
- if (ierr .ne. 0) return
+ if (ierr /= 0) return
read(string_read,*,iostat=ierr) value_to_read
end subroutine read_value_double_precision
@@ -71,7 +71,7 @@
integer ierr
call param_read(string_read, len(string_read), name, len(name), ierr)
- if (ierr .ne. 0) return
+ if (ierr /= 0) return
read(string_read,*,iostat=ierr) value_to_read
end subroutine read_value_logical
@@ -88,7 +88,7 @@
integer ierr
call param_read(string_read, len(string_read), name, len(name), ierr)
- if (ierr .ne. 0) return
+ if (ierr /= 0) return
value_to_read = string_read
end subroutine read_value_string
@@ -102,7 +102,7 @@
filename = 'DATA/Par_file'
call param_open(filename, len(filename), ierr)
- if (ierr .ne. 0) return
+ if (ierr /= 0) return
end subroutine open_parameter_file
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_Dev.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_Dev.F90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_Dev.F90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -313,7 +313,7 @@
NGLOB1D_RADIAL_IC,NCHUNKS_VAL,iphase)
!$OMP END MASTER
!$OMP BARRIER
- end if
+ endif
if(INCLUDE_CENTRAL_CUBE) then
if(iphase > 7 .and. iphase_CC <= 4) then
@@ -325,7 +325,7 @@
ibelm_bottom_inner_core,NSPEC2D_BOTTOM_IC,accel_inner_core,NDIM,iphase_CC)
!$OMP END MASTER
!$OMP BARRIER
- end if
+ endif
endif
endif
@@ -640,8 +640,8 @@
epsilondev(:,:,:,:,ispec) = epsilondev_loc(:,:,:,:)
endif
! end ispec loop
- end do
-!$OMP END DO
+ enddo
+!$OMP enddo
! end ispec_globe strided loop
enddo ! spectral element loop NSPEC_CRUST_MANTLE
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/convert_time.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/convert_time.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/convert_time.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -84,10 +84,10 @@
if (mon == 2) then
if (is_leap_year(yr) .and. (day < 1 .or. day > 29)) then
stop 'Error in convtime: February day out of range (1-29)'
- elseif (.not. is_leap_year(yr) .and. (day < 1 .or. day > 28)) then
+ else if (.not. is_leap_year(yr) .and. (day < 1 .or. day > 28)) then
stop 'Error in convtime: February day out of range (1-28)'
endif
- elseif (mon == 4 .or. mon == 6 .or. mon == 9 .or. mon == 11) then
+ else if (mon == 4 .or. mon == 6 .or. mon == 9 .or. mon == 11) then
if (day < 1 .or. day > 30) stop 'Error in convtime: day out of range (1-30)'
else
if (day < 1 .or. day > 31) stop 'Error in convtime: day out of range (1-31)'
@@ -214,7 +214,7 @@
tmon=itime-leap_mon(imon)
if (tmon > 0) then
goto 30
- elseif (tmon < 0) then
+ else if (tmon < 0) then
imon=imon-1
itime=itime-month(imon)
else
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_regular_points.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_regular_points.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_regular_points.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -485,7 +485,7 @@
if (z > 0) xi = 10
else
stop 'chunk number k < 6'
- end if
+ endif
end subroutine chunk_map
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_sources.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_sources.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_sources.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -720,7 +720,7 @@
else
if(isource < 10) then
write(plot_file,"('/plot_source_time_function',i1,'.txt')") isource
- elseif(isource < 100) then
+ else if(isource < 100) then
write(plot_file,"('/plot_source_time_function',i2,'.txt')") isource
else
write(plot_file,"('/plot_source_time_function',i3,'.txt')") isource
@@ -776,7 +776,7 @@
else
if(isource < 10) then
write(plot_file,"('/plot_source_spectrum',i1,'.txt')") isource
- elseif(isource < 100) then
+ else if(isource < 100) then
write(plot_file,"('/plot_source_spectrum',i2,'.txt')") isource
else
write(plot_file,"('/plot_source_spectrum',i3,'.txt')") isource
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/prepare_timerun.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/prepare_timerun.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/prepare_timerun.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -679,7 +679,7 @@
enddo
enddo
enddo
- enddo ! END DO CRUST MANTLE
+ enddo ! enddo CRUST MANTLE
! rescale in inner core
@@ -708,7 +708,7 @@
enddo
enddo
enddo
- enddo ! END DO INNER CORE
+ enddo ! enddo INNER CORE
! precompute Runge-Kutta coefficients
call get_attenuation_memory_values(tau_sigma_dble, deltat, alphaval_dble, betaval_dble, gammaval_dble)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -3235,7 +3235,7 @@
NSTEP,accel_crust_mantle,noise_sourcearray, &
ibool_crust_mantle,islice_selected_rec,ispec_selected_rec, &
it,irec_master_noise)
- elseif ( NOISE_TOMOGRAPHY == 2 ) then
+ else if ( NOISE_TOMOGRAPHY == 2 ) then
! second step of noise tomography, i.e., read the surface movie saved at every timestep
! use the movie to drive the ensemble forward wavefield
call noise_read_add_surface_movie(nmovie_points,accel_crust_mantle, &
@@ -3248,7 +3248,7 @@
! note the ensemble forward sources are generally distributed on the surface of the earth
! that's to say, the ensemble forward source is kind of a surface force density, not a body force density
! therefore, we must add it here, before applying the inverse of mass matrix
- elseif ( NOISE_TOMOGRAPHY == 3 ) then
+ else if ( NOISE_TOMOGRAPHY == 3 ) then
! third step of noise tomography, i.e., read the surface movie saved at every timestep
! use the movie to reconstruct the ensemble forward wavefield
! the ensemble adjoint wavefield is done as usual
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_output_SAC.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_output_SAC.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_output_SAC.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -286,10 +286,10 @@
if (NZJDAY > 365 .and. .not. is_leap_year(NZYEAR)) then
NZJDAY = mod(NZJDAY,365)
NZYEAR = yr + 1
- elseif (NZJDAY > 366 .and. is_leap_year(NZYEAR)) then
+ else if (NZJDAY > 366 .and. is_leap_year(NZYEAR)) then
NZJDAY = mod(NZJDAY,366)
NZYEAR = yr + 1
- elseif (NZJDAY == 366 .and. is_leap_year(NZYEAR)) then
+ else if (NZJDAY == 366 .and. is_leap_year(NZYEAR)) then
NZJDAY = 366
endif
endif
@@ -631,7 +631,7 @@
if (CUSTOM_REAL == SIZE_REAL) then
call write_n_real(seismogram_tmp(iorientation,1:seismo_current),seismo_current)
- elseif (CUSTOM_REAL == SIZE_DOUBLE) then
+ else if (CUSTOM_REAL == SIZE_DOUBLE) then
call write_n_real(real(seismogram_tmp(iorientation,1:seismo_current)),seismo_current)
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_seismograms.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_seismograms.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_seismograms.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -391,9 +391,9 @@
phi = backaz
if (phi>180.) then
phi = phi-180.
- elseif (phi<180.) then
+ else if (phi<180.) then
phi = phi+180.
- elseif (phi==180.) then
+ else if (phi==180.) then
phi = backaz
endif
@@ -407,7 +407,7 @@
seismogram_tmp(iorientation,isample) = &
cphi * one_seismogram(1,isample) + sphi * one_seismogram(2,isample)
enddo
- elseif (iorientation == 5) then ! transverse component
+ else if (iorientation == 5) then ! transverse component
do isample = 1,seismo_current
seismogram_tmp(iorientation,isample) = &
-1*sphi * one_seismogram(1,isample) + cphi * one_seismogram(2,isample)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/utils/Profiles/write_profile.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/utils/Profiles/write_profile.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/utils/Profiles/write_profile.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -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
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/utils/attenuation/attenuation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/utils/attenuation/attenuation.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/utils/attenuation/attenuation.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -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/trunk/utils/attenuation/attenuation_output.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/utils/attenuation/attenuation_output.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/utils/attenuation/attenuation_output.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -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/trunk/utils/attenuation/attenuation_simplex.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/utils/attenuation/attenuation_simplex.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/utils/attenuation/attenuation_simplex.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -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
-
+ enddo
+
call qsort(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)
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,17 +553,17 @@
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
@@ -593,14 +593,14 @@
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,11 +608,11 @@
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
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/utils/extract_database/extract_database.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/utils/extract_database/extract_database.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/utils/extract_database/extract_database.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -21,7 +21,7 @@
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
+ .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/trunk/utils/funaro_routines_bernhard/routines_SEM_Bernhard_diff_lagrange.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/utils/funaro_routines_bernhard/routines_SEM_Bernhard_diff_lagrange.f90 2013-04-18 19:15:40 UTC (rev 21898)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/utils/funaro_routines_bernhard/routines_SEM_Bernhard_diff_lagrange.f90 2013-04-18 19:23:17 UTC (rev 21899)
@@ -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