Skip to content

Commit

Permalink
Merge pull request #848 from danielpeter/devel
Browse files Browse the repository at this point in the history 
adds compatibility flag (for older version comparisons); updates adjacency array setup
  • Loading branch information
@danielpeter
danielpeter authored Oct 7, 2024
2 parents ee26690 + 5a25fa1 commit ed6d376
Show file tree
Hide file tree
Showing 9 changed files with 459 additions and 341 deletions.
13 changes: 12 additions & 1 deletion setup/constants.h.in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -564,6 +564,17 @@
! be careful when changing this flag when computing classical kernel
integer, parameter :: NTSTEP_BETWEEN_COMPUTE_KERNELS = 1

!!-----------------------------------------------------------
!!
!! new version compatibility
!!
!!-----------------------------------------------------------
! for comparison against older versions
! in version 7.0: tiso was constrained to below moho, new version allows also in crust
! in version 8.0: conversions between geodetic & geocentric latitudes were always applied,
! and additional differences in 1-chunk centering & source positioning
logical, parameter :: USE_OLD_VERSION_FORMAT = .false.


!!-----------------------------------------------------------
!!
Expand Down Expand Up @@ -1024,7 +1035,7 @@
! parameters for GLL model (used for iterative model inversions)
character(len=*), parameter :: PATHNAME_GLL_modeldir = 'DATA/GLL/'

! to create a reference model based on 1D_REF but with 3D crust and 410/660 topography
! to create a reference model based on the 1D reference model but with 3D crust and 410/660 topography
logical,parameter :: USE_1D_REFERENCE = .false.

