[cig-commits] [commit] devel,master: tried to use the same notations in model_crust_1_0.f90 and model_crust_2_0.f90; also cleaned useless white spaces in all source files using a script (615d1f7)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Thu Nov 6 08:16:35 PST 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branches: devel,master
Link : https://github.com/geodynamics/specfem3d_globe/compare/bc58e579b3b0838a0968725a076f5904845437ca...be63f20cbb6f462104e949894dbe205d2398cd7f
>---------------------------------------------------------------
commit 615d1f7927dcf8407b725c5be312ad1219948704
Author: Dimitri Komatitsch <komatitsch at lma.cnrs-mrs.fr>
Date: Mon May 19 00:16:17 2014 +0200
tried to use the same notations in model_crust_1_0.f90 and model_crust_2_0.f90;
also cleaned useless white spaces in all source files using a script
>---------------------------------------------------------------
615d1f7927dcf8407b725c5be312ad1219948704
src/gpu/assemble_MPI_scalar_gpu.c | 0
src/gpu/assemble_MPI_vector_gpu.c | 0
src/gpu/check_fields_gpu.c | 0
src/gpu/compute_add_sources_elastic_gpu.c | 0
src/gpu/compute_coupling_gpu.c | 0
src/gpu/compute_forces_crust_mantle_gpu.c | 0
src/gpu/compute_forces_inner_core_gpu.c | 0
src/gpu/compute_forces_outer_core_gpu.c | 0
src/gpu/compute_kernels_gpu.c | 0
src/gpu/compute_stacey_acoustic_gpu.c | 0
src/gpu/compute_stacey_elastic_gpu.c | 0
src/gpu/initialize_gpu.c | 0
src/gpu/mesh_constants_gpu.h | 0
src/gpu/mesh_constants_ocl.h | 0
src/gpu/noise_tomography_gpu.c | 0
src/gpu/prepare_mesh_constants_gpu.c | 0
src/gpu/update_displacement_gpu.c | 0
src/gpu/write_seismograms_gpu.c | 2 +-
src/meshfem3D/compute_coordinates_grid.f90 | 0
src/meshfem3D/meshfem3D_models.f90 | 14 +-
src/meshfem3D/model_crust_1_0.f90 | 129 ++++++-----
src/meshfem3D/model_crust_2_0.f90 | 348 ++++++++---------------------
src/meshfem3D/model_crustmaps.f90 | 173 ++++++++++++++
src/specfem3D/specfem3D_gpu_method_stubs.c | 0
24 files changed, 333 insertions(+), 333 deletions(-)
diff --git a/src/gpu/write_seismograms_gpu.c b/src/gpu/write_seismograms_gpu.c
index 2e0019d..78e4fdc 100644
--- a/src/gpu/write_seismograms_gpu.c
+++ b/src/gpu/write_seismograms_gpu.c
@@ -133,7 +133,7 @@ void write_seismograms_transfer_from_device (Mesh *mp,
print_CUDA_error_if_any(cudaMemcpy(mp->h_station_seismo_field,mp->d_station_seismo_field.cuda,
3*NGLL3*(mp->nrec_local)*sizeof(realw),cudaMemcpyDeviceToHost),77000);
- }
+ }
}
#endif
if (!GPU_ASYNC_COPY) {
diff --git a/src/meshfem3D/meshfem3D_models.f90 b/src/meshfem3D/meshfem3D_models.f90
index 60cb170..aa673a8 100644
--- a/src/meshfem3D/meshfem3D_models.f90
+++ b/src/meshfem3D/meshfem3D_models.f90
@@ -777,24 +777,14 @@
case (ICRUST_CRUST2)
! default
! crust 2.0
- call model_crust(lat,lon,r,vpc,vsc,rhoc,moho,found_crust,elem_in_crust)
+ call model_crust_2_0(lat,lon,r,vpc,vsc,rhoc,moho,found_crust,elem_in_crust)
case (ICRUST_CRUSTMAPS)
! general crustmaps
call model_crustmaps(lat,lon,r,vpc,vsc,rhoc,moho,found_crust,elem_in_crust)
case (ICRUST_EPCRUST)
-! call model_crust(lat,lon,r,vpc,vsc,rhoc,moho,found_crust,elem_in_crust)
- ! within EPCRUST region
-! if (lat >= EPCRUST_LAT_MIN .and. lat <= EPCRUST_LAT_MAX &
-! .and. lon >= EPCRUST_LON_MIN .and. lon<=EPCRUST_LON_MAX ) then
-! vpc=0.0d0
-! vsc=0.0d0
-! rhoc=0.0d0
-! moho=0.0d0
-! found_crust = .false.
- call model_epcrust(lat,lon,r,vpc,vsc,rhoc,moho,found_crust,elem_in_crust)
-! endif
+ call model_epcrust(lat,lon,r,vpc,vsc,rhoc,moho,found_crust,elem_in_crust)
case default
stop 'crustal model type not defined'
diff --git a/src/meshfem3D/model_crust_1_0.f90 b/src/meshfem3D/model_crust_1_0.f90
index b57112e..4e2f3d8 100644
--- a/src/meshfem3D/model_crust_1_0.f90
+++ b/src/meshfem3D/model_crust_1_0.f90
@@ -57,19 +57,16 @@
! crustal_model_constants
! crustal model parameters for crust1.0
- integer, parameter :: CRUST1_NP = 9
- integer, parameter :: CRUST1_NLO = 360
- integer, parameter :: CRUST1_NLA = 180
-
- ! cap smoothing
- integer, parameter :: NCAP_CRUST1 = 180
+ integer, parameter :: CRUST_NP = 9
+ integer, parameter :: CRUST_NLO = 360
+ integer, parameter :: CRUST_NLA = 180
! model_crust_variables
! Vp, Vs and density
- double precision, dimension(:,:,:),allocatable :: crust1_vp,crust1_vs,crust1_rho
+ double precision, dimension(:,:,:), allocatable :: crust_vp,crust_vs,crust_rho
! layer thickness
- double precision, dimension(:,:,:),allocatable :: crust1_thk
+ double precision, dimension(:,:,:), allocatable :: crust_thickness
end module model_crust_1_0_par
@@ -90,27 +87,27 @@
integer :: ier
! allocate crustal arrays
- allocate(crust1_vp(CRUST1_NP,CRUST1_NLA,CRUST1_NLO), &
- crust1_vs(CRUST1_NP,CRUST1_NLA,CRUST1_NLO), &
- crust1_rho(CRUST1_NP,CRUST1_NLA,CRUST1_NLO), &
- crust1_thk(CRUST1_NP,CRUST1_NLA,CRUST1_NLO), &
+ allocate(crust_thickness(CRUST_NP,CRUST_NLA,CRUST_NLO), &
+ crust_vp(CRUST_NP,CRUST_NLA,CRUST_NLO), &
+ crust_vs(CRUST_NP,CRUST_NLA,CRUST_NLO), &
+ crust_rho(CRUST_NP,CRUST_NLA,CRUST_NLO), &
stat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'error allocating crustal arrays')
! initializes
- crust1_vp(:,:,:) = ZERO
- crust1_vs(:,:,:) = ZERO
- crust1_rho(:,:,:) = ZERO
- crust1_thk(:,:,:) = ZERO
+ crust_vp(:,:,:) = ZERO
+ crust_vs(:,:,:) = ZERO
+ crust_rho(:,:,:) = ZERO
+ crust_thickness(:,:,:) = ZERO
! the variables read are declared and stored in structure model_crust_1_0_par
if(myrank == 0) call read_crust_1_0_model()
! broadcast the information read on the master to the nodes
- call bcast_all_dp(crust1_vp,CRUST1_NP*CRUST1_NLA*CRUST1_NLO)
- call bcast_all_dp(crust1_vs,CRUST1_NP*CRUST1_NLA*CRUST1_NLO)
- call bcast_all_dp(crust1_rho,CRUST1_NP*CRUST1_NLA*CRUST1_NLO)
- call bcast_all_dp(crust1_thk,CRUST1_NP*CRUST1_NLA*CRUST1_NLO)
+ call bcast_all_dp(crust_thickness,CRUST_NP*CRUST_NLA*CRUST_NLO)
+ call bcast_all_dp(crust_vp,CRUST_NP*CRUST_NLA*CRUST_NLO)
+ call bcast_all_dp(crust_vs,CRUST_NP*CRUST_NLA*CRUST_NLO)
+ call bcast_all_dp(crust_rho,CRUST_NP*CRUST_NLA*CRUST_NLO)
end subroutine model_crust_1_0_broadcast
@@ -128,13 +125,13 @@
double precision,intent(in) :: lat,lon,x
double precision,intent(out) :: vp,vs,rho,moho
logical,intent(out) :: found_crust
- logical,intent(in):: elem_in_crust
+ logical,intent(in) :: elem_in_crust
! local parameters
double precision :: h_sed,h_uc
double precision :: x3,x4,x5,x6,x7,x8
double precision :: scaleval
- double precision,dimension(CRUST1_NP):: vps,vss,rhos,thicks
+ double precision,dimension(CRUST_NP):: vps,vss,rhos,thicks
! initializes
vp = ZERO
@@ -173,15 +170,15 @@
! checks moho value
!moho = h_uc + thicks(6) + thicks(7)
- !if( moho /= thicks(NLAYERS_CRUST) ) then
- ! print*,'moho:',moho,thicks(NLAYERS_CRUST)
+ !if( moho /= thicks(CRUST_NP) ) then
+ ! print*,'moho:',moho,thicks(CRUST_NP)
! print*,' lat/lon/x:',lat,lon,x
!endif
- ! No matter found_crust true or false, output moho thickness
+ ! no matter found_crust true or false, output moho thickness
moho = (h_uc+thicks(7)+thicks(8)) * scaleval
- ! initializes
+ ! gets corresponding crustal velocities and density
found_crust = .true.
! gets corresponding crustal velocities and density
@@ -265,7 +262,7 @@
write(IMAIN,*)
! allocates temporary array
- allocate(bnd(CRUST1_NP,CRUST1_NLA,CRUST1_NLO),stat=ier)
+ allocate(bnd(CRUST_NP,CRUST_NLA,CRUST_NLO),stat=ier)
if( ier /= 0 ) call exit_MPI(0,'error allocating crustal arrays in read routine')
! initializes
@@ -297,12 +294,12 @@
endif
! reads in data values
- do j = 1,CRUST1_NLA
- do i = 1,CRUST1_NLO
- read(51,*)(crust1_vp(k,j,i),k=1,CRUST1_NP)
- read(52,*)(crust1_vs(k,j,i),k=1,CRUST1_NP)
- read(53,*)(crust1_rho(k,j,i),k=1,CRUST1_NP)
- read(54,*)(bnd(k,j,i),k=1,CRUST1_NP)
+ do j = 1,CRUST_NLA
+ do i = 1,CRUST_NLO
+ read(51,*)(crust_vp(k,j,i),k=1,CRUST_NP)
+ read(52,*)(crust_vs(k,j,i),k=1,CRUST_NP)
+ read(53,*)(crust_rho(k,j,i),k=1,CRUST_NP)
+ read(54,*)(bnd(k,j,i),k=1,CRUST_NP)
enddo
enddo
@@ -313,10 +310,10 @@
close(54)
! determines layer thickness
- do j = 1,CRUST1_NLA
- do i = 1,CRUST1_NLO
- do k = 1,CRUST1_NP - 1
- crust1_thk(k,j,i) = - (bnd(k+1,j,i) - bnd(k,j,i))
+ do j = 1,CRUST_NLA
+ do i = 1,CRUST_NLO
+ do k = 1,CRUST_NP - 1
+ crust_thickness(k,j,i) = - (bnd(k+1,j,i) - bnd(k,j,i))
enddo
enddo
enddo
@@ -327,8 +324,8 @@
! additional info
if( DEBUG_FILE_OUTPUT ) then
! allocates temporary arrays
- allocate(thc(CRUST1_NLA,CRUST1_NLO), &
- ths(CRUST1_NLA,CRUST1_NLO), &
+ allocate(thc(CRUST_NLA,CRUST_NLO), &
+ ths(CRUST_NLA,CRUST_NLO), &
stat=ier)
if( ier /= 0 ) call exit_MPI(0,'error allocating crustal arrays in read routine')
@@ -344,28 +341,29 @@
! crustal thickness
! thickness = ice (layer index 2) + sediment (index 3+4+5) + crystalline crust (index 6+7+8)
- do j = 1,CRUST1_NLA
- do i = 1,CRUST1_NLO
+ do j = 1,CRUST_NLA
+ do i = 1,CRUST_NLO
+
! crustal thickness with ice
- !thc(j,i) = crust1_thk(2,j,i) &
- ! + crust1_thk(3,j,i) + crust1_thk(4,j,i) + crust1_thk(5,j,i) &
- ! + crust1_thk(6,j,i) + crust1_thk(7,j,i) + crust1_thk(8,j,i)
+ !thc(j,i) = crust_thickness(2,j,i) &
+ ! + crust_thickness(3,j,i) + crust_thickness(4,j,i) + crust_thickness(5,j,i) &
+ ! + crust_thickness(6,j,i) + crust_thickness(7,j,i) + crust_thickness(8,j,i)
! sediment thickness
- ths(j,i) = crust1_thk(3,j,i) + crust1_thk(4,j,i) + crust1_thk(5,j,i)
+ ths(j,i) = crust_thickness(3,j,i) + crust_thickness(4,j,i) + crust_thickness(5,j,i)
! crustal thickness without ice
! note: etopo1 has topography including ice ("ice surface" version) and at base of ice sheets ("bedrock" version)
! see: http://www.ngdc.noaa.gov/mgg/global/global.html
- thc(j,i) = crust1_thk(3,j,i) + crust1_thk(4,j,i) + crust1_thk(5,j,i) &
- + crust1_thk(6,j,i) + crust1_thk(7,j,i) + crust1_thk(8,j,i)
+ thc(j,i) = crust_thickness(3,j,i) + crust_thickness(4,j,i) + crust_thickness(5,j,i) &
+ + crust_thickness(6,j,i) + crust_thickness(7,j,i) + crust_thickness(8,j,i)
! limit moho thickness
if(thc(j,i) > h_moho_max) h_moho_max = thc(j,i)
if(thc(j,i) < h_moho_min) h_moho_min = thc(j,i)
- write(77,*)(90.0-j+0.5),(-180.0+i-0.5),thc(j,i)+crust1_thk(2,j,i),thc(j,i),ths(j,i),crust1_thk(2,j,i)
+ write(77,*)(90.0-j+0.5),(-180.0+i-0.5),thc(j,i)+crust_thickness(2,j,i),thc(j,i),ths(j,i),crust_thickness(2,j,i)
enddo
enddo
close(77)
@@ -381,10 +379,10 @@
h_moho_max = -HUGEVAL
! smoothed version
- do j = 1,CRUST1_NLA
- lat = 90.d0 - j + 0.5
- do i = 1,CRUST1_NLO
- lon = -180.d0 + i - 0.5
+ do j = 1,CRUST_NLA
+ lat = 90.d0 - j + 0.5d0
+ do i = 1,CRUST_NLO
+ lon = -180.d0 + i - 0.5d0
x = 1.0d0
call model_crust_1_0(lat,lon,x,vp,vs,rho,moho,found_crust,.false.)
@@ -417,7 +415,7 @@
!
! crust1.0 gets smoothed with a cap of size CAP using NTHETA points
! in the theta direction and NPHI in the phi direction.
-! The cap is rotated to the North Pole.
+! The cap is first rotated to the North Pole for easier implementation.
use constants
use model_crust_1_0_par
@@ -430,7 +428,7 @@
! argument variables
double precision :: lat,lon
- double precision,dimension(CRUST1_NP) :: rho,thick,velp,vels
+ double precision,dimension(CRUST_NP) :: rho,thick,velp,vels
! work-around to avoid Jacobian problems when stretching mesh elements;
! one could also try to slightly change the shape of the doubling element bricks (which cause the problem)...
@@ -444,7 +442,7 @@
! local variables
double precision :: xlon(NTHETA*NPHI),xlat(NTHETA*NPHI),weight(NTHETA*NPHI)
- double precision :: rhol(CRUST1_NP),thickl(CRUST1_NP),velpl(CRUST1_NP),velsl(CRUST1_NP)
+ double precision :: rhol(CRUST_NP),thickl(CRUST_NP),velpl(CRUST_NP),velsl(CRUST_NP)
double precision :: weightl,cap_degree
double precision :: dist
@@ -480,7 +478,7 @@
endif
endif
- ! gets smoothing points and weights (see routine in model_crust_2_0.f90)
+ ! gets smoothing points and weights
call CAP_vardegree(lon,lat,xlon,xlat,weight,cap_degree,NTHETA,NPHI)
! initializes
@@ -548,7 +546,6 @@
end subroutine crust_1_0_CAPsmoothed
-
!
!-------------------------------------------------------------------------------------------------
!
@@ -560,19 +557,19 @@
implicit none
! argument variables
- integer,intent(in) :: icolat,ilon
- double precision,intent(out) :: rhtyp(CRUST1_NP),thtp(CRUST1_NP)
- double precision,intent(out) :: vptyp(CRUST1_NP),vstyp(CRUST1_NP)
+ integer, intent(in) :: icolat,ilon
+ double precision, intent(out) :: rhtyp(CRUST_NP),thtp(CRUST_NP)
+ double precision, intent(out) :: vptyp(CRUST_NP),vstyp(CRUST_NP)
- ! sets vp,vs and rho for all layers
- vptyp(:) = crust1_vp(:,icolat,ilon)
- vstyp(:) = crust1_vs(:,icolat,ilon)
- rhtyp(:) = crust1_rho(:,icolat,ilon)
- thtp(:) = crust1_thk(:,icolat,ilon)
+ ! set vp,vs and rho for all layers
+ vptyp(:) = crust_vp(:,icolat,ilon)
+ vstyp(:) = crust_vs(:,icolat,ilon)
+ rhtyp(:) = crust_rho(:,icolat,ilon)
+ thtp(:) = crust_thickness(:,icolat,ilon)
! get distance to Moho from the bottom of the ocean or the ice
! value could be used for checking, but is unused so far...
- thtp(CRUST1_NP) = thtp(CRUST1_NP) - thtp(1) - thtp(2)
+ thtp(CRUST_NP) = thtp(CRUST_NP) - thtp(1) - thtp(2)
end subroutine get_crust_1_0_structure
diff --git a/src/meshfem3D/model_crust_2_0.f90 b/src/meshfem3D/model_crust_2_0.f90
index 35bd78d..dbfff89 100644
--- a/src/meshfem3D/model_crust_2_0.f90
+++ b/src/meshfem3D/model_crust_2_0.f90
@@ -50,14 +50,18 @@
! crustal_model_constants
! crustal model parameters for crust2.0
- integer, parameter :: NKEYS_CRUST = 359
- integer, parameter :: NLAYERS_CRUST = 8
- integer, parameter :: NCAP_CRUST = 180
+ integer, parameter :: CRUST_NP = 8
+ integer, parameter :: CRUST_NLO = 359
+ integer, parameter :: CRUST_NLA = 180
! model_crust_variables
- double precision, dimension(:,:),allocatable :: thlr,velocp,velocs,dens
- character(len=2) :: abbreviation(NCAP_CRUST/2,NCAP_CRUST)
- character(len=2) :: code(NKEYS_CRUST)
+ ! Vp, Vs and density
+ double precision, dimension(:,:), allocatable :: crust_vp,crust_vs,crust_rho
+ character(len=2) :: abbreviation(CRUST_NLA/2,CRUST_NLA)
+ character(len=2) :: code(CRUST_NLO)
+
+ ! layer thickness
+ double precision, dimension(:,:), allocatable :: crust_thickness
end module model_crust_2_0_par
@@ -78,24 +82,30 @@
integer :: ier
! allocate crustal arrays
- allocate(thlr(NLAYERS_CRUST,NKEYS_CRUST), &
- velocp(NLAYERS_CRUST,NKEYS_CRUST), &
- velocs(NLAYERS_CRUST,NKEYS_CRUST), &
- dens(NLAYERS_CRUST,NKEYS_CRUST), &
+ allocate(crust_thickness(CRUST_NP,CRUST_NLO), &
+ crust_vp(CRUST_NP,CRUST_NLO), &
+ crust_vs(CRUST_NP,CRUST_NLO), &
+ crust_rho(CRUST_NP,CRUST_NLO), &
stat=ier)
if( ier /= 0 ) call exit_MPI(myrank,'error allocating crustal arrays')
+ ! initializes
+ crust_vp(:,:) = ZERO
+ crust_vs(:,:) = ZERO
+ crust_rho(:,:) = ZERO
+ crust_thickness(:,:) = ZERO
+
! the variables read are declared and stored in structure model_crust_2_0_par
- if(myrank == 0) call read_crust_model()
+ if(myrank == 0) call read_crust_2_0_model()
! broadcast the information read on the master to the nodes
- call bcast_all_dp(thlr,NKEYS_CRUST*NLAYERS_CRUST)
- call bcast_all_dp(velocp,NKEYS_CRUST*NLAYERS_CRUST)
- call bcast_all_dp(velocs,NKEYS_CRUST*NLAYERS_CRUST)
- call bcast_all_dp(dens,NKEYS_CRUST*NLAYERS_CRUST)
+ call bcast_all_dp(crust_thickness,CRUST_NLO*CRUST_NP)
+ call bcast_all_dp(crust_vp,CRUST_NLO*CRUST_NP)
+ call bcast_all_dp(crust_vs,CRUST_NLO*CRUST_NP)
+ call bcast_all_dp(crust_rho,CRUST_NLO*CRUST_NP)
- call bcast_all_ch_array2(abbreviation,NCAP_CRUST/2,NCAP_CRUST,2)
- call bcast_all_ch_array(code,NKEYS_CRUST,2)
+ call bcast_all_ch_array2(abbreviation,CRUST_NLA/2,CRUST_NLA,2)
+ call bcast_all_ch_array(code,CRUST_NLO,2)
end subroutine model_crust_2_0_broadcast
@@ -103,7 +113,7 @@
!-------------------------------------------------------------------------------------------------
!
- subroutine model_crust(lat,lon,x,vp,vs,rho,moho,found_crust,elem_in_crust)
+ subroutine model_crust_2_0(lat,lon,x,vp,vs,rho,moho,found_crust,elem_in_crust)
use constants
use model_crust_2_0_par
@@ -117,21 +127,25 @@
! local parameters
double precision :: h_sed,h_uc
- double precision :: x3,x4,x5,x6,x7,scaleval
- double precision,dimension(NLAYERS_CRUST):: vps,vss,rhos,thicks
+ double precision :: x3,x4,x5,x6,x7
+ double precision :: scaleval
+ double precision,dimension(CRUST_NP):: vps,vss,rhos,thicks
! initializes
vp = ZERO
vs = ZERO
rho = ZERO
+ moho = ZERO
- ! gets smoothed crust2.0 structure
- call crust_CAPsmoothed(lat,lon,vps,vss,rhos,thicks,abbreviation, &
- code,thlr,velocp,velocs,dens)
+ ! gets smoothed structure
+ call crust_2_0_CAPsmoothed(lat,lon,vps,vss,rhos,thicks,abbreviation, &
+ code,crust_thickness,crust_vp,crust_vs,crust_rho)
+ ! non-dimensionalization factor
scaleval = ONE / R_EARTH_KM
! non-dimensionalizes thickness (given in km)
+ ! upper sediment
x3 = ONE - thicks(3) * scaleval
h_sed = thicks(3) + thicks(4)
x4 = ONE - h_sed * scaleval
@@ -142,24 +156,22 @@
! checks moho value
!moho = h_uc + thicks(6) + thicks(7)
- !if( moho /= thicks(NLAYERS_CRUST) ) then
- ! print*,'moho:',moho,thicks(NLAYERS_CRUST)
+ !if( moho /= thicks(CRUST_NP) ) then
+ ! print*,'moho:',moho,thicks(CRUST_NP)
! print*,' lat/lon/x:',lat,lon,x
!endif
- ! No matter found_crust true or false, output moho thickness
+ ! no matter found_crust true or false, output moho thickness
moho = (h_uc+thicks(6)+thicks(7)) * scaleval
! gets corresponding crustal velocities and density
found_crust = .true.
-! if(x > x3 .and. INCLUDE_SEDIMENTS_CRUST &
-! .and. h_sed >= MINIMUM_SEDIMENT_THICKNESS) then
+
+ ! gets corresponding crustal velocities and density
if(x > x3 .and. INCLUDE_SEDIMENTS_CRUST ) then
vp = vps(3)
vs = vss(3)
rho = rhos(3)
-! else if(x > x4 .and. INCLUDE_SEDIMENTS_CRUST &
-! .and. h_sed >= MINIMUM_SEDIMENT_THICKNESS) then
else if(x > x4 .and. INCLUDE_SEDIMENTS_CRUST ) then
vp = vps(4)
vs = vss(4)
@@ -195,13 +207,13 @@
rho = rho * 1000.0d0 / RHOAV
endif
- end subroutine model_crust
+ end subroutine model_crust_2_0
!
!-------------------------------------------------------------------------------------------------
!
- subroutine read_crust_model()
+ subroutine read_crust_2_0_model()
use constants
use model_crust_2_0_par
@@ -226,11 +238,11 @@
write(IMAIN,*) 'error opening "', trim(CNtype2), '": ', ier
call flush_IMAIN()
! stop
- call exit_MPI(0, 'error model crust2.0')
+ call exit_MPI(0,'error model crust2.0')
endif
- do ila=1,NCAP_CRUST/2
- read(1,*) icolat,(abbreviation(ila,i),i=1,NCAP_CRUST)
+ do ila=1,CRUST_NLA/2
+ read(1,*) icolat,(abbreviation(ila,i),i=1,CRUST_NLA)
enddo
close(1)
@@ -238,44 +250,43 @@
open(unit=1,file=CNtype2_key_modif,status='old',action='read',iostat=ier)
if ( ier /= 0 ) then
write(IMAIN,*) 'error opening "', trim(CNtype2_key_modif), '": ', ier
- call exit_MPI(0, 'error model crust2.0')
+ call exit_MPI(0,'error model crust2.0')
endif
h_moho_min = HUGEVAL
h_moho_max = -HUGEVAL
- do ikey=1,NKEYS_CRUST
+ do ikey=1,CRUST_NLO
read (1,"(a2)") code(ikey)
- read (1,*) (velocp(i,ikey),i=1,NLAYERS_CRUST)
- read (1,*) (velocs(i,ikey),i=1,NLAYERS_CRUST)
- read (1,*) (dens(i,ikey),i=1,NLAYERS_CRUST)
- read (1,*) (thlr(i,ikey),i=1,NLAYERS_CRUST-1),thlr(NLAYERS_CRUST,ikey)
+ read (1,*) (crust_vp(i,ikey),i=1,CRUST_NP)
+ read (1,*) (crust_vs(i,ikey),i=1,CRUST_NP)
+ read (1,*) (crust_rho(i,ikey),i=1,CRUST_NP)
+ read (1,*) (crust_thickness(i,ikey),i=1,CRUST_NP-1),crust_thickness(CRUST_NP,ikey)
! limit moho thickness
- if(thlr(NLAYERS_CRUST,ikey) > h_moho_max) h_moho_max = thlr(NLAYERS_CRUST,ikey)
- if(thlr(NLAYERS_CRUST,ikey) < h_moho_min) h_moho_min = thlr(NLAYERS_CRUST,ikey)
+ if(crust_thickness(CRUST_NP,ikey) > h_moho_max) h_moho_max = crust_thickness(CRUST_NP,ikey)
+ if(crust_thickness(CRUST_NP,ikey) < h_moho_min) h_moho_min = crust_thickness(CRUST_NP,ikey)
enddo
close(1)
- if(h_moho_min == HUGEVAL .or. h_moho_max == -HUGEVAL) &
- stop 'incorrect moho depths in read_crust_model'
+ if(h_moho_min == HUGEVAL .or. h_moho_max == -HUGEVAL) stop 'incorrect moho depths in read_crust_2_0_model'
- end subroutine read_crust_model
+ end subroutine read_crust_2_0_model
!
!-------------------------------------------------------------------------------------------------
!
- subroutine crust_CAPsmoothed(lat,lon,velp,vels,rho,thick,abbreviation,code,thlr,velocp,velocs,dens)
+ subroutine crust_2_0_CAPsmoothed(lat,lon,velp,vels,rho,thick,abbreviation,code,crust_thickness,crust_vp,crust_vs,crust_rho)
! crustal vp and vs in km/s, layer thickness in km
!
-! crust2.0 is smoothed with a cap of size CAP using NTHETA points
+! crust2.0 gets smoothed with a cap of size CAP using NTHETA points
! in the theta direction and NPHI in the phi direction.
! The cap is first rotated to the North Pole for easier implementation.
use constants
- use model_crust_2_0_par,only: NLAYERS_CRUST,NKEYS_CRUST,NCAP_CRUST
+ use model_crust_2_0_par,only: CRUST_NP,CRUST_NLO,CRUST_NLA
implicit none
@@ -285,11 +296,11 @@
! argument variables
double precision :: lat,lon
- double precision,dimension(NLAYERS_CRUST) :: rho,thick,velp,vels
- double precision,dimension(NLAYERS_CRUST,NKEYS_CRUST) :: thlr,velocp,velocs,dens
+ double precision,dimension(CRUST_NP) :: rho,thick,velp,vels
+ double precision,dimension(CRUST_NP,CRUST_NLO) :: crust_thickness,crust_vp,crust_vs,crust_rho
- character(len=2) :: code(NKEYS_CRUST)
- character(len=2) :: abbreviation(NCAP_CRUST/2,NCAP_CRUST)
+ character(len=2) :: code(CRUST_NLO)
+ character(len=2) :: abbreviation(CRUST_NLA/2,CRUST_NLA)
!-------------------------------
! work-around to avoid Jacobian problems when stretching mesh elements;
@@ -304,9 +315,10 @@
! local variables
double precision :: xlon(NTHETA*NPHI),xlat(NTHETA*NPHI),weight(NTHETA*NPHI)
- double precision :: rhol(NLAYERS_CRUST),thickl(NLAYERS_CRUST),velpl(NLAYERS_CRUST),velsl(NLAYERS_CRUST)
+ double precision :: rhol(CRUST_NP),thickl(CRUST_NP),velpl(CRUST_NP),velsl(CRUST_NP)
- double precision :: weightl,cap_degree,dist
+ double precision :: weightl,cap_degree
+ double precision :: dist
double precision :: h_sed
integer :: i,icolat,ilon
character(len=2) :: crustaltype
@@ -317,15 +329,17 @@
! fill in the hash table
crustalhash_to_key = -1
- do i=1,NKEYS_CRUST
+ do i=1,CRUST_NLO
call hash_crustal_type(code(i), ihash)
- if (crustalhash_to_key(ihash) /= -1) stop 'error in crust_CAPsmoothed: hash table collision'
+ if (crustalhash_to_key(ihash) /= -1) stop 'error in crust_2_0_CAPsmoothed: hash table collision'
crustalhash_to_key(ihash) = i
enddo
! checks latitude/longitude
- if(lat > 90.0d0 .or. lat < -90.0d0 .or. lon > 180.0d0 .or. lon < -180.0d0) &
- stop 'error in latitude/longitude range in crust'
+ if(lat > 90.0d0 .or. lat < -90.0d0 .or. lon > 180.0d0 .or. lon < -180.0d0) then
+ print*,'error in lat/lon:',lat,lon
+ stop 'error in latitude/longitude range in crust2.0'
+ endif
! makes sure lat/lon are within crust2.0 range
if(lat==90.0d0) lat=89.9999d0
@@ -361,7 +375,7 @@
! loops over weight points
do i=1,NTHETA*NPHI
- ! gets crust values
+ ! gets lat/lon indices
call icolat_ilon(xlat(i),xlon(i),icolat,ilon)
crustaltype = abbreviation(icolat,ilon)
@@ -370,8 +384,9 @@
crustalkey = crustalhash_to_key(ihash)
if(crustalkey == -1) stop 'error in retrieving crust type key'
- call get_crust_structure(crustalkey,velpl,velsl,rhol,thickl, &
- thlr,velocp,velocs,dens)
+ ! gets crust values
+ call get_crust_2_0_structure(crustalkey,velpl,velsl,rhol,thickl, &
+ crust_thickness,crust_vp,crust_vs,crust_rho)
! sediment thickness
h_sed = thickl(3) + thickl(4)
@@ -396,214 +411,39 @@
vels(:) = vels(:) + weightl*velsl(:)
enddo
- end subroutine crust_CAPsmoothed
-
-!
-!-------------------------------------------------------------------------------------------------
-!
-
-! hash table to define the crustal type using an integer instead of characters
-! because working with integers in the rest of the routines results in much faster code
- subroutine hash_crustal_type(crustaltype, ihash)
-
- implicit none
-
- character(len=2), intent(in) :: crustaltype
- integer, intent(out) :: ihash
-
- ihash = iachar(crustaltype(1:1)) + 128*iachar(crustaltype(2:2)) + 1
-
- end subroutine hash_crustal_type
-
-!
-!-------------------------------------------------------------------------------------------------
-!
-
- subroutine icolat_ilon(xlat,xlon,icolat,ilon)
-
- implicit none
-
-! argument variables
- double precision :: xlat,xlon
- integer :: icolat,ilon
-
- if(xlat > 90.0d0 .or. xlat < -90.0d0 .or. xlon > 180.0d0 .or. xlon < -180.0d0) &
- stop 'error in latitude/longitude range in icolat_ilon'
-
- icolat = int(1+( (90.d0-xlat)*0.5d0 ))
- if(icolat == 91) icolat = 90
-
- ilon = int(1+( (180.d0+xlon)*0.5d0 ))
- if(ilon == 181) ilon = 1
-
- if(icolat>90 .or. icolat<1) stop 'error in routine icolat_ilon'
- if(ilon<1 .or. ilon>180) stop 'error in routine icolat_ilon'
-
- end subroutine icolat_ilon
+ end subroutine crust_2_0_CAPsmoothed
!
!-------------------------------------------------------------------------------------------------
!
- subroutine get_crust_structure(ikey,vptyp,vstyp,rhtyp,thtp,thlr,velocp,velocs,dens)
+ subroutine get_crust_2_0_structure(ikey,vptyp,vstyp,rhtyp,thtp,crust_thickness,crust_vp,crust_vs,crust_rho)
- use model_crust_2_0_par,only: NLAYERS_CRUST,NKEYS_CRUST
+ use model_crust_2_0_par,only: CRUST_NP,CRUST_NLO
implicit none
! argument variables
- double precision :: rhtyp(NLAYERS_CRUST),thtp(NLAYERS_CRUST)
- double precision :: vptyp(NLAYERS_CRUST),vstyp(NLAYERS_CRUST)
- double precision :: thlr(NLAYERS_CRUST,NKEYS_CRUST),velocp(NLAYERS_CRUST,NKEYS_CRUST)
- double precision :: velocs(NLAYERS_CRUST,NKEYS_CRUST),dens(NLAYERS_CRUST,NKEYS_CRUST)
- integer :: ikey
+ integer, intent(in) :: ikey
+ double precision, intent(out) :: rhtyp(CRUST_NP),thtp(CRUST_NP)
+ double precision, intent(out) :: vptyp(CRUST_NP),vstyp(CRUST_NP)
+ double precision :: crust_thickness(CRUST_NP,CRUST_NLO),crust_vp(CRUST_NP,CRUST_NLO)
+ double precision :: crust_vs(CRUST_NP,CRUST_NLO),crust_rho(CRUST_NP,CRUST_NLO)
! local variables
integer :: i
- do i=1,NLAYERS_CRUST
- vptyp(i)=velocp(i,ikey)
- vstyp(i)=velocs(i,ikey)
- rhtyp(i)=dens(i,ikey)
- enddo
-
- do i=1,NLAYERS_CRUST-1
- thtp(i)=thlr(i,ikey)
+ ! set vp,vs and rho for all layers
+ do i=1,CRUST_NP
+ vptyp(i)=crust_vp(i,ikey)
+ vstyp(i)=crust_vs(i,ikey)
+ rhtyp(i)=crust_rho(i,ikey)
+ thtp(i)=crust_thickness(i,ikey)
enddo
! get distance to Moho from the bottom of the ocean or the ice
- thtp(NLAYERS_CRUST) = thlr(NLAYERS_CRUST,ikey) - thtp(1) - thtp(2)
-
- end subroutine get_crust_structure
-
-
-!
-!-------------------------------------------------------------------------------------------------
-!
-
- subroutine CAP_vardegree(lon,lat,xlon,xlat,weight,CAP_DEGREE,NTHETA,NPHI)
-
-! calculates weighting points to smooth around lon/lat location with
-! a smoothing range of CAP_DEGREE
-!
-! The cap is rotated to the North Pole.
-!
-! returns: xlon,xlat,weight
-
- use constants
-
- implicit none
-
- ! sampling rate
- integer :: NTHETA
- integer :: NPHI
-
- ! smoothing size (in degrees)
- double precision :: CAP_DEGREE
-
- ! argument variables
- double precision :: lat,lon
- double precision,dimension(NTHETA*NPHI) :: xlon,xlat,weight
-
- ! local variables
- double precision :: CAP
- double precision :: theta,phi,sint,cost,sinp,cosp,wght,total
- double precision :: r_rot,theta_rot,phi_rot
- double precision :: rotation_matrix(3,3),x(3),xc(3)
- double precision :: dtheta,dphi,cap_area,dweight,pi_over_nphi
- integer :: i,j,k,itheta,iphi
-
- ! initializes
- xlon(:) = ZERO
- xlat(:) = ZERO
- weight(:) = ZERO
-
- ! checks cap degree size
- if( CAP_DEGREE < TINYVAL ) then
- ! no cap smoothing
- print*,'error cap:',CAP_DEGREE
- print*,' lat/lon:',lat,lon
- stop 'error cap_degree too small'
- endif
-
- ! pre-compute parameters
- CAP = CAP_DEGREE * DEGREES_TO_RADIANS
- dtheta = 0.5d0 * CAP / dble(NTHETA)
- dphi = TWO_PI / dble(NPHI)
-
- cap_area = TWO_PI * ( ONE - dcos(CAP) )
- dweight = CAP / dble(NTHETA) * dphi / cap_area
- pi_over_nphi = PI/dble(NPHI)
-
- ! colatitude/longitude in radian
- theta = (90.0d0 - lat ) * DEGREES_TO_RADIANS
- phi = lon * DEGREES_TO_RADIANS
-
- sint = dsin(theta)
- cost = dcos(theta)
- sinp = dsin(phi)
- cosp = dcos(phi)
-
- ! set up rotation matrix to go from cap at North pole to cap around point of interest
- rotation_matrix(1,1) = cosp*cost
- rotation_matrix(1,2) = -sinp
- rotation_matrix(1,3) = cosp*sint
- rotation_matrix(2,1) = sinp*cost
- rotation_matrix(2,2) = cosp
- rotation_matrix(2,3) = sinp*sint
- rotation_matrix(3,1) = -sint
- rotation_matrix(3,2) = ZERO
- rotation_matrix(3,3) = cost
-
- ! calculates points over a cap at the North pole and rotates them to specified lat/lon point
- i = 0
- total = ZERO
- do itheta = 1,NTHETA
-
- theta = dble(2*itheta-1)*dtheta
- cost = dcos(theta)
- sint = dsin(theta)
- wght = sint*dweight
-
- do iphi = 1,NPHI
-
- i = i+1
-
- ! get the weight associated with this integration point (same for all phi)
- weight(i) = wght
-
- total = total + weight(i)
- phi = dble(2*iphi-1)*pi_over_nphi
- cosp = dcos(phi)
- sinp = dsin(phi)
-
- ! x,y,z coordinates of integration point in cap at North pole
- xc(1) = sint*cosp
- xc(2) = sint*sinp
- xc(3) = cost
-
- ! get x,y,z coordinates in cap around point of interest
- do j=1,3
- x(j) = ZERO
- do k=1,3
- x(j) = x(j)+rotation_matrix(j,k)*xc(k)
- enddo
- enddo
-
- ! get latitude and longitude (degrees) of integration point
- call xyz_2_rthetaphi_dble(x(1),x(2),x(3),r_rot,theta_rot,phi_rot)
- call reduce(theta_rot,phi_rot)
- xlat(i) = (PI_OVER_TWO - theta_rot) * RADIANS_TO_DEGREES
- xlon(i) = phi_rot * RADIANS_TO_DEGREES
- if(xlon(i) > 180.0d0) xlon(i) = xlon(i) - 360.0d0
-
- enddo
-
- enddo
- if(abs(total - ONE) > 0.001d0) then
- print*,'error cap:',total,CAP_DEGREE
- stop 'error in cap integration for variable degree'
- endif
+ ! value could be used for checking, but is unused so far...
+ thtp(CRUST_NP) = thtp(CRUST_NP) - thtp(1) - thtp(2)
- end subroutine CAP_vardegree
+ end subroutine get_crust_2_0_structure
diff --git a/src/meshfem3D/model_crustmaps.f90 b/src/meshfem3D/model_crustmaps.f90
index 6935dbb..f365415 100644
--- a/src/meshfem3D/model_crustmaps.f90
+++ b/src/meshfem3D/model_crustmaps.f90
@@ -552,3 +552,176 @@
end subroutine ibilinearmap
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+! hash table to define the crustal type using an integer instead of characters
+! because working with integers in the rest of the routines results in much faster code
+ subroutine hash_crustal_type(crustaltype, ihash)
+
+ implicit none
+
+ character(len=2), intent(in) :: crustaltype
+ integer, intent(out) :: ihash
+
+ ihash = iachar(crustaltype(1:1)) + 128*iachar(crustaltype(2:2)) + 1
+
+ end subroutine hash_crustal_type
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine icolat_ilon(xlat,xlon,icolat,ilon)
+
+ implicit none
+
+! argument variables
+ double precision :: xlat,xlon
+ integer :: icolat,ilon
+
+ if(xlat > 90.0d0 .or. xlat < -90.0d0 .or. xlon > 180.0d0 .or. xlon < -180.0d0) &
+ stop 'error in latitude/longitude range in icolat_ilon'
+
+ icolat = int(1+( (90.d0-xlat)*0.5d0 ))
+ if(icolat == 91) icolat = 90
+
+ ilon = int(1+( (180.d0+xlon)*0.5d0 ))
+ if(ilon == 181) ilon = 1
+
+ if(icolat>90 .or. icolat<1) stop 'error in routine icolat_ilon'
+ if(ilon<1 .or. ilon>180) stop 'error in routine icolat_ilon'
+
+ end subroutine icolat_ilon
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine CAP_vardegree(lon,lat,xlon,xlat,weight,CAP_DEGREE,NTHETA,NPHI)
+
+! calculates weighting points to smooth around lon/lat location with
+! a smoothing range of CAP_DEGREE
+!
+! The cap is rotated to the North Pole.
+!
+! returns: xlon,xlat,weight
+
+ use constants
+
+ implicit none
+
+ ! sampling rate
+ integer :: NTHETA
+ integer :: NPHI
+
+ ! smoothing size (in degrees)
+ double precision :: CAP_DEGREE
+
+ ! argument variables
+ double precision :: lat,lon
+ double precision,dimension(NTHETA*NPHI) :: xlon,xlat,weight
+
+ ! local variables
+ double precision :: CAP
+ double precision :: theta,phi,sint,cost,sinp,cosp,wght,total
+ double precision :: r_rot,theta_rot,phi_rot
+ double precision :: rotation_matrix(3,3),x(3),xc(3)
+ double precision :: dtheta,dphi,cap_area,dweight,pi_over_nphi
+ integer :: i,j,k,itheta,iphi
+
+ ! initializes
+ xlon(:) = ZERO
+ xlat(:) = ZERO
+ weight(:) = ZERO
+
+ ! checks cap degree size
+ if( CAP_DEGREE < TINYVAL ) then
+ ! no cap smoothing
+ print*,'error cap:',CAP_DEGREE
+ print*,' lat/lon:',lat,lon
+ stop 'error cap_degree too small'
+ endif
+
+ ! pre-compute parameters
+ CAP = CAP_DEGREE * DEGREES_TO_RADIANS
+ dtheta = 0.5d0 * CAP / dble(NTHETA)
+ dphi = TWO_PI / dble(NPHI)
+
+ cap_area = TWO_PI * ( ONE - dcos(CAP) )
+ dweight = CAP / dble(NTHETA) * dphi / cap_area
+ pi_over_nphi = PI/dble(NPHI)
+
+ ! colatitude/longitude in radian
+ theta = (90.0d0 - lat ) * DEGREES_TO_RADIANS
+ phi = lon * DEGREES_TO_RADIANS
+
+ sint = dsin(theta)
+ cost = dcos(theta)
+ sinp = dsin(phi)
+ cosp = dcos(phi)
+
+ ! set up rotation matrix to go from cap at North pole to cap around point of interest
+ rotation_matrix(1,1) = cosp*cost
+ rotation_matrix(1,2) = -sinp
+ rotation_matrix(1,3) = cosp*sint
+ rotation_matrix(2,1) = sinp*cost
+ rotation_matrix(2,2) = cosp
+ rotation_matrix(2,3) = sinp*sint
+ rotation_matrix(3,1) = -sint
+ rotation_matrix(3,2) = ZERO
+ rotation_matrix(3,3) = cost
+
+ ! calculates points over a cap at the North pole and rotates them to specified lat/lon point
+ i = 0
+ total = ZERO
+ do itheta = 1,NTHETA
+
+ theta = dble(2*itheta-1)*dtheta
+ cost = dcos(theta)
+ sint = dsin(theta)
+ wght = sint*dweight
+
+ do iphi = 1,NPHI
+
+ i = i+1
+
+ ! get the weight associated with this integration point (same for all phi)
+ weight(i) = wght
+
+ total = total + weight(i)
+ phi = dble(2*iphi-1)*pi_over_nphi
+ cosp = dcos(phi)
+ sinp = dsin(phi)
+
+ ! x,y,z coordinates of integration point in cap at North pole
+ xc(1) = sint*cosp
+ xc(2) = sint*sinp
+ xc(3) = cost
+
+ ! get x,y,z coordinates in cap around point of interest
+ do j=1,3
+ x(j) = ZERO
+ do k=1,3
+ x(j) = x(j)+rotation_matrix(j,k)*xc(k)
+ enddo
+ enddo
+
+ ! get latitude and longitude (degrees) of integration point
+ call xyz_2_rthetaphi_dble(x(1),x(2),x(3),r_rot,theta_rot,phi_rot)
+ call reduce(theta_rot,phi_rot)
+ xlat(i) = (PI_OVER_TWO - theta_rot) * RADIANS_TO_DEGREES
+ xlon(i) = phi_rot * RADIANS_TO_DEGREES
+ if(xlon(i) > 180.0d0) xlon(i) = xlon(i) - 360.0d0
+
+ enddo
+
+ enddo
+ if(abs(total - ONE) > 0.001d0) then
+ print*,'error cap:',total,CAP_DEGREE
+ stop 'error in cap integration for variable degree'
+ endif
+
+ end subroutine CAP_vardegree
+
More information about the CIG-COMMITS
mailing list