[cig-commits] r22842 - seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sun Sep 22 17:41:18 PDT 2013
Author: dkomati1
Date: 2013-09-22 17:41:18 -0700 (Sun, 22 Sep 2013)
New Revision: 22842
Modified:
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/add_topography.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/add_topography_410_650.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/add_topography_cmb.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/add_topography_icb.f90
seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/compute_element_properties.f90
Log:
fixed a few (but not all) problems in src/meshfem3D/add_topography_410_650.f90, and commented the call to it in src/meshfem3D/compute_element_properties.f90
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/add_topography.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/add_topography.f90 2013-09-22 18:43:02 UTC (rev 22841)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/add_topography.f90 2013-09-23 00:41:18 UTC (rev 22842)
@@ -129,7 +129,6 @@
phi = phi + 0.0000001d0
call reduce(theta,phi)
-
! convert the geocentric colatitude to a geographic colatitude
colat = PI_OVER_TWO - datan(1.006760466d0*dcos(theta)/dmax1(TINYVAL,dsin(theta)))
@@ -145,14 +144,12 @@
! stretching topography between d220 and the surface
gamma = (r - R220/R_EARTH) / (R_UNIT_SPHERE - R220/R_EARTH)
- !
! add elevation to all the points of that element
! also make sure factor makes sense
if(gamma < -0.02 .or. gamma > 1.02) then
call exit_MPI(myrank,'incorrect value of factor for topography gll points')
endif
- !
! since not all GLL points are exactlly at R220, use a small
! tolerance for R220 detection
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/add_topography_410_650.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/add_topography_410_650.f90 2013-09-22 18:43:02 UTC (rev 22841)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/add_topography_410_650.f90 2013-09-23 00:41:18 UTC (rev 22842)
@@ -51,13 +51,29 @@
do ia = 1,NGNOD
! convert to r theta phi
+!! DK DK old
+! call xyz_2_rthetaphi_dble(xelm(ia),yelm(ia),zelm(ia),r,theta,phi)
+! call reduce(theta,phi)
+!! DK DK new
+ ! gets elevation of point
+ ! convert to r theta phi
+ ! slightly move points to avoid roundoff problem when exactly on the polar axis
call xyz_2_rthetaphi_dble(xelm(ia),yelm(ia),zelm(ia),r,theta,phi)
+ theta = theta + 0.0000001d0
+ phi = phi + 0.0000001d0
call reduce(theta,phi)
-! get colatitude and longitude in degrees
- xcolat = sngl(theta*RADIANS_TO_DEGREES)
- xlon = sngl(phi*RADIANS_TO_DEGREES)
+!! DK DK old
+! get colatitude in degrees
+! xcolat = sngl(theta*RADIANS_TO_DEGREES) !! DK DK bug here, confusion between theta and PI_OVER_TWO - theta
+!! DK DK new
+ ! convert the geocentric colatitude to a geographic colatitude in degrees
+ xcolat = sngl((PI_OVER_TWO - datan(1.006760466d0*dcos(theta)/dmax1(TINYVAL,dsin(theta)))) * RADIANS_TO_DEGREES)
+
+ ! get geographic longitude in degrees
+ xlon = sngl(phi * RADIANS_TO_DEGREES)
+
! compute topography on 410 and 650 at current point
call model_s362ani_subtopo(xcolat,xlon,topo410out,topo650out)
@@ -101,7 +117,6 @@
! JAN08, 2010
subroutine add_topography_410_650_gll(myrank,xstore,ystore,zstore,ispec,nspec)
-
use constants
use meshfem3D_par,only: R220,R400,R670,R771
@@ -125,14 +140,25 @@
do j = 1,NGLLY
do i = 1,NGLLX
+ ! gets elevation of point
! convert to r theta phi
+ ! slightly move points to avoid roundoff problem when exactly on the polar axis
call xyz_2_rthetaphi_dble(xstore(i,j,k,ispec),ystore(i,j,k,ispec),zstore(i,j,k,ispec),r,theta,phi)
+ theta = theta + 0.0000001d0
+ phi = phi + 0.0000001d0
call reduce(theta,phi)
- ! get colatitude and longitude in degrees
- xcolat = sngl(theta*RADIANS_TO_DEGREES)
- xlon = sngl(phi*RADIANS_TO_DEGREES)
+!! DK DK old
+! get colatitude in degrees
+! xcolat = sngl(theta*RADIANS_TO_DEGREES) !! DK DK bug here, confusion between theta and PI_OVER_TWO - theta
+!! DK DK new
+ ! convert the geocentric colatitude to a geographic colatitude in degrees
+ xcolat = sngl((PI_OVER_TWO - datan(1.006760466d0*dcos(theta)/dmax1(TINYVAL,dsin(theta)))) * RADIANS_TO_DEGREES)
+
+ ! get geographic longitude in degrees
+ xlon = sngl(phi * RADIANS_TO_DEGREES)
+
! compute topography on 410 and 650 at current point
call model_s362ani_subtopo(xcolat,xlon,topo410out,topo650out)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/add_topography_cmb.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/add_topography_cmb.f90 2013-09-22 18:43:02 UTC (rev 22841)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/add_topography_cmb.f90 2013-09-23 00:41:18 UTC (rev 22842)
@@ -51,10 +51,12 @@
! convert to r theta phi
call xyz_2_rthetaphi_dble(xelm(ia),yelm(ia),zelm(ia),r,theta,phi)
+ theta = theta + 0.0000001d0
+ phi = phi + 0.0000001d0
call reduce(theta,phi)
! compute topography on CMB; routine subtopo_cmb needs to be supplied by the user
-! call subtopo_cmb(theta,phi,topocmb)
+! call subtopo_cmb(theta,phi,topocmb)
topocmb = 0.0d0
! non-dimensionalize the topography, which is in km
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/add_topography_icb.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/add_topography_icb.f90 2013-09-22 18:43:02 UTC (rev 22841)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/add_topography_icb.f90 2013-09-23 00:41:18 UTC (rev 22842)
@@ -51,10 +51,12 @@
! convert to r theta phi
call xyz_2_rthetaphi_dble(xelm(ia),yelm(ia),zelm(ia),r,theta,phi)
+ theta = theta + 0.0000001d0
+ phi = phi + 0.0000001d0
call reduce(theta,phi)
! compute topography on ICB; the routine subtopo_icb needs to be supplied by the user
-! call subtopo_icb(theta,phi,topoicb)
+! call subtopo_icb(theta,phi,topoicb)
topoicb = 0.0d0
! non-dimensionalize the topography, which is in km
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/compute_element_properties.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/compute_element_properties.f90 2013-09-22 18:43:02 UTC (rev 22841)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/SPECFEM3D_GLOBE_SUNFLOWER/src/meshfem3D/compute_element_properties.f90 2013-09-23 00:41:18 UTC (rev 22842)
@@ -216,11 +216,12 @@
.or. THREE_D_MODEL == THREE_D_MODEL_S362ANI_PREM .or. THREE_D_MODEL == THREE_D_MODEL_S29EA) then
if( USE_GLL ) then
! stretches every gll point accordingly
- call add_topography_410_650_gll(myrank,xstore,ystore,zstore,ispec,nspec)
-
+!! DK DK commented this out because it makes the mesher crash
+!! DK DK 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)
+!! DK DK commented this out because it makes the mesher crash
+!! DK DK call add_topography_410_650(myrank,xelm,yelm,zelm)
endif
endif
More information about the CIG-COMMITS
mailing list