[cig-commits] [commit] devel: adds flags in constants.h to suppress mesh stretching for 3D moho/410/660 surfaces (ad2d0fe)
cig_noreply at geodynamics.org
cig_noreply at geodynamics.org
Tue Nov 25 06:56:23 PST 2014
Repository : https://github.com/geodynamics/specfem3d_globe
On branch : devel
Link : https://github.com/geodynamics/specfem3d_globe/compare/0e5b55c6f30be94583639fd325373eecd6facc6d...8be3e0b0267c8d4cf5af3bc26e8903da17bc4fd1
>---------------------------------------------------------------
commit ad2d0fe151101b8b4014a2ced93d75b29ab616c6
Author: daniel peter <peterda at ethz.ch>
Date: Thu Nov 13 09:42:33 2014 +0100
adds flags in constants.h to suppress mesh stretching for 3D moho/410/660 surfaces
>---------------------------------------------------------------
ad2d0fe151101b8b4014a2ced93d75b29ab616c6
setup/constants.h.in | 8 +-
src/meshfem3D/add_topography_410_650.f90 | 116 +++++++++++++-------------
src/meshfem3D/compute_element_properties.f90 | 55 +++++++------
src/meshfem3D/moho_stretching.f90 | 118 ++++++++++++++-------------
src/meshfem3D/setup_model.f90 | 46 +++++++++++
5 files changed, 201 insertions(+), 142 deletions(-)
diff --git a/setup/constants.h.in b/setup/constants.h.in
index 2721a43..7d048b2 100644
--- a/setup/constants.h.in
+++ b/setup/constants.h.in
@@ -299,6 +299,13 @@
! the mesher) and use them for the computation of boundary kernel (in the solver)
logical, parameter :: SAVE_BOUNDARY_MESH = .false.
+! to suppress element stretching for 3D moho surface
+ logical,parameter :: SUPPRESS_MOHO_STRETCHING = .false.
+
+! to suppress element stretching at 410/660 internal topography
+! (i.e. creates mesh without 410/660 topography for Harvard model (s362ani,..))
+ logical,parameter :: SUPPRESS_INTERNAL_TOPOGRAPHY = .false.
+
!!-----------------------------------------------------------
!!
!! GPU optimization
@@ -764,7 +771,6 @@
! integer, parameter :: GLL_REFERENCE_1D_MODEL = REFERENCE_MODEL_1DREF
! integer, parameter :: GLL_REFERENCE_MODEL = THREE_D_MODEL_S29EA
-
!!-----------------------------------------------------------
!!
!! crustal stretching
diff --git a/src/meshfem3D/add_topography_410_650.f90 b/src/meshfem3D/add_topography_410_650.f90
index c4b779e..31967f3 100644
--- a/src/meshfem3D/add_topography_410_650.f90
+++ b/src/meshfem3D/add_topography_410_650.f90
@@ -186,64 +186,64 @@
! we loop on all GLL points of the element
do k = 1,NGLLZ
- do j = 1,NGLLY
- do i = 1,NGLLX
-
- x = xstore(i,j,k,ispec)
- y = ystore(i,j,k,ispec)
- z = zstore(i,j,k,ispec)
-
- if (USE_OLD_VERSION_5_1_5_FORMAT) then
- ! convert to r theta phi
- call xyz_2_rthetaphi_dble(x,y,z,r,theta,phi)
- call reduce(theta,phi)
- ! get colatitude and longitude in degrees
- xcolat = sngl(theta*RADIANS_TO_DEGREES)
- xlon = sngl(phi*RADIANS_TO_DEGREES)
- else
- ! converts geocentric coordinates x/y/z to geographic radius/latitude/longitude (in degrees)
- call xyz_2_rlatlon_dble(x,y,z,r,lat,lon)
-
- ! converts lat/lon to (real) colat/lon in degrees
- xcolat = sngl(90.d0 - lat) ! colatitude
- xlon = sngl(lon)
- endif
-
- ! stretching occurs between 220 and 770
- if (r > R220/R_EARTH .or. r < R771/R_EARTH) cycle
-
- ! compute topography on 410 and 650 at current point
- call model_s362ani_subtopo(xcolat,xlon,topo410out,topo650out)
-
- ! non-dimensionalize the topography, which is in km
- ! positive for a depression, so change the sign for a perturbation in radius
- topo410 = -dble(topo410out) / R_EARTH_KM
- topo650 = -dble(topo650out) / R_EARTH_KM
-
- gamma = 0.d0
- if (r >= R400/R_EARTH .and. r <= R220/R_EARTH) then
- ! stretching between R220 and R400
- gamma = (R220/R_EARTH - r) / (R220/R_EARTH - R400/R_EARTH)
- xstore(i,j,k,ispec) = x*(ONE + gamma * topo410 / r)
- ystore(i,j,k,ispec) = y*(ONE + gamma * topo410 / r)
- zstore(i,j,k,ispec) = z*(ONE + gamma * topo410 / r)
- else if (r>= R771/R_EARTH .and. r <= R670/R_EARTH) then
- ! stretching between R771 and R670
- gamma = (r - R771/R_EARTH) / (R670/R_EARTH - R771/R_EARTH)
- xstore(i,j,k,ispec) = x*(ONE + gamma * topo650 / r)
- ystore(i,j,k,ispec) = y*(ONE + gamma * topo650 / r)
- zstore(i,j,k,ispec) = z*(ONE + gamma * topo650 / r)
- else if (r > R670/R_EARTH .and. r < R400/R_EARTH) then
- ! stretching between R670 and R400
- gamma = (R400/R_EARTH - r) / (R400/R_EARTH - R670/R_EARTH)
- xstore(i,j,k,ispec) = x*(ONE + (topo410 + gamma * (topo650 - topo410)) / r)
- ystore(i,j,k,ispec) = y*(ONE + (topo410 + gamma * (topo650 - topo410)) / r)
- zstore(i,j,k,ispec) = z*(ONE + (topo410 + gamma * (topo650 - topo410)) / r)
- endif
- if (gamma < -0.0001 .or. gamma > 1.0001) call exit_MPI(myrank,'incorrect value of gamma for 410-650 topography')
-
- enddo
- enddo
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
+ x = xstore(i,j,k,ispec)
+ y = ystore(i,j,k,ispec)
+ z = zstore(i,j,k,ispec)
+
+ if (USE_OLD_VERSION_5_1_5_FORMAT) then
+ ! convert to r theta phi
+ call xyz_2_rthetaphi_dble(x,y,z,r,theta,phi)
+ call reduce(theta,phi)
+ ! get colatitude and longitude in degrees
+ xcolat = sngl(theta*RADIANS_TO_DEGREES)
+ xlon = sngl(phi*RADIANS_TO_DEGREES)
+ else
+ ! converts geocentric coordinates x/y/z to geographic radius/latitude/longitude (in degrees)
+ call xyz_2_rlatlon_dble(x,y,z,r,lat,lon)
+
+ ! converts lat/lon to (real) colat/lon in degrees
+ xcolat = sngl(90.d0 - lat) ! colatitude
+ xlon = sngl(lon)
+ endif
+
+ ! stretching occurs between 220 and 770
+ if (r > R220/R_EARTH .or. r < R771/R_EARTH) cycle
+
+ ! compute topography on 410 and 650 at current point
+ call model_s362ani_subtopo(xcolat,xlon,topo410out,topo650out)
+
+ ! non-dimensionalize the topography, which is in km
+ ! positive for a depression, so change the sign for a perturbation in radius
+ topo410 = -dble(topo410out) / R_EARTH_KM
+ topo650 = -dble(topo650out) / R_EARTH_KM
+
+ gamma = 0.d0
+ if (r >= R400/R_EARTH .and. r <= R220/R_EARTH) then
+ ! stretching between R220 and R400
+ gamma = (R220/R_EARTH - r) / (R220/R_EARTH - R400/R_EARTH)
+ xstore(i,j,k,ispec) = x*(ONE + gamma * topo410 / r)
+ ystore(i,j,k,ispec) = y*(ONE + gamma * topo410 / r)
+ zstore(i,j,k,ispec) = z*(ONE + gamma * topo410 / r)
+ else if (r>= R771/R_EARTH .and. r <= R670/R_EARTH) then
+ ! stretching between R771 and R670
+ gamma = (r - R771/R_EARTH) / (R670/R_EARTH - R771/R_EARTH)
+ xstore(i,j,k,ispec) = x*(ONE + gamma * topo650 / r)
+ ystore(i,j,k,ispec) = y*(ONE + gamma * topo650 / r)
+ zstore(i,j,k,ispec) = z*(ONE + gamma * topo650 / r)
+ else if (r > R670/R_EARTH .and. r < R400/R_EARTH) then
+ ! stretching between R670 and R400
+ gamma = (R400/R_EARTH - r) / (R400/R_EARTH - R670/R_EARTH)
+ xstore(i,j,k,ispec) = x*(ONE + (topo410 + gamma * (topo650 - topo410)) / r)
+ ystore(i,j,k,ispec) = y*(ONE + (topo410 + gamma * (topo650 - topo410)) / r)
+ zstore(i,j,k,ispec) = z*(ONE + (topo410 + gamma * (topo650 - topo410)) / r)
+ endif
+ if (gamma < -0.0001 .or. gamma > 1.0001) call exit_MPI(myrank,'incorrect value of gamma for 410-650 topography')
+
+ enddo
+ enddo
enddo
end subroutine add_topography_410_650_gll
diff --git a/src/meshfem3D/compute_element_properties.f90 b/src/meshfem3D/compute_element_properties.f90
index 05f2fab..6a766ff 100644
--- a/src/meshfem3D/compute_element_properties.f90
+++ b/src/meshfem3D/compute_element_properties.f90
@@ -190,7 +190,7 @@
! interpolates and stores GLL point locations
call compute_element_GLL_locations(xelm,yelm,zelm,ispec,nspec, &
- xstore,ystore,zstore,shape3D)
+ xstore,ystore,zstore,shape3D)
! computes velocity/density/... values for the chosen Earth model
! (only needed for second meshing phase)
@@ -231,35 +231,38 @@
endif
! adds topography on 410 km and 650 km discontinuity in model S362ANI
- if (THREE_D_MODEL == THREE_D_MODEL_S362ANI .or. THREE_D_MODEL == THREE_D_MODEL_S362WMANI &
- .or. THREE_D_MODEL == THREE_D_MODEL_S362ANI_PREM .or. THREE_D_MODEL == THREE_D_MODEL_S29EA) then
- ! stretching between 220 and 770
- if (idoubling(ispec) == IFLAG_670_220 .or. &
- idoubling(ispec) == IFLAG_MANTLE_NORMAL) then
- if (USE_GLL) then
- ! stretches every GLL point accordingly
- call add_topography_410_650_gll(myrank,xstore,ystore,zstore,ispec,nspec)
- else
- ! stretches anchor points only, interpolates GLL points later on
- call add_topography_410_650(myrank,xelm,yelm,zelm)
+ if (.not. SUPPRESS_INTERNAL_TOPOGRAPHY) then
+ if (THREE_D_MODEL == THREE_D_MODEL_S362ANI .or. THREE_D_MODEL == THREE_D_MODEL_S362WMANI &
+ .or. THREE_D_MODEL == THREE_D_MODEL_S362ANI_PREM .or. THREE_D_MODEL == THREE_D_MODEL_S29EA) then
+ ! stretching between 220 and 770
+ if (idoubling(ispec) == IFLAG_670_220 .or. &
+ idoubling(ispec) == IFLAG_MANTLE_NORMAL) then
+ if (USE_GLL) then
+ ! stretches every GLL point accordingly
+ call add_topography_410_650_gll(myrank,xstore,ystore,zstore,ispec,nspec)
+ else
+ ! stretches anchor points only, interpolates GLL points later on
+ call add_topography_410_650(myrank,xelm,yelm,zelm)
+ endif
endif
endif
+
+ ! these are placeholders:
+ ! their corresponding subroutines subtopo_cmb() and subtopo_icb() are not implemented yet....
+ ! must be done/supplied by the user; uncomment in case
+ ! CMB topography
+ ! if (THREE_D_MODEL == THREE_D_MODEL_S362ANI .and. (idoubling(ispec)==IFLAG_MANTLE_NORMAL &
+ ! .or. idoubling(ispec)==IFLAG_OUTER_CORE_NORMAL)) &
+ ! call add_topography_cmb(myrank,xelm,yelm,zelm,RTOPDDOUBLEPRIME,RCMB)
+
+ ! ICB topography
+ ! if (THREE_D_MODEL == THREE_D_MODEL_S362ANI .and. (idoubling(ispec)==IFLAG_OUTER_CORE_NORMAL &
+ ! .or. idoubling(ispec)==IFLAG_INNER_CORE_NORMAL .or. idoubling(ispec)==IFLAG_MIDDLE_CENTRAL_CUBE &
+ ! .or. idoubling(ispec)==IFLAG_BOTTOM_CENTRAL_CUBE .or. idoubling(ispec)==IFLAG_TOP_CENTRAL_CUBE &
+ ! .or. idoubling(ispec)==IFLAG_IN_FICTITIOUS_CUBE)) &
+ ! call add_topography_icb(myrank,xelm,yelm,zelm,RICB,RCMB)
endif
- ! these are placeholders:
- ! their corresponding subroutines subtopo_cmb() and subtopo_icb() are not implemented yet....
- ! must be done/supplied by the user; uncomment in case
- ! CMB topography
- ! if (THREE_D_MODEL == THREE_D_MODEL_S362ANI .and. (idoubling(ispec)==IFLAG_MANTLE_NORMAL &
- ! .or. idoubling(ispec)==IFLAG_OUTER_CORE_NORMAL)) &
- ! call add_topography_cmb(myrank,xelm,yelm,zelm,RTOPDDOUBLEPRIME,RCMB)
-
- ! ICB topography
- ! if (THREE_D_MODEL == THREE_D_MODEL_S362ANI .and. (idoubling(ispec)==IFLAG_OUTER_CORE_NORMAL &
- ! .or. idoubling(ispec)==IFLAG_INNER_CORE_NORMAL .or. idoubling(ispec)==IFLAG_MIDDLE_CENTRAL_CUBE &
- ! .or. idoubling(ispec)==IFLAG_BOTTOM_CENTRAL_CUBE .or. idoubling(ispec)==IFLAG_TOP_CENTRAL_CUBE &
- ! .or. idoubling(ispec)==IFLAG_IN_FICTITIOUS_CUBE)) &
- ! call add_topography_icb(myrank,xelm,yelm,zelm,RICB,RCMB)
! make the Earth elliptical
if (ELLIPTICITY) then
diff --git a/src/meshfem3D/moho_stretching.f90 b/src/meshfem3D/moho_stretching.f90
index 7ba55f8..af68037 100644
--- a/src/meshfem3D/moho_stretching.f90
+++ b/src/meshfem3D/moho_stretching.f90
@@ -34,12 +34,13 @@
use constants,only: &
NGNOD,R_EARTH_KM,R_EARTH,R_UNIT_SPHERE, &
- PI_OVER_TWO,RADIANS_TO_DEGREES,TINYVAL,SMALLVAL,ONE,USE_OLD_VERSION_5_1_5_FORMAT
+ PI_OVER_TWO,RADIANS_TO_DEGREES,TINYVAL,SMALLVAL,ONE,USE_OLD_VERSION_5_1_5_FORMAT, &
+ SUPPRESS_MOHO_STRETCHING
use meshfem3D_par,only: &
RMOHO_FICTITIOUS_IN_MESHER,R220,RMIDDLE_CRUST
- use meshfem3D_models_par,only: &
+ use meshfem3D_par,only: &
TOPOGRAPHY
implicit none
@@ -125,69 +126,68 @@
! note: we will honor the moho only, if the moho depth is below R_moho (~35km)
! or above R_middlecrust (~15km). otherwise, the moho will be "interpolated"
! within the element
-
- if (TOPOGRAPHY) then
- ! globe surface honors topography, elements stretched for moho
+ if (.not. SUPPRESS_MOHO_STRETCHING) then
+ ! globe surface must honor topography for elements to be stretched for moho
!
! note: if no topography is honored, stretching may lead to distorted elements and invalid Jacobian
-
- if (moho < R_moho) then
- ! actual moho below fictitious moho
- ! elements in second layer will stretch down to honor moho topography
-
- elevation = moho - R_moho
-
- if (r >= R_moho) then
- ! point above fictitious moho
- ! gamma ranges from 0 (point at surface) to 1 (point at fictitious moho depth)
- gamma = (( R_UNIT_SPHERE - r )/( R_UNIT_SPHERE - R_moho ))
- else
- ! point below fictitious moho
- ! gamma ranges from 0 (point at R220) to 1 (point at fictitious moho depth)
- gamma = (( r - R220/R_EARTH)/( R_moho - R220/R_EARTH))
-
- ! since not all GLL points are exactly at R220, use a small
- ! tolerance for R220 detection, fix R220
- if (abs(gamma) < SMALLVAL) then
- gamma = 0.0d0
+ if (TOPOGRAPHY) then
+ if (moho < R_moho) then
+ ! actual moho below fictitious moho
+ ! elements in second layer will stretch down to honor moho topography
+
+ elevation = moho - R_moho
+
+ if (r >= R_moho) then
+ ! point above fictitious moho
+ ! gamma ranges from 0 (point at surface) to 1 (point at fictitious moho depth)
+ gamma = (( R_UNIT_SPHERE - r )/( R_UNIT_SPHERE - R_moho ))
+ else
+ ! point below fictitious moho
+ ! gamma ranges from 0 (point at R220) to 1 (point at fictitious moho depth)
+ gamma = (( r - R220/R_EARTH)/( R_moho - R220/R_EARTH))
+
+ ! since not all GLL points are exactly at R220, use a small
+ ! tolerance for R220 detection, fix R220
+ if (abs(gamma) < SMALLVAL) then
+ gamma = 0.0d0
+ endif
endif
- endif
- if (gamma < -0.0001d0 .or. gamma > 1.0001d0) &
- call exit_MPI(myrank,'incorrect value of gamma for moho from crust 2.0')
+ if (gamma < -0.0001d0 .or. gamma > 1.0001d0) &
+ call exit_MPI(myrank,'incorrect value of gamma for moho from crust 2.0')
- call move_point(ia,xelm,yelm,zelm,x,y,z,gamma,elevation,r)
+ call move_point(ia,xelm,yelm,zelm,x,y,z,gamma,elevation,r)
- else if (moho > R_middlecrust) then
- ! moho above middle crust
- ! elements in first layer will squeeze into crust above moho
+ else if (moho > R_middlecrust) then
+ ! moho above middle crust
+ ! elements in first layer will squeeze into crust above moho
- elevation = moho - R_middlecrust
+ elevation = moho - R_middlecrust
- if (r > R_middlecrust) then
- ! point above middle crust
- ! gamma ranges from 0 (point at surface) to 1 (point at middle crust depth)
- gamma = (R_UNIT_SPHERE-r)/(R_UNIT_SPHERE - R_middlecrust )
- else
- ! point below middle crust
- ! gamma ranges from 0 (point at R220) to 1 (point at middle crust depth)
- gamma = (r - R220/R_EARTH)/( R_middlecrust - R220/R_EARTH )
+ if (r > R_middlecrust) then
+ ! point above middle crust
+ ! gamma ranges from 0 (point at surface) to 1 (point at middle crust depth)
+ gamma = (R_UNIT_SPHERE-r)/(R_UNIT_SPHERE - R_middlecrust )
+ else
+ ! point below middle crust
+ ! gamma ranges from 0 (point at R220) to 1 (point at middle crust depth)
+ gamma = (r - R220/R_EARTH)/( R_middlecrust - R220/R_EARTH )
- ! since not all GLL points are exactly at R220, use a small
- ! tolerance for R220 detection, fix R220
- if (abs(gamma) < SMALLVAL) then
- gamma = 0.0d0
+ ! since not all GLL points are exactly at R220, use a small
+ ! tolerance for R220 detection, fix R220
+ if (abs(gamma) < SMALLVAL) then
+ gamma = 0.0d0
+ endif
endif
- endif
-
- if (gamma < -0.0001d0 .or. gamma > 1.0001d0) &
- call exit_MPI(myrank,'incorrect value of gamma for moho from crust 2.0')
- call move_point(ia,xelm,yelm,zelm,x,y,z,gamma,elevation,r)
+ if (gamma < -0.0001d0 .or. gamma > 1.0001d0) &
+ call exit_MPI(myrank,'incorrect value of gamma for moho from crust 2.0')
- endif
+ call move_point(ia,xelm,yelm,zelm,x,y,z,gamma,elevation,r)
- endif ! TOPOGRAPHY
+ endif
+ endif ! TOPOGRAPHY
+ endif
! counts corners in above moho
! note: uses a small tolerance
@@ -236,7 +236,8 @@
use constants,only: &
NGNOD,R_EARTH_KM,R_EARTH,R_UNIT_SPHERE, &
- PI_OVER_TWO,RADIANS_TO_DEGREES,TINYVAL,SMALLVAL,ONE,HONOR_DEEP_MOHO,USE_OLD_VERSION_5_1_5_FORMAT
+ PI_OVER_TWO,RADIANS_TO_DEGREES,TINYVAL,SMALLVAL,ONE,HONOR_DEEP_MOHO,USE_OLD_VERSION_5_1_5_FORMAT, &
+ SUPPRESS_MOHO_STRETCHING
use meshfem3D_par,only: &
R220
@@ -304,10 +305,13 @@
! - between 25km and 45km
! - below 60 km (in HONOR_DEEP_MOHO case)
! otherwise, the moho will be "interpolated" within the element
- if (HONOR_DEEP_MOHO) then
- call stretch_deep_moho(ia,xelm,yelm,zelm,x,y,z,r,moho)
- else
- call stretch_moho(ia,xelm,yelm,zelm,x,y,z,r,moho)
+ if (.not. SUPPRESS_MOHO_STRETCHING) then
+ ! distinguish between regions with very deep moho (e.g. Himalayan)
+ if (HONOR_DEEP_MOHO) then
+ call stretch_deep_moho(ia,xelm,yelm,zelm,x,y,z,r,moho)
+ else
+ call stretch_moho(ia,xelm,yelm,zelm,x,y,z,r,moho)
+ endif
endif
! counts corners in above moho
diff --git a/src/meshfem3D/setup_model.f90 b/src/meshfem3D/setup_model.f90
index e7d1ff2..9334e6f 100644
--- a/src/meshfem3D/setup_model.f90
+++ b/src/meshfem3D/setup_model.f90
@@ -59,6 +59,51 @@
R80,R220,R670,RCMB,RICB, &
LOCAL_PATH)
+ ! infos about additional mesh optimizations
+ if (myrank == 0) then
+ write(IMAIN,*)
+ write(IMAIN,*) 'additional mesh optimizations'
+ write(IMAIN,*)
+
+ ! moho stretching
+ write(IMAIN,*) 'moho:'
+ if (CRUSTAL .and. CASE_3D) then
+ if (REGIONAL_MOHO_MESH) then
+ write(IMAIN,*) ' enforcing regional 3-layer crust'
+ if (SUPPRESS_MOHO_STRETCHING) then
+ write(IMAIN,*) ' no element stretching for 3-D moho surface'
+ else
+ if (HONOR_DEEP_MOHO) then
+ write(IMAIN,*) ' incorporating element stretching for 3-D moho surface with deep moho'
+ else
+ write(IMAIN,*) ' incorporating element stretching for 3-D moho surface'
+ endif
+ endif
+ else
+ write(IMAIN,'(a,i1,a)') ' default ',NER_CRUST,'-layer crust'
+ if (SUPPRESS_MOHO_STRETCHING .or. (.not. TOPOGRAPHY)) then
+ write(IMAIN,*) ' no element stretching for 3-D moho surface'
+ else
+ write(IMAIN,*) ' incorporating element stretching for 3-D moho surface'
+ endif
+ endif
+ else
+ write(IMAIN,*) ' no element stretching for 3-D moho surface'
+ endif
+ write(IMAIN,*)
+
+ ! internal topography
+ write(IMAIN,*) 'internal topography 410/660:'
+ if ((.not. SUPPRESS_INTERNAL_TOPOGRAPHY) .and. &
+ (THREE_D_MODEL == THREE_D_MODEL_S362ANI .or. THREE_D_MODEL == THREE_D_MODEL_S362WMANI &
+ .or. THREE_D_MODEL == THREE_D_MODEL_S362ANI_PREM .or. THREE_D_MODEL == THREE_D_MODEL_S29EA)) then
+ write(IMAIN,*) ' incorporating element stretching for 3-D internal surfaces'
+ else
+ write(IMAIN,*) ' no element stretching for 3-D internal surfaces'
+ endif
+ call flush_IMAIN()
+ endif
+
! user output
if (myrank == 0) then
write(IMAIN,*)
@@ -176,6 +221,7 @@
write(IMAIN,*) ' no general mantle anisotropy'
endif
write(IMAIN,*)
+
write(IMAIN,*) 'Reference radius of the Earth used is ',R_EARTH_KM,' km'
write(IMAIN,*)
write(IMAIN,*) 'Central cube is at a radius of ',R_CENTRAL_CUBE/1000.d0,' km'
More information about the CIG-COMMITS
mailing list