[cig-commits] r17180 - in seismo/3D/SPECFEM3D_GLOBE/trunk: . UTILS/Profiles
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Wed Sep 8 13:38:25 PDT 2010
Author: danielpeter
Date: 2010-09-08 13:38:25 -0700 (Wed, 08 Sep 2010)
New Revision: 17180
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/Profiles/write_profile.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_1D.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_2D.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_corners_chunks.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_faces_chunks.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/combine_vol_data.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_crust_mantle.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_crust_mantle_Dev.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_outer_core.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_outer_core_Dev.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in
seismo/3D/SPECFEM3D_GLOBE/trunk/flags.guess
seismo/3D/SPECFEM3D_GLOBE/trunk/get_model_parameters.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/initialize_simulation.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/memory_eval.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D_models.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/model_attenuation.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/model_eucrust.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/model_gll.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/model_jp3d.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/model_ppm.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/model_sea99_s.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/prepare_timerun.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/read_arrays_solver.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/read_mesh_databases.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/save_arrays_solver.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/save_kernels.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/write_seismograms.f90
Log:
adds transverse isotropic kernel calculations; minor bug fixes, e.g. in SAC origin time
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/Profiles/write_profile.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/Profiles/write_profile.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/Profiles/write_profile.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -279,8 +279,6 @@
! attenuation_simplex_variables
type attenuation_simplex_variables
sequence
- integer nf ! nf = Number of Frequencies
- integer nsls ! nsls = Number of Standard Linear Solids
double precision Q ! Q = Desired Value of Attenuation or Q
double precision iQ ! iQ = 1/Q
double precision, dimension(:), pointer :: f
@@ -288,6 +286,8 @@
double precision, dimension(:), pointer :: tau_s
! tau_s = Tau_sigma defined by the frequency range and
! number of standard linear solids
+ integer nf ! nf = Number of Frequencies
+ integer nsls ! nsls = Number of Standard Linear Solids
end type attenuation_simplex_variables
type(attenuation_simplex_variables) AS_V
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_1D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_1D.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_1D.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -60,7 +60,8 @@
NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
NTSTEP_BETWEEN_READ_ADJSRC,NSTEP,NTSTEP_BETWEEN_FRAMES, &
NTSTEP_BETWEEN_OUTPUT_INFO,NUMBER_OF_RUNS, &
- NUMBER_OF_THIS_RUN,NCHUNKS,SIMULATION_TYPE,REFERENCE_1D_MODEL,THREE_D_MODEL,MOVIE_VOLUME_TYPE,MOVIE_START,MOVIE_STOP,NOISE_TOMOGRAPHY
+ NUMBER_OF_THIS_RUN,NCHUNKS,SIMULATION_TYPE,REFERENCE_1D_MODEL, &
+ THREE_D_MODEL,MOVIE_VOLUME_TYPE,MOVIE_START,MOVIE_STOP,NOISE_TOMOGRAPHY
double precision DT,ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,CENTER_LONGITUDE_IN_DEGREES, &
CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,ROCEAN,RMIDDLE_CRUST, &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_2D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_2D.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_2D.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -66,7 +66,7 @@
NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
NTSTEP_BETWEEN_READ_ADJSRC,NSTEP,NTSTEP_BETWEEN_FRAMES, &
NTSTEP_BETWEEN_OUTPUT_INFO,NUMBER_OF_RUNS,NUMBER_OF_THIS_RUN,NCHUNKS,SIMULATION_TYPE, &
- REFERENCE_1D_MODEL,THREE_D_MODEL,MOVIE_VOLUME_TYPE,MOVIE_START,MOVIE_STOP
+ REFERENCE_1D_MODEL,THREE_D_MODEL,MOVIE_VOLUME_TYPE,MOVIE_START,MOVIE_STOP,NOISE_TOMOGRAPHY
double precision DT,ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,CENTER_LONGITUDE_IN_DEGREES, &
CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,ROCEAN,RMIDDLE_CRUST, &
@@ -152,7 +152,7 @@
OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
ROTATE_SEISMOGRAMS_RT,ratio_divide_central_cube,HONOR_1D_SPHERICAL_MOHO,CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,&
DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA,&
- WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,.false.)
+ WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,.false.,NOISE_TOMOGRAPHY)
! get the base pathname for output files
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_corners_chunks.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_corners_chunks.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_corners_chunks.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -58,7 +58,7 @@
NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
NTSTEP_BETWEEN_READ_ADJSRC,NSTEP,NTSTEP_BETWEEN_FRAMES, &
NTSTEP_BETWEEN_OUTPUT_INFO,NUMBER_OF_RUNS,NUMBER_OF_THIS_RUN,NCHUNKS,SIMULATION_TYPE, &
- REFERENCE_1D_MODEL,THREE_D_MODEL,MOVIE_VOLUME_TYPE,MOVIE_START,MOVIE_STOP
+ REFERENCE_1D_MODEL,THREE_D_MODEL,MOVIE_VOLUME_TYPE,MOVIE_START,MOVIE_STOP,NOISE_TOMOGRAPHY
double precision DT,ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,CENTER_LONGITUDE_IN_DEGREES, &
CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,ROCEAN,RMIDDLE_CRUST, &
@@ -141,7 +141,7 @@
OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
ROTATE_SEISMOGRAMS_RT,ratio_divide_central_cube,HONOR_1D_SPHERICAL_MOHO,CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,&
DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA,&
- WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,.false.)
+ WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,.false.,NOISE_TOMOGRAPHY)
print *
print *,'There are ',NPROCTOT,' slices numbered from 0 to ',NPROCTOT-1
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_faces_chunks.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_faces_chunks.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/check_buffers_faces_chunks.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -62,7 +62,7 @@
NPROC_XI,NPROC_ETA,NTSTEP_BETWEEN_OUTPUT_SEISMOS, &
NTSTEP_BETWEEN_READ_ADJSRC,NSTEP,NTSTEP_BETWEEN_FRAMES, &
NTSTEP_BETWEEN_OUTPUT_INFO,NUMBER_OF_RUNS,NUMBER_OF_THIS_RUN,NCHUNKS,SIMULATION_TYPE, &
- REFERENCE_1D_MODEL,THREE_D_MODEL,MOVIE_VOLUME_TYPE,MOVIE_START,MOVIE_STOP
+ REFERENCE_1D_MODEL,THREE_D_MODEL,MOVIE_VOLUME_TYPE,MOVIE_START,MOVIE_STOP,NOISE_TOMOGRAPHY
double precision DT,ANGULAR_WIDTH_XI_IN_DEGREES,ANGULAR_WIDTH_ETA_IN_DEGREES,CENTER_LONGITUDE_IN_DEGREES, &
CENTER_LATITUDE_IN_DEGREES,GAMMA_ROTATION_AZIMUTH,ROCEAN,RMIDDLE_CRUST, &
@@ -142,7 +142,7 @@
OUTPUT_SEISMOS_ASCII_TEXT,OUTPUT_SEISMOS_SAC_ALPHANUM,OUTPUT_SEISMOS_SAC_BINARY, &
ROTATE_SEISMOGRAMS_RT,ratio_divide_central_cube,HONOR_1D_SPHERICAL_MOHO,CUT_SUPERBRICK_XI,CUT_SUPERBRICK_ETA,&
DIFF_NSPEC1D_RADIAL,DIFF_NSPEC2D_XI,DIFF_NSPEC2D_ETA,&
- WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,.false.)
+ WRITE_SEISMOGRAMS_BY_MASTER,SAVE_ALL_SEISMOS_IN_ONE_FILE,USE_BINARY_FOR_LARGE_FILE,.false.,NOISE_TOMOGRAPHY)
! get the base pathname for output files
call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/combine_vol_data.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/combine_vol_data.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/combine_vol_data.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -59,6 +59,19 @@
integer:: ibool_count(NGLOB_CRUST_MANTLE)
real(kind=CUSTOM_REAL):: ibool_dat(NGLOB_CRUST_MANTLE)
+ ! note:
+ ! if one wants to remove the topography and ellipticity distortion, you would run the mesher again
+ ! but turning the flags: TOPOGRAPHY and ELLIPTICITY to .false.
+ ! then, use those as topo files: proc***_array_dims.txt and proc***_solver_data_2.bin
+ ! of course, this would also work by just turning ELLIPTICITY to .false. so that the CORRECT_ELLIPTICITY below
+ ! becomes unneccessary
+ !
+ ! puts point locations back into a perfectly spherical shape by removing the ellipticity factor;
+ ! useful for plotting spherical cuts at certain depths
+ logical,parameter:: CORRECT_ELLIPTICITY = .false.
+ integer :: nspl
+ double precision :: rspl(NR),espl(NR),espl2(NR)
+ logical,parameter :: ONE_CRUST = .false. ! if you want to correct a model with one layer only in PREM crust
! starts here--------------------------------------------------------------------------------------------------
@@ -143,6 +156,10 @@
print *, ' low resolution ', HIGH_RESOLUTION_MESH
endif
+ ! sets up ellipticity splines in order to remove ellipticity from point coordinates
+ if( CORRECT_ELLIPTICITY ) call make_ellipticity(nspl,rspl,espl,espl2,ONE_CRUST)
+
+
do ir = irs, ire
print *, '----------- Region ', ir, '----------------'
@@ -291,6 +308,9 @@
y = ystore(iglob)
z = zstore(iglob)
+ ! remove ellipticity
+ if( CORRECT_ELLIPTICITY ) call reverse_ellipticity(x,y,z,nspl,rspl,espl,espl2)
+
!dat = data(i,j,k,ispec)
dat = ibool_dat(iglob)
@@ -308,6 +328,10 @@
x = xstore(iglob)
y = ystore(iglob)
z = zstore(iglob)
+
+ ! remove ellipticity
+ if( CORRECT_ELLIPTICITY ) call reverse_ellipticity(x,y,z,nspl,rspl,espl,espl2)
+
dat = data(i,j,k,ispec)
call write_real_fd(pfd,x)
call write_real_fd(pfd,y)
@@ -398,3 +422,578 @@
end program combine_vol_data
+!
+! ------------------------------------------------------------------------------------------------
+!
+
+
+ subroutine reverse_ellipticity(x,y,z,nspl,rspl,espl,espl2)
+
+ implicit none
+
+ include "constants.h"
+
+ real(kind=CUSTOM_REAL) :: x,y,z
+ integer nspl
+ double precision rspl(NR),espl(NR),espl2(NR)
+ double precision x1,y1,z1
+
+ double precision ell
+ double precision r,theta,phi,factor
+ double precision cost,p20
+
+ ! gets spherical coordinates
+ x1 = x
+ y1 = y
+ z1 = z
+ call xyz_2_rthetaphi_dble(x1,y1,z1,r,theta,phi)
+
+ cost=dcos(theta)
+ p20=0.5d0*(3.0d0*cost*cost-1.0d0)
+
+ ! get ellipticity using spline evaluation
+ call spline_evaluation(rspl,espl,espl2,nspl,r,ell)
+
+ factor=ONE-(TWO/3.0d0)*ell*p20
+
+ ! removes ellipticity factor
+ x = x / factor
+ y = y / factor
+ z = z / factor
+
+ end subroutine reverse_ellipticity
+
+!
+! ------------------------------------------------------------------------------------------------
+!
+
+! copy from make_ellipticity.f90 to avoid compiling issues
+
+ subroutine make_ellipticity(nspl,rspl,espl,espl2,ONE_CRUST)
+
+! creates a spline for the ellipticity profile in PREM
+! radius and density are non-dimensional
+
+ implicit none
+
+ include "constants.h"
+
+ integer nspl
+
+ logical ONE_CRUST
+
+! radius of the Earth for gravity calculation
+ double precision, parameter :: R_EARTH_ELLIPTICITY = 6371000.d0
+! radius of the ocean floor for gravity calculation
+ double precision, parameter :: ROCEAN_ELLIPTICITY = 6368000.d0
+
+ double precision rspl(NR),espl(NR),espl2(NR)
+
+ integer i
+ double precision ROCEAN,RMIDDLE_CRUST,RMOHO,R80,R220,R400,R600,R670, &
+ R771,RTOPDDOUBLEPRIME,RCMB,RICB
+ double precision r_icb,r_cmb,r_topddoubleprime,r_771,r_670,r_600
+ double precision r_400,r_220,r_80,r_moho,r_middle_crust,r_ocean,r_0
+ double precision r(NR),rho(NR),epsilonval(NR),eta(NR)
+ double precision radau(NR),z,k(NR),g_a,bom,exponentval,i_rho,i_radau
+ double precision s1(NR),s2(NR),s3(NR)
+ double precision yp1,ypn
+
+! PREM
+ ROCEAN = 6368000.d0
+ RMIDDLE_CRUST = 6356000.d0
+ RMOHO = 6346600.d0
+ R80 = 6291000.d0
+ R220 = 6151000.d0
+ R400 = 5971000.d0
+ R600 = 5771000.d0
+ R670 = 5701000.d0
+ R771 = 5600000.d0
+ RTOPDDOUBLEPRIME = 3630000.d0
+ RCMB = 3480000.d0
+ RICB = 1221000.d0
+
+! non-dimensionalize
+ r_icb = RICB/R_EARTH_ELLIPTICITY
+ r_cmb = RCMB/R_EARTH_ELLIPTICITY
+ r_topddoubleprime = RTOPDDOUBLEPRIME/R_EARTH_ELLIPTICITY
+ r_771 = R771/R_EARTH_ELLIPTICITY
+ r_670 = R670/R_EARTH_ELLIPTICITY
+ r_600 = R600/R_EARTH_ELLIPTICITY
+ r_400 = R400/R_EARTH_ELLIPTICITY
+ r_220 = R220/R_EARTH_ELLIPTICITY
+ r_80 = R80/R_EARTH_ELLIPTICITY
+ r_moho = RMOHO/R_EARTH_ELLIPTICITY
+ r_middle_crust = RMIDDLE_CRUST/R_EARTH_ELLIPTICITY
+ r_ocean = ROCEAN_ELLIPTICITY/R_EARTH_ELLIPTICITY
+ r_0 = 1.d0
+
+ do i=1,163
+ r(i) = r_icb*dble(i-1)/dble(162)
+ enddo
+ do i=164,323
+ r(i) = r_icb+(r_cmb-r_icb)*dble(i-164)/dble(159)
+ enddo
+ do i=324,336
+ r(i) = r_cmb+(r_topddoubleprime-r_cmb)*dble(i-324)/dble(12)
+ enddo
+ do i=337,517
+ r(i) = r_topddoubleprime+(r_771-r_topddoubleprime)*dble(i-337)/dble(180)
+ enddo
+ do i=518,530
+ r(i) = r_771+(r_670-r_771)*dble(i-518)/dble(12)
+ enddo
+ do i=531,540
+ r(i) = r_670+(r_600-r_670)*dble(i-531)/dble(9)
+ enddo
+ do i=541,565
+ r(i) = r_600+(r_400-r_600)*dble(i-541)/dble(24)
+ enddo
+ do i=566,590
+ r(i) = r_400+(r_220-r_400)*dble(i-566)/dble(24)
+ enddo
+ do i=591,609
+ r(i) = r_220+(r_80-r_220)*dble(i-591)/dble(18)
+ enddo
+ do i=610,619
+ r(i) = r_80+(r_moho-r_80)*dble(i-610)/dble(9)
+ enddo
+ do i=620,626
+ r(i) = r_moho+(r_middle_crust-r_moho)*dble(i-620)/dble(6)
+ enddo
+ do i=627,633
+ r(i) = r_middle_crust+(r_ocean-r_middle_crust)*dble(i-627)/dble(6)
+ enddo
+ do i=634,NR
+ r(i) = r_ocean+(r_0-r_ocean)*dble(i-634)/dble(6)
+ enddo
+
+
+! use PREM to get the density profile for ellipticity (fine for other 1D reference models)
+ do i=1,NR
+ call prem_density(r(i),rho(i),ONE_CRUST,RICB,RCMB,RTOPDDOUBLEPRIME, &
+ R600,R670,R220,R771,R400,R80,RMOHO,RMIDDLE_CRUST,ROCEAN)
+ radau(i)=rho(i)*r(i)*r(i)
+ enddo
+
+ eta(1)=0.0d0
+
+ k(1)=0.0d0
+
+ do i=2,NR
+ call intgrl(i_rho,r,1,i,rho,s1,s2,s3)
+ call intgrl(i_radau,r,1,i,radau,s1,s2,s3)
+ z=(2.0d0/3.0d0)*i_radau/(i_rho*r(i)*r(i))
+ eta(i)=(25.0d0/4.0d0)*((1.0d0-(3.0d0/2.0d0)*z)**2.0d0)-1.0d0
+ k(i)=eta(i)/(r(i)**3.0d0)
+ enddo
+
+ g_a=4.0D0*i_rho
+ bom=TWO_PI/(24.0d0*3600.0d0)
+ bom=bom/sqrt(PI*GRAV*RHOAV)
+ epsilonval(NR)=15.0d0*(bom**2.0d0)/(24.0d0*i_rho*(eta(NR)+2.0d0))
+
+ do i=1,NR-1
+ call intgrl(exponentval,r,i,NR,k,s1,s2,s3)
+ epsilonval(i)=epsilonval(NR)*exp(-exponentval)
+ enddo
+
+! get ready to spline epsilonval
+ nspl=1
+ rspl(1)=r(1)
+ espl(1)=epsilonval(1)
+ do i=2,NR
+ if(r(i) /= r(i-1)) then
+ nspl=nspl+1
+ rspl(nspl)=r(i)
+ espl(nspl)=epsilonval(i)
+ endif
+ enddo
+
+! spline epsilonval
+ yp1=0.0d0
+ ypn=(5.0d0/2.0d0)*(bom**2)/g_a-2.0d0*epsilonval(NR)
+ call spline_construction(rspl,espl,nspl,yp1,ypn,espl2)
+
+ end subroutine make_ellipticity
+
+!
+! ------------------------------------------------------------------------------------------------
+!
+
+! copy from model_prem.f90 to avoid compiling issues
+
+ subroutine prem_density(x,rho,ONE_CRUST,RICB,RCMB,RTOPDDOUBLEPRIME, &
+ R600,R670,R220,R771,R400,R80,RMOHO,RMIDDLE_CRUST,ROCEAN)
+
+ implicit none
+
+ include "constants.h"
+
+ double precision x,rho,RICB,RCMB,RTOPDDOUBLEPRIME, &
+ R600,R670,R220,R771,R400,R80,RMOHO,RMIDDLE_CRUST,ROCEAN
+
+ logical ONE_CRUST
+
+ double precision r
+
+ ! compute real physical radius in meters
+ r = x * R_EARTH
+
+ ! calculates density according to radius
+ if(r <= RICB) then
+ rho=13.0885d0-8.8381d0*x*x
+ else if(r > RICB .and. r <= RCMB) then
+ rho=12.5815d0-1.2638d0*x-3.6426d0*x*x-5.5281d0*x*x*x
+ else if(r > RCMB .and. r <= RTOPDDOUBLEPRIME) then
+ rho=7.9565d0-6.4761d0*x+5.5283d0*x*x-3.0807d0*x*x*x
+ else if(r > RTOPDDOUBLEPRIME .and. r <= R771) then
+ rho=7.9565d0-6.4761d0*x+5.5283d0*x*x-3.0807d0*x*x*x
+ else if(r > R771 .and. r <= R670) then
+ rho=7.9565d0-6.4761d0*x+5.5283d0*x*x-3.0807d0*x*x*x
+ else if(r > R670 .and. r <= R600) then
+ rho=5.3197d0-1.4836d0*x
+ else if(r > R600 .and. r <= R400) then
+ rho=11.2494d0-8.0298d0*x
+ else if(r > R400 .and. r <= R220) then
+ rho=7.1089d0-3.8045d0*x
+ else if(r > R220 .and. r <= R80) then
+ rho=2.6910d0+0.6924d0*x
+ else
+ if(r > R80 .and. r <= RMOHO) then
+ rho=2.6910d0+0.6924d0*x
+ else if(r > RMOHO .and. r <= RMIDDLE_CRUST) then
+ if(ONE_CRUST) then
+ rho=2.6d0
+ else
+ rho=2.9d0
+ endif
+ else if(r > RMIDDLE_CRUST .and. r <= ROCEAN) then
+ rho=2.6d0
+ else if(r > ROCEAN) then
+ rho=2.6d0
+ endif
+ endif
+
+ rho=rho*1000.0d0/RHOAV
+
+ end subroutine prem_density
+
+!
+! ------------------------------------------------------------------------------------------------
+!
+
+! copy from intgrl.f90 to avoid compiling issues
+
+
+ subroutine intgrl(sum,r,nir,ner,f,s1,s2,s3)
+
+! Computes the integral of f[i]*r[i]*r[i] from i=nir to i=ner for
+! radii values as in model PREM_an640
+
+ implicit none
+
+! Argument variables
+ integer ner,nir
+ double precision f(640),r(640),s1(640),s2(640)
+ double precision s3(640),sum
+
+! Local variables
+ double precision, parameter :: third = 1.0d0/3.0d0
+ double precision, parameter :: fifth = 1.0d0/5.0d0
+ double precision, parameter :: sixth = 1.0d0/6.0d0
+
+ double precision rji,yprime(640)
+ double precision s1l,s2l,s3l
+
+ integer i,j,n,kdis(28)
+ integer ndis,nir1
+
+
+
+ data kdis/163,323,336,517,530,540,565,590,609,619,626,633,16*0/
+
+ ndis = 12
+ n = 640
+
+ call deriv(f,yprime,n,r,ndis,kdis,s1,s2,s3)
+ nir1 = nir + 1
+ sum = 0.0d0
+ do i=nir1,ner
+ j = i-1
+ rji = r(i) - r(j)
+ s1l = s1(j)
+ s2l = s2(j)
+ s3l = s3(j)
+ sum = sum + r(j)*r(j)*rji*(f(j) &
+ + rji*(0.5d0*s1l + rji*(third*s2l + rji*0.25d0*s3l))) &
+ + 2.0d0*r(j)*rji*rji*(0.5d0*f(j) + rji*(third*s1l + rji*(0.25d0*s2l + rji*fifth*s3l))) &
+ + rji*rji*rji*(third*f(j) + rji*(0.25d0*s1l + rji*(fifth*s2l + rji*sixth*s3l)))
+ enddo
+
+ end subroutine intgrl
+
+! -------------------------------
+
+ subroutine deriv(y,yprime,n,r,ndis,kdis,s1,s2,s3)
+
+ implicit none
+
+! Argument variables
+ integer kdis(28),n,ndis
+ double precision r(n),s1(n),s2(n),s3(n)
+ double precision y(n),yprime(n)
+
+! Local variables
+ integer i,j,j1,j2
+ integer k,nd,ndp
+ double precision a0,b0,b1
+ double precision f(3,1000),h,h2,h2a
+ double precision h2b,h3a,ha,s13
+ double precision s21,s32,yy(3)
+
+ yy(1) = 0.d0
+ yy(2) = 0.d0
+ yy(3) = 0.d0
+
+ ndp=ndis+1
+ do 3 nd=1,ndp
+ if(nd == 1) goto 4
+ if(nd == ndp) goto 5
+ j1=kdis(nd-1)+1
+ j2=kdis(nd)-2
+ goto 6
+ 4 j1=1
+ j2=kdis(1)-2
+ goto 6
+ 5 j1=kdis(ndis)+1
+ j2=n-2
+ 6 if((j2+1-j1)>0) goto 11
+ j2=j2+2
+ yy(1)=(y(j2)-y(j1))/(r(j2)-r(j1))
+ s1(j1)=yy(1)
+ s1(j2)=yy(1)
+ s2(j1)=yy(2)
+ s2(j2)=yy(2)
+ s3(j1)=yy(3)
+ s3(j2)=yy(3)
+ goto 3
+ 11 a0=0.0d0
+ if(j1 == 1) goto 7
+ h=r(j1+1)-r(j1)
+ h2=r(j1+2)-r(j1)
+ yy(1)=h*h2*(h2-h)
+ h=h*h
+ h2=h2*h2
+ b0=(y(j1)*(h-h2)+y(j1+1)*h2-y(j1+2)*h)/yy(1)
+ goto 8
+ 7 b0=0.0d0
+ 8 b1=b0
+
+ if(j2 > 1000) stop 'error in subroutine deriv for j2'
+
+ do i=j1,j2
+ h=r(i+1)-r(i)
+ yy(1)=y(i+1)-y(i)
+ h2=h*h
+ ha=h-a0
+ h2a=h-2.0d0*a0
+ h3a=2.0d0*h-3.0d0*a0
+ h2b=h2*b0
+ s1(i)=h2/ha
+ s2(i)=-ha/(h2a*h2)
+ s3(i)=-h*h2a/h3a
+ f(1,i)=(yy(1)-h*b0)/(h*ha)
+ f(2,i)=(h2b-yy(1)*(2.0d0*h-a0))/(h*h2*h2a)
+ f(3,i)=-(h2b-3.0d0*yy(1)*ha)/(h*h3a)
+ a0=s3(i)
+ b0=f(3,i)
+ enddo
+
+ i=j2+1
+ h=r(i+1)-r(i)
+ yy(1)=y(i+1)-y(i)
+ h2=h*h
+ ha=h-a0
+ h2a=h*ha
+ h2b=h2*b0-yy(1)*(2.d0*h-a0)
+ s1(i)=h2/ha
+ f(1,i)=(yy(1)-h*b0)/h2a
+ ha=r(j2)-r(i+1)
+ yy(1)=-h*ha*(ha+h)
+ ha=ha*ha
+ yy(1)=(y(i+1)*(h2-ha)+y(i)*ha-y(j2)*h2)/yy(1)
+ s3(i)=(yy(1)*h2a+h2b)/(h*h2*(h-2.0d0*a0))
+ s13=s1(i)*s3(i)
+ s2(i)=f(1,i)-s13
+
+ do j=j1,j2
+ k=i-1
+ s32=s3(k)*s2(i)
+ s1(i)=f(3,k)-s32
+ s21=s2(k)*s1(i)
+ s3(k)=f(2,k)-s21
+ s13=s1(k)*s3(k)
+ s2(k)=f(1,k)-s13
+ i=k
+ enddo
+
+ s1(i)=b1
+ j2=j2+2
+ s1(j2)=yy(1)
+ s2(j2)=yy(2)
+ s3(j2)=yy(3)
+ 3 continue
+
+ do i=1,n
+ yprime(i)=s1(i)
+ enddo
+
+ end subroutine deriv
+
+!
+! ------------------------------------------------------------------------------------------------
+!
+
+! copy from spline_routines.f90 to avoid compiling issues
+
+! compute spline coefficients
+
+ subroutine spline_construction(xpoint,ypoint,npoint,tangent_first_point,tangent_last_point,spline_coefficients)
+
+ implicit none
+
+! tangent to the spline imposed at the first and last points
+ double precision, intent(in) :: tangent_first_point,tangent_last_point
+
+! number of input points and coordinates of the input points
+ integer, intent(in) :: npoint
+ double precision, dimension(npoint), intent(in) :: xpoint,ypoint
+
+! spline coefficients output by the routine
+ double precision, dimension(npoint), intent(out) :: spline_coefficients
+
+ integer :: i
+
+ double precision, dimension(:), allocatable :: temporary_array
+
+ allocate(temporary_array(npoint))
+
+ spline_coefficients(1) = - 1.d0 / 2.d0
+
+ temporary_array(1) = (3.d0/(xpoint(2)-xpoint(1)))*((ypoint(2)-ypoint(1))/(xpoint(2)-xpoint(1))-tangent_first_point)
+
+ do i = 2,npoint-1
+
+ spline_coefficients(i) = ((xpoint(i)-xpoint(i-1))/(xpoint(i+1)-xpoint(i-1))-1.d0) &
+ / ((xpoint(i)-xpoint(i-1))/(xpoint(i+1)-xpoint(i-1))*spline_coefficients(i-1)+2.d0)
+
+ temporary_array(i) = (6.d0*((ypoint(i+1)-ypoint(i))/(xpoint(i+1)-xpoint(i)) &
+ - (ypoint(i)-ypoint(i-1))/(xpoint(i)-xpoint(i-1)))/(xpoint(i+1)-xpoint(i-1)) &
+ - (xpoint(i)-xpoint(i-1))/(xpoint(i+1)-xpoint(i-1))*temporary_array(i-1)) &
+ / ((xpoint(i)-xpoint(i-1))/(xpoint(i+1)-xpoint(i-1))*spline_coefficients(i-1)+2.d0)
+
+ enddo
+
+ spline_coefficients(npoint) = ((3.d0/(xpoint(npoint)-xpoint(npoint-1))) &
+ * (tangent_last_point-(ypoint(npoint)-ypoint(npoint-1))/(xpoint(npoint)-xpoint(npoint-1))) &
+ - 1.d0/2.d0*temporary_array(npoint-1))/(1.d0/2.d0*spline_coefficients(npoint-1)+1.d0)
+
+ do i = npoint-1,1,-1
+ spline_coefficients(i) = spline_coefficients(i)*spline_coefficients(i+1) + temporary_array(i)
+ enddo
+
+ deallocate(temporary_array)
+
+ end subroutine spline_construction
+
+! --------------
+
+! evaluate a spline
+
+ subroutine spline_evaluation(xpoint,ypoint,spline_coefficients,npoint,x_evaluate_spline,y_spline_obtained)
+
+ implicit none
+
+! number of input points and coordinates of the input points
+ integer, intent(in) :: npoint
+ double precision, dimension(npoint), intent(in) :: xpoint,ypoint
+
+! spline coefficients to use
+ double precision, dimension(npoint), intent(in) :: spline_coefficients
+
+! abscissa at which we need to evaluate the value of the spline
+ double precision, intent(in):: x_evaluate_spline
+
+! ordinate evaluated by the routine for the spline at this abscissa
+ double precision, intent(out):: y_spline_obtained
+
+ integer :: index_loop,index_lower,index_higher
+
+ double precision :: coef1,coef2
+
+! initialize to the whole interval
+ index_lower = 1
+ index_higher = npoint
+
+! determine the right interval to use, by dichotomy
+ do while (index_higher - index_lower > 1)
+! compute the middle of the interval
+ index_loop = (index_higher + index_lower) / 2
+ if(xpoint(index_loop) > x_evaluate_spline) then
+ index_higher = index_loop
+ else
+ index_lower = index_loop
+ endif
+ enddo
+
+! test that the interval obtained does not have a size of zero
+! (this could happen for instance in the case of duplicates in the input list of points)
+ if(xpoint(index_higher) == xpoint(index_lower)) stop 'incorrect interval found in spline evaluation'
+
+ coef1 = (xpoint(index_higher) - x_evaluate_spline) / (xpoint(index_higher) - xpoint(index_lower))
+ coef2 = (x_evaluate_spline - xpoint(index_lower)) / (xpoint(index_higher) - xpoint(index_lower))
+
+ y_spline_obtained = coef1*ypoint(index_lower) + coef2*ypoint(index_higher) + &
+ ((coef1**3 - coef1)*spline_coefficients(index_lower) + &
+ (coef2**3 - coef2)*spline_coefficients(index_higher))*((xpoint(index_higher) - xpoint(index_lower))**2)/6.d0
+
+ end subroutine spline_evaluation
+
+
+!
+! ------------------------------------------------------------------------------------------------
+!
+
+! copy from rthetaphi_xyz.f90 to avoid compiling issues
+
+
+ subroutine xyz_2_rthetaphi_dble(x,y,z,r,theta,phi)
+
+! convert x y z to r theta phi, double precision call
+
+ implicit none
+
+ include "constants.h"
+
+ double precision x,y,z,r,theta,phi
+ double precision xmesh,ymesh,zmesh
+
+ xmesh = x
+ ymesh = y
+ zmesh = z
+
+ if(zmesh > -SMALL_VAL_ANGLE .and. zmesh <= ZERO) zmesh = -SMALL_VAL_ANGLE
+ if(zmesh < SMALL_VAL_ANGLE .and. zmesh >= ZERO) zmesh = SMALL_VAL_ANGLE
+
+ theta = datan2(dsqrt(xmesh*xmesh+ymesh*ymesh),zmesh)
+
+ if(xmesh > -SMALL_VAL_ANGLE .and. xmesh <= ZERO) xmesh = -SMALL_VAL_ANGLE
+ if(xmesh < SMALL_VAL_ANGLE .and. xmesh >= ZERO) xmesh = SMALL_VAL_ANGLE
+
+ phi = datan2(ymesh,xmesh)
+
+ r = dsqrt(xmesh*xmesh + ymesh*ymesh + zmesh*zmesh)
+
+ end subroutine xyz_2_rthetaphi_dble
+
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_crust_mantle.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_crust_mantle.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_crust_mantle.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -47,24 +47,24 @@
include "OUTPUT_FILES/values_from_mesher.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
+! 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
+! 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, dimension(:), pointer :: interval_Q ! Steps
+! integer :: Qn ! Number of points
+! integer dummy_pad ! padding 4 bytes to align the structure
+! end type model_attenuation_variables
! array with the local to global mapping per slice
integer, dimension(NSPEC_CRUST_MANTLE) :: idoubling
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_crust_mantle_Dev.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_crust_mantle_Dev.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_crust_mantle_Dev.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -453,6 +453,16 @@
else
+! note : mesh is built such that anisotropic elements are created first in anisotropic layers,
+! thus they are listed first ( see in create_regions_mesh.f90: perm_layer() ordering )
+! this is therefore still in bounds of 1:NSPECMAX_TISO_MANTLE even if NSPECMAX_TISO is less than NSPEC
+
+ ! uncomment to debug
+ !if ( ispec > NSPECMAX_TISO_MANTLE ) then
+ ! print*,'error tiso: ispec = ',ispec,'max = ',NSPECMAX_TISO_MANTLE
+ ! call exit_mpi(0,'error tiso ispec bounds')
+ !endif
+
! use Kappa and mu from transversely isotropic model
kappavl = kappavstore(i,j,k,ispec)
muvl = muvstore(i,j,k,ispec)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_outer_core.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_outer_core.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_outer_core.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -46,10 +46,10 @@
include "OUTPUT_FILES/values_from_mesher.h"
! displacement and acceleration
- real(kind=CUSTOM_REAL), dimension(nglob_outer_core) :: displfluid,accelfluid
+ real(kind=CUSTOM_REAL), dimension(NGLOB_OUTER_CORE) :: displfluid,accelfluid
! divergence of displacement
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec_outer_core) :: div_displfluid
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_OUTER_CORE_ADJOINT) :: div_displfluid
! arrays with mesh parameters per slice
integer, dimension(NGLLX,NGLLY,NGLLZ,nspec_outer_core) :: ibool
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_outer_core_Dev.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_outer_core_Dev.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_outer_core_Dev.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -46,10 +46,10 @@
include "OUTPUT_FILES/values_from_mesher.h"
! displacement and acceleration
- real(kind=CUSTOM_REAL), dimension(nglob_outer_core) :: displfluid,accelfluid
+ real(kind=CUSTOM_REAL), dimension(NGLOB_OUTER_CORE) :: displfluid,accelfluid
! divergence of displacement
- real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec_outer_core) :: div_displfluid
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_OUTER_CORE_ADJOINT) :: div_displfluid
! arrays with mesh parameters per slice
integer, dimension(NGLLX,NGLLY,NGLLZ,nspec_outer_core) :: ibool
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in 2010-09-08 20:38:25 UTC (rev 17180)
@@ -186,6 +186,10 @@
! default is .false. to compute isotropic kernels (related to alpha and beta)
logical, parameter :: ANISOTROPIC_KL = .false.
+! output only transverse isotropic kernels (alpha_v,alpha_h,beta_v,beta_h,eta,rho)
+! rather than fully anisotropic kernels in case ANISOTROPIC_KL is set to .true.
+ logical, parameter :: SAVE_TRANSVERSE_KL = .false.
+
! print date and time estimate of end of run in another country,
! in addition to local time.
! For instance: the code runs at Caltech in California but the person
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/flags.guess
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/flags.guess 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/flags.guess 2010-09-08 20:38:25 UTC (rev 17180)
@@ -24,6 +24,7 @@
ifort|*/ifort)
#
# Intel ifort Fortran90 for Linux
+ # check: http://software.intel.com/sites/products/documentation/hpc/compilerpro/en-us/fortran/lin/compiler_f/index.htm
#
if test x"$FLAGS_CHECK" = x; then
FLAGS_CHECK="-O3 -fpe0 -ftz -align sequence -assume byterecl -vec-report0 -std95 -implicitnone -warn truncated_source -warn argument_checking -warn unused -warn declarations -warn alignments -warn ignore_loc -warn usage" # -mcmodel=medium
@@ -69,7 +70,7 @@
if test x"$FLAGS_CHECK" = x; then
FLAGS_CHECK="\$(FLAGS_NO_CHECK)"
fi
- # useful for debugging...
+ # useful for debugging... -ffpe-trap=underflow,overflow,zero -fbacktrace
#if test x"$FLAGS_NO_CHECK" = x; then
# FLAGS_NO_CHECK="-std=f95 -fimplicit-none -frange-check -O3 -fmax-errors=10 -pedantic -pedantic-errors -Waliasing -Wampersand -Wcharacter-truncation -Wline-truncation -Wsurprising -Wno-tabs -Wunderflow -fno-trapping-math" # -mcmodel=medium
#
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/get_model_parameters.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/get_model_parameters.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/get_model_parameters.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -372,11 +372,11 @@
else if(MODEL_ROOT == 'GLL') then
! model will be given on local basis, at all GLL points,
! as from meshfem3d output from routine save_arrays_solver()
+ ! based on model s29ea
CASE_3D = .true.
CRUSTAL = .true.
ISOTROPIC_3D_MANTLE = .true.
ONE_CRUST = .true.
- ! based on model s29ea
REFERENCE_1D_MODEL = REFERENCE_MODEL_1DREF
THREE_D_MODEL = THREE_D_MODEL_GLL
TRANSVERSE_ISOTROPY = .true.
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/initialize_simulation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/initialize_simulation.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/initialize_simulation.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -454,6 +454,17 @@
if (SIMULATION_TYPE == 3 .and. (ANISOTROPIC_3D_MANTLE_VAL .or. ANISOTROPIC_INNER_CORE_VAL)) &
call exit_MPI(myrank, 'anisotropic model is not implemented for kernel simulations yet')
+ ! checks model for transverse isotropic kernel computation
+ if( SAVE_TRANSVERSE_KL ) then
+ if( ANISOTROPIC_3D_MANTLE_VAL ) then
+ call exit_mpi(myrank,'error SAVE_TRANSVERSE_KL: Earth model not supported yet')
+ endif
+ if( SIMULATION_TYPE == 3 ) then
+ if( .not. ANISOTROPIC_KL ) then
+ call exit_mpi(myrank,'error SAVE_TRANSVERSE_KL: needs anisotropic kernel calculations')
+ endif
+ endif
+ endif
! make ellipticity
if(ELLIPTICITY_VAL) call make_ellipticity(nspl,rspl,espl,espl2,ONE_CRUST)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/memory_eval.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/memory_eval.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/memory_eval.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -117,7 +117,11 @@
NSPECMAX_ISO_MANTLE = NSPEC(IREGION_CRUST_MANTLE)
if(TRANSVERSE_ISOTROPY) then
- NSPECMAX_TISO_MANTLE = ispec_aniso
+! note: the number of transverse isotropic elements is ispec_aniso
+! however for transverse isotropic kernels, the arrays muhstore,kappahstore,eta_anisostore,
+! will be needed for the crust_mantle region everywhere still...
+! originally: NSPECMAX_TISO_MANTLE = ispec_aniso
+ NSPECMAX_TISO_MANTLE = NSPEC(IREGION_CRUST_MANTLE)
else
NSPECMAX_TISO_MANTLE = 1
endif
@@ -142,7 +146,8 @@
NSPEC_INNER_CORE_STR_OR_ATT = 1
endif
- if(ATTENUATION .and. SIMULATION_TYPE == 3) then
+ if(ATTENUATION .and. &
+ ( SIMULATION_TYPE == 3 .or. (SIMULATION_TYPE == 1 .and. SAVE_FORWARD)) ) then
NSPEC_CRUST_MANTLE_STR_AND_ATT = NSPEC(IREGION_CRUST_MANTLE)
NSPEC_INNER_CORE_STR_AND_ATT = NSPEC(IREGION_INNER_CORE)
else
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D_models.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D_models.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D_models.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -57,7 +57,6 @@
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
@@ -65,6 +64,7 @@
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, dimension(:), pointer :: interval_Q ! Steps
integer :: Qn ! Number of points
integer dummy_pad ! padding 4 bytes to align the structure
end type model_attenuation_variables
@@ -176,13 +176,7 @@
! model_jp3d_variables
type model_jp3d_variables
sequence
- ! vmod3d
- integer :: NPA
- integer :: NRA
- integer :: NHA
- integer :: NPB
- integer :: NRB
- integer :: NHB
+ ! vmod3d
double precision :: PNA(MPA)
double precision :: RNA(MRA)
double precision :: HNA(MHA)
@@ -191,34 +185,22 @@
double precision :: HNB(MHB)
double precision :: VELAP(MPA,MRA,MHA)
double precision :: VELBP(MPB,MRB,MHB)
- ! discon
+ ! discon
double precision :: PN(51)
double precision :: RRN(63)
double precision :: DEPA(51,63)
double precision :: DEPB(51,63)
double precision :: DEPC(51,63)
- ! locate
- integer :: IPLOCA(MKA)
- integer :: IRLOCA(MKA)
- integer :: IHLOCA(MKA)
- integer :: IPLOCB(MKB)
- integer :: IRLOCB(MKB)
- integer :: IHLOCB(MKB)
+ ! locate
double precision :: PLA
double precision :: RLA
double precision :: HLA
double precision :: PLB
double precision :: RLB
double precision :: HLB
- ! weight
- integer :: IP
- integer :: JP
- integer :: KP
- integer :: IP1
- integer :: JP1
- integer :: KP1
+ ! weight
double precision :: WV(8)
- ! prhfd
+ ! prhfd
double precision :: P
double precision :: R
double precision :: H
@@ -231,11 +213,32 @@
double precision :: PD
double precision :: RD
double precision :: HD
- ! jpmodv
+ ! jpmodv
double precision :: VP(29)
double precision :: VS(29)
double precision :: RA(29)
double precision :: DEPJ(29)
+ ! locate integers
+ integer :: IPLOCA(MKA)
+ integer :: IRLOCA(MKA)
+ integer :: IHLOCA(MKA)
+ integer :: IPLOCB(MKB)
+ integer :: IRLOCB(MKB)
+ integer :: IHLOCB(MKB)
+ ! vmod3D integers
+ integer :: NPA
+ integer :: NRA
+ integer :: NHA
+ integer :: NPB
+ integer :: NRB
+ integer :: NHB
+ ! weight integers
+ integer :: IP
+ integer :: JP
+ integer :: KP
+ integer :: IP1
+ integer :: JP1
+ integer :: KP1
end type model_jp3d_variables
type (model_jp3d_variables) JP3DM_V
! model_jp3d_variables
@@ -243,17 +246,17 @@
! model_sea99_s_variables
type model_sea99_s_variables
sequence
- integer :: sea99_ndep
- integer :: sea99_nlat
- integer :: sea99_nlon
- integer :: dummy_pad ! padding 4 bytes to align the structure
+ double precision :: sea99_vs(100,100,100)
+ double precision :: sea99_depth(100)
double precision :: sea99_ddeg
double precision :: alatmin
double precision :: alatmax
double precision :: alonmin
double precision :: alonmax
- double precision :: sea99_vs(100,100,100)
- double precision :: sea99_depth(100)
+ integer :: sea99_ndep
+ integer :: sea99_nlat
+ integer :: sea99_nlon
+ integer :: dummy_pad ! padding 4 bytes to align the structure
end type model_sea99_s_variables
type (model_sea99_s_variables) SEA99M_V
! model_sea99_s_variables
@@ -279,7 +282,7 @@
eucrust_vp_uppercrust,eucrust_vp_lowercrust,eucrust_mohodepth,&
eucrust_basement,eucrust_ucdepth
integer :: num_eucrust
- integer :: dummy ! padding 4 bytes to align the structure
+ integer :: dummy_pad ! padding 4 bytes to align the structure
end type model_eucrust_variables
type (model_eucrust_variables) EUCM_V
@@ -305,10 +308,10 @@
! model_attenuation_storage_var
type model_attenuation_storage_var
sequence
- integer Q_resolution
- integer Q_max
double precision, dimension(:,:), pointer :: tau_e_storage
double precision, dimension(:), pointer :: Qmu_storage
+ integer Q_resolution
+ integer Q_max
end type model_attenuation_storage_var
type (model_attenuation_storage_var) AM_S
! model_attenuation_storage_var
@@ -316,8 +319,6 @@
! attenuation_simplex_variables
type attenuation_simplex_variables
sequence
- integer nf ! nf = Number of Frequencies
- integer nsls ! nsls = Number of Standard Linear Solids
double precision Q ! Q = Desired Value of Attenuation or Q
double precision iQ ! iQ = 1/Q
double precision, dimension(:), pointer :: f
@@ -325,8 +326,9 @@
double precision, dimension(:), pointer :: tau_s
! tau_s = Tau_sigma defined by the frequency range and
! number of standard linear solids
+ integer nf ! nf = Number of Frequencies
+ integer nsls ! nsls = Number of Standard Linear Solids
end type attenuation_simplex_variables
-
type(attenuation_simplex_variables) AS_V
! attenuation_simplex_variables
@@ -337,7 +339,7 @@
double precision :: maxlat,maxlon,minlat,minlon,maxdepth,mindepth
double precision :: dlat,dlon,ddepth,max_dvs,min_dvs
integer :: num_v,num_latperlon,num_lonperdepth
- integer :: dummy ! padding 4 bytes to align the structure
+ integer :: dummy_pad ! padding 4 bytes to align the structure
end type model_ppm_variables
type (model_ppm_variables) PPM_V
@@ -346,7 +348,11 @@
sequence
! tomographic iteration model on GLL points
double precision :: scale_velocity,scale_density
+ ! isotropic model
real(kind=CUSTOM_REAL),dimension(:,:,:,:),pointer :: vs_new,vp_new,rho_new
+ ! transverse isotropic model
+ real(kind=CUSTOM_REAL),dimension(:,:,:,:),pointer :: vsv_new,vpv_new, &
+ vsh_new,vph_new,eta_new
logical :: MODEL_GLL
logical,dimension(3) :: dummy_pad ! padding 3 bytes to align the structure
end type model_gll_variables
@@ -362,26 +368,11 @@
integer, parameter :: maxver=1000
integer, parameter :: maxhpa=2
- integer numker
- integer numhpa,numcof
- integer ihpa,lmax,nylm
- integer lmxhpa(maxhpa)
- integer itypehpa(maxhpa)
- integer ihpakern(maxker)
- integer numcoe(maxhpa)
- integer ivarkern(maxker)
- integer itpspl(maxcoe,maxhpa)
-
- integer nconpt(maxhpa),iver
- integer iconpt(maxver,maxhpa)
real(kind=4) conpt(maxver,maxhpa)
-
real(kind=4) xlaspl(maxcoe,maxhpa)
real(kind=4) xlospl(maxcoe,maxhpa)
real(kind=4) radspl(maxcoe,maxhpa)
real(kind=4) coe(maxcoe,maxker)
- character(len=80) hsplfl(maxhpa)
- character(len=40) dskker(maxker)
real(kind=4) vercof(maxker)
real(kind=4) vercofd(maxker)
@@ -390,14 +381,29 @@
real(kind=4) wk2(maxl+1)
real(kind=4) wk3(maxl+1)
+ integer lmxhpa(maxhpa)
+ integer itypehpa(maxhpa)
+ integer ihpakern(maxker)
+ integer numcoe(maxhpa)
+ integer ivarkern(maxker)
+ integer itpspl(maxcoe,maxhpa)
+
+ integer nconpt(maxhpa),iver
+ integer iconpt(maxver,maxhpa)
+ integer numker
+ integer numhpa,numcof
+ integer ihpa,lmax,nylm
+
character(len=80) kerstr
character(len=80) refmdl
character(len=40) varstr(maxker)
+ character(len=80) hsplfl(maxhpa)
+ character(len=40) dskker(maxker)
! for ellipticity
+ double precision rspl(NR),espl(NR),espl2(NR)
integer nspl
- double precision rspl(NR),espl(NR),espl2(NR)
! model parameter and flags
integer REFERENCE_1D_MODEL,THREE_D_MODEL
@@ -1297,32 +1303,72 @@
! model GLL
if( MGLL_V%MODEL_GLL .and. iregion_code == IREGION_CRUST_MANTLE ) then
- !check
- if( ispec > size(MGLL_V%vp_new(1,1,1,:)) ) then
- call exit_MPI(myrank,'model gll: ispec too big')
- endif
- ! takes stored gll values from file
- if(CUSTOM_REAL == SIZE_REAL) then
- vp = dble( MGLL_V%vp_new(i,j,k,ispec) )
- vs = dble( MGLL_V%vs_new(i,j,k,ispec) )
- rho = dble( MGLL_V%rho_new(i,j,k,ispec) )
+
+ ! isotropic model
+ if( .not. TRANSVERSE_ISOTROPY ) then
+
+ !check
+ if( ispec > size(MGLL_V%vp_new(1,1,1,:)) ) then
+ call exit_MPI(myrank,'model gll: ispec too big')
+ endif
+
+ ! takes stored gll values from file
+ if(CUSTOM_REAL == SIZE_REAL) then
+ vp = dble( MGLL_V%vp_new(i,j,k,ispec) )
+ vs = dble( MGLL_V%vs_new(i,j,k,ispec) )
+ rho = dble( MGLL_V%rho_new(i,j,k,ispec) )
+ else
+ vp = MGLL_V%vp_new(i,j,k,ispec)
+ vs = MGLL_V%vs_new(i,j,k,ispec)
+ rho = MGLL_V%rho_new(i,j,k,ispec)
+ endif
+ ! non-dimensionalize
+ vp = vp * MGLL_V%scale_velocity
+ vs = vs * MGLL_V%scale_velocity
+ rho = rho * MGLL_V%scale_density
+ ! isotropic model
+ vpv = vp
+ vph = vp
+ vsv = vs
+ vsh = vs
+ rho = rho
+ dvp = 0.0d0
+ eta_aniso = 1.0d0
+
+ ! transverse isotropic model
else
- vp = MGLL_V%vp_new(i,j,k,ispec)
- vs = MGLL_V%vs_new(i,j,k,ispec)
- rho = MGLL_V%rho_new(i,j,k,ispec)
+
+ !check
+ if( ispec > size(MGLL_V%vpv_new(1,1,1,:)) ) then
+ call exit_MPI(myrank,'model gll: ispec too big')
+ endif
+
+ ! takes stored gll values from file
+ if(CUSTOM_REAL == SIZE_REAL) then
+ vph = dble( MGLL_V%vph_new(i,j,k,ispec) )
+ vsh = dble( MGLL_V%vsh_new(i,j,k,ispec) )
+ vpv = dble( MGLL_V%vpv_new(i,j,k,ispec) )
+ vsv = dble( MGLL_V%vsv_new(i,j,k,ispec) )
+ rho = dble( MGLL_V%rho_new(i,j,k,ispec) )
+ eta_aniso = dble( MGLL_V%eta_new(i,j,k,ispec) )
+ else
+ vph = MGLL_V%vph_new(i,j,k,ispec)
+ vsh = MGLL_V%vsh_new(i,j,k,ispec)
+ vpv = MGLL_V%vpv_new(i,j,k,ispec)
+ vsv = MGLL_V%vsv_new(i,j,k,ispec)
+ rho = MGLL_V%rho_new(i,j,k,ispec)
+ eta_aniso = MGLL_V%eta_new(i,j,k,ispec)
+ endif
+ ! non-dimensionalize
+ ! transverse isotropic model
+ vpv = vpv * MGLL_V%scale_velocity
+ vph = vph * MGLL_V%scale_velocity
+ vsv = vsv * MGLL_V%scale_velocity
+ vsh = vsh * MGLL_V%scale_velocity
+ rho = rho * MGLL_V%scale_density
+ dvp = 0.0d0
endif
- ! non-dimensionalize
- vp = vp * MGLL_V%scale_velocity
- vs = vs * MGLL_V%scale_velocity
- rho = rho * MGLL_V%scale_density
- ! isotropic model
- vpv = vp
- vph = vp
- vsv = vs
- vsh = vs
- rho = rho
- dvp = 0.0d0
- eta_aniso = 1.0d0
+
endif ! MODEL_GLL
end subroutine meshfem3D_models_impose_val
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/model_attenuation.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/model_attenuation.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/model_attenuation.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -69,7 +69,6 @@
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
@@ -77,6 +76,7 @@
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, dimension(:), pointer :: interval_Q ! Steps
integer :: Qn ! Number of points
integer dummy_pad ! padding 4 bytes to align the structure
end type model_attenuation_variables
@@ -119,7 +119,6 @@
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
@@ -127,6 +126,7 @@
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, dimension(:), pointer :: interval_Q ! Steps
integer :: Qn ! Number of points
integer dummy_pad ! padding 4 bytes to align the structure
end type model_attenuation_variables
@@ -173,7 +173,6 @@
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
@@ -181,6 +180,7 @@
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, dimension(:), pointer :: interval_Q ! Steps
integer :: Qn ! Number of points
integer dummy_pad ! padding 4 bytes to align the structure
end type model_attenuation_variables
@@ -250,10 +250,10 @@
! model_attenuation_storage_var
type model_attenuation_storage_var
sequence
- integer Q_resolution
- integer Q_max
double precision, dimension(:,:), pointer :: tau_e_storage
double precision, dimension(:), pointer :: Qmu_storage
+ integer Q_resolution
+ integer Q_max
end type model_attenuation_storage_var
type (model_attenuation_storage_var) AM_S
@@ -262,8 +262,6 @@
! attenuation_simplex_variables
type attenuation_simplex_variables
sequence
- integer nf ! nf = Number of Frequencies
- integer nsls ! nsls = Number of Standard Linear Solids
double precision Q ! Q = Desired Value of Attenuation or Q
double precision iQ ! iQ = 1/Q
double precision, dimension(:), pointer :: f
@@ -271,6 +269,8 @@
double precision, dimension(:), pointer :: tau_s
! tau_s = Tau_sigma defined by the frequency range and
! number of standard linear solids
+ integer nf ! nf = Number of Frequencies
+ integer nsls ! nsls = Number of Standard Linear Solids
end type attenuation_simplex_variables
type(attenuation_simplex_variables) AS_V
@@ -371,7 +371,6 @@
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
@@ -379,6 +378,7 @@
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, dimension(:), pointer :: interval_Q ! Steps
integer :: Qn ! Number of points
integer dummy_pad ! padding 4 bytes to align the structure
end type model_attenuation_variables
@@ -389,10 +389,10 @@
! model_attenuation_storage_var
type model_attenuation_storage_var
sequence
- integer Q_resolution
- integer Q_max
double precision, dimension(:,:), pointer :: tau_e_storage
double precision, dimension(:), pointer :: Qmu_storage
+ integer Q_resolution
+ integer Q_max
end type model_attenuation_storage_var
type (model_attenuation_storage_var) AM_S
@@ -401,8 +401,6 @@
! attenuation_simplex_variables
type attenuation_simplex_variables
sequence
- integer nf ! nf = Number of Frequencies
- integer nsls ! nsls = Number of Standard Linear Solids
double precision Q ! Q = Desired Value of Attenuation or Q
double precision iQ ! iQ = 1/Q
double precision, dimension(:), pointer :: f
@@ -410,6 +408,8 @@
double precision, dimension(:), pointer :: tau_s
! tau_s = Tau_sigma defined by the frequency range and
! number of standard linear solids
+ integer nf ! nf = Number of Frequencies
+ integer nsls ! nsls = Number of Standard Linear Solids
end type attenuation_simplex_variables
type(attenuation_simplex_variables) AS_V
@@ -446,10 +446,10 @@
! model_attenuation_storage_var
type model_attenuation_storage_var
sequence
- integer Q_resolution
- integer Q_max
double precision, dimension(:,:), pointer :: tau_e_storage
double precision, dimension(:), pointer :: Qmu_storage
+ integer Q_resolution
+ integer Q_max
end type model_attenuation_storage_var
type (model_attenuation_storage_var) AM_S
@@ -587,8 +587,6 @@
! attenuation_simplex_variables
type attenuation_simplex_variables
sequence
- integer nf ! nf = Number of Frequencies
- integer nsls ! nsls = Number of Standard Linear Solids
double precision Q ! Q = Desired Value of Attenuation or Q
double precision iQ ! iQ = 1/Q
double precision, dimension(:), pointer :: f
@@ -596,6 +594,8 @@
double precision, dimension(:), pointer :: tau_s
! tau_s = Tau_sigma defined by the frequency range and
! number of standard linear solids
+ integer nf ! nf = Number of Frequencies
+ integer nsls ! nsls = Number of Standard Linear Solids
end type attenuation_simplex_variables
type(attenuation_simplex_variables) AS_V
@@ -690,8 +690,6 @@
! attenuation_simplex_variables
type attenuation_simplex_variables
sequence
- integer nf ! nf = Number of Frequencies
- integer nsls ! nsls = Number of Standard Linear Solids
double precision Q ! Q = Desired Value of Attenuation or Q
double precision iQ ! iQ = 1/Q
double precision, dimension(:), pointer :: f
@@ -699,6 +697,8 @@
double precision, dimension(:), pointer :: tau_s
! tau_s = Tau_sigma defined by the frequency range and
! number of standard linear solids
+ integer nf ! nf = Number of Frequencies
+ integer nsls ! nsls = Number of Standard Linear Solids
end type attenuation_simplex_variables
type(attenuation_simplex_variables) AS_V
@@ -718,8 +718,6 @@
! attenuation_simplex_variables
type attenuation_simplex_variables
sequence
- integer nf ! nf = Number of Frequencies
- integer nsls ! nsls = Number of Standard Linear Solids
double precision Q ! Q = Desired Value of Attenuation or Q
double precision iQ ! iQ = 1/Q
double precision, dimension(:), pointer :: f
@@ -727,6 +725,8 @@
double precision, dimension(:), pointer :: tau_s
! tau_s = Tau_sigma defined by the frequency range and
! number of standard linear solids
+ integer nf ! nf = Number of Frequencies
+ integer nsls ! nsls = Number of Standard Linear Solids
end type attenuation_simplex_variables
type(attenuation_simplex_variables) AS_V
@@ -843,8 +843,6 @@
! attenuation_simplex_variables
type attenuation_simplex_variables
sequence
- integer nf ! nf = Number of Frequencies
- integer nsls ! nsls = Number of Standard Linear Solids
double precision Q ! Q = Desired Value of Attenuation or Q
double precision iQ ! iQ = 1/Q
double precision, dimension(:), pointer :: f
@@ -852,6 +850,8 @@
double precision, dimension(:), pointer :: tau_s
! tau_s = Tau_sigma defined by the frequency range and
! number of standard linear solids
+ integer nf ! nf = Number of Frequencies
+ integer nsls ! nsls = Number of Standard Linear Solids
end type attenuation_simplex_variables
type(attenuation_simplex_variables) AS_V
@@ -921,8 +921,6 @@
! attenuation_simplex_variables
type attenuation_simplex_variables
sequence
- integer nf ! nf = Number of Frequencies
- integer nsls ! nsls = Number of Standard Linear Solids
double precision Q ! Q = Desired Value of Attenuation or Q
double precision iQ ! iQ = 1/Q
double precision, dimension(:), pointer :: f
@@ -930,6 +928,8 @@
double precision, dimension(:), pointer :: tau_s
! tau_s = Tau_sigma defined by the frequency range and
! number of standard linear solids
+ integer nf ! nf = Number of Frequencies
+ integer nsls ! nsls = Number of Standard Linear Solids
end type attenuation_simplex_variables
type(attenuation_simplex_variables) AS_V
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/model_eucrust.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/model_eucrust.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/model_eucrust.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -50,7 +50,7 @@
eucrust_vp_uppercrust,eucrust_vp_lowercrust,eucrust_mohodepth,&
eucrust_basement,eucrust_ucdepth
integer :: num_eucrust
- integer :: dummy ! padding 4 bytes to align the structure
+ integer :: dummy_pad ! padding 4 bytes to align the structure
end type model_eucrust_variables
type (model_eucrust_variables) EUCM_V
@@ -94,7 +94,7 @@
eucrust_vp_uppercrust,eucrust_vp_lowercrust,eucrust_mohodepth,&
eucrust_basement,eucrust_ucdepth
integer :: num_eucrust
- integer :: dummy ! padding 4 bytes to align the structure
+ integer :: dummy_pad ! padding 4 bytes to align the structure
end type model_eucrust_variables
type (model_eucrust_variables) EUCM_V
@@ -163,7 +163,7 @@
eucrust_vp_uppercrust,eucrust_vp_lowercrust,eucrust_mohodepth,&
eucrust_basement,eucrust_ucdepth
integer :: num_eucrust
- integer :: dummy ! padding 4 bytes to align the structure
+ integer :: dummy_pad ! padding 4 bytes to align the structure
end type model_eucrust_variables
type (model_eucrust_variables) EUCM_V
@@ -212,7 +212,7 @@
eucrust_vp_uppercrust,eucrust_vp_lowercrust,eucrust_mohodepth,&
eucrust_basement,eucrust_ucdepth
integer :: num_eucrust
- integer :: dummy ! padding 4 bytes to align the structure
+ integer :: dummy_pad ! padding 4 bytes to align the structure
end type model_eucrust_variables
type (model_eucrust_variables) EUCM_V
@@ -306,7 +306,7 @@
eucrust_vp_uppercrust,eucrust_vp_lowercrust,eucrust_mohodepth,&
eucrust_basement,eucrust_ucdepth
integer :: num_eucrust
- integer :: dummy ! padding 4 bytes to align the structure
+ integer :: dummy_pad ! padding 4 bytes to align the structure
end type model_eucrust_variables
type (model_eucrust_variables) EUCM_V
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/model_gll.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/model_gll.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/model_gll.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -38,6 +38,8 @@
! standard routine to setup model
+ use meshfem3D_models_par,only: TRANSVERSE_ISOTROPY
+
implicit none
include "constants.h"
@@ -49,7 +51,11 @@
sequence
! tomographic iteration model on GLL points
double precision :: scale_velocity,scale_density
+ ! isotropic model
real(kind=CUSTOM_REAL),dimension(:,:,:,:),pointer :: vs_new,vp_new,rho_new
+ ! transverse isotropic model
+ real(kind=CUSTOM_REAL),dimension(:,:,:,:),pointer :: vsv_new,vpv_new, &
+ vsh_new,vph_new,eta_new
logical :: MODEL_GLL
logical,dimension(3) :: dummy_pad ! padding 3 bytes to align the structure
end type model_gll_variables
@@ -62,8 +68,19 @@
double precision :: min_dvs,max_dvs,min_dvs_all,max_dvs_all,scaleval
integer :: ier
- allocate( MGLL_V%vp_new(NGLLX,NGLLY,NGLLZ,NSPEC(IREGION_CRUST_MANTLE)) )
- allocate( MGLL_V%vs_new(NGLLX,NGLLY,NGLLZ,NSPEC(IREGION_CRUST_MANTLE)) )
+ ! differs for isotropic model or transverse isotropic models
+ if( .not. TRANSVERSE_ISOTROPY ) then
+ ! isotropic model
+ allocate( MGLL_V%vp_new(NGLLX,NGLLY,NGLLZ,NSPEC(IREGION_CRUST_MANTLE)) )
+ allocate( MGLL_V%vs_new(NGLLX,NGLLY,NGLLZ,NSPEC(IREGION_CRUST_MANTLE)) )
+ else
+ ! transverse isotropic model
+ allocate( MGLL_V%vpv_new(NGLLX,NGLLY,NGLLZ,NSPEC(IREGION_CRUST_MANTLE)) )
+ allocate( MGLL_V%vsv_new(NGLLX,NGLLY,NGLLZ,NSPEC(IREGION_CRUST_MANTLE)) )
+ allocate( MGLL_V%vph_new(NGLLX,NGLLY,NGLLZ,NSPEC(IREGION_CRUST_MANTLE)) )
+ allocate( MGLL_V%vsh_new(NGLLX,NGLLY,NGLLZ,NSPEC(IREGION_CRUST_MANTLE)) )
+ allocate( MGLL_V%eta_new(NGLLX,NGLLY,NGLLZ,NSPEC(IREGION_CRUST_MANTLE)) )
+ endif
allocate( MGLL_V%rho_new(NGLLX,NGLLY,NGLLZ,NSPEC(IREGION_CRUST_MANTLE)) )
! non-dimensionalize scaling values
@@ -75,8 +92,13 @@
call read_gll_model(myrank,MGLL_V,NSPEC)
! checks velocity range
- max_dvs = maxval( MGLL_V%vs_new )
- min_dvs = minval( MGLL_V%vs_new )
+ if( .not. TRANSVERSE_ISOTROPY ) then
+ max_dvs = maxval( MGLL_V%vs_new )
+ min_dvs = minval( MGLL_V%vs_new )
+ else
+ max_dvs = maxval( MGLL_V%vsv_new + MGLL_V%vsh_new )
+ min_dvs = minval( MGLL_V%vsv_new + MGLL_V%vsh_new )
+ endif
call mpi_reduce(max_dvs, max_dvs_all, 1, MPI_DOUBLE_PRECISION, MPI_MAX, 0, MPI_COMM_WORLD,ier)
call mpi_reduce(min_dvs, min_dvs_all, 1, MPI_DOUBLE_PRECISION, MPI_MIN, 0, MPI_COMM_WORLD,ier)
if( myrank == 0 ) then
@@ -94,6 +116,8 @@
subroutine read_gll_model(myrank,MGLL_V,NSPEC)
+ use meshfem3D_models_par,only: TRANSVERSE_ISOTROPY
+
implicit none
include "constants.h"
@@ -103,7 +127,11 @@
sequence
! tomographic iteration model on GLL points
double precision :: scale_velocity,scale_density
+ ! isotropic model
real(kind=CUSTOM_REAL),dimension(:,:,:,:),pointer :: vs_new,vp_new,rho_new
+ ! transverse isotropic model
+ real(kind=CUSTOM_REAL),dimension(:,:,:,:),pointer :: vsv_new,vpv_new, &
+ vsh_new,vph_new,eta_new
logical :: MODEL_GLL
logical,dimension(3) :: dummy_pad ! padding 3 bytes to align the structure
end type model_gll_variables
@@ -130,25 +158,81 @@
! only crust and mantle
write(prname,'(a,i6.6,a)') MGLL_path(1:len_trim(MGLL_path))//'proc',myrank,'_reg1_'
- ! vp mesh
- open(unit=27,file=prname(1:len_trim(prname))//'vp_new.bin',&
- status='old',action='read',form='unformatted',iostat=ier)
- if( ier /= 0 ) then
- write(IMAIN,*) 'error opening: ',prname(1:len_trim(prname))//'vp_new.bin'
- call exit_MPI(myrank,'error model gll')
- endif
- read(27) MGLL_V%vp_new(:,:,:,1:nspec(IREGION_CRUST_MANTLE))
- close(27)
+ ! reads in model for each partition
+ if( .not. TRANSVERSE_ISOTROPY ) then
+ ! isotropic model
+ ! vp mesh
+ open(unit=27,file=prname(1:len_trim(prname))//'vp_new.bin',&
+ status='old',action='read',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ write(IMAIN,*) 'error opening: ',prname(1:len_trim(prname))//'vp_new.bin'
+ call exit_MPI(myrank,'error model gll')
+ endif
+ read(27) MGLL_V%vp_new(:,:,:,1:nspec(IREGION_CRUST_MANTLE))
+ close(27)
- ! vs mesh
- open(unit=27,file=prname(1:len_trim(prname))//'vs_new.bin', &
- status='old',action='read',form='unformatted',iostat=ier)
- if( ier /= 0 ) then
- print*,'error opening: ',prname(1:len_trim(prname))//'vs_new.bin'
- call exit_MPI(myrank,'error model gll')
+ ! vs mesh
+ open(unit=27,file=prname(1:len_trim(prname))//'vs_new.bin', &
+ status='old',action='read',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',prname(1:len_trim(prname))//'vs_new.bin'
+ call exit_MPI(myrank,'error model gll')
+ endif
+ read(27) MGLL_V%vs_new(:,:,:,1:nspec(IREGION_CRUST_MANTLE))
+ close(27)
+
+ else
+
+ ! transverse isotropic model
+ ! vp mesh
+ open(unit=27,file=prname(1:len_trim(prname))//'vpv_new.bin',&
+ status='old',action='read',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ write(IMAIN,*) 'error opening: ',prname(1:len_trim(prname))//'vpv_new.bin'
+ call exit_MPI(myrank,'error model gll')
+ endif
+ read(27) MGLL_V%vpv_new(:,:,:,1:nspec(IREGION_CRUST_MANTLE))
+ close(27)
+
+ open(unit=27,file=prname(1:len_trim(prname))//'vph_new.bin',&
+ status='old',action='read',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ write(IMAIN,*) 'error opening: ',prname(1:len_trim(prname))//'vph_new.bin'
+ call exit_MPI(myrank,'error model gll')
+ endif
+ read(27) MGLL_V%vph_new(:,:,:,1:nspec(IREGION_CRUST_MANTLE))
+ close(27)
+
+ ! vs mesh
+ open(unit=27,file=prname(1:len_trim(prname))//'vsv_new.bin', &
+ status='old',action='read',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',prname(1:len_trim(prname))//'vsv_new.bin'
+ call exit_MPI(myrank,'error model gll')
+ endif
+ read(27) MGLL_V%vsv_new(:,:,:,1:nspec(IREGION_CRUST_MANTLE))
+ close(27)
+
+ open(unit=27,file=prname(1:len_trim(prname))//'vsh_new.bin', &
+ status='old',action='read',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',prname(1:len_trim(prname))//'vsh_new.bin'
+ call exit_MPI(myrank,'error model gll')
+ endif
+ read(27) MGLL_V%vsh_new(:,:,:,1:nspec(IREGION_CRUST_MANTLE))
+ close(27)
+
+ ! eta mesh
+ open(unit=27,file=prname(1:len_trim(prname))//'eta_new.bin', &
+ status='old',action='read',form='unformatted',iostat=ier)
+ if( ier /= 0 ) then
+ print*,'error opening: ',prname(1:len_trim(prname))//'eta_new.bin'
+ call exit_MPI(myrank,'error model gll')
+ endif
+ read(27) MGLL_V%eta_new(:,:,:,1:nspec(IREGION_CRUST_MANTLE))
+ close(27)
+
endif
- read(27) MGLL_V%vs_new(:,:,:,1:nspec(IREGION_CRUST_MANTLE))
- close(27)
! rho mesh
open(unit=27,file=prname(1:len_trim(prname))//'rho_new.bin', &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/model_jp3d.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/model_jp3d.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/model_jp3d.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -72,66 +72,69 @@
! model_jp3d_variables
type model_jp3d_variables
sequence
-! vmod3d
- integer :: NPA
- integer :: NRA
- integer :: NHA
- integer :: NPB
- integer :: NRB
- integer :: NHB
- double precision :: PNA(MPA)
- double precision :: RNA(MRA)
- double precision :: HNA(MHA)
- double precision :: PNB(MPB)
- double precision :: RNB(MRB)
- double precision :: HNB(MHB)
- double precision :: VELAP(MPA,MRA,MHA)
- double precision :: VELBP(MPB,MRB,MHB)
-! discon
- double precision :: PN(51)
- double precision :: RRN(63)
- double precision :: DEPA(51,63)
- double precision :: DEPB(51,63)
- double precision :: DEPC(51,63)
-! locate
- integer :: IPLOCA(MKA)
- integer :: IRLOCA(MKA)
- integer :: IHLOCA(MKA)
- integer :: IPLOCB(MKB)
- integer :: IRLOCB(MKB)
- integer :: IHLOCB(MKB)
- double precision :: PLA
- double precision :: RLA
- double precision :: HLA
- double precision :: PLB
- double precision :: RLB
- double precision :: HLB
-! weight
- integer :: IP
- integer :: JP
- integer :: KP
- integer :: IP1
- integer :: JP1
- integer :: KP1
- double precision :: WV(8)
-! prhfd
- double precision :: P
- double precision :: R
- double precision :: H
- double precision :: PF
- double precision :: RF
- double precision :: HF
- double precision :: PF1
- double precision :: RF1
- double precision :: HF1
- double precision :: PD
- double precision :: RD
- double precision :: HD
-! jpmodv
- double precision :: VP(29)
- double precision :: VS(29)
- double precision :: RA(29)
- double precision :: DEPJ(29)
+ ! vmod3d
+ double precision :: PNA(MPA)
+ double precision :: RNA(MRA)
+ double precision :: HNA(MHA)
+ double precision :: PNB(MPB)
+ double precision :: RNB(MRB)
+ double precision :: HNB(MHB)
+ double precision :: VELAP(MPA,MRA,MHA)
+ double precision :: VELBP(MPB,MRB,MHB)
+ ! discon
+ double precision :: PN(51)
+ double precision :: RRN(63)
+ double precision :: DEPA(51,63)
+ double precision :: DEPB(51,63)
+ double precision :: DEPC(51,63)
+ ! locate
+ double precision :: PLA
+ double precision :: RLA
+ double precision :: HLA
+ double precision :: PLB
+ double precision :: RLB
+ double precision :: HLB
+ ! weight
+ double precision :: WV(8)
+ ! prhfd
+ double precision :: P
+ double precision :: R
+ double precision :: H
+ double precision :: PF
+ double precision :: RF
+ double precision :: HF
+ double precision :: PF1
+ double precision :: RF1
+ double precision :: HF1
+ double precision :: PD
+ double precision :: RD
+ double precision :: HD
+ ! jpmodv
+ double precision :: VP(29)
+ double precision :: VS(29)
+ double precision :: RA(29)
+ double precision :: DEPJ(29)
+ ! locate integers
+ integer :: IPLOCA(MKA)
+ integer :: IRLOCA(MKA)
+ integer :: IHLOCA(MKA)
+ integer :: IPLOCB(MKB)
+ integer :: IRLOCB(MKB)
+ integer :: IHLOCB(MKB)
+ ! vmod3D integers
+ integer :: NPA
+ integer :: NRA
+ integer :: NHA
+ integer :: NPB
+ integer :: NRB
+ integer :: NHB
+ ! weight integers
+ integer :: IP
+ integer :: JP
+ integer :: KP
+ integer :: IP1
+ integer :: JP1
+ integer :: KP1
end type model_jp3d_variables
type (model_jp3d_variables) JP3DM_V
@@ -213,66 +216,69 @@
! model_jp3d_variables
type model_jp3d_variables
sequence
-! vmod3d
- integer :: NPA
- integer :: NRA
- integer :: NHA
- integer :: NPB
- integer :: NRB
- integer :: NHB
- double precision :: PNA(MPA)
- double precision :: RNA(MRA)
- double precision :: HNA(MHA)
- double precision :: PNB(MPB)
- double precision :: RNB(MRB)
- double precision :: HNB(MHB)
- double precision :: VELAP(MPA,MRA,MHA)
- double precision :: VELBP(MPB,MRB,MHB)
-! discon
- double precision :: PN(51)
- double precision :: RRN(63)
- double precision :: DEPA(51,63)
- double precision :: DEPB(51,63)
- double precision :: DEPC(51,63)
-! locate
- integer :: IPLOCA(MKA)
- integer :: IRLOCA(MKA)
- integer :: IHLOCA(MKA)
- integer :: IPLOCB(MKB)
- integer :: IRLOCB(MKB)
- integer :: IHLOCB(MKB)
- double precision :: PLA
- double precision :: RLA
- double precision :: HLA
- double precision :: PLB
- double precision :: RLB
- double precision :: HLB
-! weight
- integer :: IP
- integer :: JP
- integer :: KP
- integer :: IP1
- integer :: JP1
- integer :: KP1
- double precision :: WV(8)
-! prhfd
- double precision :: P
- double precision :: R
- double precision :: H
- double precision :: PF
- double precision :: RF
- double precision :: HF
- double precision :: PF1
- double precision :: RF1
- double precision :: HF1
- double precision :: PD
- double precision :: RD
- double precision :: HD
-! jpmodv
- double precision :: VP(29)
- double precision :: VS(29)
- double precision :: RA(29)
- double precision :: DEPJ(29)
+ ! vmod3d
+ double precision :: PNA(MPA)
+ double precision :: RNA(MRA)
+ double precision :: HNA(MHA)
+ double precision :: PNB(MPB)
+ double precision :: RNB(MRB)
+ double precision :: HNB(MHB)
+ double precision :: VELAP(MPA,MRA,MHA)
+ double precision :: VELBP(MPB,MRB,MHB)
+ ! discon
+ double precision :: PN(51)
+ double precision :: RRN(63)
+ double precision :: DEPA(51,63)
+ double precision :: DEPB(51,63)
+ double precision :: DEPC(51,63)
+ ! locate
+ double precision :: PLA
+ double precision :: RLA
+ double precision :: HLA
+ double precision :: PLB
+ double precision :: RLB
+ double precision :: HLB
+ ! weight
+ double precision :: WV(8)
+ ! prhfd
+ double precision :: P
+ double precision :: R
+ double precision :: H
+ double precision :: PF
+ double precision :: RF
+ double precision :: HF
+ double precision :: PF1
+ double precision :: RF1
+ double precision :: HF1
+ double precision :: PD
+ double precision :: RD
+ double precision :: HD
+ ! jpmodv
+ double precision :: VP(29)
+ double precision :: VS(29)
+ double precision :: RA(29)
+ double precision :: DEPJ(29)
+ ! locate integers
+ integer :: IPLOCA(MKA)
+ integer :: IRLOCA(MKA)
+ integer :: IHLOCA(MKA)
+ integer :: IPLOCB(MKB)
+ integer :: IRLOCB(MKB)
+ integer :: IHLOCB(MKB)
+ ! vmod3D integers
+ integer :: NPA
+ integer :: NRA
+ integer :: NHA
+ integer :: NPB
+ integer :: NRB
+ integer :: NHB
+ ! weight integers
+ integer :: IP
+ integer :: JP
+ integer :: KP
+ integer :: IP1
+ integer :: JP1
+ integer :: KP1
end type model_jp3d_variables
type (model_jp3d_variables) JP3DM_V
@@ -299,66 +305,69 @@
! model_jp3d_variables
type model_jp3d_variables
sequence
-! vmod3d
- integer :: NPA
- integer :: NRA
- integer :: NHA
- integer :: NPB
- integer :: NRB
- integer :: NHB
- double precision :: PNA(MPA)
- double precision :: RNA(MRA)
- double precision :: HNA(MHA)
- double precision :: PNB(MPB)
- double precision :: RNB(MRB)
- double precision :: HNB(MHB)
- double precision :: VELAP(MPA,MRA,MHA)
- double precision :: VELBP(MPB,MRB,MHB)
-! discon
- double precision :: PN(51)
- double precision :: RRN(63)
- double precision :: DEPA(51,63)
- double precision :: DEPB(51,63)
- double precision :: DEPC(51,63)
-! locate
- integer :: IPLOCA(MKA)
- integer :: IRLOCA(MKA)
- integer :: IHLOCA(MKA)
- integer :: IPLOCB(MKB)
- integer :: IRLOCB(MKB)
- integer :: IHLOCB(MKB)
- double precision :: PLA
- double precision :: RLA
- double precision :: HLA
- double precision :: PLB
- double precision :: RLB
- double precision :: HLB
-! weight
- integer :: IP
- integer :: JP
- integer :: KP
- integer :: IP1
- integer :: JP1
- integer :: KP1
- double precision :: WV(8)
-! prhfd
- double precision :: P
- double precision :: R
- double precision :: H
- double precision :: PF
- double precision :: RF
- double precision :: HF
- double precision :: PF1
- double precision :: RF1
- double precision :: HF1
- double precision :: PD
- double precision :: RD
- double precision :: HD
-! jpmodv
- double precision :: VP(29)
- double precision :: VS(29)
- double precision :: RA(29)
- double precision :: DEPJ(29)
+ ! vmod3d
+ double precision :: PNA(MPA)
+ double precision :: RNA(MRA)
+ double precision :: HNA(MHA)
+ double precision :: PNB(MPB)
+ double precision :: RNB(MRB)
+ double precision :: HNB(MHB)
+ double precision :: VELAP(MPA,MRA,MHA)
+ double precision :: VELBP(MPB,MRB,MHB)
+ ! discon
+ double precision :: PN(51)
+ double precision :: RRN(63)
+ double precision :: DEPA(51,63)
+ double precision :: DEPB(51,63)
+ double precision :: DEPC(51,63)
+ ! locate
+ double precision :: PLA
+ double precision :: RLA
+ double precision :: HLA
+ double precision :: PLB
+ double precision :: RLB
+ double precision :: HLB
+ ! weight
+ double precision :: WV(8)
+ ! prhfd
+ double precision :: P
+ double precision :: R
+ double precision :: H
+ double precision :: PF
+ double precision :: RF
+ double precision :: HF
+ double precision :: PF1
+ double precision :: RF1
+ double precision :: HF1
+ double precision :: PD
+ double precision :: RD
+ double precision :: HD
+ ! jpmodv
+ double precision :: VP(29)
+ double precision :: VS(29)
+ double precision :: RA(29)
+ double precision :: DEPJ(29)
+ ! locate integers
+ integer :: IPLOCA(MKA)
+ integer :: IRLOCA(MKA)
+ integer :: IHLOCA(MKA)
+ integer :: IPLOCB(MKB)
+ integer :: IRLOCB(MKB)
+ integer :: IHLOCB(MKB)
+ ! vmod3D integers
+ integer :: NPA
+ integer :: NRA
+ integer :: NHA
+ integer :: NPB
+ integer :: NRB
+ integer :: NHB
+ ! weight integers
+ integer :: IP
+ integer :: JP
+ integer :: KP
+ integer :: IP1
+ integer :: JP1
+ integer :: KP1
end type model_jp3d_variables
type (model_jp3d_variables) JP3DM_V
@@ -435,66 +444,69 @@
! model_jp3d_variables
type model_jp3d_variables
sequence
-! vmod3d
- integer :: NPA
- integer :: NRA
- integer :: NHA
- integer :: NPB
- integer :: NRB
- integer :: NHB
- double precision :: PNA(MPA)
- double precision :: RNA(MRA)
- double precision :: HNA(MHA)
- double precision :: PNB(MPB)
- double precision :: RNB(MRB)
- double precision :: HNB(MHB)
- double precision :: VELAP(MPA,MRA,MHA)
- double precision :: VELBP(MPB,MRB,MHB)
-! discon
- double precision :: PN(51)
- double precision :: RRN(63)
- double precision :: DEPA(51,63)
- double precision :: DEPB(51,63)
- double precision :: DEPC(51,63)
-! locate
- integer :: IPLOCA(MKA)
- integer :: IRLOCA(MKA)
- integer :: IHLOCA(MKA)
- integer :: IPLOCB(MKB)
- integer :: IRLOCB(MKB)
- integer :: IHLOCB(MKB)
- double precision :: PLA
- double precision :: RLA
- double precision :: HLA
- double precision :: PLB
- double precision :: RLB
- double precision :: HLB
-! weight
- integer :: IP
- integer :: JP
- integer :: KP
- integer :: IP1
- integer :: JP1
- integer :: KP1
- double precision :: WV(8)
-! prhfd
- double precision :: P
- double precision :: R
- double precision :: H
- double precision :: PF
- double precision :: RF
- double precision :: HF
- double precision :: PF1
- double precision :: RF1
- double precision :: HF1
- double precision :: PD
- double precision :: RD
- double precision :: HD
-! jpmodv
- double precision :: VP(29)
- double precision :: VS(29)
- double precision :: RA(29)
- double precision :: DEPJ(29)
+ ! vmod3d
+ double precision :: PNA(MPA)
+ double precision :: RNA(MRA)
+ double precision :: HNA(MHA)
+ double precision :: PNB(MPB)
+ double precision :: RNB(MRB)
+ double precision :: HNB(MHB)
+ double precision :: VELAP(MPA,MRA,MHA)
+ double precision :: VELBP(MPB,MRB,MHB)
+ ! discon
+ double precision :: PN(51)
+ double precision :: RRN(63)
+ double precision :: DEPA(51,63)
+ double precision :: DEPB(51,63)
+ double precision :: DEPC(51,63)
+ ! locate
+ double precision :: PLA
+ double precision :: RLA
+ double precision :: HLA
+ double precision :: PLB
+ double precision :: RLB
+ double precision :: HLB
+ ! weight
+ double precision :: WV(8)
+ ! prhfd
+ double precision :: P
+ double precision :: R
+ double precision :: H
+ double precision :: PF
+ double precision :: RF
+ double precision :: HF
+ double precision :: PF1
+ double precision :: RF1
+ double precision :: HF1
+ double precision :: PD
+ double precision :: RD
+ double precision :: HD
+ ! jpmodv
+ double precision :: VP(29)
+ double precision :: VS(29)
+ double precision :: RA(29)
+ double precision :: DEPJ(29)
+ ! locate integers
+ integer :: IPLOCA(MKA)
+ integer :: IRLOCA(MKA)
+ integer :: IHLOCA(MKA)
+ integer :: IPLOCB(MKB)
+ integer :: IRLOCB(MKB)
+ integer :: IHLOCB(MKB)
+ ! vmod3D integers
+ integer :: NPA
+ integer :: NRA
+ integer :: NHA
+ integer :: NPB
+ integer :: NRB
+ integer :: NHB
+ ! weight integers
+ integer :: IP
+ integer :: JP
+ integer :: KP
+ integer :: IP1
+ integer :: JP1
+ integer :: KP1
end type model_jp3d_variables
type (model_jp3d_variables) JP3DM_V
@@ -537,66 +549,69 @@
! model_jp3d_variables
type model_jp3d_variables
sequence
-! vmod3d
- integer :: NPA
- integer :: NRA
- integer :: NHA
- integer :: NPB
- integer :: NRB
- integer :: NHB
- double precision :: PNA(MPA)
- double precision :: RNA(MRA)
- double precision :: HNA(MHA)
- double precision :: PNB(MPB)
- double precision :: RNB(MRB)
- double precision :: HNB(MHB)
- double precision :: VELAP(MPA,MRA,MHA)
- double precision :: VELBP(MPB,MRB,MHB)
-! discon
- double precision :: PN(51)
- double precision :: RRN(63)
- double precision :: DEPA(51,63)
- double precision :: DEPB(51,63)
- double precision :: DEPC(51,63)
-! locate
- integer :: IPLOCA(MKA)
- integer :: IRLOCA(MKA)
- integer :: IHLOCA(MKA)
- integer :: IPLOCB(MKB)
- integer :: IRLOCB(MKB)
- integer :: IHLOCB(MKB)
- double precision :: PLA
- double precision :: RLA
- double precision :: HLA
- double precision :: PLB
- double precision :: RLB
- double precision :: HLB
-! weight
- integer :: IP
- integer :: JP
- integer :: KP
- integer :: IP1
- integer :: JP1
- integer :: KP1
- double precision :: WV(8)
-! prhfd
- double precision :: P
- double precision :: R
- double precision :: H
- double precision :: PF
- double precision :: RF
- double precision :: HF
- double precision :: PF1
- double precision :: RF1
- double precision :: HF1
- double precision :: PD
- double precision :: RD
- double precision :: HD
-! jpmodv
- double precision :: VP(29)
- double precision :: VS(29)
- double precision :: RA(29)
- double precision :: DEPJ(29)
+ ! vmod3d
+ double precision :: PNA(MPA)
+ double precision :: RNA(MRA)
+ double precision :: HNA(MHA)
+ double precision :: PNB(MPB)
+ double precision :: RNB(MRB)
+ double precision :: HNB(MHB)
+ double precision :: VELAP(MPA,MRA,MHA)
+ double precision :: VELBP(MPB,MRB,MHB)
+ ! discon
+ double precision :: PN(51)
+ double precision :: RRN(63)
+ double precision :: DEPA(51,63)
+ double precision :: DEPB(51,63)
+ double precision :: DEPC(51,63)
+ ! locate
+ double precision :: PLA
+ double precision :: RLA
+ double precision :: HLA
+ double precision :: PLB
+ double precision :: RLB
+ double precision :: HLB
+ ! weight
+ double precision :: WV(8)
+ ! prhfd
+ double precision :: P
+ double precision :: R
+ double precision :: H
+ double precision :: PF
+ double precision :: RF
+ double precision :: HF
+ double precision :: PF1
+ double precision :: RF1
+ double precision :: HF1
+ double precision :: PD
+ double precision :: RD
+ double precision :: HD
+ ! jpmodv
+ double precision :: VP(29)
+ double precision :: VS(29)
+ double precision :: RA(29)
+ double precision :: DEPJ(29)
+ ! locate integers
+ integer :: IPLOCA(MKA)
+ integer :: IRLOCA(MKA)
+ integer :: IHLOCA(MKA)
+ integer :: IPLOCB(MKB)
+ integer :: IRLOCB(MKB)
+ integer :: IHLOCB(MKB)
+ ! vmod3D integers
+ integer :: NPA
+ integer :: NRA
+ integer :: NHA
+ integer :: NPB
+ integer :: NRB
+ integer :: NHB
+ ! weight integers
+ integer :: IP
+ integer :: JP
+ integer :: KP
+ integer :: IP1
+ integer :: JP1
+ integer :: KP1
end type model_jp3d_variables
type (model_jp3d_variables) JP3DM_V
@@ -633,66 +648,69 @@
! model_jp3d_variables
type model_jp3d_variables
sequence
-! vmod3d
- integer :: NPA
- integer :: NRA
- integer :: NHA
- integer :: NPB
- integer :: NRB
- integer :: NHB
- double precision :: PNA(MPA)
- double precision :: RNA(MRA)
- double precision :: HNA(MHA)
- double precision :: PNB(MPB)
- double precision :: RNB(MRB)
- double precision :: HNB(MHB)
- double precision :: VELAP(MPA,MRA,MHA)
- double precision :: VELBP(MPB,MRB,MHB)
-! discon
- double precision :: PN(51)
- double precision :: RRN(63)
- double precision :: DEPA(51,63)
- double precision :: DEPB(51,63)
- double precision :: DEPC(51,63)
-! locate
- integer :: IPLOCA(MKA)
- integer :: IRLOCA(MKA)
- integer :: IHLOCA(MKA)
- integer :: IPLOCB(MKB)
- integer :: IRLOCB(MKB)
- integer :: IHLOCB(MKB)
- double precision :: PLA
- double precision :: RLA
- double precision :: HLA
- double precision :: PLB
- double precision :: RLB
- double precision :: HLB
-! weight
- integer :: IP
- integer :: JP
- integer :: KP
- integer :: IP1
- integer :: JP1
- integer :: KP1
- double precision :: WV(8)
-! prhfd
- double precision :: P
- double precision :: R
- double precision :: H
- double precision :: PF
- double precision :: RF
- double precision :: HF
- double precision :: PF1
- double precision :: RF1
- double precision :: HF1
- double precision :: PD
- double precision :: RD
- double precision :: HD
-! jpmodv
- double precision :: VP(29)
- double precision :: VS(29)
- double precision :: RA(29)
- double precision :: DEPJ(29)
+ ! vmod3d
+ double precision :: PNA(MPA)
+ double precision :: RNA(MRA)
+ double precision :: HNA(MHA)
+ double precision :: PNB(MPB)
+ double precision :: RNB(MRB)
+ double precision :: HNB(MHB)
+ double precision :: VELAP(MPA,MRA,MHA)
+ double precision :: VELBP(MPB,MRB,MHB)
+ ! discon
+ double precision :: PN(51)
+ double precision :: RRN(63)
+ double precision :: DEPA(51,63)
+ double precision :: DEPB(51,63)
+ double precision :: DEPC(51,63)
+ ! locate
+ double precision :: PLA
+ double precision :: RLA
+ double precision :: HLA
+ double precision :: PLB
+ double precision :: RLB
+ double precision :: HLB
+ ! weight
+ double precision :: WV(8)
+ ! prhfd
+ double precision :: P
+ double precision :: R
+ double precision :: H
+ double precision :: PF
+ double precision :: RF
+ double precision :: HF
+ double precision :: PF1
+ double precision :: RF1
+ double precision :: HF1
+ double precision :: PD
+ double precision :: RD
+ double precision :: HD
+ ! jpmodv
+ double precision :: VP(29)
+ double precision :: VS(29)
+ double precision :: RA(29)
+ double precision :: DEPJ(29)
+ ! locate integers
+ integer :: IPLOCA(MKA)
+ integer :: IRLOCA(MKA)
+ integer :: IHLOCA(MKA)
+ integer :: IPLOCB(MKB)
+ integer :: IRLOCB(MKB)
+ integer :: IHLOCB(MKB)
+ ! vmod3D integers
+ integer :: NPA
+ integer :: NRA
+ integer :: NHA
+ integer :: NPB
+ integer :: NRB
+ integer :: NHB
+ ! weight integers
+ integer :: IP
+ integer :: JP
+ integer :: KP
+ integer :: IP1
+ integer :: JP1
+ integer :: KP1
end type model_jp3d_variables
type (model_jp3d_variables) JP3DM_V
@@ -752,66 +770,69 @@
! model_jp3d_variables
type model_jp3d_variables
sequence
-! vmod3d
- integer :: NPA
- integer :: NRA
- integer :: NHA
- integer :: NPB
- integer :: NRB
- integer :: NHB
- double precision :: PNA(MPA)
- double precision :: RNA(MRA)
- double precision :: HNA(MHA)
- double precision :: PNB(MPB)
- double precision :: RNB(MRB)
- double precision :: HNB(MHB)
- double precision :: VELAP(MPA,MRA,MHA)
- double precision :: VELBP(MPB,MRB,MHB)
-! discon
- double precision :: PN(51)
- double precision :: RRN(63)
- double precision :: DEPA(51,63)
- double precision :: DEPB(51,63)
- double precision :: DEPC(51,63)
-! locate
- integer :: IPLOCA(MKA)
- integer :: IRLOCA(MKA)
- integer :: IHLOCA(MKA)
- integer :: IPLOCB(MKB)
- integer :: IRLOCB(MKB)
- integer :: IHLOCB(MKB)
- double precision :: PLA
- double precision :: RLA
- double precision :: HLA
- double precision :: PLB
- double precision :: RLB
- double precision :: HLB
-! weight
- integer :: IP
- integer :: JP
- integer :: KP
- integer :: IP1
- integer :: JP1
- integer :: KP1
- double precision :: WV(8)
-! prhfd
- double precision :: P
- double precision :: R
- double precision :: H
- double precision :: PF
- double precision :: RF
- double precision :: HF
- double precision :: PF1
- double precision :: RF1
- double precision :: HF1
- double precision :: PD
- double precision :: RD
- double precision :: HD
-! jpmodv
- double precision :: VP(29)
- double precision :: VS(29)
- double precision :: RA(29)
- double precision :: DEPJ(29)
+ ! vmod3d
+ double precision :: PNA(MPA)
+ double precision :: RNA(MRA)
+ double precision :: HNA(MHA)
+ double precision :: PNB(MPB)
+ double precision :: RNB(MRB)
+ double precision :: HNB(MHB)
+ double precision :: VELAP(MPA,MRA,MHA)
+ double precision :: VELBP(MPB,MRB,MHB)
+ ! discon
+ double precision :: PN(51)
+ double precision :: RRN(63)
+ double precision :: DEPA(51,63)
+ double precision :: DEPB(51,63)
+ double precision :: DEPC(51,63)
+ ! locate
+ double precision :: PLA
+ double precision :: RLA
+ double precision :: HLA
+ double precision :: PLB
+ double precision :: RLB
+ double precision :: HLB
+ ! weight
+ double precision :: WV(8)
+ ! prhfd
+ double precision :: P
+ double precision :: R
+ double precision :: H
+ double precision :: PF
+ double precision :: RF
+ double precision :: HF
+ double precision :: PF1
+ double precision :: RF1
+ double precision :: HF1
+ double precision :: PD
+ double precision :: RD
+ double precision :: HD
+ ! jpmodv
+ double precision :: VP(29)
+ double precision :: VS(29)
+ double precision :: RA(29)
+ double precision :: DEPJ(29)
+ ! locate integers
+ integer :: IPLOCA(MKA)
+ integer :: IRLOCA(MKA)
+ integer :: IHLOCA(MKA)
+ integer :: IPLOCB(MKB)
+ integer :: IRLOCB(MKB)
+ integer :: IHLOCB(MKB)
+ ! vmod3D integers
+ integer :: NPA
+ integer :: NRA
+ integer :: NHA
+ integer :: NPB
+ integer :: NRB
+ integer :: NHB
+ ! weight integers
+ integer :: IP
+ integer :: JP
+ integer :: KP
+ integer :: IP1
+ integer :: JP1
+ integer :: KP1
end type model_jp3d_variables
type (model_jp3d_variables) JP3DM_V
@@ -864,66 +885,69 @@
! model_jp3d_variables
type model_jp3d_variables
sequence
-! vmod3d
- integer :: NPA
- integer :: NRA
- integer :: NHA
- integer :: NPB
- integer :: NRB
- integer :: NHB
- double precision :: PNA(MPA)
- double precision :: RNA(MRA)
- double precision :: HNA(MHA)
- double precision :: PNB(MPB)
- double precision :: RNB(MRB)
- double precision :: HNB(MHB)
- double precision :: VELAP(MPA,MRA,MHA)
- double precision :: VELBP(MPB,MRB,MHB)
-! discon
- double precision :: PN(51)
- double precision :: RRN(63)
- double precision :: DEPA(51,63)
- double precision :: DEPB(51,63)
- double precision :: DEPC(51,63)
-! locate
- integer :: IPLOCA(MKA)
- integer :: IRLOCA(MKA)
- integer :: IHLOCA(MKA)
- integer :: IPLOCB(MKB)
- integer :: IRLOCB(MKB)
- integer :: IHLOCB(MKB)
- double precision :: PLA
- double precision :: RLA
- double precision :: HLA
- double precision :: PLB
- double precision :: RLB
- double precision :: HLB
-! weight
- integer :: IP
- integer :: JP
- integer :: KP
- integer :: IP1
- integer :: JP1
- integer :: KP1
- double precision :: WV(8)
-! prhfd
- double precision :: P
- double precision :: R
- double precision :: H
- double precision :: PF
- double precision :: RF
- double precision :: HF
- double precision :: PF1
- double precision :: RF1
- double precision :: HF1
- double precision :: PD
- double precision :: RD
- double precision :: HD
-! jpmodv
- double precision :: VP(29)
- double precision :: VS(29)
- double precision :: RA(29)
- double precision :: DEPJ(29)
+ ! vmod3d
+ double precision :: PNA(MPA)
+ double precision :: RNA(MRA)
+ double precision :: HNA(MHA)
+ double precision :: PNB(MPB)
+ double precision :: RNB(MRB)
+ double precision :: HNB(MHB)
+ double precision :: VELAP(MPA,MRA,MHA)
+ double precision :: VELBP(MPB,MRB,MHB)
+ ! discon
+ double precision :: PN(51)
+ double precision :: RRN(63)
+ double precision :: DEPA(51,63)
+ double precision :: DEPB(51,63)
+ double precision :: DEPC(51,63)
+ ! locate
+ double precision :: PLA
+ double precision :: RLA
+ double precision :: HLA
+ double precision :: PLB
+ double precision :: RLB
+ double precision :: HLB
+ ! weight
+ double precision :: WV(8)
+ ! prhfd
+ double precision :: P
+ double precision :: R
+ double precision :: H
+ double precision :: PF
+ double precision :: RF
+ double precision :: HF
+ double precision :: PF1
+ double precision :: RF1
+ double precision :: HF1
+ double precision :: PD
+ double precision :: RD
+ double precision :: HD
+ ! jpmodv
+ double precision :: VP(29)
+ double precision :: VS(29)
+ double precision :: RA(29)
+ double precision :: DEPJ(29)
+ ! locate integers
+ integer :: IPLOCA(MKA)
+ integer :: IRLOCA(MKA)
+ integer :: IHLOCA(MKA)
+ integer :: IPLOCB(MKB)
+ integer :: IRLOCB(MKB)
+ integer :: IHLOCB(MKB)
+ ! vmod3D integers
+ integer :: NPA
+ integer :: NRA
+ integer :: NHA
+ integer :: NPB
+ integer :: NRB
+ integer :: NHB
+ ! weight integers
+ integer :: IP
+ integer :: JP
+ integer :: KP
+ integer :: IP1
+ integer :: JP1
+ integer :: KP1
end type model_jp3d_variables
type (model_jp3d_variables) JP3DM_V
@@ -959,66 +983,69 @@
! model_jp3d_variables
type model_jp3d_variables
sequence
-! vmod3d
- integer :: NPA
- integer :: NRA
- integer :: NHA
- integer :: NPB
- integer :: NRB
- integer :: NHB
- double precision :: PNA(MPA)
- double precision :: RNA(MRA)
- double precision :: HNA(MHA)
- double precision :: PNB(MPB)
- double precision :: RNB(MRB)
- double precision :: HNB(MHB)
- double precision :: VELAP(MPA,MRA,MHA)
- double precision :: VELBP(MPB,MRB,MHB)
-! discon
- double precision :: PN(51)
- double precision :: RRN(63)
- double precision :: DEPA(51,63)
- double precision :: DEPB(51,63)
- double precision :: DEPC(51,63)
-! locate
- integer :: IPLOCA(MKA)
- integer :: IRLOCA(MKA)
- integer :: IHLOCA(MKA)
- integer :: IPLOCB(MKB)
- integer :: IRLOCB(MKB)
- integer :: IHLOCB(MKB)
- double precision :: PLA
- double precision :: RLA
- double precision :: HLA
- double precision :: PLB
- double precision :: RLB
- double precision :: HLB
-! weight
- integer :: IP
- integer :: JP
- integer :: KP
- integer :: IP1
- integer :: JP1
- integer :: KP1
- double precision :: WV(8)
-! prhfd
- double precision :: P
- double precision :: R
- double precision :: H
- double precision :: PF
- double precision :: RF
- double precision :: HF
- double precision :: PF1
- double precision :: RF1
- double precision :: HF1
- double precision :: PD
- double precision :: RD
- double precision :: HD
-! jpmodv
- double precision :: VP(29)
- double precision :: VS(29)
- double precision :: RA(29)
- double precision :: DEPJ(29)
+ ! vmod3d
+ double precision :: PNA(MPA)
+ double precision :: RNA(MRA)
+ double precision :: HNA(MHA)
+ double precision :: PNB(MPB)
+ double precision :: RNB(MRB)
+ double precision :: HNB(MHB)
+ double precision :: VELAP(MPA,MRA,MHA)
+ double precision :: VELBP(MPB,MRB,MHB)
+ ! discon
+ double precision :: PN(51)
+ double precision :: RRN(63)
+ double precision :: DEPA(51,63)
+ double precision :: DEPB(51,63)
+ double precision :: DEPC(51,63)
+ ! locate
+ double precision :: PLA
+ double precision :: RLA
+ double precision :: HLA
+ double precision :: PLB
+ double precision :: RLB
+ double precision :: HLB
+ ! weight
+ double precision :: WV(8)
+ ! prhfd
+ double precision :: P
+ double precision :: R
+ double precision :: H
+ double precision :: PF
+ double precision :: RF
+ double precision :: HF
+ double precision :: PF1
+ double precision :: RF1
+ double precision :: HF1
+ double precision :: PD
+ double precision :: RD
+ double precision :: HD
+ ! jpmodv
+ double precision :: VP(29)
+ double precision :: VS(29)
+ double precision :: RA(29)
+ double precision :: DEPJ(29)
+ ! locate integers
+ integer :: IPLOCA(MKA)
+ integer :: IRLOCA(MKA)
+ integer :: IHLOCA(MKA)
+ integer :: IPLOCB(MKB)
+ integer :: IRLOCB(MKB)
+ integer :: IHLOCB(MKB)
+ ! vmod3D integers
+ integer :: NPA
+ integer :: NRA
+ integer :: NHA
+ integer :: NPB
+ integer :: NRB
+ integer :: NHB
+ ! weight integers
+ integer :: IP
+ integer :: JP
+ integer :: KP
+ integer :: IP1
+ integer :: JP1
+ integer :: KP1
end type model_jp3d_variables
type (model_jp3d_variables) JP3DM_V
@@ -1060,66 +1087,69 @@
! model_jp3d_variables
type model_jp3d_variables
sequence
-! vmod3d
- integer :: NPA
- integer :: NRA
- integer :: NHA
- integer :: NPB
- integer :: NRB
- integer :: NHB
- double precision :: PNA(MPA)
- double precision :: RNA(MRA)
- double precision :: HNA(MHA)
- double precision :: PNB(MPB)
- double precision :: RNB(MRB)
- double precision :: HNB(MHB)
- double precision :: VELAP(MPA,MRA,MHA)
- double precision :: VELBP(MPB,MRB,MHB)
-! discon
- double precision :: PN(51)
- double precision :: RRN(63)
- double precision :: DEPA(51,63)
- double precision :: DEPB(51,63)
- double precision :: DEPC(51,63)
-! locate
- integer :: IPLOCA(MKA)
- integer :: IRLOCA(MKA)
- integer :: IHLOCA(MKA)
- integer :: IPLOCB(MKB)
- integer :: IRLOCB(MKB)
- integer :: IHLOCB(MKB)
- double precision :: PLA
- double precision :: RLA
- double precision :: HLA
- double precision :: PLB
- double precision :: RLB
- double precision :: HLB
-! weight
- integer :: IP
- integer :: JP
- integer :: KP
- integer :: IP1
- integer :: JP1
- integer :: KP1
- double precision :: WV(8)
-! prhfd
- double precision :: P
- double precision :: R
- double precision :: H
- double precision :: PF
- double precision :: RF
- double precision :: HF
- double precision :: PF1
- double precision :: RF1
- double precision :: HF1
- double precision :: PD
- double precision :: RD
- double precision :: HD
-! jpmodv
- double precision :: VP(29)
- double precision :: VS(29)
- double precision :: RA(29)
- double precision :: DEPJ(29)
+ ! vmod3d
+ double precision :: PNA(MPA)
+ double precision :: RNA(MRA)
+ double precision :: HNA(MHA)
+ double precision :: PNB(MPB)
+ double precision :: RNB(MRB)
+ double precision :: HNB(MHB)
+ double precision :: VELAP(MPA,MRA,MHA)
+ double precision :: VELBP(MPB,MRB,MHB)
+ ! discon
+ double precision :: PN(51)
+ double precision :: RRN(63)
+ double precision :: DEPA(51,63)
+ double precision :: DEPB(51,63)
+ double precision :: DEPC(51,63)
+ ! locate
+ double precision :: PLA
+ double precision :: RLA
+ double precision :: HLA
+ double precision :: PLB
+ double precision :: RLB
+ double precision :: HLB
+ ! weight
+ double precision :: WV(8)
+ ! prhfd
+ double precision :: P
+ double precision :: R
+ double precision :: H
+ double precision :: PF
+ double precision :: RF
+ double precision :: HF
+ double precision :: PF1
+ double precision :: RF1
+ double precision :: HF1
+ double precision :: PD
+ double precision :: RD
+ double precision :: HD
+ ! jpmodv
+ double precision :: VP(29)
+ double precision :: VS(29)
+ double precision :: RA(29)
+ double precision :: DEPJ(29)
+ ! locate integers
+ integer :: IPLOCA(MKA)
+ integer :: IRLOCA(MKA)
+ integer :: IHLOCA(MKA)
+ integer :: IPLOCB(MKB)
+ integer :: IRLOCB(MKB)
+ integer :: IHLOCB(MKB)
+ ! vmod3D integers
+ integer :: NPA
+ integer :: NRA
+ integer :: NHA
+ integer :: NPB
+ integer :: NRB
+ integer :: NHB
+ ! weight integers
+ integer :: IP
+ integer :: JP
+ integer :: KP
+ integer :: IP1
+ integer :: JP1
+ integer :: KP1
end type model_jp3d_variables
type (model_jp3d_variables) JP3DM_V
@@ -1180,66 +1210,69 @@
! model_jp3d_variables
type model_jp3d_variables
sequence
-! vmod3d
- integer :: NPA
- integer :: NRA
- integer :: NHA
- integer :: NPB
- integer :: NRB
- integer :: NHB
- double precision :: PNA(MPA)
- double precision :: RNA(MRA)
- double precision :: HNA(MHA)
- double precision :: PNB(MPB)
- double precision :: RNB(MRB)
- double precision :: HNB(MHB)
- double precision :: VELAP(MPA,MRA,MHA)
- double precision :: VELBP(MPB,MRB,MHB)
-! discon
- double precision :: PN(51)
- double precision :: RRN(63)
- double precision :: DEPA(51,63)
- double precision :: DEPB(51,63)
- double precision :: DEPC(51,63)
-! locate
- integer :: IPLOCA(MKA)
- integer :: IRLOCA(MKA)
- integer :: IHLOCA(MKA)
- integer :: IPLOCB(MKB)
- integer :: IRLOCB(MKB)
- integer :: IHLOCB(MKB)
- double precision :: PLA
- double precision :: RLA
- double precision :: HLA
- double precision :: PLB
- double precision :: RLB
- double precision :: HLB
-! weight
- integer :: IP
- integer :: JP
- integer :: KP
- integer :: IP1
- integer :: JP1
- integer :: KP1
- double precision :: WV(8)
-! prhfd
- double precision :: P
- double precision :: R
- double precision :: H
- double precision :: PF
- double precision :: RF
- double precision :: HF
- double precision :: PF1
- double precision :: RF1
- double precision :: HF1
- double precision :: PD
- double precision :: RD
- double precision :: HD
-! jpmodv
- double precision :: VP(29)
- double precision :: VS(29)
- double precision :: RA(29)
- double precision :: DEPJ(29)
+ ! vmod3d
+ double precision :: PNA(MPA)
+ double precision :: RNA(MRA)
+ double precision :: HNA(MHA)
+ double precision :: PNB(MPB)
+ double precision :: RNB(MRB)
+ double precision :: HNB(MHB)
+ double precision :: VELAP(MPA,MRA,MHA)
+ double precision :: VELBP(MPB,MRB,MHB)
+ ! discon
+ double precision :: PN(51)
+ double precision :: RRN(63)
+ double precision :: DEPA(51,63)
+ double precision :: DEPB(51,63)
+ double precision :: DEPC(51,63)
+ ! locate
+ double precision :: PLA
+ double precision :: RLA
+ double precision :: HLA
+ double precision :: PLB
+ double precision :: RLB
+ double precision :: HLB
+ ! weight
+ double precision :: WV(8)
+ ! prhfd
+ double precision :: P
+ double precision :: R
+ double precision :: H
+ double precision :: PF
+ double precision :: RF
+ double precision :: HF
+ double precision :: PF1
+ double precision :: RF1
+ double precision :: HF1
+ double precision :: PD
+ double precision :: RD
+ double precision :: HD
+ ! jpmodv
+ double precision :: VP(29)
+ double precision :: VS(29)
+ double precision :: RA(29)
+ double precision :: DEPJ(29)
+ ! locate integers
+ integer :: IPLOCA(MKA)
+ integer :: IRLOCA(MKA)
+ integer :: IHLOCA(MKA)
+ integer :: IPLOCB(MKB)
+ integer :: IRLOCB(MKB)
+ integer :: IHLOCB(MKB)
+ ! vmod3D integers
+ integer :: NPA
+ integer :: NRA
+ integer :: NHA
+ integer :: NPB
+ integer :: NRB
+ integer :: NHB
+ ! weight integers
+ integer :: IP
+ integer :: JP
+ integer :: KP
+ integer :: IP1
+ integer :: JP1
+ integer :: KP1
end type model_jp3d_variables
type (model_jp3d_variables) JP3DM_V
@@ -1273,66 +1306,69 @@
! model_jp3d_variables
type model_jp3d_variables
sequence
-! vmod3d
- integer :: NPA
- integer :: NRA
- integer :: NHA
- integer :: NPB
- integer :: NRB
- integer :: NHB
- double precision :: PNA(MPA)
- double precision :: RNA(MRA)
- double precision :: HNA(MHA)
- double precision :: PNB(MPB)
- double precision :: RNB(MRB)
- double precision :: HNB(MHB)
- double precision :: VELAP(MPA,MRA,MHA)
- double precision :: VELBP(MPB,MRB,MHB)
-! discon
- double precision :: PN(51)
- double precision :: RRN(63)
- double precision :: DEPA(51,63)
- double precision :: DEPB(51,63)
- double precision :: DEPC(51,63)
-! locate
- integer :: IPLOCA(MKA)
- integer :: IRLOCA(MKA)
- integer :: IHLOCA(MKA)
- integer :: IPLOCB(MKB)
- integer :: IRLOCB(MKB)
- integer :: IHLOCB(MKB)
- double precision :: PLA
- double precision :: RLA
- double precision :: HLA
- double precision :: PLB
- double precision :: RLB
- double precision :: HLB
-! weight
- integer :: IP
- integer :: JP
- integer :: KP
- integer :: IP1
- integer :: JP1
- integer :: KP1
- double precision :: WV(8)
-! prhfd
- double precision :: P
- double precision :: R
- double precision :: H
- double precision :: PF
- double precision :: RF
- double precision :: HF
- double precision :: PF1
- double precision :: RF1
- double precision :: HF1
- double precision :: PD
- double precision :: RD
- double precision :: HD
-! jpmodv
- double precision :: VP(29)
- double precision :: VS(29)
- double precision :: RA(29)
- double precision :: DEPJ(29)
+ ! vmod3d
+ double precision :: PNA(MPA)
+ double precision :: RNA(MRA)
+ double precision :: HNA(MHA)
+ double precision :: PNB(MPB)
+ double precision :: RNB(MRB)
+ double precision :: HNB(MHB)
+ double precision :: VELAP(MPA,MRA,MHA)
+ double precision :: VELBP(MPB,MRB,MHB)
+ ! discon
+ double precision :: PN(51)
+ double precision :: RRN(63)
+ double precision :: DEPA(51,63)
+ double precision :: DEPB(51,63)
+ double precision :: DEPC(51,63)
+ ! locate
+ double precision :: PLA
+ double precision :: RLA
+ double precision :: HLA
+ double precision :: PLB
+ double precision :: RLB
+ double precision :: HLB
+ ! weight
+ double precision :: WV(8)
+ ! prhfd
+ double precision :: P
+ double precision :: R
+ double precision :: H
+ double precision :: PF
+ double precision :: RF
+ double precision :: HF
+ double precision :: PF1
+ double precision :: RF1
+ double precision :: HF1
+ double precision :: PD
+ double precision :: RD
+ double precision :: HD
+ ! jpmodv
+ double precision :: VP(29)
+ double precision :: VS(29)
+ double precision :: RA(29)
+ double precision :: DEPJ(29)
+ ! locate integers
+ integer :: IPLOCA(MKA)
+ integer :: IRLOCA(MKA)
+ integer :: IHLOCA(MKA)
+ integer :: IPLOCB(MKB)
+ integer :: IRLOCB(MKB)
+ integer :: IHLOCB(MKB)
+ ! vmod3D integers
+ integer :: NPA
+ integer :: NRA
+ integer :: NHA
+ integer :: NPB
+ integer :: NRB
+ integer :: NHB
+ ! weight integers
+ integer :: IP
+ integer :: JP
+ integer :: KP
+ integer :: IP1
+ integer :: JP1
+ integer :: KP1
end type model_jp3d_variables
type (model_jp3d_variables) JP3DM_V
@@ -1370,66 +1406,69 @@
! model_jp3d_variables
type model_jp3d_variables
sequence
-! vmod3d
- integer :: NPA
- integer :: NRA
- integer :: NHA
- integer :: NPB
- integer :: NRB
- integer :: NHB
- double precision :: PNA(MPA)
- double precision :: RNA(MRA)
- double precision :: HNA(MHA)
- double precision :: PNB(MPB)
- double precision :: RNB(MRB)
- double precision :: HNB(MHB)
- double precision :: VELAP(MPA,MRA,MHA)
- double precision :: VELBP(MPB,MRB,MHB)
-! discon
- double precision :: PN(51)
- double precision :: RRN(63)
- double precision :: DEPA(51,63)
- double precision :: DEPB(51,63)
- double precision :: DEPC(51,63)
-! locate
- integer :: IPLOCA(MKA)
- integer :: IRLOCA(MKA)
- integer :: IHLOCA(MKA)
- integer :: IPLOCB(MKB)
- integer :: IRLOCB(MKB)
- integer :: IHLOCB(MKB)
- double precision :: PLA
- double precision :: RLA
- double precision :: HLA
- double precision :: PLB
- double precision :: RLB
- double precision :: HLB
-! weight
- integer :: IP
- integer :: JP
- integer :: KP
- integer :: IP1
- integer :: JP1
- integer :: KP1
- double precision :: WV(8)
-! prhfd
- double precision :: P
- double precision :: R
- double precision :: H
- double precision :: PF
- double precision :: RF
- double precision :: HF
- double precision :: PF1
- double precision :: RF1
- double precision :: HF1
- double precision :: PD
- double precision :: RD
- double precision :: HD
-! jpmodv
- double precision :: VP(29)
- double precision :: VS(29)
- double precision :: RA(29)
- double precision :: DEPJ(29)
+ ! vmod3d
+ double precision :: PNA(MPA)
+ double precision :: RNA(MRA)
+ double precision :: HNA(MHA)
+ double precision :: PNB(MPB)
+ double precision :: RNB(MRB)
+ double precision :: HNB(MHB)
+ double precision :: VELAP(MPA,MRA,MHA)
+ double precision :: VELBP(MPB,MRB,MHB)
+ ! discon
+ double precision :: PN(51)
+ double precision :: RRN(63)
+ double precision :: DEPA(51,63)
+ double precision :: DEPB(51,63)
+ double precision :: DEPC(51,63)
+ ! locate
+ double precision :: PLA
+ double precision :: RLA
+ double precision :: HLA
+ double precision :: PLB
+ double precision :: RLB
+ double precision :: HLB
+ ! weight
+ double precision :: WV(8)
+ ! prhfd
+ double precision :: P
+ double precision :: R
+ double precision :: H
+ double precision :: PF
+ double precision :: RF
+ double precision :: HF
+ double precision :: PF1
+ double precision :: RF1
+ double precision :: HF1
+ double precision :: PD
+ double precision :: RD
+ double precision :: HD
+ ! jpmodv
+ double precision :: VP(29)
+ double precision :: VS(29)
+ double precision :: RA(29)
+ double precision :: DEPJ(29)
+ ! locate integers
+ integer :: IPLOCA(MKA)
+ integer :: IRLOCA(MKA)
+ integer :: IHLOCA(MKA)
+ integer :: IPLOCB(MKB)
+ integer :: IRLOCB(MKB)
+ integer :: IHLOCB(MKB)
+ ! vmod3D integers
+ integer :: NPA
+ integer :: NRA
+ integer :: NHA
+ integer :: NPB
+ integer :: NRB
+ integer :: NHB
+ ! weight integers
+ integer :: IP
+ integer :: JP
+ integer :: KP
+ integer :: IP1
+ integer :: JP1
+ integer :: KP1
end type model_jp3d_variables
type (model_jp3d_variables) JP3DM_V
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/model_ppm.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/model_ppm.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/model_ppm.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -107,7 +107,7 @@
double precision :: maxlat,maxlon,minlat,minlon,maxdepth,mindepth
double precision :: dlat,dlon,ddepth,max_dvs,min_dvs
integer :: num_v,num_latperlon,num_lonperdepth
- integer :: dummy ! padding 4 bytes to align the structure
+ integer :: dummy_pad ! padding 4 bytes to align the structure
end type model_ppm_variables
type (model_ppm_variables) PPM_V
@@ -158,7 +158,7 @@
double precision :: maxlat,maxlon,minlat,minlon,maxdepth,mindepth
double precision :: dlat,dlon,ddepth,max_dvs,min_dvs
integer :: num_v,num_latperlon,num_lonperdepth
- integer :: dummy ! padding 4 bytes to align the structure
+ integer :: dummy_pad ! padding 4 bytes to align the structure
end type model_ppm_variables
type (model_ppm_variables) PPM_V
@@ -341,7 +341,7 @@
double precision :: maxlat,maxlon,minlat,minlon,maxdepth,mindepth
double precision :: dlat,dlon,ddepth,max_dvs,min_dvs
integer :: num_v,num_latperlon,num_lonperdepth
- integer :: dummy ! padding 4 bytes to align the structure
+ integer :: dummy_pad ! padding 4 bytes to align the structure
end type model_ppm_variables
type (model_ppm_variables) PPM_V
@@ -455,7 +455,7 @@
double precision :: maxlat,maxlon,minlat,minlon,maxdepth,mindepth
double precision :: dlat,dlon,ddepth,max_dvs,min_dvs
integer :: num_v,num_latperlon,num_lonperdepth
- integer :: dummy ! padding 4 bytes to align the structure
+ integer :: dummy_pad ! padding 4 bytes to align the structure
end type model_ppm_variables
type (model_ppm_variables) PPM_V
@@ -558,7 +558,7 @@
double precision :: maxlat,maxlon,minlat,minlon,maxdepth,mindepth
double precision :: dlat,dlon,ddepth,max_dvs,min_dvs
integer :: num_v,num_latperlon,num_lonperdepth
- integer :: dummy ! padding 4 bytes to align the structure
+ integer :: dummy_pad ! padding 4 bytes to align the structure
end type model_ppm_variables
type (model_ppm_variables) PPM_V
@@ -588,54 +588,11 @@
real(kind=CUSTOM_REAL) rho_vp(NGLLX,NGLLY,NGLLZ,nspec_stacey)
real(kind=CUSTOM_REAL) rho_vs(NGLLX,NGLLY,NGLLZ,nspec_stacey)
-! 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
-
-! attenuation
- !logical ATTENUATION,ATTENUATION_3D
- !integer vx, vy, vz, vnspec
- !double precision T_c_source
- !double precision, dimension(N_SLS) :: tau_s
- !double precision, dimension(vx, vy, vz, vnspec) :: Qmu_store
- !double precision, dimension(N_SLS, vx, vy, vz, vnspec) :: tau_e_store
-
- !integer NEX_PER_PROC_XI,NEX_PER_PROC_ETA,NEX_XI,ichunk
- !integer nglob
-
- !integer nspec_ani
- !real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec_ani) :: &
- ! c11store,c12store,c13store,c14store,c15store,c16store, &
- ! c22store,c23store,c24store,c25store,c26store,c33store,c34store, &
- ! c35store,c36store,c44store,c45store,c46store,c55store,c56store,c66store
-
- !logical TRANSVERSE_ISOTROPY,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE
- !logical OCEANS
-
! local parameters
integer i,j,k,ispec
integer iregion_code
-! only include the neighboring 3 x 3 slices
+ ! only include the neighboring 3 x 3 slices
integer, parameter :: NSLICES = 3
integer ,parameter :: NSLICES2 = NSLICES * NSLICES
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/model_sea99_s.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/model_sea99_s.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/model_sea99_s.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -54,18 +54,18 @@
! model_sea99_s_variables
type model_sea99_s_variables
sequence
- integer :: sea99_ndep
- integer :: sea99_nlat
- integer :: sea99_nlon
- integer :: dummy_pad ! padding 4 bytes to align the structure
+ double precision :: sea99_vs(100,100,100)
+ double precision :: sea99_depth(100)
double precision :: sea99_ddeg
double precision :: alatmin
double precision :: alatmax
double precision :: alonmin
double precision :: alonmax
- double precision :: sea99_vs(100,100,100)
- double precision :: sea99_depth(100)
- end type model_sea99_s_variables
+ integer :: sea99_ndep
+ integer :: sea99_nlat
+ integer :: sea99_nlon
+ integer :: dummy_pad ! padding 4 bytes to align the structure
+ end type model_sea99_s_variables
type (model_sea99_s_variables) SEA99M_V
! model_sea99_s_variables
@@ -103,18 +103,18 @@
! model_sea99_s_variables
type model_sea99_s_variables
sequence
- integer :: sea99_ndep
- integer :: sea99_nlat
- integer :: sea99_nlon
- integer :: dummy_pad ! padding 4 bytes to align the structure
+ double precision :: sea99_vs(100,100,100)
+ double precision :: sea99_depth(100)
double precision :: sea99_ddeg
double precision :: alatmin
double precision :: alatmax
double precision :: alonmin
double precision :: alonmax
- double precision :: sea99_vs(100,100,100)
- double precision :: sea99_depth(100)
- end type model_sea99_s_variables
+ integer :: sea99_ndep
+ integer :: sea99_nlat
+ integer :: sea99_nlon
+ integer :: dummy_pad ! padding 4 bytes to align the structure
+ end type model_sea99_s_variables
type (model_sea99_s_variables) SEA99M_V
! model_sea99_s_variables
@@ -167,19 +167,19 @@
! model_sea99_s_variables
type model_sea99_s_variables
- sequence
- integer :: sea99_ndep
- integer :: sea99_nlat
- integer :: sea99_nlon
- integer :: dummy_pad ! padding 4 bytes to align the structure
- double precision :: sea99_ddeg
- double precision :: alatmin
- double precision :: alatmax
- double precision :: alonmin
- double precision :: alonmax
- double precision :: sea99_vs(100,100,100)
- double precision :: sea99_depth(100)
- end type model_sea99_s_variables
+ sequence
+ double precision :: sea99_vs(100,100,100)
+ double precision :: sea99_depth(100)
+ double precision :: sea99_ddeg
+ double precision :: alatmin
+ double precision :: alatmax
+ double precision :: alonmin
+ double precision :: alonmax
+ integer :: sea99_ndep
+ integer :: sea99_nlat
+ integer :: sea99_nlon
+ integer :: dummy_pad ! padding 4 bytes to align the structure
+ end type model_sea99_s_variables
type (model_sea99_s_variables) SEA99M_V
! model_sea99_s_variables
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/prepare_timerun.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/prepare_timerun.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/prepare_timerun.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -364,8 +364,9 @@
two_omega_earth = - 2.d0 * TWO_PI / (HOURS_PER_DAY * 3600.d0 * scale_t_inv)
endif
endif
- A_array_rotation = 0.
- B_array_rotation = 0.
+
+ A_array_rotation = 0._CUSTOM_REAL
+ B_array_rotation = 0._CUSTOM_REAL
if (SIMULATION_TYPE == 3) then
if(CUSTOM_REAL == SIZE_REAL) then
@@ -709,6 +710,3 @@
endif
end subroutine prepare_timerun_attenuation
-
-
-
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/read_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/read_arrays_solver.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/read_arrays_solver.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -95,7 +95,8 @@
! create the name for the database of the current slide and region
call create_name_database(prname,myrank,iregion_code,LOCAL_PATH)
- open(unit=IIN,file=prname(1:len_trim(prname))//'solver_data_1.bin',status='old',action='read',form='unformatted')
+ open(unit=IIN,file=prname(1:len_trim(prname))//'solver_data_1.bin', &
+ status='old',action='read',form='unformatted')
read(IIN) xix
read(IIN) xiy
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/read_mesh_databases.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/read_mesh_databases.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/read_mesh_databases.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -602,7 +602,7 @@
! Stacey put back
open(unit=27,file=prname(1:len_trim(prname))//'boundary.bin', &
- status='unknown',form='unformatted',action='read')
+ status='old',form='unformatted',action='read')
read(27) nspec2D_xmin_crust_mantle
read(27) nspec2D_xmax_crust_mantle
read(27) nspec2D_ymin_crust_mantle
@@ -645,7 +645,7 @@
! Stacey put back
open(unit=27,file=prname(1:len_trim(prname))//'boundary.bin', &
- status='unknown',form='unformatted',action='read')
+ status='old',form='unformatted',action='read')
read(27) nspec2D_xmin_outer_core
read(27) nspec2D_xmax_outer_core
read(27) nspec2D_ymin_outer_core
@@ -685,7 +685,7 @@
! read info for vertical edges for central cube matching in inner core
open(unit=27,file=prname(1:len_trim(prname))//'boundary.bin', &
- status='old',action='read',form='unformatted')
+ status='old',form='unformatted',action='read')
read(27) nspec2D_xmin_inner_core
read(27) nspec2D_xmax_inner_core
read(27) nspec2D_ymin_inner_core
@@ -709,7 +709,7 @@
call create_name_database(prname,myrank,IREGION_CRUST_MANTLE,LOCAL_PATH)
open(unit=27,file=prname(1:len_trim(prname))//'boundary_disc.bin', &
- status='old',form='unformatted')
+ status='old',form='unformatted',action='read')
read(27) njunk1,njunk2,njunk3
if (njunk1 /= NSPEC2D_MOHO .and. njunk2 /= NSPEC2D_400 .and. njunk3 /= NSPEC2D_670) &
call exit_mpi(myrank, 'Error reading ibelm_disc.bin file')
@@ -801,7 +801,7 @@
! read arrays for Stacey conditions
open(unit=27,file=prname(1:len_trim(prname))//'stacey.bin', &
- status='unknown',form='unformatted',action='read')
+ status='old',form='unformatted',action='read')
read(27) nimin_crust_mantle
read(27) nimax_crust_mantle
read(27) njmin_crust_mantle
@@ -897,7 +897,7 @@
! read arrays for Stacey conditions
open(unit=27,file=prname(1:len_trim(prname))//'stacey.bin', &
- status='unknown',form='unformatted',action='read')
+ status='old',form='unformatted',action='read')
read(27) nimin_outer_core
read(27) nimax_outer_core
read(27) njmin_outer_core
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/save_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/save_arrays_solver.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/save_arrays_solver.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -58,7 +58,6 @@
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
@@ -66,6 +65,7 @@
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, dimension(:), pointer :: interval_Q ! Steps
integer :: Qn ! Number of points
integer dummy_pad ! padding 4 bytes to align the structure
end type model_attenuation_variables
@@ -327,7 +327,6 @@
enddo
enddo
enddo
-
write(27) rmass
write(27) ibool
@@ -384,15 +383,49 @@
if( SAVE_MESH_FILES ) then
scaleval1 = sngl( sqrt(PI*GRAV*RHOAV)*(R_EARTH/1000.0d0) )
scaleval2 = sngl( RHOAV/1000.0d0 )
+
+ ! isotropic model
+ ! vp
open(unit=27,file=prname(1:len_trim(prname))//'vp.bin',status='unknown',form='unformatted',action='write')
write(27) sqrt( (kappavstore+4.*muvstore/3.)/rhostore )*scaleval1
close(27)
+ ! vs
open(unit=27,file=prname(1:len_trim(prname))//'vs.bin',status='unknown',form='unformatted',action='write')
write(27) sqrt( muvstore/rhostore )*scaleval1
close(27)
+ ! rho
open(unit=27,file=prname(1:len_trim(prname))//'rho.bin',status='unknown',form='unformatted',action='write')
write(27) rhostore*scaleval2
close(27)
+
+ ! transverse isotropic model
+ if( TRANSVERSE_ISOTROPY ) then
+ ! vpv
+ open(unit=27,file=prname(1:len_trim(prname))//'vpv.bin',status='unknown',form='unformatted',action='write')
+ write(27) sqrt( (kappavstore+4.*muvstore/3.)/rhostore )*scaleval1
+ close(27)
+ ! vph
+ open(unit=27,file=prname(1:len_trim(prname))//'vph.bin',status='unknown',form='unformatted',action='write')
+ write(27) sqrt( (kappahstore+4.*muhstore/3.)/rhostore )*scaleval1
+ close(27)
+ ! vsv
+ open(unit=27,file=prname(1:len_trim(prname))//'vsv.bin',status='unknown',form='unformatted',action='write')
+ write(27) sqrt( muvstore/rhostore )*scaleval1
+ close(27)
+ ! vsh
+ open(unit=27,file=prname(1:len_trim(prname))//'vsh.bin',status='unknown',form='unformatted',action='write')
+ write(27) sqrt( muhstore/rhostore )*scaleval1
+ close(27)
+ ! rho
+ open(unit=27,file=prname(1:len_trim(prname))//'rho.bin',status='unknown',form='unformatted',action='write')
+ write(27) rhostore*scaleval2
+ close(27)
+ ! eta
+ open(unit=27,file=prname(1:len_trim(prname))//'eta.bin',status='unknown',form='unformatted',action='write')
+ write(27) eta_anisostore
+ close(27)
+ endif
+
endif
end subroutine save_arrays_solver
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/save_kernels.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/save_kernels.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/save_kernels.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -32,6 +32,8 @@
ystore_crust_mantle,zstore_crust_mantle, &
rhostore_crust_mantle,muvstore_crust_mantle, &
kappavstore_crust_mantle,ibool_crust_mantle, &
+ kappahstore_crust_mantle,muhstore_crust_mantle, &
+ eta_anisostore_crust_mantle,idoubling_crust_mantle, &
LOCAL_PATH)
implicit none
@@ -55,6 +57,11 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_ISO_MANTLE) :: &
rhostore_crust_mantle,kappavstore_crust_mantle,muvstore_crust_mantle
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_TISO_MANTLE) :: &
+ kappahstore_crust_mantle,muhstore_crust_mantle,eta_anisostore_crust_mantle
+
+ integer, dimension(NSPEC_CRUST_MANTLE) :: idoubling_crust_mantle
+
integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: ibool_crust_mantle
character(len=150) LOCAL_PATH
@@ -67,8 +74,23 @@
real(kind=CUSTOM_REAL) :: rhol,mul,kappal,rho_kl,alpha_kl,beta_kl
integer :: ispec,i,j,k,iglob
character(len=150) prname
+
+ ! transverse isotropic parameters
+ real(kind=CUSTOM_REAL), dimension(21) :: an_kl
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:),allocatable :: &
+ alphav_kl_crust_mantle,alphah_kl_crust_mantle, &
+ betav_kl_crust_mantle,betah_kl_crust_mantle, &
+ eta_kl_crust_mantle
-
+ ! bulk parameterization
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:),allocatable :: &
+ bulk_c_kl_crust_mantle,bulk_beta_kl_crust_mantle, &
+ bulk_betav_kl_crust_mantle,bulk_betah_kl_crust_mantle
+ real(kind=CUSTOM_REAL) :: A,C,F,L,N,eta
+ real(kind=CUSTOM_REAL) :: muvl,kappavl,muhl,kappahl
+ real(kind=CUSTOM_REAL) :: alphav_sq,alphah_sq,betav_sq,betah_sq,bulk_sq
+
+ ! scaling factors
scale_kl = scale_t/scale_displ * 1.d9
! For anisotropic kernels
! final unit : [s km^(-3) GPa^(-1)]
@@ -76,6 +98,28 @@
! final unit : [s km^(-3) (kg/m^3)^(-1)]
scale_kl_rho = scale_t / scale_displ / RHOAV * 1.d9
+ ! allocates temporary arrays
+ if( SAVE_TRANSVERSE_KL ) then
+ ! transverse isotropic kernel arrays for file output
+ allocate(alphav_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT), &
+ alphah_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT), &
+ betav_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT), &
+ betah_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT), &
+ eta_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT))
+
+ ! isotropic kernel arrays for file output
+ allocate(bulk_c_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT), &
+ bulk_betav_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT), &
+ bulk_betah_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT), &
+ bulk_beta_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT))
+ endif
+
+ if( .not. ANISOTROPIC_KL ) then
+ ! allocates temporary isotropic kernel arrays for file output
+ allocate(bulk_c_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT), &
+ bulk_beta_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT))
+ endif
+
! crust_mantle
do ispec = 1, NSPEC_CRUST_MANTLE
do k = 1, NGLLZ
@@ -96,26 +140,210 @@
cijkl_kl_crust_mantle(:,i,j,k,ispec) = cijkl_kl_local * scale_kl_ani
rho_kl_crust_mantle(i,j,k,ispec) = rho_kl_crust_mantle(i,j,k,ispec) * scale_kl_rho
+ ! transverse isotropic kernel calculations
+ if( SAVE_TRANSVERSE_KL ) then
+ ! note: transverse isotropic kernels are calculated for all elements
+ !
+ ! however, the factors A,C,L,N,F are based only on transverse elements
+ ! in between Moho and 220 km layer, otherwise they are taken from isotropic values
+
+ rhol = rhostore_crust_mantle(i,j,k,ispec)
+
+ ! transverse isotropic parameters from compute_force_crust_mantle.f90
+ ! C=rhovpvsq A=rhovphsq L=rhovsvsq N=rhovshsq eta=F/(A - 2 L)
+
+ ! Get A,C,F,L,N,eta from kappa,mu
+ ! element can have transverse isotropy if between d220 and Moho
+ if( .not. (TRANSVERSE_ISOTROPY_VAL .and. &
+ (idoubling_crust_mantle(ispec) == IFLAG_80_MOHO .or. &
+ idoubling_crust_mantle(ispec) == IFLAG_220_80))) then
+
+ ! layer with no transverse isotropy
+ ! A,C,L,N,F from isotropic model
+
+ mul = muvstore_crust_mantle(i,j,k,ispec)
+ kappal = kappavstore_crust_mantle(i,j,k,ispec)
+ muvl = mul
+ muhl = mul
+ kappavl = kappal
+ kappahl = kappal
+
+ A = kappal + FOUR_THIRDS * mul
+ C = A
+ L = mul
+ N = mul
+ F = kappal - 2._CUSTOM_REAL/3._CUSTOM_REAL * mul
+ eta = 1._CUSTOM_REAL
+
+ else
+
+ ! A,C,L,N,F from transverse isotropic model
+ kappavl = kappavstore_crust_mantle(i,j,k,ispec)
+ kappahl = kappahstore_crust_mantle(i,j,k,ispec)
+ muvl = muvstore_crust_mantle(i,j,k,ispec)
+ muhl = muhstore_crust_mantle(i,j,k,ispec)
+ kappal = kappavl
+
+ A = kappahl + FOUR_THIRDS * muhl
+ C = kappavl + FOUR_THIRDS * muvl
+ L = muvl
+ N = muhl
+ eta = eta_anisostore_crust_mantle(i,j,k,ispec) ! that is F / (A - 2 L)
+ F = eta * ( A - 2._CUSTOM_REAL * L )
+
+ endif
+
+ ! note: cijkl_kl_local() is fully anisotropic C_ij kernel components (non-dimensionalized)
+ ! for GLL point at (i,j,k,ispec)
+
+ ! Purpose : compute the kernels for the An coeffs (an_kl)
+ ! from the kernels for Cij (cijkl_kl_local)
+ ! At r,theta,phi fixed
+ ! kernel def : dx = kij * dcij + krho * drho
+ ! = kAn * dAn + krho * drho
+
+ ! Definition of the input array cij_kl :
+ ! cij_kl(1) = C11 ; cij_kl(2) = C12 ; cij_kl(3) = C13
+ ! cij_kl(4) = C14 ; cij_kl(5) = C15 ; cij_kl(6) = C16
+ ! cij_kl(7) = C22 ; cij_kl(8) = C23 ; cij_kl(9) = C24
+ ! cij_kl(10) = C25 ; cij_kl(11) = C26 ; cij_kl(12) = C33
+ ! cij_kl(13) = C34 ; cij_kl(14) = C35 ; cij_kl(15) = C36
+ ! cij_kl(16) = C44 ; cij_kl(17) = C45 ; cij_kl(18) = C46
+ ! cij_kl(19) = C55 ; cij_kl(20) = C56 ; cij_kl(21) = C66
+ ! where the Cij (Voigt's notation) are defined as function of
+ ! the components of the elastic tensor in spherical coordinates
+ ! by eq. (A.1) of Chen & Tromp, GJI 168 (2007)
+
+ ! From the relations giving Cij in function of An
+ ! Checked with Min Chen's results (routine build_cij)
+
+ an_kl(1) = cijkl_kl_local(1)+cijkl_kl_local(2)+cijkl_kl_local(7) !A
+ an_kl(2) = cijkl_kl_local(12) !C
+ an_kl(3) = -2*cijkl_kl_local(2)+cijkl_kl_local(21) !N
+ an_kl(4) = cijkl_kl_local(16)+cijkl_kl_local(19) !L
+ an_kl(5) = cijkl_kl_local(3)+cijkl_kl_local(8) !F
+
+ ! not used yet
+ !an_kl(6)=2*cijkl_kl_local(5)+2*cijkl_kl_local(10)+2*cijkl_kl_local(14) !Jc
+ !an_kl(7)=2*cijkl_kl_local(4)+2*cijkl_kl_local(9)+2*cijkl_kl_local(13) !Js
+ !an_kl(8)=-2*cijkl_kl_local(14) !Kc
+ !an_kl(9)=-2*cijkl_kl_local(13) !Ks
+ !an_kl(10)=-2*cijkl_kl_local(10)+cijkl_kl_local(18) !Mc
+ !an_kl(11)=2*cijkl_kl_local(4)-cijkl_kl_local(20) !Ms
+ !an_kl(12)=cijkl_kl_local(1)-cijkl_kl_local(7) !Bc
+ !an_kl(13)=-1./2.*(cijkl_kl_local(6)+cijkl_kl_local(11)) !Bs
+ !an_kl(14)=cijkl_kl_local(3)-cijkl_kl_local(8) !Hc
+ !an_kl(15)=-cijkl_kl_local(15) !Hs
+ !an_kl(16)=-cijkl_kl_local(16)+cijkl_kl_local(19) !Gc
+ !an_kl(17)=-cijkl_kl_local(17) !Gs
+ !an_kl(18)=cijkl_kl_local(5)-cijkl_kl_local(10)-cijkl_kl_local(18) !Dc
+ !an_kl(19)=cijkl_kl_local(4)-cijkl_kl_local(9)+cijkl_kl_local(20) !Ds
+ !an_kl(20)=cijkl_kl_local(1)-cijkl_kl_local(2)+cijkl_kl_local(7)-cijkl_kl_local(21) !Ec
+ !an_kl(21)=-cijkl_kl_local(6)+cijkl_kl_local(11) !Es
+
+ ! K_rho (primary kernel, for a parameterization (A,C,L,N,F,rho) )
+ rhonotprime_kl_crust_mantle(i,j,k,ispec) = rhol * rho_kl_crust_mantle(i,j,k,ispec) / scale_kl_rho
+
+ ! note: transverse isotropic kernels are calculated for ALL elements,
+ ! and not just transverse elements
+ !
+ ! note: the kernels are for relative perturbations (delta ln (m_i) = (m_i - m_0)/m_i )
+ !
+ ! Gets transverse isotropic kernels
+ ! (see Appendix B of Sieminski et al., GJI 171, 2007)
+
+ ! for parameterization: ( alpha_v, alpha_h, beta_v, beta_h, eta, rho )
+ ! K_alpha_v
+ alphav_kl_crust_mantle(i,j,k,ispec) = 2*C*an_kl(2)
+ ! K_alpha_h
+ alphah_kl_crust_mantle(i,j,k,ispec) = 2*A*an_kl(1) + 2*A*eta*an_kl(5)
+ ! K_beta_v
+ betav_kl_crust_mantle(i,j,k,ispec) = 2*L*an_kl(4) - 4*L*eta*an_kl(5)
+ ! K_beta_h
+ betah_kl_crust_mantle(i,j,k,ispec) = 2*N*an_kl(3)
+ ! K_eta
+ eta_kl_crust_mantle(i,j,k,ispec) = F*an_kl(5)
+ ! K_rhoprime (for a parameterization (alpha_v, ..., rho) )
+ rho_kl_crust_mantle(i,j,k,ispec) = A*an_kl(1) + C*an_kl(2) &
+ + N*an_kl(3) + L*an_kl(4) + F*an_kl(5) &
+ + rhonotprime_kl_crust_mantle(i,j,k,ispec)
+
+ ! write the kernel in physical units (01/05/2006)
+ rhonotprime_kl_crust_mantle(i,j,k,ispec) = - rhonotprime_kl_crust_mantle(i,j,k,ispec) * scale_kl
+
+ alphav_kl_crust_mantle(i,j,k,ispec) = - alphav_kl_crust_mantle(i,j,k,ispec) * scale_kl
+ alphah_kl_crust_mantle(i,j,k,ispec) = - alphah_kl_crust_mantle(i,j,k,ispec) * scale_kl
+ betav_kl_crust_mantle(i,j,k,ispec) = - betav_kl_crust_mantle(i,j,k,ispec) * scale_kl
+ betah_kl_crust_mantle(i,j,k,ispec) = - betah_kl_crust_mantle(i,j,k,ispec) * scale_kl
+ eta_kl_crust_mantle(i,j,k,ispec) = - eta_kl_crust_mantle(i,j,k,ispec) * scale_kl
+ rho_kl_crust_mantle(i,j,k,ispec) = - rho_kl_crust_mantle(i,j,k,ispec) * scale_kl
+
+ ! for parameterization: ( bulk, beta_v, beta_h, eta, rho )
+ ! where kappa_v = kappa_h = kappa and bulk c = sqrt( kappa / rho )
+ betav_sq = muvl / rhol
+ betah_sq = muhl / rhol
+ alphav_sq = ( kappal + FOUR_THIRDS * muvl ) / rhol
+ alphah_sq = ( kappal + FOUR_THIRDS * muhl ) / rhol
+ bulk_sq = kappal / rhol
+
+ bulk_c_kl_crust_mantle(i,j,k,ispec) = &
+ bulk_sq / alphav_sq * alphav_kl_crust_mantle(i,j,k,ispec) + &
+ bulk_sq / alphah_sq * alphah_kl_crust_mantle(i,j,k,ispec)
+
+ bulk_betah_kl_crust_mantle(i,j,k,ispec ) = &
+ betah_kl_crust_mantle(i,j,k,ispec) + &
+ FOUR_THIRDS * betah_sq / alphah_sq * alphah_kl_crust_mantle(i,j,k,ispec)
+
+ bulk_betav_kl_crust_mantle(i,j,k,ispec ) = &
+ betav_kl_crust_mantle(i,j,k,ispec) + &
+ FOUR_THIRDS * betav_sq / alphav_sq * alphav_kl_crust_mantle(i,j,k,ispec)
+ ! the rest, K_eta and K_rho are the same as above
+
+ ! to check: isotropic kernels from transverse isotropic ones
+ alpha_kl_crust_mantle(i,j,k,ispec) = alphav_kl_crust_mantle(i,j,k,ispec) &
+ + alphah_kl_crust_mantle(i,j,k,ispec)
+ beta_kl_crust_mantle(i,j,k,ispec) = betav_kl_crust_mantle(i,j,k,ispec) &
+ + betah_kl_crust_mantle(i,j,k,ispec)
+ !rho_kl_crust_mantle(i,j,k,ispec) = rhonotprime_kl_crust_mantle(i,j,k,ispec) &
+ ! + alpha_kl_crust_mantle(i,j,k,ispec) &
+ ! + beta_kl_crust_mantle(i,j,k,ispec)
+ bulk_beta_kl_crust_mantle(i,j,k,ispec) = bulk_betah_kl_crust_mantle(i,j,k,ispec ) &
+ + bulk_betav_kl_crust_mantle(i,j,k,ispec )
+
+ endif ! SAVE_TRANSVERSE_KL
+
else
+ ! isotropic kernels
+
rhol = rhostore_crust_mantle(i,j,k,ispec)
mul = muvstore_crust_mantle(i,j,k,ispec)
kappal = kappavstore_crust_mantle(i,j,k,ispec)
- ! kernel values for rho, kappa, mu
+ ! kernel values for rho, kappa, mu (primary kernel values)
rho_kl = - rhol * rho_kl_crust_mantle(i,j,k,ispec)
- alpha_kl = - kappal * alpha_kl_crust_mantle(i,j,k,ispec)
- beta_kl = - 2 * mul * beta_kl_crust_mantle(i,j,k,ispec)
+ alpha_kl = - kappal * alpha_kl_crust_mantle(i,j,k,ispec) ! note: alpha_kl corresponds to K_kappa
+ beta_kl = - 2 * mul * beta_kl_crust_mantle(i,j,k,ispec) ! note: beta_kl corresponds to K_mu
+ ! for a parameterization: (rho,mu,kappa) "primary" kernels
rhonotprime_kl_crust_mantle(i,j,k,ispec) = rho_kl * scale_kl
mu_kl_crust_mantle(i,j,k,ispec) = beta_kl * scale_kl
kappa_kl_crust_mantle(i,j,k,ispec) = alpha_kl * scale_kl
+ ! for a parameterization: (rho,alpha,beta)
! kernels rho^prime, beta, alpha
rho_kl_crust_mantle(i,j,k,ispec) = (rho_kl + alpha_kl + beta_kl) * scale_kl
- beta_kl_crust_mantle(i,j,k,ispec) = 2 * (beta_kl - FOUR_THIRDS * mul * alpha_kl / kappal) * scale_kl
- alpha_kl_crust_mantle(i,j,k,ispec) = 2 * (1 + FOUR_THIRDS * mul / kappal) * alpha_kl * scale_kl
+ beta_kl_crust_mantle(i,j,k,ispec) = &
+ 2._CUSTOM_REAL * (beta_kl - FOUR_THIRDS * mul * alpha_kl / kappal) * scale_kl
+ alpha_kl_crust_mantle(i,j,k,ispec) = &
+ 2._CUSTOM_REAL * (1 + FOUR_THIRDS * mul / kappal) * alpha_kl * scale_kl
+ ! for a parameterization: (rho,bulk, beta)
+ ! where bulk wave speed is c = sqrt( kappa / rho)
+ ! note: rhoprime is the same as for (rho,alpha,beta) parameterization
+ bulk_c_kl_crust_mantle(i,j,k,ispec) = 2._CUSTOM_REAL * alpha_kl * scale_kl
+ bulk_beta_kl_crust_mantle(i,j,k,ispec ) = 2._CUSTOM_REAL * beta_kl * scale_kl
+
endif
enddo
@@ -128,15 +356,72 @@
! For anisotropic kernels
if (ANISOTROPIC_KL) then
- open(unit=27,file=trim(prname)//'rho_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) -rho_kl_crust_mantle
- close(27)
- open(unit=27,file=trim(prname)//'cijkl_kernel.bin',status='unknown',form='unformatted',action='write')
- write(27) -cijkl_kl_crust_mantle
- close(27)
+ ! outputs transverse isotropic kernels only
+ if( SAVE_TRANSVERSE_KL ) then
+ ! transverse isotropic kernels
+ ! (alpha_v, alpha_h, beta_v, beta_h, eta, rho ) parameterization
+ open(unit=27,file=trim(prname)//'alphav_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) alphav_kl_crust_mantle
+ close(27)
+ open(unit=27,file=trim(prname)//'alphah_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) alphah_kl_crust_mantle
+ close(27)
+ open(unit=27,file=trim(prname)//'betav_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) betav_kl_crust_mantle
+ close(27)
+ open(unit=27,file=trim(prname)//'betah_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) betah_kl_crust_mantle
+ close(27)
+ open(unit=27,file=trim(prname)//'eta_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) eta_kl_crust_mantle
+ close(27)
+ open(unit=27,file=trim(prname)//'rho_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) rho_kl_crust_mantle
+ close(27)
+ ! in case one is interested in primary kernel K_rho
+ !open(unit=27,file=trim(prname)//'rhonotprime_kernel.bin',status='unknown',form='unformatted',action='write')
+ !write(27) rhonotprime_kl_crust_mantle
+ !close(27)
+
+ ! (bulk, beta_v, beta_h, eta, rho ) parameterization: K_eta and K_rho same as above
+ open(unit=27,file=trim(prname)//'bulk_c_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) bulk_c_kl_crust_mantle
+ close(27)
+ open(unit=27,file=trim(prname)//'bulk_betav_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) bulk_betav_kl_crust_mantle
+ close(27)
+ open(unit=27,file=trim(prname)//'bulk_betah_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) bulk_betah_kl_crust_mantle
+ close(27)
+
+ ! to check: isotropic kernels
+ open(unit=27,file=trim(prname)//'alpha_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) alpha_kl_crust_mantle
+ close(27)
+ open(unit=27,file=trim(prname)//'beta_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) beta_kl_crust_mantle
+ close(27)
+ open(unit=27,file=trim(prname)//'bulk_beta_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) bulk_beta_kl_crust_mantle
+ close(27)
+
+ else
+
+ ! fully anisotropic kernels
+ ! note: the C_ij and density kernels are not for relative perturbations (delta ln( m_i) = delta m_i / m_i),
+ ! but absolute perturbations (delta m_i = m_i - m_0)
+ open(unit=27,file=trim(prname)//'rho_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) - rho_kl_crust_mantle
+ close(27)
+ open(unit=27,file=trim(prname)//'cijkl_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) - cijkl_kl_crust_mantle
+ close(27)
+
+ endif
+
else
-
+ ! primary kernels: (rho,kappa,mu) parameterization
open(unit=27,file=trim(prname)//'rhonotprime_kernel.bin',status='unknown',form='unformatted',action='write')
write(27) rhonotprime_kl_crust_mantle
close(27)
@@ -147,7 +432,7 @@
write(27) mu_kl_crust_mantle
close(27)
-
+ ! (rho, alpha, beta ) parameterization
open(unit=27,file=trim(prname)//'rho_kernel.bin',status='unknown',form='unformatted',action='write')
write(27) rho_kl_crust_mantle
close(27)
@@ -158,9 +443,29 @@
write(27) beta_kl_crust_mantle
close(27)
+ ! (rho, bulk, beta ) parameterization, K_rho same as above
+ open(unit=27,file=trim(prname)//'bulk_c_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) bulk_c_kl_crust_mantle
+ close(27)
+ open(unit=27,file=trim(prname)//'bulk_beta_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) bulk_beta_kl_crust_mantle
+ close(27)
+
+
endif
-
+ ! cleans up temporary kernel arrays
+ if( SAVE_TRANSVERSE_KL ) then
+ deallocate(alphav_kl_crust_mantle,alphah_kl_crust_mantle, &
+ betav_kl_crust_mantle,betah_kl_crust_mantle, &
+ eta_kl_crust_mantle)
+ deallocate(bulk_c_kl_crust_mantle,bulk_betah_kl_crust_mantle, &
+ bulk_betav_kl_crust_mantle,bulk_beta_kl_crust_mantle)
+ endif
+ if( .not. ANISOTROPIC_KL ) then
+ deallocate(bulk_c_kl_crust_mantle,bulk_beta_kl_crust_mantle)
+ endif
+
end subroutine save_kernels_crust_mantle
!
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -1917,7 +1917,8 @@
!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*,' 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))
@@ -3321,6 +3322,8 @@
ystore_crust_mantle,zstore_crust_mantle, &
rhostore_crust_mantle,muvstore_crust_mantle, &
kappavstore_crust_mantle,ibool_crust_mantle, &
+ kappahstore_crust_mantle,muhstore_crust_mantle, &
+ eta_anisostore_crust_mantle,idoubling_crust_mantle, &
LOCAL_PATH)
!<YANGL
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/write_seismograms.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/write_seismograms.f90 2010-09-08 02:47:37 UTC (rev 17179)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/write_seismograms.f90 2010-09-08 20:38:25 UTC (rev 17180)
@@ -635,7 +635,7 @@
if (NZSEC >= 60) then
time_sec = jda*24*3600 + ho*3600 + mi*60 + int(sec+t_shift)
NZJDAY = int(time_sec/(24*3600))
- NZHOUR = int(mod(time_sec,24)/3600)
+ NZHOUR = int(mod(time_sec,24*3600)/3600)
NZMIN = int(mod(time_sec,3600)/60)
NZSEC = mod(time_sec,60)
if (NZJDAY > 365 .and. .not. is_leap_year(NZYEAR)) then
More information about the CIG-COMMITS
mailing list