!!-- uncomment for using PREM as reference model (used in CEM inversion)
Expand Down
74 changes: 52 additions & 22 deletions src/meshfem3D/compute_element_properties.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -480,7 +480,8 @@ subroutine compute_element_tiso_flag(elem_is_tiso,elem_in_mantle,iregion_code,is
use constants, only: IMAIN,myrank, &
IFLAG_CRUST,IFLAG_220_80,IFLAG_80_MOHO,IFLAG_670_220,IFLAG_MANTLE_NORMAL,IREGION_CRUST_MANTLE, &
REFERENCE_MODEL_1DREF,REFERENCE_MODEL_1DREF, &
THREE_D_MODEL_S362WMANI,THREE_D_MODEL_SGLOBE
THREE_D_MODEL_S362WMANI,THREE_D_MODEL_SGLOBE, &
USE_OLD_VERSION_FORMAT

use meshfem_models_par, only: &
TRANSVERSE_ISOTROPY,USE_FULL_TISO_MANTLE,REFERENCE_1D_MODEL,THREE_D_MODEL
Expand Down Expand Up @@ -511,6 +512,8 @@ subroutine compute_element_tiso_flag(elem_is_tiso,elem_in_mantle,iregion_code,is
write(IMAIN,*) ' setting tiso flags in mantle model'
if (USE_FULL_TISO_MANTLE) &
write(IMAIN,*) ' using fully transverse isotopic mantle'
if (USE_OLD_VERSION_FORMAT) &
write(IMAIN,*) ' using formatting from old versions (7.0 to 8.0)'
call flush_IMAIN()
endif
endif
Expand All @@ -521,12 +524,16 @@ subroutine compute_element_tiso_flag(elem_is_tiso,elem_in_mantle,iregion_code,is
if (USE_FULL_TISO_MANTLE) then
! all elements below the actual moho will be used for transverse isotropy
! note: this will increase the computation time by ~ 45 %
if (idoubling(ispec) == IFLAG_MANTLE_NORMAL &
.or. idoubling(ispec) == IFLAG_670_220 &
.or. idoubling(ispec) == IFLAG_220_80 &
.or. idoubling(ispec) == IFLAG_80_MOHO &
.or. idoubling(ispec) == IFLAG_CRUST) then
elem_is_tiso = .true.
if (USE_OLD_VERSION_FORMAT) then
if (elem_in_mantle) elem_is_tiso = .true.
else
if (idoubling(ispec) == IFLAG_MANTLE_NORMAL &
.or. idoubling(ispec) == IFLAG_670_220 &
.or. idoubling(ispec) == IFLAG_220_80 &
.or. idoubling(ispec) == IFLAG_80_MOHO &
.or. idoubling(ispec) == IFLAG_CRUST) then
elem_is_tiso = .true.
endif
endif

! all done
Expand All @@ -545,26 +552,49 @@ subroutine compute_element_tiso_flag(elem_is_tiso,elem_in_mantle,iregion_code,is
! THREE_D_MODEL_S29EA
! THREE_D_MODEL_GLL
! which show significant transverse isotropy also below 220km depth
! assigns TI to elements in crust and mantle down to 670
if (idoubling(ispec) == IFLAG_670_220 &
.or. idoubling(ispec) == IFLAG_220_80 &
.or. idoubling(ispec) == IFLAG_80_MOHO &
.or. idoubling(ispec) == IFLAG_CRUST) then
elem_is_tiso = .true.
if (USE_OLD_VERSION_FORMAT) then
! assigns TI to elements in mantle elements just below actual moho down to 670
if (idoubling(ispec) == IFLAG_670_220 &
.or. idoubling(ispec) == IFLAG_220_80 &
.or. idoubling(ispec) == IFLAG_80_MOHO &
.or. (idoubling(ispec) == IFLAG_CRUST &
.and. elem_in_mantle) ) then
elem_is_tiso = .true.
endif
else
! assigns TI to elements in crust and mantle down to 670
if (idoubling(ispec) == IFLAG_670_220 &
.or. idoubling(ispec) == IFLAG_220_80 &
.or. idoubling(ispec) == IFLAG_80_MOHO &
.or. idoubling(ispec) == IFLAG_CRUST) then
elem_is_tiso = .true.
endif
endif

case default
! default reference models
! for example, PREM assigns transverse isotropy between Moho and 220km
! assigns TI to elements in crust and mantle elements (down to 220),
! to allow for tiso in crust and below actual moho (especially for oceanic crusts);
! the crustal models will decide if model parameters are tiso or iso
if (idoubling(ispec) == IFLAG_220_80 &
.or. idoubling(ispec) == IFLAG_80_MOHO &
.or. idoubling(ispec) == IFLAG_CRUST) then
! default case for PREM reference models:
! models use only transverse isotropy between moho and 220 km depth
elem_is_tiso = .true.
if (USE_OLD_VERSION_FORMAT) then
! assigns TI to elements in mantle elements just below actual moho
if (idoubling(ispec) == IFLAG_220_80 &
.or. idoubling(ispec) == IFLAG_80_MOHO &
.or. (idoubling(ispec) == IFLAG_CRUST &
.and. elem_in_mantle )) then
! default case for PREM reference models:
! models use only transverse isotropy between moho and 220 km depth
elem_is_tiso = .true.
endif
else
! assigns TI to elements in crust and mantle elements (down to 220),
! to allow for tiso in crust and below actual moho (especially for oceanic crusts);
! the crustal models will decide if model parameters are tiso or iso
if (idoubling(ispec) == IFLAG_220_80 &
.or. idoubling(ispec) == IFLAG_80_MOHO &
.or. idoubling(ispec) == IFLAG_CRUST) then
! default case for PREM reference models:
! models use only transverse isotropy between moho and 220 km depth
elem_is_tiso = .true.
endif
endif
end select

Expand Down
94 changes: 66 additions & 28 deletions src/meshfem3D/meshfem3D_models.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -703,7 +703,10 @@ subroutine meshfem3D_models_get3Dmntl_val(iregion_code,r_prem,rho, &
!
! note: in general, models here make use of perturbation values with respect to their
! corresponding 1-D reference models
if (MODEL_3D_MANTLE_PERTUBATIONS .and. r_prem > RCMB/R_PLANET .and. .not. suppress_mantle_extension) then
if (MODEL_3D_MANTLE_PERTUBATIONS &
.and. r_prem > RCMB/R_PLANET &
.and. .not. suppress_mantle_extension &
.and. .not. USE_1D_REFERENCE) then

! extend 3-D mantle model above the Moho to the surface before adding the crust
if (r_prem > RCMB/R_PLANET .and. r_prem < RMOHO/R_PLANET) then
Expand Down Expand Up @@ -749,6 +752,7 @@ subroutine meshfem3D_models_get3Dmntl_val(iregion_code,r_prem,rho, &

case(THREE_D_MODEL_MANTLE_SH)
! full_sh model
! gets lat/lon coordinates from geocentric theta/phi
lat = (PI/2.0d0-theta)*180.0d0/PI
lon = phi*180.0d0/PI
if (lon > 180.0d0) lon = lon - 360.0d0
Expand Down Expand Up @@ -794,20 +798,12 @@ subroutine meshfem3D_models_get3Dmntl_val(iregion_code,r_prem,rho, &
case (THREE_D_MODEL_S362ANI,THREE_D_MODEL_S362WMANI, &
THREE_D_MODEL_S362ANI_PREM,THREE_D_MODEL_S29EA)
! 3D Harvard models s362ani, s362wmani, s362ani_prem and s2.9ea
! gets colat/lon coordinates from geocentric theta/phi
xcolat = sngl(theta*180.0d0/PI)
xlon = sngl(phi*180.0d0/PI)
xrad = sngl(r_used*R_PLANET_KM)
call model_s362ani_subshsv(xcolat,xlon,xrad,xdvsh,xdvsv,xdvph,xdvpv)

! to use speed values from the 1D reference model but with 3D mesh variations
if (USE_1D_REFERENCE) then
! sets all 3D variations in the mantle to zero
xdvpv = 0.d0
xdvph = 0.d0
xdvsv = 0.d0
xdvsh = 0.d0
endif

if (TRANSVERSE_ISOTROPY) then
! tiso perturbation
vpv = vpv*(1.0d0+dble(xdvpv))
Expand Down Expand Up @@ -852,7 +848,6 @@ subroutine meshfem3D_models_get3Dmntl_val(iregion_code,r_prem,rho, &

case (THREE_D_MODEL_SGLOBE,THREE_D_MODEL_SGLOBE_ISO)
! 3D SGLOBE-rani model (Chang)

! normally mantle perturbations are taken from 24.4km (R_MOHO) up.
! we need to add the if statement for sgloberani_iso or sgloberani_aniso to take from 50km up:
if (r_prem > RCMB/R_PLANET .and. r_prem < 6321000.d0/R_PLANET) then
Expand All @@ -861,7 +856,6 @@ subroutine meshfem3D_models_get3Dmntl_val(iregion_code,r_prem,rho, &
! this will then "extend the mantle up to the surface" from 50km depth
r_used = 6321000.d0/R_PLANET
endif

call mantle_sglobe(r_used,theta,phi,vsv,vsh,dvsv,dvsh,dvp,drho)

! updates velocities from reference model
Expand Down Expand Up @@ -902,6 +896,7 @@ subroutine meshfem3D_models_get3Dmntl_val(iregion_code,r_prem,rho, &

case (THREE_D_MODEL_SPIRAL)
! SPiRaL v1.4 anisotropic model
! gets lat/lon coordinates from geocentric theta/phi
lat = (PI/2.0d0-theta)*180.0d0/PI
lon = phi*180.0d0/PI
! input: lat in range [-90,90]; lon in range [-180,180]
Expand Down Expand Up @@ -938,7 +933,8 @@ subroutine meshfem3D_models_get3Dmntl_val(iregion_code,r_prem,rho, &
! (based on density variations) on top of reference 3D model
if (HETEROGEN_3D_MANTLE &
.and. .not. suppress_mantle_extension &
.and. THREE_D_MODEL /= THREE_D_MODEL_HETEROGEN_PREM) then
.and. THREE_D_MODEL /= THREE_D_MODEL_HETEROGEN_PREM &
.and. .not. USE_1D_REFERENCE) then
! gets spherical coordinates of actual point location
r_used = r
! adds hetergeneous perturbations (isotropic)
Expand Down Expand Up @@ -1199,7 +1195,6 @@ subroutine meshfem3D_models_get3Dcrust_val(iregion_code,r,theta,phi, &

! returns velocities and density for points in 3D crustal region

use shared_parameters, only: ELLIPTICITY
use meshfem_models_par

implicit none
Expand All @@ -1218,7 +1213,6 @@ subroutine meshfem3D_models_get3Dcrust_val(iregion_code,r,theta,phi, &
double precision,intent(inout) :: moho,sediment

! local parameters
double precision :: lat,lon
double precision :: vpvc,vphc,vsvc,vshc,etac
double precision :: vpc,vsc,rhoc !vpc_eu
double precision :: c11c,c12c,c13c,c14c,c15c,c16c,c22c,c23c,c24c,c25c,c26c, &
Expand All @@ -1231,13 +1225,6 @@ subroutine meshfem3D_models_get3Dcrust_val(iregion_code,r,theta,phi, &
! for point radius smaller than deepest possible crust radius (~80 km depth)
if (r < R_DEEPEST_CRUST) return

! converts geocentric colatitude/lon (theta/phi) to geographic latitude/lon (lat/lon)
! lat/lon in degrees (range lat/lon = [-90,90] / [-180,180]
call thetaphi_2_geographic_latlon_dble(theta,phi,lat,lon,ELLIPTICITY)

! double-check lon in range [-180,180]
if (lon > 180.0d0 ) lon = lon - 360.0d0

!---
!
! ADD YOUR MODEL HERE
Expand All @@ -1260,14 +1247,14 @@ subroutine meshfem3D_models_get3Dcrust_val(iregion_code,r,theta,phi, &

if (.not. is_inside_region) then
! uses default crust outside of model region
call meshfem3D_model_crust(lat,lon,r,vpvc,vphc,vsvc,vshc,etac,rhoc,moho,sediment,found_crust,elem_in_crust,moho_only, &
call meshfem3D_model_crust(r,theta,phi,vpvc,vphc,vsvc,vshc,etac,rhoc,moho,sediment,found_crust,elem_in_crust,moho_only, &
c11c,c12c,c13c,c14c,c15c,c16c,c22c,c23c,c24c,c25c,c26c, &
c33c,c34c,c35c,c36c,c44c,c45c,c46c,c55c,c56c,c66c)
endif

case default
! default crust
call meshfem3D_model_crust(lat,lon,r,vpvc,vphc,vsvc,vshc,etac,rhoc,moho,sediment,found_crust,elem_in_crust,moho_only, &
call meshfem3D_model_crust(r,theta,phi,vpvc,vphc,vsvc,vshc,etac,rhoc,moho,sediment,found_crust,elem_in_crust,moho_only, &
c11c,c12c,c13c,c14c,c15c,c16c,c22c,c23c,c24c,c25c,c26c, &
c33c,c34c,c35c,c36c,c44c,c45c,c46c,c55c,c56c,c66c)

Expand Down Expand Up @@ -1358,20 +1345,21 @@ end subroutine meshfem3D_models_get3Dcrust_val
!


subroutine meshfem3D_model_crust(lat,lon,r,vpvc,vphc,vsvc,vshc,etac,rhoc, &
subroutine meshfem3D_model_crust(r,theta,phi,vpvc,vphc,vsvc,vshc,etac,rhoc, &
moho,sediment,found_crust,elem_in_crust,moho_only, &
c11c,c12c,c13c,c14c,c15c,c16c,c22c,c23c,c24c,c25c,c26c, &
c33c,c34c,c35c,c36c,c44c,c45c,c46c,c55c,c56c,c66c)

! returns velocity/density for default crust

use constants, only: USE_OLD_VERSION_FORMAT
use meshfem_models_par

implicit none

! lat/lon - in degrees (range lat/lon = [-90,90] / [-180,180]
! radius r - normalized by globe radius [0,1.x]
double precision,intent(in) :: lat,lon,r
! theta/phi - geocentric coordinates
double precision,intent(in) :: r,theta,phi

double precision,intent(out) :: vpvc,vphc,vsvc,vshc,etac,rhoc
double precision,intent(out) :: moho,sediment
Expand All @@ -1380,6 +1368,7 @@ subroutine meshfem3D_model_crust(lat,lon,r,vpvc,vphc,vsvc,vshc,etac,rhoc, &
double precision,intent(out) :: c11c,c12c,c13c,c14c,c15c,c16c,c22c,c23c,c24c,c25c,c26c, &
c33c,c34c,c35c,c36c,c44c,c45c,c46c,c55c,c56c,c66c
! local parameters
double precision :: lat,lon
! for isotropic crust
double precision :: vpc,vsc
double precision :: vpc_area,vsc_area,rhoc_area,moho_area,sediment_area
Expand Down Expand Up @@ -1420,7 +1409,56 @@ subroutine meshfem3D_model_crust(lat,lon,r,vpvc,vphc,vsvc,vshc,etac,rhoc, &
! ADD YOUR MODEL HERE
!
!---
! lat/lon range: [-90,90] / [-180,180]

! note: Here we assume that the crustal models are defined with respect to geographic lat/lon positions.
! This is different to mantle models, where we assume that those are referenced to a spherical Earth.
!
! Geocentric and geographic colatitude (theta) are only identical for a spherical Earth.
! Depending on the Par_file 'ELLIPTICITY' flag however, we will stretch out the spherical mesh after assigining
! the model velocities.
!
! Given we want to assign the crustal velocities at point (x/y/z geocentric) which then corresponds to the
! coordinates (lat/lon geographic) in the final, stretched out mesh, we need to calculate the (lat/lon) position
! from (theta/phi) with the ellipticity correction.
!
! again, in other words (see comment in moho_stretching.f90):
! "
! the mesh here by default starts as a spherical mesh. the position (x/y/z) is thus a geocentric position
! and we want to assign the velocity & depth of the resulting geographic (lat/lon) position after stretching.
!
! we would have two options to achieve this:
! (a) we convert first all the crustal model (lat/lon) positions to the geocentric (lat'/lon') ones
! and then search for the corresponding geocentric (lat'/lon') position here.
! (b) we calculate the geographic (lat/lon) for this geocentric point location (x/y/z) and then
! search for the geographic (lat/lon) position in the original crustal model positions.
! the implementation here follows option (b) to leave the crustal model in its original form.
!
! thus, having the geocentric (x/y/z) position and its geocentric (lat'/lon') coordinate, we account here for the
! ellipticity stretching afterwards to find the proper resulting geographic (lat/lon) coordinate
! for which the crustal model is defined.
! in case the mesh stays spherical (no ELLIPTICITY), the geocentric (lat'/lon') coordinates are identical
! to geographic (lat/lon) coordinates and no correction is applied.
!
! this assumes that the crustal models are given for geographic (lat/lon) positions and
! not corrected geocentric (lat'/lon') positions.
! "
!
! TODO: we might want to re-visit this assumption.

if (USE_OLD_VERSION_FORMAT) then
! lat/lon from geocentric position w/out geodetic correction
lat = (PI_OVER_TWO - theta) * RADIANS_TO_DEGREES
lon = phi * RADIANS_TO_DEGREES
else
! converts geocentric colatitude/lon (theta/phi) to geographic latitude/lon (lat/lon)
! lat/lon in degrees (range lat/lon = [-90,90] / [-180,180]
call thetaphi_2_geographic_latlon_dble(theta,phi,lat,lon,ELLIPTICITY)
endif

! lat/lon in degrees (range lat/lon = [-90,90] / [-180,180]
! double-check lon in range [-180,180]
if (lon > 180.d0) lon = lon - 360.0d0
if (lon < -180.d0) lon = lon + 360.0d0

select case (REFERENCE_CRUSTAL_MODEL)

Expand Down
Loading

0 comments on commit ed6d376

Please sign in to comment.