[cig-commits] r16535 - seismo/3D/SPECFEM3D_GLOBE/tags/v5.0.1
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Mon Apr 12 15:36:44 PDT 2010
Author: danielpeter
Date: 2010-04-12 15:36:44 -0700 (Mon, 12 Apr 2010)
New Revision: 16535
Modified:
seismo/3D/SPECFEM3D_GLOBE/tags/v5.0.1/check_simulation_stability.f90
seismo/3D/SPECFEM3D_GLOBE/tags/v5.0.1/exit_mpi.f90
seismo/3D/SPECFEM3D_GLOBE/tags/v5.0.1/model_crust.f90
seismo/3D/SPECFEM3D_GLOBE/tags/v5.0.1/model_crustmaps.f90
seismo/3D/SPECFEM3D_GLOBE/tags/v5.0.1/specfem3D.f90
Log:
copying changes to tags version 5.0.1: smoothes sediment basin edges
Modified: seismo/3D/SPECFEM3D_GLOBE/tags/v5.0.1/check_simulation_stability.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/tags/v5.0.1/check_simulation_stability.f90 2010-04-12 22:34:21 UTC (rev 16534)
+++ seismo/3D/SPECFEM3D_GLOBE/tags/v5.0.1/check_simulation_stability.f90 2010-04-12 22:36:44 UTC (rev 16535)
@@ -87,7 +87,7 @@
timestamp_remote,year_remote,mon_remote,day_remote,hr_remote,minutes_remote,day_of_week_remote
integer :: ier
integer, external :: idaywk
-
+
double precision,parameter :: scale_displ = R_EARTH
@@ -98,7 +98,7 @@
maxval(sqrt(displ_inner_core(1,:)**2 + displ_inner_core(2,:)**2 + displ_inner_core(3,:)**2)))
Ufluidnorm = maxval(abs(displ_outer_core))
-
+
! compute the maximum of the maxima for all the slices using an MPI reduction
call MPI_REDUCE(Usolidnorm,Usolidnorm_all,1,CUSTOM_MPI_TYPE,MPI_MAX,0, &
MPI_COMM_WORLD,ier)
Modified: seismo/3D/SPECFEM3D_GLOBE/tags/v5.0.1/exit_mpi.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/tags/v5.0.1/exit_mpi.f90 2010-04-12 22:34:21 UTC (rev 16534)
+++ seismo/3D/SPECFEM3D_GLOBE/tags/v5.0.1/exit_mpi.f90 2010-04-12 22:36:44 UTC (rev 16535)
@@ -63,9 +63,18 @@
if(myrank == 0 .and. IMAIN /= ISTANDARD_OUTPUT) close(IMAIN)
! stop all the MPI processes, and exit
+! note: MPI_ABORT does not return, and does exit the
+! program with an error code of 30
call MPI_ABORT(MPI_COMM_WORLD,30,ier)
- stop 'error, program ended in exit_MPI'
+! otherwise: there is no standard behaviour to exit with an error code in fortran,
+! however most compilers do recognize this as an error code stop statement;
+! to check stop code in terminal: > echo $?
+ stop 30
+
+ ! or just exit with message:
+ !stop 'error, program ended in exit_MPI'
+
end subroutine exit_MPI
!
Modified: seismo/3D/SPECFEM3D_GLOBE/tags/v5.0.1/model_crust.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/tags/v5.0.1/model_crust.f90 2010-04-12 22:34:21 UTC (rev 16534)
+++ seismo/3D/SPECFEM3D_GLOBE/tags/v5.0.1/model_crust.f90 2010-04-12 22:36:44 UTC (rev 16535)
@@ -130,13 +130,15 @@
found_crust = .true.
- if(x > x3 .and. INCLUDE_SEDIMENTS_CRUST &
- .and. h_sed >= MINIMUM_SEDIMENT_THICKNESS) then
+! if(x > x3 .and. INCLUDE_SEDIMENTS_CRUST &
+! .and. h_sed >= MINIMUM_SEDIMENT_THICKNESS) then
+ if(x > x3 .and. INCLUDE_SEDIMENTS_CRUST ) then
vp = vps(3)
vs = vss(3)
rho = rhos(3)
- else if(x > x4 .and. INCLUDE_SEDIMENTS_CRUST &
- .and. h_sed >= MINIMUM_SEDIMENT_THICKNESS) then
+! else if(x > x4 .and. INCLUDE_SEDIMENTS_CRUST &
+! .and. h_sed >= MINIMUM_SEDIMENT_THICKNESS) then
+ else if(x > x4 .and. INCLUDE_SEDIMENTS_CRUST ) then
vp = vps(4)
vs = vss(4)
rho = rhos(4)
@@ -282,6 +284,7 @@
double precision xlon(NTHETA*NPHI),xlat(NTHETA*NPHI),weight(NTHETA*NPHI)
double precision rhol(NLAYERS_CRUST),thickl(NLAYERS_CRUST),velpl(NLAYERS_CRUST),velsl(NLAYERS_CRUST)
double precision weightl,cap_degree,dist
+ double precision h_sed
integer i,icolat,ilon,ierr
character(len=2) crustaltype
@@ -331,6 +334,19 @@
code,thlr,velocp,velocs,dens,ierr)
if(ierr /= 0) stop 'error in routine get_crust_structure'
+ ! sediment thickness
+ h_sed = thickl(3) + thickl(4)
+
+ ! takes upper crust value if sediment too thin
+ if( h_sed < MINIMUM_SEDIMENT_THICKNESS ) then
+ velpl(3) = velpl(5)
+ velpl(4) = velpl(5)
+ velsl(3) = velsl(5)
+ velsl(4) = velsl(5)
+ rhol(3) = rhol(5)
+ rhol(4) = rhol(5)
+ endif
+
! weighting value
weightl = weight(i)
Modified: seismo/3D/SPECFEM3D_GLOBE/tags/v5.0.1/model_crustmaps.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/tags/v5.0.1/model_crustmaps.f90 2010-04-12 22:34:21 UTC (rev 16534)
+++ seismo/3D/SPECFEM3D_GLOBE/tags/v5.0.1/model_crustmaps.f90 2010-04-12 22:36:44 UTC (rev 16535)
@@ -403,13 +403,15 @@
x7 = (R_EARTH-(h_uc+thicks(4)+thicks(5))*1000.0d0)/R_EARTH
found_crust = .true.
- if(x > x3 .and. INCLUDE_SEDIMENTS_CRUST &
- .and. h_sed > MINIMUM_SEDIMENT_THICKNESS) then
+! if(x > x3 .and. INCLUDE_SEDIMENTS_CRUST &
+! .and. h_sed > MINIMUM_SEDIMENT_THICKNESS) then
+ if(x > x3 .and. INCLUDE_SEDIMENTS_CRUST ) then
vp = vps(1)
vs = vss(1)
rho = rhos(1)
- else if(x > x4 .and. INCLUDE_SEDIMENTS_CRUST &
- .and. h_sed > MINIMUM_SEDIMENT_THICKNESS) then
+! else if(x > x4 .and. INCLUDE_SEDIMENTS_CRUST &
+! .and. h_sed > MINIMUM_SEDIMENT_THICKNESS) then
+ else if(x > x4 .and. INCLUDE_SEDIMENTS_CRUST ) then
vp = vps(2)
vs = vss(2)
rho = rhos(2)
@@ -516,6 +518,7 @@
double precision rhol(NLAYERS_CRUSTMAP),thickl(NLAYERS_CRUSTMAP), &
velpl(NLAYERS_CRUSTMAP),velsl(NLAYERS_CRUSTMAP)
double precision weightl,cap_degree,dist
+ double precision h_sed
integer num_points
integer i,ipoin,iupcolat,ileftlng,irightlng
@@ -652,6 +655,19 @@
enddo
endif
+ ! sediment thickness
+ h_sed = thickl(1) + thickl(2)
+
+ ! takes upper crust value if sediment too thin
+ if( h_sed < MINIMUM_SEDIMENT_THICKNESS ) then
+ velpl(1) = velpl(3)
+ velpl(2) = velpl(3)
+ velsl(1) = velsl(3)
+ velsl(2) = velsl(3)
+ rhol(1) = rhol(3)
+ rhol(2) = rhol(3)
+ endif
+
! total, smoothed values
rhos(:) = rhos(:) + weightl*rhol(:)
thicks(:) = thicks(:) + weightl*thickl(:)
Modified: seismo/3D/SPECFEM3D_GLOBE/tags/v5.0.1/specfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/tags/v5.0.1/specfem3D.f90 2010-04-12 22:34:21 UTC (rev 16534)
+++ seismo/3D/SPECFEM3D_GLOBE/tags/v5.0.1/specfem3D.f90 2010-04-12 22:36:44 UTC (rev 16535)
@@ -750,8 +750,11 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NSPEC2D_ICB) :: icb_kl, icb_kl_top, icb_kl_bot
logical :: fluid_solid_boundary
- integer i,ier
+ integer :: i,ier
+ !daniel: debugging
+ !integer:: indx(1)
+
! ************** PROGRAM STARTS HERE **************
!
!-------------------------------------------------------------------------------------------------
@@ -1874,6 +1877,23 @@
+ deltat*eps_trace_over_3_crust_mantle(:,:,:,:)
endif
+ ! daniel: debugging
+ !if( maxval(displ_crust_mantle(1,:)**2 + &
+ ! displ_crust_mantle(2,:)**2 + displ_crust_mantle(3,:)**2) > 1.e4 ) then
+ ! print*,'slice',myrank
+ ! print*,' crust_mantle displ:', maxval(displ_crust_mantle(1,:)),maxval(displ_crust_mantle(2,:)),maxval(displ_crust_mantle(3,:))
+ ! print*,' indxs: ',maxloc( displ_crust_mantle(1,:)),maxloc( displ_crust_mantle(2,:)),maxloc( displ_crust_mantle(3,:))
+ ! indx = maxloc( displ_crust_mantle(3,:) )
+ ! rval = xstore_crust_mantle(indx(1))
+ ! thetaval = ystore_crust_mantle(indx(1))
+ ! phival = zstore_crust_mantle(indx(1))
+ ! !thetaval = PI/2.0d0-datan(1.006760466d0*dcos(dble(thetaval))/dmax1(TINYVAL,dsin(dble(thetaval))))
+ ! print*,'r/lat/lon:',rval*R_EARTH_KM,90.0-thetaval*180./PI,phival*180./PI
+ ! call rthetaphi_2_xyz(rval,thetaval,phival,xstore_crust_mantle(indx(1)),&
+ ! ystore_crust_mantle(indx(1)),zstore_crust_mantle(indx(1)))
+ ! print*,'x/y/z:',rval,thetaval,phival
+ ! call exit_MPI(myrank,'error stability')
+ !endif
! compute the maximum of the norm of the displacement
More information about the CIG-COMMITS
mailing list