[cig-commits] r12588 - in seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta: . OUTPUT_FILES src
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Thu Aug 7 18:55:06 PDT 2008
Author: dkomati1
Date: 2008-08-07 18:55:06 -0700 (Thu, 07 Aug 2008)
New Revision: 12588
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/OUTPUT_FILES/values_from_mesher.h
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/attenuation_model.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/declarations_mesher.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_jacobian_discontinuities.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/meshfem3D.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/read_compute_parameters.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/recompute_jacobian.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/s362ani.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.f90
Log:
converted create_regions_mesh.f90 from heap to stack and moved many of its reusable arrays to meshfem3D.f90; this final version works (but anisotropy is not tested, and attenuation is not implemented yet).
Also converted the Makefile from Fortran95 to Fortran2003 for Intel ifort and GNU gfortran.
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile 2008-08-07 23:32:02 UTC (rev 12587)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/Makefile 2008-08-08 01:55:06 UTC (rev 12588)
@@ -36,16 +36,15 @@
#
FC = ifort
MPIFC = mpif90
-#FLAGS_NO_CHECK = -O1 -vec-report0 -e95 -implicitnone -warn truncated_source -warn argument_checking -warn unused -warn declarations -std95 -check nobounds -align sequence -assume byterecl -ftrapuv -fpe0 -ftz -traceback
-FLAGS_NO_CHECK = -O1 -vec-report0 -e95 -implicitnone -warn truncated_source -warn argument_checking -warn unused -warn declarations -std95 -check all -align sequence -assume byterecl -ftrapuv -fpe0 -ftz -traceback
-#FLAGS_NO_CHECK = -O3 -xP -vec-report0 -e95 -implicitnone -warn truncated_source -warn argument_checking -warn unused -warn declarations -std95 -check nobounds -align sequence -assume byterecl -fpe3 -ftz
+FLAGS_NO_CHECK = -O1 -vec-report0 -e03 -std03 -implicitnone -warn truncated_source -warn argument_checking -warn unused -warn declarations -warn alignments -warn ignore_loc -warn usage -check all -align sequence -assume byterecl -fpe0 -ftz -traceback -ftrapuv
+#FLAGS_NO_CHECK = -O3 -xP -vec-report0 -e03 -std03 -implicitnone -warn truncated_source -warn argument_checking -warn unused -warn declarations -warn alignments -warn ignore_loc -warn usage -check nobounds -align sequence -assume byterecl -fpe3 -ftz
#
# GNU gfortran
#
#FC = gfortran
#MPIFC = /opt/mpich2_gfortran/bin/mpif90
-#FLAGS_NO_CHECK = -std=gnu -fimplicit-none -frange-check -O3 -fmax-errors=10 -pedantic -pedantic-errors -Waliasing -Wampersand -Wcharacter-truncation -Wline-truncation -Wsurprising -Wno-tabs -Wunderflow -fno-trapping-math -fbounds-check
+#FLAGS_NO_CHECK = -std=f2003 -fimplicit-none -frange-check -O3 -fmax-errors=10 -pedantic -pedantic-errors -Waliasing -Wampersand -Wcharacter-truncation -Wline-truncation -Wsurprising -Wno-tabs -Wunderflow -fno-trapping-math # -fbounds-check
#
# Portland pgf90
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/OUTPUT_FILES/values_from_mesher.h
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/OUTPUT_FILES/values_from_mesher.h 2008-08-07 23:32:02 UTC (rev 12587)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/OUTPUT_FILES/values_from_mesher.h 2008-08-08 01:55:06 UTC (rev 12588)
@@ -6,11 +6,11 @@
! ---------------
!
!
- ! number of chunks = 1
+ ! number of chunks = 3
!
! these statistics do not include the central cube
!
- ! number of processors = 4
+ ! number of processors = 12
!
! maximum number of points per region = 576013
!
@@ -83,8 +83,8 @@
! (if significantly more, the job will not run by lack of memory)
! (if significantly less, you waste a significant amount of memory)
!
- ! size of static arrays for all slices = 0.325386047363281 GB
- ! = 3.177598118782043E-004 TB
+ ! size of static arrays for all slices = 0.976158142089844 GB
+ ! = 9.532794356346130E-004 TB
!
integer, parameter :: NEX_XI_VAL = 64
@@ -157,12 +157,12 @@
integer, parameter :: NGLOB2DMAX_YMIN_YMAX_IC = 178
integer, parameter :: NPROC_XI_VAL = 2
integer, parameter :: NPROC_ETA_VAL = 2
- integer, parameter :: NCHUNKS_VAL = 1
- integer, parameter :: NPROCTOT_VAL = 4
+ integer, parameter :: NCHUNKS_VAL = 3
+ integer, parameter :: NPROCTOT_VAL = 12
integer, parameter :: NGLOB2DMAX_XY_VAL_CM = 8574
integer, parameter :: NGLOB2DMAX_XY_VAL_OC = 2134
integer, parameter :: NGLOB2DMAX_XY_VAL_IC = 178
- integer, parameter :: NUMMSGS_FACES_VAL = 2
+ integer, parameter :: NUMMSGS_FACES_VAL = 6
integer, parameter :: NCORNERSCHUNKS_VAL = 1
integer, parameter :: ATT1 = 1
integer, parameter :: ATT2 = 1
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/attenuation_model.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/attenuation_model.f90 2008-08-07 23:32:02 UTC (rev 12587)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/attenuation_model.f90 2008-08-08 01:55:06 UTC (rev 12588)
@@ -245,12 +245,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 attenuation_conversion(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 attenuation_model_setup
@@ -1119,7 +1119,7 @@
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,AS_V)
+ call fminsearch(attenuation_eval, tau_e, n, iterations, min_value, err,AS_V)
if(err > 0) then
write(*,*)'Search did not converge for an attenuation of ', Q_real
write(*,*)' Iterations: ', iterations
@@ -1249,12 +1249,10 @@
do i = 1,nf
w = 2.0d0 * PI * 10**f(i)
do j = 1,nsls
-! write(*,*)j,tau_s(j),tau_e(j)
demon = 1.0d0 + w**2 * tau_s(j)**2
A(i) = A(i) + ((1.0d0 + (w**2 * tau_e(j) * tau_s(j)))/ demon)
B(i) = B(i) + ((w * (tau_e(j) - tau_s(j))) / demon)
- end do
-! write(*,*)A(i),B(i),10**f(i)
+ enddo
enddo
end subroutine attenuation_maxwell
@@ -1326,8 +1324,8 @@
! subroutine fminsearch
! - Computes the minimization of funk(x(n)) using the simplex method
-! - This subroutine is copied from Matlab fminsearch.m
-! and modified to suit my nefarious needs
+! - This subroutine is similar to Matlab's fminsearch.m
+! and modified to suit our needs
! Input
! funk = double precision function with one input parameter
! double precision function the_funk(x)
@@ -1344,16 +1342,12 @@
! 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
-! 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, AS_V)
+subroutine fminsearch(funk, x, n, itercount, tolf, err, AS_V)
implicit none
@@ -1379,7 +1373,7 @@
integer n
double precision x(n) ! Also Output
- integer itercount, prnt, err
+ integer itercount, err
double precision tolf
!Internal
@@ -1462,22 +1456,12 @@
how = initial
itercount = 1
func_evals = n+1
- if(prnt == 3) then
- write(*,*)'Iterations Funk Evals Value How'
- write(*,*)itercount, func_evals, fv(1), how
- endif
- if(prnt == 4) then
- write(*,*)'How: ',how
- write(*,*)'V: ', v
- write(*,*)'fv: ',fv
- write(*,*)'evals: ',func_evals
- endif
do while (func_evals < maxfun .AND. itercount < maxiter)
if(max_size_simplex(v,n) <= tolx .AND. &
max_value(fv,n+1) <= tolf) then
- goto 666
+ goto 888
endif
how = none
@@ -1565,27 +1549,18 @@
v = vtmp
itercount = itercount + 1
- if (prnt == 3) then
- write(*,*)itercount, func_evals, fv(1), how
- elseif (prnt == 4) then
- write(*,*)
- write(*,*)'How: ',how
- write(*,*)'v: ',v
- write(*,*)'fv: ',fv
- write(*,*)'evals: ',func_evals
- endif
enddo
if(func_evals > maxfun) then
write(*,*)'function evaluations exceeded prescribed limit', maxfun
- err = 1
+ stop 'err = 1'
endif
if(itercount > maxiter) then
write(*,*)'iterations exceeded prescribed limit', maxiter
- err = 2
+ stop 'err = 2'
endif
-666 continue
+888 continue
x = v(:,1)
tolf = fv(1)
@@ -1771,7 +1746,7 @@
steps(r) = i-1
r = r + 1
endif
- end do
+ enddo
steps(r) = n
! Run spline for each piece
@@ -1894,7 +1869,6 @@
Qmu=600.0d0
Qkappa = 57827.0d0
else
- write(*,*)'iflag:',iflag
call exit_MPI_without_rank('Invalid idoubling flag in attenuation_model_1D_prem from get_model()')
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.f90 2008-08-07 23:32:02 UTC (rev 12587)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/create_regions_mesh.f90 2008-08-08 01:55:06 UTC (rev 12588)
@@ -59,7 +59,15 @@
normal_xmin,normal_xmax,normal_ymin, &
normal_ymax,normal_bottom,normal_top, &
kappavstore,kappahstore,muvstore,muhstore,eta_anisostore,rmass,xelm_store,yelm_store,zelm_store, &
- npoin2D_xi,npoin2D_eta)
+ npoin2D_xi,npoin2D_eta, &
+ xigll,wxgll, yigll,wygll, zigll,wzgll, shape3D, dershape3D, shape2D_x, shape2D_y, shape2D_bottom, shape2D_top, &
+ dershape2D_x, dershape2D_y, dershape2D_bottom, dershape2D_top, rhostore_local,kappavstore_local, &
+ c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
+ c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
+ c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
+ iboun, locval, ifseg, xp,yp,zp, rmass_ocean_load, mask_ibool, copy_ibool_ori, iMPIcut_xi,iMPIcut_eta, &
+ rho_vp,rho_vs, Qmu_store, tau_e_store, ibelm_moho_top, ibelm_moho_bot, ibelm_400_top, ibelm_400_bot, &
+ ibelm_670_top, ibelm_670_bot, normal_moho, normal_400, normal_670, jacobian2D_moho, jacobian2D_400, jacobian2D_670)
! create the different regions of the mesh
@@ -119,9 +127,11 @@
logical, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: this_region_has_a_doubling
integer :: ignod,ner_without_doubling,ispec_superbrick,ilayer,ilayer_loop,ix_elem,iy_elem,iz_elem, &
- ifirst_region,ilast_region,ratio_divide_central_cube
- integer, dimension(:), allocatable :: perm_layer
+ ifirst_layer,ilast_layer,ratio_divide_central_cube
+! allocate this automatic array in the memory stack to avoid memory fragmentation with "allocate()"
+ integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: permutation_layer
+
! mesh doubling superbrick
integer, dimension(NGNOD_EIGHT_CORNERS,NSPEC_DOUBLING_SUPERBRICK) :: ibool_superbrick
@@ -412,17 +422,26 @@
! code for the four regions of the mesh
integer iregion_code
-! Gauss-Lobatto-Legendre points and weights of integration
- double precision, dimension(:), allocatable :: xigll,yigll,zigll,wxgll,wygll,wzgll
+! Gauss-Lobatto-Legendre points of integration and weights
+ double precision, dimension(NGLLX) :: xigll,wxgll
+ double precision, dimension(NGLLY) :: yigll,wygll
+ double precision, dimension(NGLLZ) :: zigll,wzgll
! 3D shape functions and their derivatives
- double precision, dimension(:,:,:,:), allocatable :: shape3D
- double precision, dimension(:,:,:,:,:), allocatable :: dershape3D
+ double precision, dimension(NGNOD,NGLLX,NGLLY,NGLLZ) :: shape3D
+ double precision, dimension(NDIM,NGNOD,NGLLX,NGLLY,NGLLZ) :: dershape3D
! 2D shape functions and their derivatives
- double precision, dimension(:,:,:), allocatable :: shape2D_x,shape2D_y,shape2D_bottom,shape2D_top
- double precision, dimension(:,:,:,:), allocatable :: dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top
+ double precision shape2D_x(NGNOD2D,NGLLY,NGLLZ)
+ double precision shape2D_y(NGNOD2D,NGLLX,NGLLZ)
+ double precision shape2D_bottom(NGNOD2D,NGLLX,NGLLY)
+ double precision shape2D_top(NGNOD2D,NGLLX,NGLLY)
+ double precision dershape2D_x(NDIM2D,NGNOD2D,NGLLY,NGLLZ)
+ double precision dershape2D_y(NDIM2D,NGNOD2D,NGLLX,NGLLZ)
+ double precision dershape2D_bottom(NDIM2D,NGNOD2D,NGLLX,NGLLY)
+ double precision dershape2D_top(NDIM2D,NGNOD2D,NGLLX,NGLLY)
+
! for ellipticity
integer nspl
double precision rspl(NR),espl(NR),espl2(NR)
@@ -440,19 +459,21 @@
! for model density and anisotropy
integer nspec_ani
+
!! DK DK changed this for the merged version
- real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: rhostore_local,kappavstore_local
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: rhostore_local,kappavstore_local
+
!! DK DK added this for merged version
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: kappavstore,kappahstore,muvstore,muhstore,eta_anisostore
! the 21 coefficients for an anisotropic medium in reduced notation
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_ANISO_MANTLE) :: &
c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
c36store,c44store,c45store,c46store,c55store,c56store,c66store
! boundary locator
- logical, dimension(:,:), allocatable :: iboun
+ logical, dimension(6,nspec) :: iboun
! proc numbers for MPI
integer myrank
@@ -463,9 +484,9 @@
double precision volume_local
! variables for creating array ibool (some arrays also used for AVS or DX files)
- integer, dimension(:), allocatable :: locval
- logical, dimension(:), allocatable :: ifseg
- double precision, dimension(:), allocatable :: xp,yp,zp
+ integer, dimension(npointot) :: locval
+ logical, dimension(npointot) :: ifseg
+ double precision, dimension(npointot) :: xp,yp,zp
integer nglob,nglob_theor,ieoff,ilocnum,ier,errorcode
@@ -476,12 +497,12 @@
double precision xval,yval,zval,rval,thetaval,phival
double precision lat,lon,colat
double precision elevation,height_oceans
- real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_ocean_load
+ real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE_OCEANS) :: rmass_ocean_load
! mask to sort ibool
- integer, dimension(:), allocatable :: mask_ibool
- integer, dimension(:,:,:,:), allocatable :: copy_ibool_ori
integer :: inumber
+ integer, dimension(nglob_theor) :: mask_ibool
+ integer, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: copy_ibool_ori
! boundary parameters locator
integer, dimension(NSPEC2DMAX_XMIN_XMAX) :: ibelm_xmin,ibelm_xmax
@@ -490,11 +511,10 @@
integer, dimension(NSPEC2D_TOP) :: ibelm_top
! MPI cut-planes parameters along xi and along eta
- logical, dimension(:,:), allocatable :: iMPIcut_xi,iMPIcut_eta
+ logical, dimension(2,nspec) :: iMPIcut_xi,iMPIcut_eta
-! Stacey, indices for Clayton-Engquist absorbing conditions
- integer, dimension(:,:), allocatable :: nimin,nimax,njmin,njmax,nkmin_xi,nkmin_eta
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: rho_vp,rho_vs
+! Stacey indices for Clayton-Engquist absorbing conditions
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STACEY) :: rho_vp,rho_vs
! number of elements on the boundaries
integer nspec2D_xmin,nspec2D_xmax,nspec2D_ymin,nspec2D_ymax
@@ -508,9 +528,9 @@
double precision, dimension(NDIM,NDIM) :: rotation_matrix
! attenuation
- double precision, dimension(:,:,:,:), allocatable :: Qmu_store
- double precision, dimension(:,:,:,:,:), allocatable :: tau_e_store
- double precision, dimension(N_SLS) :: tau_s
+ double precision, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: Qmu_store
+ double precision, dimension(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: tau_e_store
+ double precision, dimension(N_SLS) :: tau_s
double precision T_c_source
! **************
@@ -561,14 +581,25 @@
! now perform two passes in this part to be able to save memory
integer :: ipass
-! Boundary Mesh
+! boundary mesh
integer nex_eta_moho
- integer, dimension(:), allocatable :: ibelm_moho_top,ibelm_moho_bot,ibelm_400_top,ibelm_400_bot, &
- ibelm_670_top,ibelm_670_bot
- real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: normal_moho,normal_400,normal_670
- real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: jacobian2D_moho,jacobian2D_400,jacobian2D_670
integer ispec2D_moho_top,ispec2D_moho_bot,ispec2D_400_top,ispec2D_400_bot,ispec2D_670_top,ispec2D_670_bot
double precision r_moho,r_400,r_670
+
+ integer ibelm_moho_top(NSPEC2D_MOHO)
+ integer ibelm_moho_bot(NSPEC2D_MOHO)
+ integer ibelm_400_top(NSPEC2D_400)
+ integer ibelm_400_bot(NSPEC2D_400)
+ integer ibelm_670_top(NSPEC2D_670)
+ integer ibelm_670_bot(NSPEC2D_670)
+
+ real(kind=CUSTOM_REAL) normal_moho(NDIM,NGLLX,NGLLY,NSPEC2D_MOHO)
+ real(kind=CUSTOM_REAL) normal_400(NDIM,NGLLX,NGLLY,NSPEC2D_400)
+ real(kind=CUSTOM_REAL) normal_670(NDIM,NGLLX,NGLLY,NSPEC2D_670)
+ real(kind=CUSTOM_REAL) jacobian2D_moho(NGLLX,NGLLY,NSPEC2D_MOHO)
+ real(kind=CUSTOM_REAL) jacobian2D_400(NGLLX,NGLLY,NSPEC2D_400)
+ real(kind=CUSTOM_REAL) jacobian2D_670(NGLLX,NGLLY,NSPEC2D_670)
+
logical :: is_superbrick
! the height at which the central cube is cut
@@ -594,362 +625,23 @@
if(ATTENUATION .and. ATTENUATION_3D) then
T_c_source = AM_V%QT_c_source
tau_s(:) = AM_V%Qtau_s(:)
- allocate(Qmu_store(NGLLX,NGLLY,NGLLZ,nspec),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(tau_e_store(N_SLS,NGLLX,NGLLY,NGLLZ,nspec),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
else
- allocate(Qmu_store(1,1,1,1),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(tau_e_store(N_SLS,1,1,1,1),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
Qmu_store(1,1,1,1) = 0.0d0
tau_e_store(:,1,1,1,1) = 0.0d0
endif
-! Gauss-Lobatto-Legendre points of integration
- allocate(xigll(NGLLX),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(yigll(NGLLY),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(zigll(NGLLZ),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
-! Gauss-Lobatto-Legendre weights of integration
- allocate(wxgll(NGLLX),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(wygll(NGLLY),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(wzgll(NGLLZ),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
-! 3D shape functions and their derivatives
- allocate(shape3D(NGNOD,NGLLX,NGLLY,NGLLZ),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(dershape3D(NDIM,NGNOD,NGLLX,NGLLY,NGLLZ),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
-
-! 2D shape functions and their derivatives
- allocate(shape2D_x(NGNOD2D,NGLLY,NGLLZ),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(shape2D_y(NGNOD2D,NGLLX,NGLLZ),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(shape2D_bottom(NGNOD2D,NGLLX,NGLLY),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(shape2D_top(NGNOD2D,NGLLX,NGLLY),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(dershape2D_x(NDIM2D,NGNOD2D,NGLLY,NGLLZ),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(dershape2D_y(NDIM2D,NGNOD2D,NGLLX,NGLLZ),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(dershape2D_bottom(NDIM2D,NGNOD2D,NGLLX,NGLLY),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(dershape2D_top(NDIM2D,NGNOD2D,NGLLX,NGLLY),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
-
-! array with model density
-!! DK DK changed this for the merged version
- allocate(rhostore_local(NGLLX,NGLLY,NGLLZ),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
-!! DK DK added this for the merged version
- allocate(kappavstore_local(NGLLX,NGLLY,NGLLZ),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
-
! Stacey
if(NCHUNKS /= 6) then
nspec_stacey = nspec
else
nspec_stacey = 1
endif
- allocate(rho_vp(NGLLX,NGLLY,NGLLZ,nspec_stacey),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
- allocate(rho_vs(NGLLX,NGLLY,NGLLZ,nspec_stacey),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
-
+! anisotropy
nspec_ani = 1
if((ANISOTROPIC_INNER_CORE .and. iregion_code == IREGION_INNER_CORE) .or. &
(ANISOTROPIC_3D_MANTLE .and. iregion_code == IREGION_CRUST_MANTLE)) nspec_ani = nspec
- allocate(c11store(NGLLX,NGLLY,NGLLZ,nspec_ani),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(c12store(NGLLX,NGLLY,NGLLZ,nspec_ani),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(c13store(NGLLX,NGLLY,NGLLZ,nspec_ani),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(c14store(NGLLX,NGLLY,NGLLZ,nspec_ani),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(c15store(NGLLX,NGLLY,NGLLZ,nspec_ani),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(c16store(NGLLX,NGLLY,NGLLZ,nspec_ani),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(c22store(NGLLX,NGLLY,NGLLZ,nspec_ani),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(c23store(NGLLX,NGLLY,NGLLZ,nspec_ani),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(c24store(NGLLX,NGLLY,NGLLZ,nspec_ani),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(c25store(NGLLX,NGLLY,NGLLZ,nspec_ani),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(c26store(NGLLX,NGLLY,NGLLZ,nspec_ani),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(c33store(NGLLX,NGLLY,NGLLZ,nspec_ani),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(c34store(NGLLX,NGLLY,NGLLZ,nspec_ani),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(c35store(NGLLX,NGLLY,NGLLZ,nspec_ani),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(c36store(NGLLX,NGLLY,NGLLZ,nspec_ani),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(c44store(NGLLX,NGLLY,NGLLZ,nspec_ani),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(c45store(NGLLX,NGLLY,NGLLZ,nspec_ani),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(c46store(NGLLX,NGLLY,NGLLZ,nspec_ani),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(c55store(NGLLX,NGLLY,NGLLZ,nspec_ani),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(c56store(NGLLX,NGLLY,NGLLZ,nspec_ani),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(c66store(NGLLX,NGLLY,NGLLZ,nspec_ani),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
-
-! boundary locator
- allocate(iboun(6,nspec),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
-
-! Stacey
- allocate(nimin(2,NSPEC2DMAX_YMIN_YMAX),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(nimax(2,NSPEC2DMAX_YMIN_YMAX),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(njmin(2,NSPEC2DMAX_XMIN_XMAX),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(njmax(2,NSPEC2DMAX_XMIN_XMAX),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(nkmin_xi(2,NSPEC2DMAX_XMIN_XMAX),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(nkmin_eta(2,NSPEC2DMAX_YMIN_YMAX),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
-
-! MPI cut-planes parameters along xi and along eta
- allocate(iMPIcut_xi(2,nspec),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(iMPIcut_eta(2,nspec),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
-
! set up coordinates of the Gauss-Lobatto-Legendre points
call zwgljd(xigll,wxgll,NGLLX,GAUSSALPHA,GAUSSBETA)
call zwgljd(yigll,wygll,NGLLY,GAUSSALPHA,GAUSSBETA)
@@ -1002,16 +694,16 @@
! define the first and last layers that define this region
if(iregion_code == IREGION_CRUST_MANTLE) then
- ifirst_region = 1
- ilast_region = 10 + layer_shift
+ ifirst_layer = 1
+ ilast_layer = 10 + layer_shift
else if(iregion_code == IREGION_OUTER_CORE) then
- ifirst_region = 11 + layer_shift
- ilast_region = NUMBER_OF_MESH_LAYERS - 1
+ ifirst_layer = 11 + layer_shift
+ ilast_layer = NUMBER_OF_MESH_LAYERS - 1
else if(iregion_code == IREGION_INNER_CORE) then
- ifirst_region = NUMBER_OF_MESH_LAYERS
- ilast_region = NUMBER_OF_MESH_LAYERS
+ ifirst_layer = NUMBER_OF_MESH_LAYERS
+ ilast_layer = NUMBER_OF_MESH_LAYERS
else
call exit_MPI(myrank,'incorrect region code detected')
@@ -1028,19 +720,14 @@
last_layer_aniso=4
nb_layer_above_aniso = 2
endif
- allocate(perm_layer(ifirst_region:ilast_region),STAT=ier)
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
- perm_layer = (/ (i, i=ilast_region,ifirst_region,-1) /)
+ permutation_layer(ifirst_layer:ilast_layer) = (/ (i, i=ilast_layer,ifirst_layer,-1) /)
if(iregion_code == IREGION_CRUST_MANTLE) then
cpt=3
- perm_layer(1)=first_layer_aniso
- perm_layer(2)=last_layer_aniso
- do i = ilast_region,ifirst_region,-1
+ permutation_layer(1)=first_layer_aniso
+ permutation_layer(2)=last_layer_aniso
+ do i = ilast_layer,ifirst_layer,-1
if (i/=first_layer_aniso .and. i/=last_layer_aniso) then
- perm_layer(cpt) = i
+ permutation_layer(cpt) = i
cpt=cpt+1
endif
enddo
@@ -1077,61 +764,6 @@
! boundary mesh
if (ipass == 2 .and. SAVE_BOUNDARY_MESH .and. iregion_code == IREGION_CRUST_MANTLE) then
- allocate(ibelm_moho_top(NSPEC2D_MOHO),ibelm_moho_bot(NSPEC2D_MOHO),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(ibelm_400_top(NSPEC2D_400),ibelm_400_bot(NSPEC2D_400),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(ibelm_670_top(NSPEC2D_670),ibelm_670_bot(NSPEC2D_670),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(normal_moho(NDIM,NGLLX,NGLLY,NSPEC2D_MOHO),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(normal_400(NDIM,NGLLX,NGLLY,NSPEC2D_400),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(normal_670(NDIM,NGLLX,NGLLY,NSPEC2D_670),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(jacobian2D_moho(NGLLX,NGLLY,NSPEC2D_MOHO),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(jacobian2D_400(NGLLX,NGLLY,NSPEC2D_400),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(jacobian2D_670(NGLLX,NGLLY,NSPEC2D_670),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
-
ispec2D_moho_top = 0; ispec2D_moho_bot = 0
ispec2D_400_top = 0; ispec2D_400_bot = 0
ispec2D_670_top = 0; ispec2D_670_bot = 0
@@ -1139,16 +771,15 @@
nex_eta_moho = NEX_PER_PROC_ETA
r_moho = RMOHO/R_EARTH; r_400 = R400 / R_EARTH; r_670 = R670/R_EARTH
-
endif
! generate and count all the elements in this region of the mesh
ispec = 0
! loop on all the layers in this region of the mesh
- do ilayer_loop = ifirst_region,ilast_region
+ do ilayer_loop = ifirst_layer,ilast_layer
- ilayer = perm_layer(ilayer_loop)
+ ilayer = permutation_layer(ilayer_loop)
! determine the radii that define the shell
rmin = rmins(ilayer)
@@ -1157,7 +788,7 @@
if(iregion_code == IREGION_CRUST_MANTLE .and. ilayer_loop==3) then
FIRST_ELT_NON_ANISO = ispec+1
endif
- if(iregion_code == IREGION_CRUST_MANTLE .and. ilayer_loop==(ilast_region-nb_layer_above_aniso+1)) then
+ if(iregion_code == IREGION_CRUST_MANTLE .and. ilayer_loop==(ilast_layer-nb_layer_above_aniso+1)) then
FIRST_ELT_ABOVE_ANISO = ispec+1
endif
@@ -1239,10 +870,10 @@
if (iproc_eta == NPROC_ETA-1) iboun(4,ispec)= .true.
endif
! zmin & zmax
- if (iz_elem == ner(ilayer) .and. ilayer == ifirst_region) then
+ if (iz_elem == ner(ilayer) .and. ilayer == ifirst_layer) then
iboun(6,ispec)= .true.
endif
- if (iz_elem == 1 .and. ilayer == ilast_region) then ! defined if no doubling in this layer
+ if (iz_elem == 1 .and. ilayer == ilast_layer) then ! defined if no doubling in this layer
iboun(5,ispec)= .true.
endif
@@ -1320,7 +951,7 @@
iz_elem = ner(ilayer)
step_mult = 2
else
- if(iregion_code==IREGION_OUTER_CORE .and. ilayer==ilast_region .and. (CUT_SUPERBRICK_XI .or. CUT_SUPERBRICK_ETA)) then
+ if(iregion_code==IREGION_OUTER_CORE .and. ilayer==ilast_layer .and. (CUT_SUPERBRICK_XI .or. CUT_SUPERBRICK_ETA)) then
nspec_sb = NSPEC_DOUBLING_BASICBRICK
step_mult = 1
else
@@ -1428,7 +1059,8 @@
if(ispec > nspec) call exit_MPI(myrank,'ispec greater than nspec in mesh creation')
! new get_flag_boundaries
-! xmin & xmax
+
+! xmin and xmax
if (ix_elem == 1) then
iMPIcut_xi(1,ispec) = iboun_sb(ispec_superbrick,1)
if (iproc_xi == 0) iboun(1,ispec)= iboun_sb(ispec_superbrick,1)
@@ -1437,7 +1069,8 @@
iMPIcut_xi(2,ispec) = iboun_sb(ispec_superbrick,2)
if (iproc_xi == NPROC_XI-1) iboun(2,ispec)= iboun_sb(ispec_superbrick,2)
endif
-!! ymin & ymax
+
+! ymin and ymax
if (iy_elem == 1) then
iMPIcut_eta(1,ispec) = iboun_sb(ispec_superbrick,3)
if (iproc_eta == 0) iboun(3,ispec)= iboun_sb(ispec_superbrick,3)
@@ -1446,13 +1079,10 @@
iMPIcut_eta(2,ispec) = iboun_sb(ispec_superbrick,4)
if (iproc_eta == NPROC_ETA-1) iboun(4,ispec)= iboun_sb(ispec_superbrick,4)
endif
+
! zmax only
- if (ilayer==ifirst_region) then
- iboun(6,ispec)= iboun_sb(ispec_superbrick,6)
- endif
- if (ilayer==ilast_region .and. iz_elem==1) then
- iboun(5,ispec)= iboun_sb(ispec_superbrick,5)
- endif
+ if (ilayer==ifirst_layer) iboun(6,ispec)= iboun_sb(ispec_superbrick,6)
+ if (ilayer==ilast_layer .and. iz_elem==1) iboun(5,ispec)= iboun_sb(ispec_superbrick,5)
! define the doubling flag of this element
if(iregion_code /= IREGION_OUTER_CORE) idoubling(ispec) = doubling_index(ilayer)
@@ -1520,13 +1150,9 @@
deallocate(stretch_tab,STAT=ier )
if (ier /= 0) then
print *,"ERROR can not deallocate stretch_tab in create_regions_mesh ier=",ier
+ stop 'error in deallocate'
endif
-
endif
- deallocate(perm_layer,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate perm_layer in create_regions_mesh ier=",ier
- endif
!---
@@ -1693,37 +1319,6 @@
! only create global addressing and the MPI buffers in the first pass
if(ipass == 1) then
- ! allocate memory for arrays
- allocate(locval(npointot),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(ifseg(npointot),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(xp(npointot),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(yp(npointot),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(zp(npointot),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
locval = 0
ifseg = .false.
xp = 0.d0
@@ -1749,22 +1344,6 @@
call get_global(nspec,xp,yp,zp,ibool,locval,ifseg,nglob,npointot)
- deallocate(xp,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate xp in create_regions_mesh ier=",ier
- endif
-
- deallocate(yp,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate yp in create_regions_mesh ier=",ier
- endif
-
- deallocate(zp,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate zp in create_regions_mesh ier=",ier
- endif
-
-
! check that number of points found equals theoretical value
if(nglob /= nglob_theor) then
write(errmsg,*) 'incorrect total number of points found: myrank,nglob,nglob_theor,ipass,iregion_code = ',&
@@ -1774,21 +1353,8 @@
if(minval(ibool) /= 1 .or. maxval(ibool) /= nglob_theor) call exit_MPI(myrank,'incorrect global numbering')
- ! create a new indirect addressing to reduce cache misses in memory access in the solver
- ! this is *critical* to improve performance in the solver
- allocate(copy_ibool_ori(NGLLX,NGLLY,NGLLZ,nspec),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate copy_ibool_ori in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
- allocate(mask_ibool(nglob),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate mask_ibool in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
-
-
+! create a new indirect addressing to reduce cache misses in memory access in the solver
+! this is *critical* to improve performance in the solver
mask_ibool(:) = -1
copy_ibool_ori(:,:,:,:) = ibool(:,:,:,:)
@@ -1811,17 +1377,6 @@
enddo
enddo
- deallocate(copy_ibool_ori,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate copy_ibool_ori in create_regions_mesh ier=",ier
- endif
-
- deallocate(mask_ibool,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate mask_ibool in create_regions_mesh ier=",ier
- endif
-
-
if(minval(ibool) /= 1 .or. maxval(ibool) /= nglob_theor) call exit_MPI(myrank,'incorrect global numbering after sorting')
! create MPI buffers
@@ -1842,16 +1397,6 @@
zread1D_leftxi_lefteta, zread1D_rightxi_lefteta, zread1D_leftxi_righteta, zread1D_rightxi_righteta, &
iregion_code)
- deallocate(locval,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate locval in create_regions_mesh ier=",ier
- endif
-
- deallocate(ifseg,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate ifseg in create_regions_mesh ier=",ier
- endif
-
! only create mass matrix and save all the final arrays in the second pass
else if(ipass == 2) then
@@ -1885,11 +1430,6 @@
! adding ocean load mass matrix at the top of the crust for oceans
nglob_oceans = nglob
- allocate(rmass_ocean_load(nglob_oceans),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate rmass_ocean_load in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
! create ocean load mass matrix for degrees of freedom at ocean bottom
rmass_ocean_load(:) = 0._CUSTOM_REAL
@@ -1965,11 +1505,6 @@
! allocate dummy array if no oceans
nglob_oceans = 1
- allocate(rmass_ocean_load(nglob_oceans),STAT=ier )
- if (ier /= 0) then
- print *,"ABORTING can not allocate rmass_ocean_load in create_regions_mesh ier=",ier
- call MPI_Abort(MPI_COMM_WORLD,errorcode,ier)
- endif
endif
@@ -2000,31 +1535,6 @@
write(27) normal_670
close(27)
- deallocate(ibelm_moho_top,ibelm_moho_bot,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate ibelm_moho_top in create_regions_mesh ier=",ier
- endif
-
- deallocate(ibelm_400_top,ibelm_400_bot,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate ibelm_400_top in create_regions_mesh ier=",ier
- endif
-
- deallocate(ibelm_670_top,ibelm_670_bot,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate ibelm_670_top,ibelm_670_bot in create_regions_mesh ier=",ier
- endif
-
- deallocate(normal_moho,normal_400,normal_670,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate normal_moho,normal_400,normal_670 in create_regions_mesh ier=",ier
- endif
-
- deallocate(jacobian2D_moho,jacobian2D_400,jacobian2D_670,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate jacobian2D_moho,jacobian2D_400,jacobian2D_670 in create_regions_mesh ier=",ier
- endif
-
endif
! compute volume, bottom and top area of that part of the slice
@@ -2040,7 +1550,7 @@
weight = wxgll(i)*wygll(j)*wzgll(k)
! compute the jacobian
-!! DK DK for merged version the jacobian is not stored anymore and therefore not valid anymore
+!! DK DK in merged version the jacobian is not stored anymore and therefore not valid anymore
goto 777
xixl = xixstore(i,j,k)
xiyl = xiystore(i,j,k)
@@ -2057,7 +1567,7 @@
+ xizl*(etaxl*gammayl-etayl*gammaxl))
volume_local = volume_local + dble(jacobianl)*weight
-!! DK DK for merged version the jacobian is not stored anymore and therefore not valid anymore
+!! DK DK in merged version the jacobian is not stored anymore and therefore not valid anymore
777 continue
enddo
@@ -2088,190 +1598,5 @@
endif ! end of test if first or second pass
-! deallocate arrays
- deallocate(rhostore_local,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate rhostore_local in create_regions_mesh ier=",ier
- endif
-
- deallocate(kappavstore_local,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate kappavstore_local in create_regions_mesh ier=",ier
- endif
-
-
- deallocate(c11store,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate c11store in create_regions_mesh ier=",ier
- endif
-
- deallocate(c12store,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate c12store in create_regions_mesh ier=",ier
- endif
-
- deallocate(c13store,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate c13store in create_regions_mesh ier=",ier
- endif
-
- deallocate(c14store,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate c14store in create_regions_mesh ier=",ier
- endif
-
- deallocate(c15store,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate c15store in create_regions_mesh ier=",ier
- endif
-
- deallocate(c16store,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate c16store in create_regions_mesh ier=",ier
- endif
-
- deallocate(c22store,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate c22store in create_regions_mesh ier=",ier
- endif
-
- deallocate(c23store,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate c23store in create_regions_mesh ier=",ier
- endif
-
- deallocate(c24store,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate c24store in create_regions_mesh ier=",ier
- endif
-
- deallocate(c25store,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate c25store in create_regions_mesh ier=",ier
- endif
-
- deallocate(c26store,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate c26store in create_regions_mesh ier=",ier
- endif
-
- deallocate(c33store,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate c33store in create_regions_mesh ier=",ier
- endif
-
- deallocate(c34store,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate c34store in create_regions_mesh ier=",ier
- endif
-
- deallocate(c35store,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate c35store in create_regions_mesh ier=",ier
- endif
-
- deallocate(c36store,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate c36store in create_regions_mesh ier=",ier
- endif
-
- deallocate(c44store,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate c44store in create_regions_mesh ier=",ier
- endif
-
- deallocate(c45store,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate c45store in create_regions_mesh ier=",ier
- endif
-
- deallocate(c46store,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate c46store in create_regions_mesh ier=",ier
- endif
-
- deallocate(c55store,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate c55store in create_regions_mesh ier=",ier
- endif
-
- deallocate(c56store,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate c56store in create_regions_mesh ier=",ier
- endif
-
- deallocate(c66store,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate c66store in create_regions_mesh ier=",ier
- endif
-
-
- deallocate(iboun,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate iboun in create_regions_mesh ier=",ier
- endif
-
- deallocate(xigll,yigll,zigll,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate xigll,yigll,zigll in create_regions_mesh ier=",ier
- endif
-
- deallocate(wxgll,wygll,wzgll,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate wxgll,wygll,wzgll in create_regions_mesh ier=",ier
- endif
-
- deallocate(shape3D,dershape3D,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate shape3D,dershape3D in create_regions_mesh ier=",ier
- endif
-
- deallocate(shape2D_x,shape2D_y,shape2D_bottom,shape2D_top,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate shape2D_x,shape2D_y,shape2D_bottom,shape2D_top in create_regions_mesh ier=",ier
- endif
-
- deallocate(dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate dershape2D_x,dershape2D_y,dershape2D_bottom,dershape2D_top" &
- ," in create_regions_mesh ier=",ier
- endif
-
- deallocate(iMPIcut_xi,iMPIcut_eta,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate iMPIcut_xi,iMPIcut_eta in create_regions_mesh ier=",ier
- endif
-
-
- deallocate(nimin,nimax,njmin,njmax,nkmin_xi,nkmin_eta,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate nimin,nimax,njmin,njmax,nkmin_xi,nkmin_eta ",&
- "in create_regions_mesh ier=",ier
- endif
-
- deallocate(rho_vp,rho_vs,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate rho_vp,rho_vs in create_regions_mesh ier=",ier
- endif
-
-
- deallocate(Qmu_store,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate Qmu_store in create_regions_mesh ier=",ier
- endif
-
- deallocate(tau_e_store,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate tau_e_store in create_regions_mesh ier=",ier
- endif
-
- if (allocated(rmass_ocean_load) ) then
- deallocate(rmass_ocean_load,STAT=ier )
- if (ier /= 0) then
- print *,"ERROR can not deallocate rmass_ocean_load in create_regions_mesh ier=",ier
- endif
- endif
-
-
end subroutine create_regions_mesh
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/declarations_mesher.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/declarations_mesher.f90 2008-08-07 23:32:02 UTC (rev 12587)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/declarations_mesher.f90 2008-08-08 01:55:06 UTC (rev 12588)
@@ -1,4 +1,66 @@
+! allocate these automatic arrays in the memory stack to avoid memory fragmentation with "allocate()"
+
+! Gauss-Lobatto-Legendre points of integration and weights
+ double precision, dimension(NGLLX) :: xigll,wxgll
+ double precision, dimension(NGLLY) :: yigll,wygll
+ double precision, dimension(NGLLZ) :: zigll,wzgll
+
+! 3D shape functions and their derivatives
+ double precision, dimension(NGNOD,NGLLX,NGLLY,NGLLZ) :: shape3D
+ double precision, dimension(NDIM,NGNOD,NGLLX,NGLLY,NGLLZ) :: dershape3D
+
+! 2D shape functions and their derivatives
+ double precision shape2D_x(NGNOD2D,NGLLY,NGLLZ)
+ double precision shape2D_y(NGNOD2D,NGLLX,NGLLZ)
+ double precision shape2D_bottom(NGNOD2D,NGLLX,NGLLY)
+ double precision shape2D_top(NGNOD2D,NGLLX,NGLLY)
+
+ double precision dershape2D_x(NDIM2D,NGNOD2D,NGLLY,NGLLZ)
+ double precision dershape2D_y(NDIM2D,NGNOD2D,NGLLX,NGLLZ)
+ double precision dershape2D_bottom(NDIM2D,NGNOD2D,NGLLX,NGLLY)
+ double precision dershape2D_top(NDIM2D,NGNOD2D,NGLLX,NGLLY)
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: rhostore_local,kappavstore_local
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_ANISO_MANTLE) :: &
+ c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
+ c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
+ c36store,c44store,c45store,c46store,c55store,c56store,c66store
+
+ logical, dimension(6,NSPEC_CRUST_MANTLE) :: iboun
+
+! very large arrays used for the sorting routine
+ integer, dimension(NSPEC_CRUST_MANTLE * NGLLX * NGLLY * NGLLZ) :: locval
+ logical, dimension(NSPEC_CRUST_MANTLE * NGLLX * NGLLY * NGLLZ) :: ifseg
+ double precision, dimension(NSPEC_CRUST_MANTLE * NGLLX * NGLLY * NGLLZ) :: xp,yp,zp
+
+ real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE_OCEANS) :: rmass_ocean_load
+
+ integer, dimension(NGLOB_CRUST_MANTLE) :: mask_ibool
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: copy_ibool_ori
+
+ logical, dimension(2,NSPEC_CRUST_MANTLE) :: iMPIcut_xi,iMPIcut_eta
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STACEY) :: rho_vp,rho_vs
+
+ double precision, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: Qmu_store
+ double precision, dimension(N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: tau_e_store
+
+ integer ibelm_moho_top(NSPEC2D_MOHO)
+ integer ibelm_moho_bot(NSPEC2D_MOHO)
+ integer ibelm_400_top(NSPEC2D_400)
+ integer ibelm_400_bot(NSPEC2D_400)
+ integer ibelm_670_top(NSPEC2D_670)
+ integer ibelm_670_bot(NSPEC2D_670)
+
+ real(kind=CUSTOM_REAL) normal_moho(NDIM,NGLLX,NGLLY,NSPEC2D_MOHO)
+ real(kind=CUSTOM_REAL) normal_400(NDIM,NGLLX,NGLLY,NSPEC2D_400)
+ real(kind=CUSTOM_REAL) normal_670(NDIM,NGLLX,NGLLY,NSPEC2D_670)
+ real(kind=CUSTOM_REAL) jacobian2D_moho(NGLLX,NGLLY,NSPEC2D_MOHO)
+ real(kind=CUSTOM_REAL) jacobian2D_400(NGLLX,NGLLY,NSPEC2D_400)
+ real(kind=CUSTOM_REAL) jacobian2D_670(NGLLX,NGLLY,NSPEC2D_670)
+
!!!!!!!!!!!!!!!! DK DK for merged version, all the arrays below are allocated statically instead
!!!!!!!!!!!!!!!! DK DK for merged version, all the arrays below are allocated statically instead
!!!!!!!!!!!!!!!! DK DK for merged version, all the arrays below are allocated statically instead
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_jacobian_discontinuities.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_jacobian_discontinuities.f90 2008-08-07 23:32:02 UTC (rev 12587)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/get_jacobian_discontinuities.f90 2008-08-08 01:55:06 UTC (rev 12588)
@@ -141,8 +141,6 @@
endif
map_ix(ispec_list(i)) = (map_isub_ix(isub_block) - 1) * 2 + irem_ix
map_iy(ispec_list(i)) = (map_isub_iy(isub_block) - 1) * 2 + irem_iy
-! if (ispec_superbrick == 1 .and. myrank == 0) &
-! write(*,'(10i4)') i, ispec_list(i), map_ix(ispec_list(i)), map_iy(ispec_list(i))
enddo
endif
@@ -197,11 +195,9 @@
ix_top = (ix_elem - 1) + ix
iy_top = (iy_elem - 1) + iy
ispec2D_moho_bot_map = (ix_top - 1) * nex_eta_moho + iy_top
-! if (myrank == 0) write(*,'(10i6)') ix_elem, iy_elem, ispec_superbrick, ix, iy, ix_top, iy_top, ispec2D_moho_bot_map
ibelm_moho_bot(ispec2D_moho_bot_map) = ispec
endif
endif
-
end subroutine get_jacobian_discontinuities
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/meshfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/meshfem3D.f90 2008-08-07 23:32:02 UTC (rev 12587)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/meshfem3D.f90 2008-08-08 01:55:06 UTC (rev 12588)
@@ -1217,7 +1217,15 @@
kappavstore_crust_mantle,kappahstore_crust_mantle,muvstore_crust_mantle,muhstore_crust_mantle,eta_anisostore_crust_mantle, &
rmass_crust_mantle,xelm_store_crust_mantle,yelm_store_crust_mantle,zelm_store_crust_mantle, &
!! DK DK this will have to change to fully support David's code to cut the superbrick
- npoin2D_xi_crust_mantle(1),npoin2D_eta_crust_mantle(1))
+ npoin2D_xi_crust_mantle(1),npoin2D_eta_crust_mantle(1), &
+ xigll,wxgll, yigll,wygll, zigll,wzgll, shape3D, dershape3D, shape2D_x, shape2D_y, shape2D_bottom, shape2D_top, &
+ dershape2D_x, dershape2D_y, dershape2D_bottom, dershape2D_top, rhostore_local,kappavstore_local, &
+ c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
+ c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
+ c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
+ iboun, locval, ifseg, xp,yp,zp, rmass_ocean_load, mask_ibool, copy_ibool_ori, iMPIcut_xi,iMPIcut_eta, &
+ rho_vp,rho_vs, Qmu_store, tau_e_store, ibelm_moho_top, ibelm_moho_bot, ibelm_400_top, ibelm_400_bot, &
+ ibelm_670_top, ibelm_670_bot, normal_moho, normal_400, normal_670, jacobian2D_moho, jacobian2D_400, jacobian2D_670)
else if(iregion_code == IREGION_OUTER_CORE) then
! outer_core
@@ -1255,7 +1263,15 @@
kappavstore_outer_core,kappahstore_outer_core,muvstore_outer_core,muhstore_outer_core,eta_anisostore_outer_core, &
rmass_outer_core,xelm_store_outer_core,yelm_store_outer_core,zelm_store_outer_core, &
!! DK DK this will have to change to fully support David's code to cut the superbrick
- npoin2D_xi_outer_core(1),npoin2D_eta_outer_core(1))
+ npoin2D_xi_outer_core(1),npoin2D_eta_outer_core(1), &
+ xigll,wxgll, yigll,wygll, zigll,wzgll, shape3D, dershape3D, shape2D_x, shape2D_y, shape2D_bottom, shape2D_top, &
+ dershape2D_x, dershape2D_y, dershape2D_bottom, dershape2D_top, rhostore_local,kappavstore_local, &
+ c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
+ c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
+ c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
+ iboun, locval, ifseg, xp,yp,zp, rmass_ocean_load, mask_ibool, copy_ibool_ori, iMPIcut_xi,iMPIcut_eta, &
+ rho_vp,rho_vs, Qmu_store, tau_e_store, ibelm_moho_top, ibelm_moho_bot, ibelm_400_top, ibelm_400_bot, &
+ ibelm_670_top, ibelm_670_bot, normal_moho, normal_400, normal_670, jacobian2D_moho, jacobian2D_400, jacobian2D_670)
else if(iregion_code == IREGION_INNER_CORE) then
! inner_core
@@ -1292,7 +1308,15 @@
kappavstore_inner_core,kappahstore_inner_core,muvstore_inner_core,muhstore_inner_core,eta_anisostore_inner_core, &
rmass_inner_core,xelm_store_inner_core,yelm_store_inner_core,zelm_store_inner_core, &
!! DK DK this will have to change to fully support David's code to cut the superbrick
- npoin2D_xi_inner_core(1),npoin2D_eta_inner_core(1))
+ npoin2D_xi_inner_core(1),npoin2D_eta_inner_core(1), &
+ xigll,wxgll, yigll,wygll, zigll,wzgll, shape3D, dershape3D, shape2D_x, shape2D_y, shape2D_bottom, shape2D_top, &
+ dershape2D_x, dershape2D_y, dershape2D_bottom, dershape2D_top, rhostore_local,kappavstore_local, &
+ c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
+ c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
+ c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
+ iboun, locval, ifseg, xp,yp,zp, rmass_ocean_load, mask_ibool, copy_ibool_ori, iMPIcut_xi,iMPIcut_eta, &
+ rho_vp,rho_vs, Qmu_store, tau_e_store, ibelm_moho_top, ibelm_moho_bot, ibelm_400_top, ibelm_400_bot, &
+ ibelm_670_top, ibelm_670_bot, normal_moho, normal_400, normal_670, jacobian2D_moho, jacobian2D_400, jacobian2D_670)
else
stop 'DK DK incorrect region in merged code'
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/read_compute_parameters.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/read_compute_parameters.f90 2008-08-07 23:32:02 UTC (rev 12587)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/read_compute_parameters.f90 2008-08-08 01:55:06 UTC (rev 12588)
@@ -752,36 +752,6 @@
call auto_time_stepping(ANGULAR_WIDTH_XI_IN_DEGREES, NEX_MAX, DT)
-!! DK DK suppressed because this routine should not write anything to the screen
-! write(*,*)'##############################################################'
-! write(*,*)
-! write(*,*)' Auto Radial Meshing Code '
-! write(*,*)' Consult read_compute_parameters.f90 and auto_ner.f90 '
-! write(*,*)' This should only be invoked for chunks less than 90 degrees'
-! write(*,*)' and for chunks greater than 1248 elements wide'
-! write(*,*)
-! write(*,*)'CHUNK WIDTH: ', ANGULAR_WIDTH_XI_IN_DEGREES
-! write(*,*)'NEX: ', NEX_MAX
-! write(*,*)'NER_CRUST: ', NER_CRUST
-! write(*,*)'NER_80_MOHO: ', NER_80_MOHO
-! write(*,*)'NER_220_80: ', NER_220_80
-! write(*,*)'NER_400_220: ', NER_400_220
-! write(*,*)'NER_600_400: ', NER_600_400
-! write(*,*)'NER_670_600: ', NER_670_600
-! write(*,*)'NER_771_670: ', NER_771_670
-! write(*,*)'NER_TOPDDOUBLEPRIME_771: ', NER_TOPDDOUBLEPRIME_771
-! write(*,*)'NER_CMB_TOPDDOUBLEPRIME: ', NER_CMB_TOPDDOUBLEPRIME
-! write(*,*)'NER_OUTER_CORE: ', NER_OUTER_CORE
-! write(*,*)'NER_TOP_CENTRAL_CUBE_ICB: ', NER_TOP_CENTRAL_CUBE_ICB
-! write(*,*)'R_CENTRAL_CUBE: ', R_CENTRAL_CUBE
-! write(*,*)'multiplication factor: ', multiplication_factor
-! write(*,*)
-! write(*,*)'DT: ',DT
-! write(*,*)'MIN_ATTENUATION_PERIOD ',MIN_ATTENUATION_PERIOD
-! write(*,*)'MAX_ATTENUATION_PERIOD ',MAX_ATTENUATION_PERIOD
-! write(*,*)
-! write(*,*)'##############################################################'
-
if (HONOR_1D_SPHERICAL_MOHO) then
if (.not. ONE_CRUST) then
! case 1D + two crustal layers
@@ -793,7 +763,6 @@
endif
endif
-
! take a 5% safety margin on the maximum stable time step
! which was obtained by trial and error
DT = DT * (1.d0 - 0.05d0)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/recompute_jacobian.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/recompute_jacobian.f90 2008-08-07 23:32:02 UTC (rev 12587)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/recompute_jacobian.f90 2008-08-08 01:55:06 UTC (rev 12588)
@@ -38,11 +38,11 @@
double precision xi,eta,gamma,jacobian
! coordinates of the control points of the surface element
- double precision xelm(NGNOD),yelm(NGNOD),zelm(NGNOD)
+ double precision, dimension(NGNOD) :: xelm,yelm,zelm
! 3D shape functions and their derivatives at receiver
- double precision shape3D(NGNOD)
- double precision dershape3D(NDIM,NGNOD)
+ double precision, dimension(NGNOD) :: shape3D
+ double precision, dimension(NDIM,NGNOD) :: dershape3D
double precision l1xi,l2xi,l3xi
double precision l1eta,l2eta,l3eta
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/s362ani.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/s362ani.f90 2008-08-07 23:32:02 UTC (rev 12587)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/s362ani.f90 2008-08-08 01:55:06 UTC (rev 12588)
@@ -563,10 +563,7 @@
1.00196657023780,1.0015515913133,1.0012554932754,1.0010368069141, &
1.00087070107920,1.0007415648034 /
- if(kmax > 13)then
- write(*,"(' kmax exceeds the limit in chebyfun')")
- stop
- endif
+ if(kmax > 13) stop 'kmax exceeds the limit in chebyfun'
f(0)=1.0
f(1)=u
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.f90 2008-08-07 23:32:02 UTC (rev 12587)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/version41_beta/src/specfem3D.f90 2008-08-08 01:55:06 UTC (rev 12588)
@@ -490,9 +490,6 @@
ibool_inner_core,NSPEC_INNER_CORE,NGLOB_INNER_CORE, &
xigll,yigll,zigll)
-!print *,yelm_store_inner_core
-!print *,etay_inner_core
-
!! DK DK for merged version, deallocate arrays that have become useless
deallocate(xelm_store_crust_mantle)
deallocate(yelm_store_crust_mantle)
More information about the cig-commits
mailing list