[cig-commits] r19881 - in seismo/3D/SPECFEM3D_GLOBE/trunk: setup src/create_header_file src/meshfem3D src/shared src/specfem3D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Tue Mar 27 06:15:30 PDT 2012
Author: dkomati1
Date: 2012-03-27 06:15:30 -0700 (Tue, 27 Mar 2012)
New Revision: 19881
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in
seismo/3D/SPECFEM3D_GLOBE/trunk/src/create_header_file/create_header_file.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/meshfem3D.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/memory_eval.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/save_header_file.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_Dev.F90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_attenuation.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_movie_volume.f90
Log:
Removed some large arrays that were useless. Adjusted the size of other large arrays when some options are off (MOVIE_VOLUME, TOPOGRAPHY etc).
Also fixed the memory size estimate in memory_eval.f90.
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in 2012-03-27 05:48:15 UTC (rev 19880)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in 2012-03-27 13:15:30 UTC (rev 19881)
@@ -456,16 +456,6 @@
integer, parameter :: THREE_D_MODEL_S40RTS = 11
integer, parameter :: THREE_D_MODEL_GAPP2 = 12
-! define flag for regions of the global Earth for attenuation
- integer, parameter :: NUM_REGIONS_ATTENUATION = 5
-
- integer, parameter :: IREGION_ATTENUATION_INNER_CORE = 1
- integer, parameter :: IREGION_ATTENUATION_CMB_670 = 2
- integer, parameter :: IREGION_ATTENUATION_670_220 = 3
- integer, parameter :: IREGION_ATTENUATION_220_80 = 4
- integer, parameter :: IREGION_ATTENUATION_80_SURFACE = 5
- integer, parameter :: IREGION_ATTENUATION_UNDEFINED = 6
-
! number of standard linear solids for attenuation
integer, parameter :: N_SLS = 3
@@ -475,10 +465,6 @@
integer, parameter :: ATTENUATION_COMP_RESOLUTION = 1
integer, parameter :: ATTENUATION_COMP_MAXIMUM = 5000
-! for lookup table for attenuation every 100 m in radial direction of Earth model
- integer, parameter :: NRAD_ATTENUATION = 70000
- double precision, parameter :: TABLE_ATTENUATION = R_EARTH_KM * 10.0d0
-
! for determination of the attenuation period range
! if this is set to .true. then the hardcoded values will be used
! otherwise they are computed automatically from the Number of elements
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/create_header_file/create_header_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/create_header_file/create_header_file.f90 2012-03-27 05:48:15 UTC (rev 19880)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/create_header_file/create_header_file.f90 2012-03-27 13:15:30 UTC (rev 19881)
@@ -177,7 +177,7 @@
! create include file for the solver
call save_header_file(NSPEC,nglob,NEX_XI,NEX_ETA,NPROC,NPROCTOT, &
TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
- ELLIPTICITY,GRAVITY,ROTATION,OCEANS,ATTENUATION,ATTENUATION_3D, &
+ ELLIPTICITY,GRAVITY,ROTATION,TOPOGRAPHY,OCEANS,ATTENUATION,ATTENUATION_3D, &
ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,NCHUNKS, &
INCLUDE_CENTRAL_CUBE,CENTER_LONGITUDE_IN_DEGREES, &
CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,NSOURCES,NSTEP,&
@@ -228,21 +228,29 @@
print *,'approximate static memory needed by the solver:'
print *,'----------------------------------------------'
print *
- print *,'size of static arrays per slice = ',static_memory_size/1073741824.d0,' GB'
+ print *,'(lower bound, usually the real amount used is 5% to 10% higher)'
print *
-! note: using less memory becomes only an issue if the code scaling is bad.
-! most users will run simulations with an executable using far less than 80% RAM per core
-! since they prefer having a faster computational time (and use a higher number of cores).
-!! DK DK reply to note: yes, but that is not "green computing" at all!
-!! DK DK Thus one day people will probably need to worry about not wasting resources...
-! print *,' (should be below and typically equal to 80% or 90% of the memory installed per core)'
+ print *,'size of static arrays per slice = ',static_memory_size/1.d6,' MB'
+ print *,' = ',static_memory_size/1048576.d0,' MiB'
+ print *,' = ',static_memory_size/1.d9,' GB'
+ print *,' = ',static_memory_size/1073741824.d0,' GiB'
+ print *
+
+! note: using less memory becomes an issue only if the strong scaling of the code is poor.
+! Some users will run simulations with an executable using far less than 80% RAM per core
+! if they prefer having a faster computational time (and use a higher number of cores).
+
print *,' (should be below 80% or 90% of the memory installed per core)'
print *,' (if significantly more, the job will not run by lack of memory)'
- print *,' (note that if significantly less, you waste a significant amount of memory per processor core)'
- print *,' (but that can be perfectly acceptable if you can afford it and want faster results by using more cores)'
+ print *,' (note that if significantly less, you waste a significant amount'
+ print *,' of memory per processor core)'
+ print *,' (but that can be perfectly acceptable if you can afford it and'
+ print *,' want faster results by using more cores)'
print *
- print *,'size of static arrays for all slices = ',static_memory_size*dble(NPROCTOT)/1073741824.d0,' GB'
- print *,' = ',static_memory_size*dble(NPROCTOT)/1099511627776.d0,' TB'
+ print *,'size of static arrays for all slices = ',static_memory_size*dble(NPROCTOT)/1.d9,' GB'
+ print *,' = ',static_memory_size*dble(NPROCTOT)/1073741824.d0,' GiB'
+ print *,' = ',static_memory_size*dble(NPROCTOT)/1.d12,' TB'
+ print *,' = ',static_memory_size*dble(NPROCTOT)/1099511627776.d0,' TiB'
print *
end program xcreate_header_file
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/meshfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/meshfem3D.f90 2012-03-27 05:48:15 UTC (rev 19880)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/meshfem3D/meshfem3D.f90 2012-03-27 13:15:30 UTC (rev 19881)
@@ -818,7 +818,7 @@
! create include file for the solver
call save_header_file(NSPEC,nglob,NEX_XI,NEX_ETA,NPROC,NPROCTOT, &
TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
- ELLIPTICITY,GRAVITY,ROTATION,OCEANS,ATTENUATION,ATTENUATION_3D, &
+ ELLIPTICITY,GRAVITY,ROTATION,TOPOGRAPHY,OCEANS,ATTENUATION,ATTENUATION_3D, &
ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,NCHUNKS, &
INCLUDE_CENTRAL_CUBE,CENTER_LONGITUDE_IN_DEGREES,&
CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,NSOURCES,NSTEP, &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/memory_eval.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/memory_eval.f90 2012-03-27 05:48:15 UTC (rev 19880)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/memory_eval.f90 2012-03-27 13:15:30 UTC (rev 19881)
@@ -220,10 +220,26 @@
static_memory_size = static_memory_size + 5.d0*dble(N_SLS)*dble(NGLLX)* &
dble(NGLLY)*dble(NGLLZ)*NSPEC_CRUST_MANTLE_ATTENUAT*dble(CUSTOM_REAL)
+! one_minus_sum_beta_crust_mantle, factor_scale_crust_mantle
+ static_memory_size = static_memory_size + 2.d0*dble(NGLLX)* &
+ dble(NGLLY)*dble(NGLLZ)*NSPEC_CRUST_MANTLE_ATTENUAT*dble(CUSTOM_REAL)
+
+! factor_common_crust_mantle
+ static_memory_size = static_memory_size + dble(N_SLS)*dble(NGLLX)* &
+ dble(NGLLY)*dble(NGLLZ)*NSPEC_CRUST_MANTLE_ATTENUAT*dble(CUSTOM_REAL)
+
! R_memory_inner_core
static_memory_size = static_memory_size + 5.d0*dble(N_SLS)*dble(NGLLX)* &
dble(NGLLY)*dble(NGLLZ)*NSPEC_INNER_CORE_ATTENUATION*dble(CUSTOM_REAL)
+! one_minus_sum_beta_inner_core, factor_scale_inner_core
+ static_memory_size = static_memory_size + 2.d0*dble(NGLLX)* &
+ dble(NGLLY)*dble(NGLLZ)*NSPEC_INNER_CORE_ATTENUATION*dble(CUSTOM_REAL)
+
+! factor_common_inner_core
+ static_memory_size = static_memory_size + dble(N_SLS)*dble(NGLLX)* &
+ dble(NGLLY)*dble(NGLLZ)*NSPEC_INNER_CORE_ATTENUATION*dble(CUSTOM_REAL)
+
! xix_crust_mantle,xiy_crust_mantle,xiz_crust_mantle
! etax_crust_mantle,etay_crust_mantle,etaz_crust_mantle,
! gammax_crust_mantle,gammay_crust_mantle,gammaz_crust_mantle
@@ -301,23 +317,19 @@
static_memory_size = static_memory_size + dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC_OUTER_CORE_ROTATION*2.d0*dble(CUSTOM_REAL)
if(ABSORBING_CONDITIONS) then
-
! rho_vp_crust_mantle,rho_vs_crust_mantle
static_memory_size = static_memory_size + dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC(IREGION_CRUST_MANTLE)*2.d0*dble(CUSTOM_REAL)
! vp_outer_core
static_memory_size = static_memory_size + dble(NGLLX)*dble(NGLLY)*dble(NGLLZ)*NSPEC(IREGION_OUTER_CORE)*dble(CUSTOM_REAL)
-
endif
if(OCEANS) then
-
! rmass_ocean_load
static_memory_size = static_memory_size + NGLOB(IREGION_CRUST_MANTLE)*dble(CUSTOM_REAL)
! updated_dof_ocean_load
static_memory_size = static_memory_size + NGLOB(IREGION_CRUST_MANTLE)*dble(SIZE_LOGICAL)
-
endif
! add arrays used to save strain for attenuation or for adjoint runs
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/save_header_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/save_header_file.f90 2012-03-27 05:48:15 UTC (rev 19880)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/shared/save_header_file.f90 2012-03-27 13:15:30 UTC (rev 19881)
@@ -29,7 +29,7 @@
subroutine save_header_file(NSPEC,nglob,NEX_XI,NEX_ETA,NPROC,NPROCTOT, &
TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
- ELLIPTICITY,GRAVITY,ROTATION,OCEANS,ATTENUATION,ATTENUATION_3D, &
+ ELLIPTICITY,GRAVITY,ROTATION,TOPOGRAPHY,OCEANS,ATTENUATION,ATTENUATION_3D, &
ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,NCHUNKS, &
INCLUDE_CENTRAL_CUBE,CENTER_LONGITUDE_IN_DEGREES, &
CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,NSOURCES,NSTEP,&
@@ -60,7 +60,7 @@
integer NEX_XI,NEX_ETA,NPROC,NPROCTOT,NCHUNKS,NSOURCES,NSTEP
logical TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE, &
- ELLIPTICITY,GRAVITY,ROTATION,OCEANS,ATTENUATION,ATTENUATION_3D,INCLUDE_CENTRAL_CUBE
+ ELLIPTICITY,GRAVITY,ROTATION,TOPOGRAPHY,OCEANS,ATTENUATION,ATTENUATION_3D,INCLUDE_CENTRAL_CUBE
double precision ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES, &
CENTER_LONGITUDE_IN_DEGREES,CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH
@@ -259,20 +259,30 @@
write(IOUT,*) '! approximate static memory needed by the solver:'
write(IOUT,*) '! ----------------------------------------------'
write(IOUT,*) '!'
- write(IOUT,*) '! size of static arrays per slice = ',static_memory_size/1073741824.d0,' GB'
+ write(IOUT,*) '! (lower bound, usually the real amount used is 5% to 10% higher)'
write(IOUT,*) '!'
- ! note: using less memory becomes only an issue if the code scaling is bad.
- ! most users will run simulations with an executable using far less than 80% RAM per core
- ! since they prefer having a faster computational time (and use a higher number of cores).
- !write(IOUT,*) '! (should be below and typically equal to 80% or 90%'
- !write(IOUT,*) '! of the memory installed per core)'
- write(IOUT,*) '! (should be below to 80% or 90% of the memory installed per core)'
- write(IOUT,*) '! (if significantly more, the job will not run by lack of memory )'
- !write(IOUT,*) '! (if significantly less, you waste a significant amount of memory)'
+ write(IOUT,*) '! size of static arrays per slice = ',static_memory_size/1.d6,' MB'
+ write(IOUT,*) '! = ',static_memory_size/1048576.d0,' MiB'
+ write(IOUT,*) '! = ',static_memory_size/1.d9,' GB'
+ write(IOUT,*) '! = ',static_memory_size/1073741824.d0,' GiB'
write(IOUT,*) '!'
- write(IOUT,*) '! size of static arrays for all slices = ',static_memory_size*dble(NPROCTOT)/1073741824.d0,' GB'
- write(IOUT,*) '! = ',static_memory_size*dble(NPROCTOT)/1099511627776.d0,' TB'
+
+ ! note: using less memory becomes an issue only if the strong scaling of the code is poor.
+ ! Some users will run simulations with an executable using far less than 80% RAM per core
+ ! if they prefer having a faster computational time (and use a higher number of cores).
+
+ write(IOUT,*) '! (should be below to 80% or 90% of the memory installed per core)'
+ write(IOUT,*) '! (if significantly more, the job will not run by lack of memory )'
+ write(IOUT,*) '! (note that if significantly less, you waste a significant amount'
+ write(IOUT,*) '! of memory per processor core)'
+ write(IOUT,*) '! (but that can be perfectly acceptable if you can afford it and'
+ write(IOUT,*) '! want faster results by using more cores)'
write(IOUT,*) '!'
+ write(IOUT,*) '! size of static arrays for all slices = ',static_memory_size*dble(NPROCTOT)/1.d9,' GB'
+ write(IOUT,*) '! = ',static_memory_size*dble(NPROCTOT)/1073741824.d0,' GiB'
+ write(IOUT,*) '! = ',static_memory_size*dble(NPROCTOT)/1.d12,' TB'
+ write(IOUT,*) '! = ',static_memory_size*dble(NPROCTOT)/1099511627776.d0,' TiB'
+ write(IOUT,*) '!'
write(IOUT,*)
write(IOUT,*) 'integer, parameter :: NEX_XI_VAL = ',NEX_XI
@@ -387,6 +397,15 @@
endif
write(IOUT,*)
+ if(TOPOGRAPHY .or. OCEANS) then
+ write(IOUT,*) 'integer, parameter :: NX_BATHY_VAL = NX_BATHY'
+ write(IOUT,*) 'integer, parameter :: NY_BATHY_VAL = NY_BATHY'
+ else
+ write(IOUT,*) 'integer, parameter :: NX_BATHY_VAL = 1'
+ write(IOUT,*) 'integer, parameter :: NY_BATHY_VAL = 1'
+ endif
+ write(IOUT,*)
+
if(ROTATION) then
write(IOUT,*) 'logical, parameter :: ROTATION_VAL = .true.'
else
@@ -525,8 +544,14 @@
write(IOUT,*) 'logical, parameter :: COMPUTE_AND_STORE_STRAIN = .false.'
endif
+ if (MOVIE_VOLUME) then
+ write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_3DMOVIE = NSPEC_CRUST_MANTLE'
+ write(IOUT,*) 'integer, parameter :: NGLOB_CRUST_MANTLE_3DMOVIE = NGLOB_CRUST_MANTLE'
+ else
+ write(IOUT,*) 'integer, parameter :: NSPEC_CRUST_MANTLE_3DMOVIE = 1'
+ write(IOUT,*) 'integer, parameter :: NGLOB_CRUST_MANTLE_3DMOVIE = 1'
+ endif
-
close(IOUT)
end subroutine save_header_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 2012-03-27 05:48:15 UTC (rev 19880)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/compute_forces_crust_mantle_Dev.F90 2012-03-27 13:15:30 UTC (rev 19881)
@@ -189,8 +189,6 @@
! this for non blocking MPI
integer :: iphase,icall
- integer :: computed_elements
-
logical, dimension(NSPEC_CRUST_MANTLE) :: is_on_a_slice_edge_crust_mantle
real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_INNER_CORE) :: accel_inner_core
@@ -264,9 +262,12 @@
!$OMP displ_crust_mantle,wgll_cube,accel_inner_core,hprime_xxt,hprime_xx,idoubling_inner_core, &
!$OMP addressing,iboolleft_xi_crust_mantle,iboolright_xi_crust_mantle,vx,vy,vz,vnspec, &
!$OMP iboolleft_eta_crust_mantle,iboolright_eta_crust_mantle,npoin2D_faces_crust_mantle, &
-!$OMP npoin2D_xi_crust_mantle,npoin2D_eta_crust_mantle,iboolfaces_crust_mantle,iboolcorner_crust_mantle,iboolleft_xi_inner_core,iboolright_xi_inner_core,ibool_inner_core, &
-!$OMP iboolleft_eta_inner_core,iboolright_eta_inner_core,npoin2D_faces_inner_core,npoin2D_xi_inner_core,npoin2D_eta_inner_core,iboolfaces_inner_core,accel_crust_mantle, &
-!$OMP iboolcorner_inner_core,iprocfrom_faces,iprocto_faces,iproc_master_corners,iproc_worker1_corners,iproc_worker2_corners,buffer_send_faces,buffer_received_faces, &
+!$OMP npoin2D_xi_crust_mantle,npoin2D_eta_crust_mantle,iboolfaces_crust_mantle,iboolcorner_crust_mantle,iboolleft_xi_inner_core, &
+!$OMP iboolright_xi_inner_core,ibool_inner_core, &
+!$OMP iboolleft_eta_inner_core,iboolright_eta_inner_core,npoin2D_faces_inner_core,npoin2D_xi_inner_core,npoin2D_eta_inner_core, &
+!$OMP iboolfaces_inner_core,accel_crust_mantle, &
+!$OMP iboolcorner_inner_core,iprocfrom_faces,iprocto_faces,iproc_master_corners,iproc_worker1_corners,iproc_worker2_corners, &
+!$OMP buffer_send_faces,buffer_received_faces, &
!$OMP buffer_send_chunkcorn_vector,buffer_recv_chunkcorn_vector, &
!$OMP sender_from_slices_to_cube,buffer_all_cube_from_slices,buffer_slices,ibool_central_cube, &
!$OMP ibelm_bottom_inner_core,hprimewgll_xx,hprimewgll_xxt,wgllwgll_xy,wgllwgll_xz, &
@@ -276,9 +277,12 @@
!$OMP npoin2D_max_all_CM_IC ) &
!$OMP PRIVATE(k,j,ispec,fac1,fac2,fac3,sum_terms,iend,ispec_glob, &
!$OMP C1_mxm_m2_m1_5points,A1_mxm_m2_m1_5points,B2_m1_m2_5points,C3_m1_m2_5points, &
-!$OMP B3_m1_m2_5points,C2_mxm_m2_m1_5points,E1_m1_m2_5points,E2_m1_m2_5points,A2_mxm_m2_m1_5points,C3_mxm_m2_m1_5points,A3_mxm_m2_m1_5points,B1_m1_m2_5points, &
-!$OMP C2_m1_m2_5points,C1_m1_m2_5points,E3_m1_m2_5points,E1_mxm_m2_m1_5points,E2_mxm_m2_m1_5points,E3_mxm_m2_m1_5points,tempx1,tempx2,tempx3, &
-!$OMP newtempx1,newtempx2,newtempx3,newtempy1,newtempy2,newtempy3,newtempz1,newtempz2,newtempz3,dummyx_loc,dummyy_loc,dummyz_loc,rho_s_H,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
+!$OMP B3_m1_m2_5points,C2_mxm_m2_m1_5points,E1_m1_m2_5points,E2_m1_m2_5points,A2_mxm_m2_m1_5points,C3_mxm_m2_m1_5points, &
+!$OMP A3_mxm_m2_m1_5points,B1_m1_m2_5points, &
+!$OMP C2_m1_m2_5points,C1_m1_m2_5points,E3_m1_m2_5points,E1_mxm_m2_m1_5points,E2_mxm_m2_m1_5points,E3_mxm_m2_m1_5points, &
+!$OMP tempx1,tempx2,tempx3, &
+!$OMP newtempx1,newtempx2,newtempx3,newtempy1,newtempy2,newtempy3,newtempz1,newtempz2,newtempz3, &
+!$OMP dummyx_loc,dummyy_loc,dummyz_loc,rho_s_H,tempy1,tempy2,tempy3,tempz1,tempz2,tempz3, &
!$OMP iglobv5,iglob1,epsilondev_loc)
do ispec_glob = 1,NSPEC_CRUST_MANTLE,ELEMENTS_NONBLOCKING_CM_IC
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_attenuation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_attenuation.f90 2012-03-27 05:48:15 UTC (rev 19880)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/get_attenuation.f90 2012-03-27 13:15:30 UTC (rev 19881)
@@ -35,8 +35,8 @@
integer myrank, vnspec
character(len=150) prname
- double precision, dimension(NGLLX,NGLLY,NGLLZ,vnspec) :: one_minus_sum_beta, scale_factor
- double precision, dimension(N_SLS,NGLLX,NGLLY,NGLLZ,vnspec) :: factor_common
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,vnspec) :: one_minus_sum_beta, scale_factor
+ real(kind=CUSTOM_REAL), dimension(N_SLS,NGLLX,NGLLY,NGLLZ,vnspec) :: factor_common
double precision, dimension(N_SLS) :: tau_s
integer i,j,k,ispec
@@ -206,516 +206,3 @@
end subroutine get_attenuation_memory_values
-
-!
-!-------------------------------------------------------------------------------------------------
-!
-! not used anymore...
-!
-! subroutine get_attenuation_model_1D(myrank, prname, iregion_code, tau_s, one_minus_sum_beta, &
-! factor_common, scale_factor, vn,vx,vy,vz, AM_V)
-!
-! implicit none
-!
-! include 'mpif.h'
-! include 'constants.h'
-!
-!! model_attenuation_variables
-! type model_attenuation_variables
-! sequence
-! double precision min_period, max_period
-! double precision :: QT_c_source ! Source Frequency
-! double precision, dimension(:), pointer :: Qtau_s ! tau_sigma
-! double precision, dimension(:), pointer :: QrDisc ! Discontinutitues Defined
-! double precision, dimension(:), pointer :: Qr ! Radius
-! integer, dimension(:), pointer :: interval_Q ! Steps
-! double precision, dimension(:), pointer :: Qmu ! Shear Attenuation
-! double precision, dimension(:,:), pointer :: Qtau_e ! tau_epsilon
-! double precision, dimension(:), pointer :: Qomsb, Qomsb2 ! one_minus_sum_beta
-! double precision, dimension(:,:), pointer :: Qfc, Qfc2 ! factor_common
-! double precision, dimension(:), pointer :: Qsf, Qsf2 ! scale_factor
-! integer, dimension(:), pointer :: Qrmin ! Max and Mins of idoubling
-! integer, dimension(:), pointer :: Qrmax ! Max and Mins of idoubling
-! integer :: Qn ! Number of points
-! integer dummy_pad ! padding 4 bytes to align the structure
-! end type model_attenuation_variables
-!
-! type (model_attenuation_variables) AM_V
-!! model_attenuation_variables
-!
-! integer myrank, iregion_code
-! character(len=150) prname
-! integer vn, vx,vy,vz
-! double precision, dimension(N_SLS) :: tau_s
-! double precision, dimension(vx,vy,vz,vn) :: scale_factor, one_minus_sum_beta
-! double precision, dimension(N_SLS, vx,vy,vz,vn) :: factor_common
-!
-! integer i,j,ier,rmax
-! double precision scale_t
-! double precision Qp1, Qpn, radius, fctmp
-! double precision, dimension(:), allocatable :: Qfctmp, Qfc2tmp
-!
-! integer, save :: first_time_called = 1
-!
-! if(myrank == 0 .AND. iregion_code == IREGION_CRUST_MANTLE .AND. first_time_called == 1) then
-! first_time_called = 0
-! open(unit=27, file=prname(1:len_trim(prname))//'1D_Q.bin', status='unknown', form='unformatted')
-! read(27) AM_V%QT_c_source
-! read(27) tau_s
-! read(27) AM_V%Qn
-!
-! allocate(AM_V%Qr(AM_V%Qn))
-! allocate(AM_V%Qmu(AM_V%Qn))
-! allocate(AM_V%Qtau_e(N_SLS,AM_V%Qn))
-!
-! read(27) AM_V%Qr
-! read(27) AM_V%Qmu
-! read(27) AM_V%Qtau_e
-! close(27)
-! endif
-!
-! ! Synch up after the Read
-! call MPI_BCAST(AM_V%QT_c_source,1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
-! call MPI_BCAST(tau_s,N_SLS,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
-! call MPI_BCAST(AM_V%Qn,1,MPI_INTEGER,0,MPI_COMM_WORLD,ier)
-!
-! if(myrank /= 0) then
-! allocate(AM_V%Qr(AM_V%Qn))
-! allocate(AM_V%Qmu(AM_V%Qn))
-! allocate(AM_V%Qtau_e(N_SLS,AM_V%Qn))
-! endif
-!
-! call MPI_BCAST(AM_V%Qr,AM_V%Qn,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
-! call MPI_BCAST(AM_V%Qmu,AM_V%Qn,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
-! call MPI_BCAST(AM_V%Qtau_e,AM_V%Qn*N_SLS,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
-!
-! scale_t = ONE/dsqrt(PI*GRAV*RHOAV)
-!
-! ! Scale the Attenuation Values
-! tau_s(:) = tau_s(:) / scale_t
-! AM_V%Qtau_e(:,:) = AM_V%Qtau_e(:,:) / scale_t
-! AM_V%QT_c_source = 1000.0d0 / AM_V%QT_c_source / scale_t
-! AM_V%Qr(:) = AM_V%Qr(:) / R_EARTH
-!
-! allocate(AM_V%Qsf(AM_V%Qn))
-! allocate(AM_V%Qomsb(AM_V%Qn))
-! allocate(AM_V%Qfc(N_SLS,AM_V%Qn))
-!
-! allocate(AM_V%Qsf2(AM_V%Qn))
-! allocate(AM_V%Qomsb2(AM_V%Qn))
-! allocate(AM_V%Qfc2(N_SLS,AM_V%Qn))
-!
-! allocate(AM_V%interval_Q(AM_V%Qn))
-!
-! allocate(Qfctmp(AM_V%Qn))
-! allocate(Qfc2tmp(AM_V%Qn))
-!
-! do i = 1,AM_V%Qn
-! if(AM_V%Qmu(i) == 0.0d0) then
-! AM_V%Qomsb(i) = 0.0d0
-! AM_V%Qfc(:,i) = 0.0d0
-! AM_V%Qsf(i) = 0.0d0
-! else
-! call attenuation_property_values(tau_s, AM_V%Qtau_e(:,i), AM_V%Qfc(:,i), AM_V%Qomsb(i))
-! call attenuation_scale_factor(myrank, AM_V%QT_c_source, AM_V%Qtau_e(:,i), tau_s, AM_V%Qmu(i), AM_V%Qsf(i))
-! endif
-! enddo
-!
-! ! Determine the Spline Coefficients or Second Derivatives
-! call pspline_construction(AM_V%Qr, AM_V%Qsf, AM_V%Qn, Qp1, Qpn, AM_V%Qsf2, AM_V%interval_Q)
-! call pspline_construction(AM_V%Qr, AM_V%Qomsb, AM_V%Qn, Qp1, Qpn, AM_V%Qomsb2, AM_V%interval_Q)
-! do i = 1,N_SLS
-!! copy the sub-arrays to temporary arrays to avoid a warning by some compilers
-!! about temporary arrays being created automatically when using this expression
-!! directly in the call to the subroutine
-! Qfctmp(:) = AM_V%Qfc(i,:)
-! Qfc2tmp(:) = AM_V%Qfc2(i,:)
-! call pspline_construction(AM_V%Qr, Qfctmp, AM_V%Qn, Qp1, Qpn, Qfc2tmp, AM_V%interval_Q)
-!! copy the arrays back to the sub-arrays, since these sub-arrays are used
-!! as input and output
-! AM_V%Qfc(i,:) = Qfctmp(:)
-! AM_V%Qfc2(i,:) = Qfc2tmp(:)
-! enddo
-!
-! radius = 0.0d0
-! rmax = nint(TABLE_ATTENUATION)
-! do i = 1,rmax
-! call attenuation_lookup_value(i, radius)
-! call pspline_evaluation(AM_V%Qr, AM_V%Qsf, AM_V%Qsf2, AM_V%Qn, radius, scale_factor(1,1,1,i), AM_V%interval_Q)
-! call pspline_evaluation(AM_V%Qr, AM_V%Qomsb, AM_V%Qomsb2, AM_V%Qn, radius, one_minus_sum_beta(1,1,1,i), AM_V%interval_Q)
-! do j = 1,N_SLS
-! Qfctmp = AM_V%Qfc(j,:)
-! Qfc2tmp = AM_V%Qfc2(j,:)
-! call pspline_evaluation(AM_V%Qr, Qfctmp, Qfc2tmp, AM_V%Qn, radius, fctmp, AM_V%interval_Q)
-! factor_common(j,1,1,1,i) = fctmp
-! enddo
-! enddo
-! do i = rmax+1,NRAD_ATTENUATION
-! scale_factor(1,1,1,i) = scale_factor(1,1,1,rmax)
-! one_minus_sum_beta(1,1,1,i) = one_minus_sum_beta(1,1,1,rmax)
-! factor_common(1,1,1,1,i) = factor_common(1,1,1,1,rmax)
-! factor_common(2,1,1,1,i) = factor_common(2,1,1,1,rmax)
-! factor_common(3,1,1,1,i) = factor_common(3,1,1,1,rmax)
-! enddo
-!
-! deallocate(AM_V%Qfc2)
-! deallocate(AM_V%Qsf2)
-! deallocate(AM_V%Qomsb2)
-! deallocate(AM_V%Qfc)
-! deallocate(AM_V%Qsf)
-! deallocate(AM_V%Qomsb)
-! deallocate(AM_V%Qtau_e)
-! deallocate(Qfctmp)
-! deallocate(Qfc2tmp)
-!
-! call MPI_BARRIER(MPI_COMM_WORLD, ier)
-!
-! end subroutine get_attenuation_model_1D
-!
-!
-!-------------------------------------------------------------------------------------------------
-!
-!
-!-------------------------------------------------------------------------------------------------
-!
-! not used anymore...
-!
-! Piecewise Continuous Splines
-! - Added Steps which describes the discontinuities
-! - Steps must be repeats in the dependent variable, X
-! - Derivates at the steps are computed using the point
-! at the derivate and the closest point within that piece
-! - A point lying directly on the discontinuity will recieve the
-! value of the first or smallest piece in terms of X
-! - Beginning and Ending points of the Function become beginning
-! and ending points of the first and last splines
-! - A Step with a value of zero is undefined
-! - Works with functions with steps or no steps
-! See the comment below about the ScS bug
-! subroutine pspline_evaluation(xa, ya, y2a, n, x, y, steps)
-!
-! implicit none
-!
-! integer n
-! double precision xa(n),ya(n),y2a(n)
-! integer steps(n)
-! double precision x, y
-!
-! integer i, l, n1, n2
-!
-! do i = 1,n-1,1
-! if(steps(i+1) == 0) return
-! if(x >= xa(steps(i)) .and. x <= xa(steps(i+1))) then
-! call pspline_piece(i,n1,n2,l,n,steps)
-! call spline_evaluation(xa(n1), ya(n1), y2a(n1), l, x, y)
-!! return <-- Commented out to fix ScS bug
-! endif
-! enddo
-!
-! end subroutine pspline_evaluation
-!
-!
-!-------------------------------------------------------------------------------------------------
-!
-! not used anymore...
-!
-! subroutine pspline_piece(i,n1,n2,l,n,s)
-!
-! implicit none
-!
-! integer i, n1, n2, l, n, s(n)
-! n1 = s(i)+1
-! if(i == 1) n1 = s(i)
-! n2 = s(i+1)
-! l = n2 - n1 + 1
-!
-! end subroutine pspline_piece
-!
-!
-!-------------------------------------------------------------------------------------------------
-!
-! not used anymore...
-!
-! subroutine pspline_construction(x, y, n, yp1, ypn, y2, steps)
-!
-! implicit none
-!
-! integer n
-! double precision x(n),y(n),y2(n)
-! double precision yp1, ypn
-! integer steps(n)
-!
-! integer i,r, l, n1,n2
-!
-! steps(:) = 0
-!
-! ! Find steps in x, defining pieces
-! steps(1) = 1
-! r = 2
-! do i = 2,n
-! if(x(i) == x(i-1)) then
-! steps(r) = i-1
-! r = r + 1
-! endif
-! end do
-! steps(r) = n
-!
-! ! Run spline for each piece
-! do i = 1,r-1
-! call pspline_piece(i,n1,n2,l,n,steps)
-! ! Determine the First Derivates at Begin/End Points
-! yp1 = ( y(n1+1) - y(n1) ) / ( x(n1+1) - x(n1))
-! ypn = ( y(n2) - y(n2-1) ) / ( x(n2) - x(n2-1))
-! call spline_construction(x(n1),y(n1),l,yp1,ypn,y2(n1))
-! enddo
-!
-! end subroutine pspline_construction
-!
-!
-!-------------------------------------------------------------------------------------------------
-!
-!
-! not used anymore...
-!
-! subroutine attenuation_lookup_value(i, r)
-!
-! implicit none
-!
-! include 'constants.h'
-!
-! integer i
-! double precision r
-!
-! r = dble(i) / TABLE_ATTENUATION
-!
-! end subroutine attenuation_lookup_value
-!
-!
-!-------------------------------------------------------------------------------------------------
-!
-! not used anymore...
-!
-! subroutine attenuation_save_arrays(prname, iregion_code, AM_V)
-!
-! implicit none
-!
-! include 'mpif.h'
-! include 'constants.h'
-!
-!! model_attenuation_variables
-! type model_attenuation_variables
-! sequence
-! double precision min_period, max_period
-! double precision :: QT_c_source ! Source Frequency
-! double precision, dimension(:), pointer :: Qtau_s ! tau_sigma
-! double precision, dimension(:), pointer :: QrDisc ! Discontinutitues Defined
-! double precision, dimension(:), pointer :: Qr ! Radius
-! integer, dimension(:), pointer :: interval_Q ! Steps
-! double precision, dimension(:), pointer :: Qmu ! Shear Attenuation
-! double precision, dimension(:,:), pointer :: Qtau_e ! tau_epsilon
-! double precision, dimension(:), pointer :: Qomsb, Qomsb2 ! one_minus_sum_beta
-! double precision, dimension(:,:), pointer :: Qfc, Qfc2 ! factor_common
-! double precision, dimension(:), pointer :: Qsf, Qsf2 ! scale_factor
-! integer, dimension(:), pointer :: Qrmin ! Max and Mins of idoubling
-! integer, dimension(:), pointer :: Qrmax ! Max and Mins of idoubling
-! integer :: Qn ! Number of points
-! end type model_attenuation_variables
-!
-! type (model_attenuation_variables) AM_V
-!! model_attenuation_variables
-!
-! integer iregion_code
-! character(len=150) prname
-! integer ier
-! integer myrank
-! integer, save :: first_time_called = 1
-!
-! call MPI_COMM_RANK(MPI_COMM_WORLD, myrank, ier)
-! if(myrank == 0 .AND. iregion_code == IREGION_CRUST_MANTLE .AND. first_time_called == 1) then
-! first_time_called = 0
-! open(unit=27,file=prname(1:len_trim(prname))//'1D_Q.bin',status='unknown',form='unformatted')
-! write(27) AM_V%QT_c_source
-! write(27) AM_V%Qtau_s
-! write(27) AM_V%Qn
-! write(27) AM_V%Qr
-! write(27) AM_V%Qmu
-! write(27) AM_V%Qtau_e
-! close(27)
-! endif
-!
-! end subroutine attenuation_save_arrays
-!
-!
-!-------------------------------------------------------------------------------------------------
-!
-! not used anymore...
-!
-! subroutine get_attenuation_index(iflag, radius, index, inner_core, AM_V)
-!
-! implicit none
-!
-! include 'constants.h'
-!
-!! model_attenuation_variables
-! type model_attenuation_variables
-! sequence
-! double precision min_period, max_period
-! double precision :: QT_c_source ! Source Frequency
-! double precision, dimension(:), pointer :: Qtau_s ! tau_sigma
-! double precision, dimension(:), pointer :: QrDisc ! Discontinutitues Defined
-! double precision, dimension(:), pointer :: Qr ! Radius
-! integer, dimension(:), pointer :: interval_Q ! Steps
-! double precision, dimension(:), pointer :: Qmu ! Shear Attenuation
-! double precision, dimension(:,:), pointer :: Qtau_e ! tau_epsilon
-! double precision, dimension(:), pointer :: Qomsb, Qomsb2 ! one_minus_sum_beta
-! double precision, dimension(:,:), pointer :: Qfc, Qfc2 ! factor_common
-! double precision, dimension(:), pointer :: Qsf, Qsf2 ! scale_factor
-! integer, dimension(:), pointer :: Qrmin ! Max and Mins of idoubling
-! integer, dimension(:), pointer :: Qrmax ! Max and Mins of idoubling
-! integer :: Qn ! Number of points
-! end type model_attenuation_variables
-!
-! type (model_attenuation_variables) AM_V
-!! model_attenuation_variables
-!
-! integer iflag, iregion, index
-! double precision radius
-!
-! ! Inner Core or not
-! logical inner_core
-!
-! index = nint(radius * TABLE_ATTENUATION)
-!
-!!! DK DK this seems incorrect and is difficult to read anyway
-!!! DK DK therefore let me rewrite it better
-!! if(inner_core) then
-!! if(iflag >= IFLAG_INNER_CORE_NORMAL) then
-!! iregion = IREGION_ATTENUATION_INNER_CORE
-!! else if(iflag >= IFLAG_OUTER_CORE_NORMAL) then
-!! iregion = 6
-!! endif
-!! else
-!! if(iflag >= IFLAG_MANTLE_NORMAL) then
-!! iregion = IREGION_ATTENUATION_CMB_670
-!! else if(iflag == IFLAG_670_220) then
-!! iregion = IREGION_ATTENUATION_670_220
-!! else if(iflag <= IFLAG_220_80) then
-!! iregion = IREGION_ATTENUATION_220_80
-!! else
-!! iregion = IREGION_ATTENUATION_80_SURFACE
-!! endif
-!! endif
-! if(inner_core) then
-!
-! if(iflag == IFLAG_INNER_CORE_NORMAL .or. iflag == IFLAG_MIDDLE_CENTRAL_CUBE .or. &
-! iflag == IFLAG_BOTTOM_CENTRAL_CUBE .or. iflag == IFLAG_TOP_CENTRAL_CUBE .or. &
-! iflag == IFLAG_IN_FICTITIOUS_CUBE) then
-! iregion = IREGION_ATTENUATION_INNER_CORE
-! else
-!! this is fictitious for the outer core, which has no Qmu attenuation since it is fluid
-!! iregion = IREGION_ATTENUATION_80_SURFACE + 1
-! iregion = IREGION_ATTENUATION_UNDEFINED
-! endif
-!
-! else
-!
-! if(iflag == IFLAG_MANTLE_NORMAL) then
-! iregion = IREGION_ATTENUATION_CMB_670
-! else if(iflag == IFLAG_670_220) then
-! iregion = IREGION_ATTENUATION_670_220
-! else if(iflag == IFLAG_220_80) then
-! iregion = IREGION_ATTENUATION_220_80
-! else if(iflag == IFLAG_CRUST .or. iflag == IFLAG_80_MOHO) then
-! iregion = IREGION_ATTENUATION_80_SURFACE
-! else
-!! this is fictitious for the outer core, which has no Qmu attenuation since it is fluid
-!! iregion = IREGION_ATTENUATION_80_SURFACE + 1
-! iregion = IREGION_ATTENUATION_UNDEFINED
-! endif
-!
-! endif
-!
-!! Clamp regions
-! if(index < AM_V%Qrmin(iregion)) index = AM_V%Qrmin(iregion)
-! if(index > AM_V%Qrmax(iregion)) index = AM_V%Qrmax(iregion)
-!
-! end subroutine get_attenuation_index
-!
-!
-!-------------------------------------------------------------------------------------------------
-!
-! not used anymore...
-!
-! subroutine set_attenuation_regions_1D(RICB, RCMB, R670, R220, R80, AM_V)
-!
-! implicit none
-!
-! include 'constants.h'
-!
-!! model_attenuation_variables
-! type model_attenuation_variables
-! sequence
-! double precision min_period, max_period
-! double precision :: QT_c_source ! Source Frequency
-! double precision, dimension(:), pointer :: Qtau_s ! tau_sigma
-! double precision, dimension(:), pointer :: QrDisc ! Discontinutitues Defined
-! double precision, dimension(:), pointer :: Qr ! Radius
-! integer, dimension(:), pointer :: interval_Q ! Steps
-! double precision, dimension(:), pointer :: Qmu ! Shear Attenuation
-! double precision, dimension(:,:), pointer :: Qtau_e ! tau_epsilon
-! double precision, dimension(:), pointer :: Qomsb, Qomsb2 ! one_minus_sum_beta
-! double precision, dimension(:,:), pointer :: Qfc, Qfc2 ! factor_common
-! double precision, dimension(:), pointer :: Qsf, Qsf2 ! scale_factor
-! integer, dimension(:), pointer :: Qrmin ! Max and Mins of idoubling
-! integer, dimension(:), pointer :: Qrmax ! Max and Mins of idoubling
-! integer :: Qn ! Number of points
-! end type model_attenuation_variables
-!
-! type (model_attenuation_variables) AM_V
-!! model_attenuation_variables
-!
-! double precision RICB, RCMB, R670, R220, R80
-! integer i
-!
-! allocate(AM_V%Qrmin(6))
-! allocate(AM_V%Qrmax(6))
-! allocate(AM_V%QrDisc(5))
-!
-! AM_V%QrDisc(1) = RICB
-! AM_V%QrDisc(2) = RCMB
-! AM_V%QrDisc(3) = R670
-! AM_V%QrDisc(4) = R220
-! AM_V%QrDisc(5) = R80
-!
-! ! INNER CORE
-! AM_V%Qrmin(IREGION_ATTENUATION_INNER_CORE) = 1 ! Center of the Earth
-! i = nint(RICB / 100.d0) ! === BOUNDARY === INNER CORE / OUTER CORE
-! AM_V%Qrmax(IREGION_ATTENUATION_INNER_CORE) = i - 1 ! Inner Core Boundary (Inner)
-!
-! ! OUTER_CORE
-! AM_V%Qrmin(6) = i ! Inner Core Boundary (Outer)
-! i = nint(RCMB / 100.d0) ! === BOUNDARY === INNER CORE / OUTER CORE
-! AM_V%Qrmax(6) = i - 1
-!
-! ! LOWER MANTLE
-! AM_V%Qrmin(IREGION_ATTENUATION_CMB_670) = i
-! i = nint(R670 / 100.d0) ! === BOUNDARY === 670 km
-! AM_V%Qrmax(IREGION_ATTENUATION_CMB_670) = i - 1
-!
-! ! UPPER MANTLE
-! AM_V%Qrmin(IREGION_ATTENUATION_670_220) = i
-! i = nint(R220 / 100.d0) ! === BOUNDARY === 220 km
-! AM_V%Qrmax(IREGION_ATTENUATION_670_220) = i - 1
-!
-! ! MANTLE ISH LITHOSPHERE
-! AM_V%Qrmin(IREGION_ATTENUATION_220_80) = i
-! i = nint(R80 / 100.d0) ! === BOUNDARY === 80 km
-! AM_V%Qrmax(IREGION_ATTENUATION_220_80) = i - 1
-!
-! ! CRUST ISH LITHOSPHERE
-! AM_V%Qrmin(IREGION_ATTENUATION_80_SURFACE) = i
-! AM_V%Qrmax(IREGION_ATTENUATION_80_SURFACE) = NRAD_ATTENUATION
-!
-! end subroutine set_attenuation_regions_1D
-!
-
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/prepare_timerun.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/prepare_timerun.f90 2012-03-27 05:48:15 UTC (rev 19880)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/prepare_timerun.f90 2012-03-27 13:15:30 UTC (rev 19881)
@@ -560,10 +560,6 @@
character(len=150) LOCAL_PATH
! local parameters
- double precision, dimension(ATT1,ATT2,ATT3,ATT4) :: omsb_crust_mantle_dble, factor_scale_crust_mantle_dble
- double precision, dimension(ATT1,ATT2,ATT3,ATT5) :: omsb_inner_core_dble, factor_scale_inner_core_dble
- double precision, dimension(N_SLS,ATT1,ATT2,ATT3,ATT4) :: factor_common_crust_mantle_dble
- double precision, dimension(N_SLS,ATT1,ATT2,ATT3,ATT5) :: factor_common_inner_core_dble
double precision, dimension(N_SLS) :: alphaval_dble, betaval_dble, gammaval_dble
double precision, dimension(N_SLS) :: tau_sigma_dble
@@ -576,32 +572,14 @@
! CRUST_MANTLE ATTENUATION
call create_name_database(prname, myrank, IREGION_CRUST_MANTLE, LOCAL_PATH)
- call get_attenuation_model_3D(myrank, prname, omsb_crust_mantle_dble, &
- factor_common_crust_mantle_dble,factor_scale_crust_mantle_dble,tau_sigma_dble,NSPEC_CRUST_MANTLE)
+ call get_attenuation_model_3D(myrank, prname, one_minus_sum_beta_crust_mantle, &
+ factor_common_crust_mantle,factor_scale_crust_mantle,tau_sigma_dble,NSPEC_CRUST_MANTLE)
! INNER_CORE ATTENUATION
call create_name_database(prname, myrank, IREGION_INNER_CORE, LOCAL_PATH)
- call get_attenuation_model_3D(myrank, prname, omsb_inner_core_dble, &
- factor_common_inner_core_dble,factor_scale_inner_core_dble,tau_sigma_dble,NSPEC_INNER_CORE)
+ call get_attenuation_model_3D(myrank, prname, one_minus_sum_beta_inner_core, &
+ factor_common_inner_core,factor_scale_inner_core,tau_sigma_dble,NSPEC_INNER_CORE)
- if(CUSTOM_REAL == SIZE_REAL) then
- factor_scale_crust_mantle = sngl(factor_scale_crust_mantle_dble)
- one_minus_sum_beta_crust_mantle = sngl(omsb_crust_mantle_dble)
- factor_common_crust_mantle = sngl(factor_common_crust_mantle_dble)
-
- factor_scale_inner_core = sngl(factor_scale_inner_core_dble)
- one_minus_sum_beta_inner_core = sngl(omsb_inner_core_dble)
- factor_common_inner_core = sngl(factor_common_inner_core_dble)
- else
- factor_scale_crust_mantle = factor_scale_crust_mantle_dble
- one_minus_sum_beta_crust_mantle = omsb_crust_mantle_dble
- factor_common_crust_mantle = factor_common_crust_mantle_dble
-
- factor_scale_inner_core = factor_scale_inner_core_dble
- one_minus_sum_beta_inner_core = omsb_inner_core_dble
- factor_common_inner_core = factor_common_inner_core_dble
- endif
-
! if attenuation is on, shift PREM to right frequency
! rescale mu in PREM to average frequency for attenuation
! the formulas to implement the scaling can be found for instance in
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90 2012-03-27 05:48:15 UTC (rev 19880)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90 2012-03-27 13:15:30 UTC (rev 19881)
@@ -423,18 +423,17 @@
store_val_ux_all,store_val_uy_all,store_val_uz_all
! to save movie volume
- integer :: npoints_3dmovie,nspecel_3dmovie
- integer, dimension(NGLOB_CRUST_MANTLE) :: num_ibool_3dmovie
+ integer :: npoints_3dmovie,nspecel_3dmovie,ispec
double precision :: scalingval
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_OR_ATT) :: muvstore_crust_mantle_3dmovie
+ integer, dimension(NGLOB_CRUST_MANTLE_3DMOVIE) :: num_ibool_3dmovie
+ logical, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_3DMOVIE) :: mask_3dmovie
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_3DMOVIE) :: muvstore_crust_mantle_3dmovie
real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: nu_3dmovie
- logical, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STRAIN_ONLY) :: mask_3dmovie
+ real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_3DMOVIE) :: Iepsilondev_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_3DMOVIE) :: Ieps_trace_over_3_crust_mantle
- real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_OR_ATT) :: Iepsilondev_crust_mantle
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STRAIN_ONLY) :: Ieps_trace_over_3_crust_mantle
-
! use integer array to store values
- integer, dimension(NX_BATHY,NY_BATHY) :: ibathy_topo
+ integer, dimension(NX_BATHY_VAL,NY_BATHY_VAL) :: ibathy_topo
! for crust/oceans coupling
integer, dimension(NSPEC2DMAX_XMIN_XMAX_CM) :: ibelm_xmin_crust_mantle,ibelm_xmax_crust_mantle
@@ -501,43 +500,6 @@
! for conversion from x y z to r theta phi
real(kind=CUSTOM_REAL) rval,thetaval,phival
-! ---- arrays to assemble between chunks
-
-! communication pattern for faces between chunks
- integer, dimension(NUMMSGS_FACES_VAL) :: iprocfrom_faces,iprocto_faces,imsg_type
-
-! communication pattern for corners between chunks
- integer, dimension(NCORNERSCHUNKS_VAL) :: iproc_master_corners,iproc_worker1_corners,iproc_worker2_corners
-
-! indirect addressing for each message for faces and corners of the chunks
-! a given slice can belong to at most one corner and at most two faces
- integer NGLOB2DMAX_XY
- integer, dimension(NGLOB2DMAX_XY_VAL,NUMFACES_SHARED) :: iboolfaces_crust_mantle, &
- iboolfaces_outer_core,iboolfaces_inner_core
-
-! this for non blocking MPI
-
-! buffers for send and receive between faces of the slices and the chunks
-! we use the same buffers to assemble scalars and vectors because vectors are
-! always three times bigger and therefore scalars can use the first part
-! of the vector buffer in memory even if it has an additional index here
- integer :: npoin2D_max_all_CM_IC
- real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: buffer_send_faces,buffer_received_faces, &
- b_buffer_send_faces,b_buffer_received_faces
-
-! for non blocking communications
- logical, dimension(NSPEC_CRUST_MANTLE) :: is_on_a_slice_edge_crust_mantle
- logical, dimension(NSPEC_OUTER_CORE) :: is_on_a_slice_edge_outer_core
- logical, dimension(NSPEC_INNER_CORE) :: is_on_a_slice_edge_inner_core
- logical, dimension(NGLOB_CRUST_MANTLE) :: mask_ibool
- real :: percentage_edge
-
-! assembling phase number for non blocking MPI
-! iphase is for the crust_mantle, outer_core and inner_core regions
-! iphase_CC is for the central cube
- integer :: iphase,iphase_CC,icall
- integer :: b_iphase,b_iphase_CC,b_icall
-
! -------- arrays specific to each region here -----------
! ----------------- crust, mantle and oceans ---------------------
@@ -682,16 +644,55 @@
integer nabs_xmin_oc,nabs_xmax_oc,nabs_ymin_oc,nabs_ymax_oc,nabs_zmin_oc
integer reclen_xmin_crust_mantle, reclen_xmax_crust_mantle, reclen_ymin_crust_mantle, &
- reclen_ymax_crust_mantle, reclen_xmin_outer_core, reclen_xmax_outer_core,&
+ reclen_ymax_crust_mantle, reclen_xmin_outer_core, reclen_xmax_outer_core, &
reclen_ymin_outer_core, reclen_ymax_outer_core, reclen_zmin
- real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_OUTER_CORE) :: vector_accel_outer_core,&
+ real(kind=CUSTOM_REAL), dimension(NDIM,NGLOB_OUTER_CORE_ADJOINT) :: vector_accel_outer_core, &
vector_displ_outer_core, b_vector_displ_outer_core
integer npoin2D_faces_crust_mantle(NUMFACES_SHARED)
integer npoin2D_faces_outer_core(NUMFACES_SHARED)
integer npoin2D_faces_inner_core(NUMFACES_SHARED)
+! ---- arrays to assemble between chunks
+
+! communication pattern for faces between chunks
+ integer, dimension(NUMMSGS_FACES_VAL) :: iprocfrom_faces,iprocto_faces,imsg_type
+
+! communication pattern for corners between chunks
+ integer, dimension(NCORNERSCHUNKS_VAL) :: iproc_master_corners,iproc_worker1_corners,iproc_worker2_corners
+
+! indirect addressing for each message for faces and corners of the chunks
+! a given slice can belong to at most one corner and at most two faces
+ integer NGLOB2DMAX_XY
+ integer, dimension(NGLOB2DMAX_XY_VAL,NUMFACES_SHARED) :: iboolfaces_crust_mantle, &
+ iboolfaces_outer_core,iboolfaces_inner_core
+
+! this for non blocking MPI
+
+! buffers for send and receive between faces of the slices and the chunks
+! we use the same buffers to assemble scalars and vectors because vectors are
+! always three times bigger and therefore scalars can use the first part
+! of the vector buffer in memory even if it has an additional index here
+ integer :: npoin2D_max_all_CM_IC
+ real(kind=CUSTOM_REAL), dimension(:,:,:), allocatable :: buffer_send_faces,buffer_received_faces, &
+ b_buffer_send_faces,b_buffer_received_faces
+
+! for non blocking communications
+ logical, dimension(NSPEC_CRUST_MANTLE) :: is_on_a_slice_edge_crust_mantle
+ logical, dimension(NSPEC_OUTER_CORE) :: is_on_a_slice_edge_outer_core
+ logical, dimension(NSPEC_INNER_CORE) :: is_on_a_slice_edge_inner_core
+ logical, dimension(NGLOB_CRUST_MANTLE) :: mask_ibool
+! added this to save memory, since this logical mask is not used inside the time loop
+ equivalence(mask_ibool,accel_crust_mantle)
+ real :: percentage_edge
+
+! assembling phase number for non blocking MPI
+! iphase is for the crust_mantle, outer_core and inner_core regions
+! iphase_CC is for the central cube
+ integer :: iphase,iphase_CC,icall
+ integer :: b_iphase,b_iphase_CC,b_icall
+
! parameters for the source
integer it
integer, dimension(:), allocatable :: islice_selected_source,ispec_selected_source
@@ -970,7 +971,7 @@
!! DK DK when turning OpenMP on, use this instead:
!! DK DK from http://mpi.deino.net/mpi_functions/MPI_Init_thread.html
!! DK DK MPI_THREAD_FUNNELED: the process may be multi-threaded, but only the main thread will make MPI calls
-!! DK DK (all MPI calls are funneled to the main thread).
+!! DK DK (all MPI calls are funneled to the main thread).
#ifdef USE_OPENMP
integer :: iprovided
call MPI_INIT_THREAD(MPI_THREAD_FUNNELED,iprovided,ier)
@@ -2209,32 +2210,16 @@
! integral of strain for adjoint movie volume
if(MOVIE_VOLUME .and. (MOVIE_VOLUME_TYPE == 2 .or. MOVIE_VOLUME_TYPE == 3) ) then
- Iepsilondev_crust_mantle(:,:,:,:,:) = Iepsilondev_crust_mantle(:,:,:,:,:) &
- + deltat*epsilondev_crust_mantle(:,:,:,:,:)
- Ieps_trace_over_3_crust_mantle(:,:,:,:) = Ieps_trace_over_3_crust_mantle(:,:,:,:) &
- + deltat*eps_trace_over_3_crust_mantle(:,:,:,:)
+! do *NOT* use array syntax for that loop, otherwise you will get a compiler error when MOVIE_VOLUME is off
+! because the shape of the arrays will not match (due to some arrays purposely declared with a dummy size of 1)
+ do ispec = 1,NSPEC_CRUST_MANTLE
+ Iepsilondev_crust_mantle(:,:,:,:,ispec) = Iepsilondev_crust_mantle(:,:,:,:,ispec) &
+ + deltat*epsilondev_crust_mantle(:,:,:,:,ispec)
+ Ieps_trace_over_3_crust_mantle(:,:,:,ispec) = Ieps_trace_over_3_crust_mantle(:,:,:,ispec) &
+ + deltat*eps_trace_over_3_crust_mantle(:,:,:,ispec)
+ enddo
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
! in all the slices using an MPI reduction
! and output timestamp file to check that simulation is running fine
@@ -4450,7 +4435,7 @@
absorb_xmax_outer_core, &
absorb_ymin_outer_core, &
absorb_ymax_outer_core, &
- absorb_zmin_outer_core)
+ absorb_zmin_outer_core)
endif
! save/read the surface movie using the same c routine as we do for absorbing boundaries (file ID is 9)
@@ -4574,7 +4559,7 @@
deallocate(iadjsrc,iadjsrc_len)
endif
endif
-
+
! receivers
deallocate(islice_selected_rec, &
ispec_selected_rec, &
@@ -4605,7 +4590,7 @@
endif
deallocate(beta_kl_outer_core)
endif
-
+
! movies
if(MOVIE_SURFACE .or. NOISE_TOMOGRAPHY /= 0 ) then
deallocate(store_val_x, &
@@ -4626,7 +4611,7 @@
if(MOVIE_VOLUME) then
deallocate(nu_3dmovie)
endif
-
+
! noise simulations
if ( NOISE_TOMOGRAPHY /= 0 ) then
deallocate(noise_sourcearray, &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_movie_volume.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_movie_volume.f90 2012-03-27 05:48:15 UTC (rev 19880)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/write_movie_volume.f90 2012-03-27 13:15:30 UTC (rev 19881)
@@ -388,9 +388,10 @@
! input
integer :: myrank,npoints_3dmovie,MOVIE_VOLUME_TYPE,it
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: ibool_crust_mantle
- real(kind=CUSTOM_REAL), dimension(3,NGLOB_CRUST_MANTLE) :: vector_crust_mantle,vector_scaled
+ real(kind=CUSTOM_REAL), dimension(3,NGLOB_CRUST_MANTLE) :: vector_crust_mantle
real(kind=CUSTOM_REAL), dimension(3,3,npoints_3dmovie) :: nu_3dmovie
double precision :: scalingval
+ real(kind=CUSTOM_REAL) :: scalingval_to_use
real(kind=CUSTOM_REAL), dimension(3) :: vector_local,vector_local_new
logical, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: mask_3dmovie
logical :: MOVIE_COARSE
@@ -420,12 +421,13 @@
endif
if(CUSTOM_REAL == SIZE_REAL) then
- vector_scaled = vector_crust_mantle*sngl(scalingval)
+ scalingval_to_use = sngl(scalingval)
else
- vector_scaled = vector_crust_mantle*scalingval
+ scalingval_to_use = scalingval
endif
- ipoints_3dmovie=0
+ ipoints_3dmovie = 0
+
do ispec=1,NSPEC_CRUST_MANTLE
do k=1,NGLLZ,NIT
do j=1,NGLLY,NIT
@@ -433,7 +435,7 @@
if(mask_3dmovie(i,j,k,ispec)) then
ipoints_3dmovie=ipoints_3dmovie+1
iglob = ibool_crust_mantle(i,j,k,ispec)
- vector_local(:) = vector_scaled(:,iglob)
+ vector_local(:) = vector_crust_mantle(:,iglob) * scalingval_to_use
! rotate eps_loc to spherical coordinates
vector_local_new(:) = matmul(nu_3dmovie(:,:,ipoints_3dmovie), vector_local(:))
More information about the CIG-COMMITS
mailing list