[cig-commits] r16156 - seismo/3D/SPECFEM3D_GLOBE/trunk
danielpeter at geodynamics.org
danielpeter at geodynamics.org
Thu Jan 21 08:14:53 PST 2010
Author: danielpeter
Date: 2010-01-21 08:14:52 -0800 (Thu, 21 Jan 2010)
New Revision: 16156
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/add_topography.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/add_topography_410_650.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/attenuation_model.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/calc_jacobian.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/compute_element_properties.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_crust_mantle.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_inner_core.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in
seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/crustal_model.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/get_ellipticity.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/get_jacobian_boundaries.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/get_model.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/moho_stretching.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/read_compute_parameters.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/save_arrays_solver.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/save_header_file.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/topo_bathy.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/write_AVS_DX_global_data.f90
Log:
This version changes how attenuation is assigned on to the GLL points and
it incorporates a new element stretching to honor surface topography and Moho topography
for 3D models.
A. attenuation model
- get_model.f90:
new attenuation assignment, changed the position.
- crustal_model.f90:
added elem_in_crust to decide whether in crust or not. uses new CAP function to smooth crustal model.
- attenuation_model.f90:
added 1D attenuation assignment subroutine; attenuation_model_1D_PREM and attenuation_model_1D_REF
- create_regions_mesh.f90:
changed attenuation allocation and deallocation
- meshfem3D.f90:
new 1D attenuation allocation
- save_arrays_solver.f90:
new 1D attenuation saver.
- specfem3D.f90:
new 1D attenuation reader
- save_header_file.f90:
new 1D header file attenuation
- compute_forces_crust_mantle.f90:
new 1D attenuation
- compute_forces_inner_core.f90:
new 1D attenuation
B. stretch R80-> 120
- read_compute_parameters.f90:
R80_FICTITIOUS_IN_MESHER
C. ETOPO1
- constants.h:
ETOPO1 data, position and points
- topo_bathy.f90:
uses bilinear interpolation rather than nearest point interpolation in get_topo_bathy
D. MOHO stretch
- compute_element_properities.f90:
added subroutine hornor_moho_crust2
- moho_stretching.f90:
added subroutine hornor_moho_crust2
E. GLL point
- compute_element_properities.f90:
uses USE_GLL to decide whether to use GLL points or anchor points only for stretching
- constants.h:
added USE_GLL parameter
- calc_jacobian.f90:
added two subroutines, recalc_jacobian_gll3D and recalc_jacobian_gll2D to recalculate
jacobian based upon GLL points
- get_jacobian_boundaries.f90:
added recalc_jacobian_gll2D to calculate 2D Jacobian based upon 25 GLL points.
- add_topography.f90:
added subroutine add_topography_gll to stretch GLL points
- add_topography_410_650.f90:
added subroutine add_topography_410_650_gll to stretch GLL points
- get_ellipticity.f90:
added subroutine get_ellipticity_gll to stretch GLL points
- write_AVS_DX_global_data.f90
added subroutine write_AVS_DX_global_data_gll to output material information for GLL points
- create_region_mesh.f90
added " call write_AVS_DX_global_data_gll " to output material information for GLL points,
can be commented out for 6 chunks simulation.
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/add_topography.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/add_topography.f90 2010-01-21 16:01:30 UTC (rev 16155)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/add_topography.f90 2010-01-21 16:14:52 UTC (rev 16156)
@@ -85,3 +85,85 @@
end subroutine add_topography
+ !> Hejun
+ ! This subroutine uses GLL points to capture topography variation rather
+ ! than using control nodes
+ ! Hejun Zhu, OCT16, 2009
+
+ ! input parameters: myrank,
+ ! xstore,ystore,zstore,
+ ! ispec,nspec,
+ ! ibathy_topo
+ ! R220
+
+ subroutine add_topography_gll(myrank,xstore,ystore,zstore,ispec,nspec,&
+ ibathy_topo,R220)
+
+ implicit none
+
+ include "constants.h"
+
+ ! input parameters
+ integer:: myrank
+ integer:: ispec,nspec
+ double precision,dimension(NGLLX,NGLLY,NGLLZ,nspec):: xstore,ystore,zstore
+ integer, dimension(NX_BATHY,NY_BATHY) :: ibathy_topo
+ double precision:: R220
+
+ ! local parameters used in this subroutine
+ integer:: i,j,k
+ double precision:: r,theta,phi,colat
+ double precision:: lat,lon,elevation,gamma
+
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
+ ! 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)
+
+
+ ! convert the geocentric colatitude to a geographic colatitude
+ colat = PI/2.0d0 - datan(1.006760466d0*dcos(theta)/dmax1(TINYVAL,dsin(theta)))
+
+ ! get geographic latitude and longitude in degrees
+ lat = 90.0d0 - colat*180.0d0/PI
+ lon = phi*180.0d0/PI
+ elevation = 0.d0
+
+ ! compute elevation at current point
+ call get_topo_bathy(lat,lon,elevation,ibathy_topo)
+ ! non-dimensionalize the elevation, which is in meters
+
+ elevation = elevation / R_EARTH
+
+ ! 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')
+ end if
+ !
+
+ ! since not all GLL points are exactlly at R220, use a small
+ ! tolerance for R220 detection
+ if (abs(gamma) < SMALLVAL) then
+ gamma = 0.0
+ end if
+ xstore(i,j,k,ispec) = xstore(i,j,k,ispec)*(ONE + gamma * elevation / r)
+ ystore(i,j,k,ispec) = ystore(i,j,k,ispec)*(ONE + gamma * elevation / r)
+ zstore(i,j,k,ispec) = zstore(i,j,k,ispec)*(ONE + gamma * elevation / r)
+
+ end do
+ end do
+ end do
+ end subroutine add_topography_gll
+ !< Hejun
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/add_topography_410_650.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/add_topography_410_650.f90 2010-01-21 16:01:30 UTC (rev 16155)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/add_topography_410_650.f90 2010-01-21 16:14:52 UTC (rev 16156)
@@ -132,3 +132,115 @@
end subroutine add_topography_410_650
+ !> Hejun
+ ! use GLL points to capture 410_650 topography
+ ! JAN08, 2010
+ subroutine add_topography_410_650_gll(myrank,xstore,ystore,zstore,ispec,nspec,R220,R400,R670,R771, &
+ numker,numhpa,numcof,ihpa,lmax,nylm, &
+ lmxhpa,itypehpa,ihpakern,numcoe,ivarkern, &
+ nconpt,iver,iconpt,conpt,xlaspl,xlospl,radspl, &
+ coe,ylmcof,wk1,wk2,wk3,varstr)
+
+ implicit none
+
+ include "constants.h"
+
+ integer myrank
+ integer:: ispec,nspec
+ double precision,dimension(NGLLX,NGLLY,NGLLZ,nspec):: xstore,ystore,zstore
+
+ double precision R220,R400,R670,R771
+
+ integer i,j,k
+
+ real(kind=4) xcolat,xlon
+ real(kind=4) topo410out,topo650out
+ double precision topo410,topo650
+
+ double precision r,theta,phi
+ double precision gamma
+
+ integer, parameter :: maxker=200
+ integer, parameter :: maxl=72
+ integer, parameter :: maxcoe=2000
+ 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 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)
+
+ real(kind=4) ylmcof((maxl+1)**2,maxhpa)
+ real(kind=4) wk1(maxl+1)
+ real(kind=4) wk2(maxl+1)
+ real(kind=4) wk3(maxl+1)
+
+ character(len=40) varstr(maxker)
+
+ ! we loop on all GLL points of the element
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
+ ! convert to r theta phi
+ call xyz_2_rthetaphi_dble(xstore(i,j,k,ispec),ystore(i,j,k,ispec),zstore(i,j,k,ispec),r,theta,phi)
+ call reduce(theta,phi)
+
+ ! get colatitude and longitude in degrees
+ xcolat = sngl(theta*180.0d0/PI)
+ xlon = sngl(phi*180.0d0/PI)
+
+ ! compute topography on 410 and 650 at current point
+ call subtopo(xcolat,xlon,topo410out,topo650out, &
+ numker,numhpa,numcof,ihpa,lmax,nylm, &
+ lmxhpa,itypehpa,ihpakern,numcoe,ivarkern, &
+ nconpt,iver,iconpt,conpt,xlaspl,xlospl,radspl, &
+ coe,ylmcof,wk1,wk2,wk3,varstr)
+
+ ! 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) = xstore(i,j,k,ispec)*(ONE + gamma * topo410 / r)
+ ystore(i,j,k,ispec) = ystore(i,j,k,ispec)*(ONE + gamma * topo410 / r)
+ zstore(i,j,k,ispec) = zstore(i,j,k,ispec)*(ONE + gamma * topo410 / r)
+ elseif(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) = xstore(i,j,k,ispec)*(ONE + gamma * topo650 / r)
+ ystore(i,j,k,ispec) = ystore(i,j,k,ispec)*(ONE + gamma * topo650 / r)
+ zstore(i,j,k,ispec) = zstore(i,j,k,ispec)*(ONE + gamma * topo650 / r)
+ elseif(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) = xstore(i,j,k,ispec)*(ONE + (topo410 + gamma * (topo650 - topo410)) / r)
+ ystore(i,j,k,ispec) = ystore(i,j,k,ispec)*(ONE + (topo410 + gamma * (topo650 - topo410)) / r)
+ zstore(i,j,k,ispec) = zstore(i,j,k,ispec)*(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
+ end do
+ end do
+
+ end subroutine add_topography_410_650_gll
+ !< Hejun
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/attenuation_model.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/attenuation_model.f90 2010-01-21 16:01:30 UTC (rev 16155)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/attenuation_model.f90 2010-01-21 16:14:52 UTC (rev 16156)
@@ -955,8 +955,10 @@
! All of the following reads use the output parameters as their temporary arrays
! use the filename to determine the actual contents of the read
-
- open(unit=27, file=prname(1:len_trim(prname))//'attenuation3D.bin',status='old',action='read',form='unformatted')
+!> Hejun
+! open(unit=27, file=prname(1:len_trim(prname))//'attenuation3D.bin',status='old',action='read',form='unformatted')
+ open(unit=27, file=prname(1:len_trim(prname))//'attenuation.bin',status='old',action='read',form='unformatted')
+!< Hejun
read(27) tau_s
read(27) factor_common
read(27) scale_factor
@@ -1787,7 +1789,7 @@
end subroutine pspline_construction
-subroutine attenuation_model_1D_PREM(x, Qmu, iflag)
+subroutine attenuation_model_1D_PREM(x, Qmu)
! x in the radius from 0 to 1 where 0 is the center and 1 is the surface
! This version is for 1D PREM.
@@ -1795,8 +1797,9 @@
implicit none
include 'constants.h'
-
- integer iflag
+!> Hejun
+! integer iflag
+!< Hejun
double precision r, x, Qmu,RICB,RCMB, &
RTOPDDOUBLEPRIME,R600,R670,R220,R771,R400,R80, ROCEAN, RMOHO, RMIDDLE_CRUST
double precision Qkappa
@@ -1869,36 +1872,122 @@
Qkappa=57827.0d0
endif
+!> Hejun
+! Since R80 may be changed, we use radius to decide the attenuation region
+! rather than doubling flag
+
! We determine the attenuation value here dependent on the doubling flag and
! which region we are sitting in. The radius reported is not accurate for
! determination of which region we are actually in, whereas the idoubling flag is
- if(iflag == IFLAG_INNER_CORE_NORMAL .or. iflag == IFLAG_MIDDLE_CENTRAL_CUBE .or. &
- iflag == IFLAG_BOTTOM_CENTRAL_CUBE .or. iflag == IFLAG_TOP_CENTRAL_CUBE .or. &
- iflag == IFLAG_IN_FICTITIOUS_CUBE) then
- Qmu = 84.6d0
- Qkappa = 1327.7d0
- else if(iflag == IFLAG_OUTER_CORE_NORMAL) then
- Qmu = 0.0d0
- Qkappa = 57827.0d0
- else if(iflag == IFLAG_MANTLE_NORMAL) then ! D'' to 670 km
- Qmu = 312.0d0
- Qkappa = 57827.0d0
- else if(iflag == IFLAG_670_220) then
- Qmu=143.0d0
- Qkappa = 57827.0d0
- else if(iflag == IFLAG_220_80) then
- Qmu=80.0d0
- Qkappa = 57827.0d0
- else if(iflag == IFLAG_80_MOHO) then
- Qmu=600.0d0
- Qkappa = 57827.0d0
- else if(iflag == IFLAG_CRUST) then
- Qmu=600.0d0
- Qkappa = 57827.0d0
- else
- write(*,*)'iflag:',iflag
- call exit_MPI_without_rank('Invalid idoubling flag in attenuation_model_1D_prem from get_model()')
- endif
+! if(iflag == IFLAG_INNER_CORE_NORMAL .or. iflag == IFLAG_MIDDLE_CENTRAL_CUBE .or. &
+! iflag == IFLAG_BOTTOM_CENTRAL_CUBE .or. iflag == IFLAG_TOP_CENTRAL_CUBE .or. &
+! iflag == IFLAG_IN_FICTITIOUS_CUBE) then
+! Qmu = 84.6d0
+! Qkappa = 1327.7d0
+! else if(iflag == IFLAG_OUTER_CORE_NORMAL) then
+! Qmu = 0.0d0
+! Qkappa = 57827.0d0
+! else if(iflag == IFLAG_MANTLE_NORMAL) then ! D'' to 670 km
+! Qmu = 312.0d0
+! Qkappa = 57827.0d0
+! else if(iflag == IFLAG_670_220) then
+! Qmu=143.0d0
+! Qkappa = 57827.0d0
+! else if(iflag == IFLAG_220_80) then
+! Qmu=80.0d0
+! Qkappa = 57827.0d0
+! else if(iflag == IFLAG_80_MOHO) then
+! Qmu=600.0d0
+! Qkappa = 57827.0d0
+! else if(iflag == IFLAG_CRUST) then
+! Qmu=600.0d0
+! Qkappa = 57827.0d0
+! else
+! write(*,*)'iflag:',iflag
+! call exit_MPI_without_rank('Invalid idoubling flag in attenuation_model_1D_prem from get_model()')
+! endif
+!< Hejun
end subroutine attenuation_model_1D_PREM
+
+!> Hejun
+! get 1D REF attenuation model according to radius
+subroutine attenuation_model_1D_REF(x, Qmu)
+
+! x in the radius from 0 to 1 where 0 is the center and 1 is the surface
+! This version is for 1D REF.
+
+ implicit none
+
+ include 'constants.h'
+
+ double precision r, x, Qmu,RICB,RCMB, &
+ RTOPDDOUBLEPRIME,R600,R670,R220,R771,R400,R80, ROCEAN, RMOHO, RMIDDLE_CRUST
+ double precision Qkappa
+
+ r = x * R_EARTH
+
+ ROCEAN = 6368000.d0
+ RMIDDLE_CRUST = 6356000.d0
+ RMOHO = 6346600.d0
+ R80 = 6291000.d0
+ R220 = 6151000.d0
+ R400 = 5961000.d0
+ R600 = 5771000.d0
+ R670 = 5721000.d0
+ R771 = 5600000.d0
+ RTOPDDOUBLEPRIME = 3630000.d0
+ RCMB = 3479958.d0
+ RICB = 1221491.d0
+
+! REF model
+!
+!--- inner core
+!
+ if(r >= 0.d0 .and. r <= RICB) then
+ Qmu=104.0d0
+ Qkappa=1327.6d0
+
+!--- outer core
+!
+ else if(r > RICB .and. r <= RCMB) then
+ Qmu=0.0d0
+ Qkappa=57822.5d0
+ if(RCMB - r < r - RICB) then
+ Qmu = 355.0d0 ! CMB
+ else
+ Qmu = 104.0d0 ! ICB
+ endif
+
+!--- D" at the base of the mantle
+!
+ else if(r > RCMB .and. r <= RTOPDDOUBLEPRIME) then
+ Qmu=355.0d0
+ Qkappa=57822.5d0
+
+!--- mantle: from top of D" to d670
+!
+ else if(r > RTOPDDOUBLEPRIME .and. r <= R670) then
+ Qmu=355.0d0
+ Qkappa=57822.5d0
+
+!--- mantle: above d670
+!
+ else if(r > R670 .and. r <= R220) then
+ Qmu=165.0d0
+ Qkappa=943.0d0
+ else if(r > R220 .and. r <= R80) then
+ Qmu=70.0d0
+ Qkappa=943.0d0
+ else if(r > R80.and. r<=RMOHO) then
+ Qmu=191.0d0
+ Qkappa=943.0d0
+ else if (r > RMOHO) then
+ Qmu=300.0d0
+ Qkappa=57822.5d0
+ endif
+
+end subroutine attenuation_model_1D_REF
+!< Hejun
+
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/calc_jacobian.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/calc_jacobian.f90 2010-01-21 16:01:30 UTC (rev 16155)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/calc_jacobian.f90 2010-01-21 16:14:52 UTC (rev 16156)
@@ -153,3 +153,322 @@
end subroutine calc_jacobian
+
+!> Hejun
+! This subroutine recompute the 3D jacobian for one element
+! based upon 125 GLL points
+! Hejun Zhu OCT16,2009
+
+! input: myrank,
+! xstore,ystore,zstore ----- input GLL point coordinate
+! xigll,yigll,zigll ----- gll points position
+! ispec,nspec ----- element number
+! ACTUALLY_STORE_ARRAYS ------ save array or not
+
+! output: xixstore,xiystore,xizstore,
+! etaxstore,etaystore,etazstore,
+! gammaxstore,gammaystore,gammazstore ------ parameters used to calculate jacobian
+
+
+ subroutine recalc_jacobian_gll3D(myrank,xstore,ystore,zstore,xigll,yigll,zigll,&
+ ispec,nspec,ACTUALLY_STORE_ARRAYS,&
+ xixstore,xiystore,xizstore, &
+ etaxstore,etaystore,etazstore, &
+ gammaxstore,gammaystore,gammazstore)
+
+ implicit none
+
+ include "constants.h"
+
+ ! input parameter
+ integer::myrank,ispec,nspec
+ double precision, dimension(NGLLX,NGLLY,NGLLZ,nspec) :: xstore,ystore,zstore
+ double precision, dimension(NGLLX):: xigll
+ double precision, dimension(NGLLY):: yigll
+ double precision, dimension(NGLLZ):: zigll
+ logical::ACTUALLY_STORE_ARRAYS
+
+
+ ! output results
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,nspec) :: &
+ xixstore,xiystore,xizstore,&
+ etaxstore,etaystore,etazstore,&
+ gammaxstore,gammaystore,gammazstore
+
+
+ ! local parameters for this subroutine
+ integer:: i,j,k,i1,j1,k1
+ double precision:: xxi,xeta,xgamma,yxi,yeta,ygamma,zxi,zeta,zgamma
+ double precision:: xi,eta,gamma
+ double precision,dimension(NGLLX):: hxir,hpxir
+ double precision,dimension(NGLLY):: hetar,hpetar
+ double precision,dimension(NGLLZ):: hgammar,hpgammar
+ double precision:: hlagrange,hlagrange_xi,hlagrange_eta,hlagrange_gamma
+ double precision:: jacobian
+ double precision:: xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz
+
+
+
+ ! test parameters which can be deleted
+ double precision:: xmesh,ymesh,zmesh
+ double precision:: sumshape,sumdershapexi,sumdershapeeta,sumdershapegamma
+
+ ! first go over all 125 GLL points
+ do k=1,NGLLZ
+ do j=1,NGLLY
+ do i=1,NGLLX
+
+ xxi = 0.0
+ xeta = 0.0
+ xgamma = 0.0
+ yxi = 0.0
+ yeta = 0.0
+ ygamma = 0.0
+ zxi = 0.0
+ zeta = 0.0
+ zgamma = 0.0
+
+ xi = xigll(i)
+ eta = yigll(j)
+ gamma = zigll(k)
+
+ ! calculate lagrange polynomial and its derivative
+ call lagrange_any(xi,NGLLX,xigll,hxir,hpxir)
+ call lagrange_any(eta,NGLLY,yigll,hetar,hpetar)
+ call lagrange_any(gamma,NGLLZ,zigll,hgammar,hpgammar)
+
+ ! test parameters
+ sumshape = 0.0
+ sumdershapexi = 0.0
+ sumdershapeeta = 0.0
+ sumdershapegamma = 0.0
+ xmesh = 0.0
+ ymesh = 0.0
+ zmesh = 0.0
+
+
+ do k1 = 1,NGLLZ
+ do j1 = 1,NGLLY
+ do i1 = 1,NGLLX
+ hlagrange = hxir(i1)*hetar(j1)*hgammar(k1)
+ hlagrange_xi = hpxir(i1)*hetar(j1)*hgammar(k1)
+ hlagrange_eta = hxir(i1)*hpetar(j1)*hgammar(k1)
+ hlagrange_gamma = hxir(i1)*hetar(j1)*hpgammar(k1)
+
+
+ xxi = xxi + xstore(i1,j1,k1,ispec)*hlagrange_xi
+ xeta = xeta + xstore(i1,j1,k1,ispec)*hlagrange_eta
+ xgamma = xgamma + xstore(i1,j1,k1,ispec)*hlagrange_gamma
+
+ yxi = yxi + ystore(i1,j1,k1,ispec)*hlagrange_xi
+ yeta = yeta + ystore(i1,j1,k1,ispec)*hlagrange_eta
+ ygamma = ygamma + ystore(i1,j1,k1,ispec)*hlagrange_gamma
+
+ zxi = zxi + zstore(i1,j1,k1,ispec)*hlagrange_xi
+ zeta = zeta + zstore(i1,j1,k1,ispec)*hlagrange_eta
+ zgamma = zgamma + zstore(i1,j1,k1,ispec)*hlagrange_gamma
+
+ ! test the lagrange polynomial and its derivate
+ xmesh = xmesh + xstore(i1,j1,k1,ispec)*hlagrange
+ ymesh = ymesh + ystore(i1,j1,k1,ispec)*hlagrange
+ zmesh = zmesh + zstore(i1,j1,k1,ispec)*hlagrange
+ sumshape = sumshape + hlagrange
+ sumdershapexi = sumdershapexi + hlagrange_xi
+ sumdershapeeta = sumdershapeeta + hlagrange_eta
+ sumdershapegamma = sumdershapegamma + hlagrange_gamma
+
+ end do
+ end do
+ end do
+
+ ! Check the lagrange polynomial and its derivative
+ if (xmesh /=xstore(i,j,k,ispec).or.ymesh/=ystore(i,j,k,ispec)&
+ .or.zmesh/=zstore(i,j,k,ispec)) then
+ call exit_MPI(myrank,'new mesh are wrong in recalc_jacobian_gall3D.f90')
+ end if
+ if(abs(sumshape-one) > TINYVAL) then
+ call exit_MPI(myrank,'error shape functions in recalc_jacobian_gll3D.f90')
+ end if
+ if(abs(sumdershapexi) > TINYVAL) then
+ call exit_MPI(myrank,'error derivative xi in recalc_jacobian_gll3D.f90')
+ end if
+ if(abs(sumdershapeeta) > TINYVAL) then
+ call exit_MPI(myrank,'error derivative eta in recalc_jacobian_gll3D.f90')
+ end if
+ if(abs(sumdershapegamma) > TINYVAL) then
+ call exit_MPI(myrank,'error derivative gamma in recalc_jacobian_gll3D.f90')
+ end if
+
+
+ jacobian = xxi*(yeta*zgamma-ygamma*zeta) - &
+ xeta*(yxi*zgamma-ygamma*zxi) + &
+ xgamma*(yxi*zeta-yeta*zxi)
+
+ ! Check the jacobian
+ if(jacobian <= ZERO) then
+ call exit_MPI(myrank,'3D Jacobian undefined in recalc_jacobian_gll3D.f90')
+ end if
+
+ ! invert the relation (Fletcher p. 50 vol. 2)
+ xix = (yeta*zgamma-ygamma*zeta) / jacobian
+ xiy = (xgamma*zeta-xeta*zgamma) / jacobian
+ xiz = (xeta*ygamma-xgamma*yeta) / jacobian
+ etax = (ygamma*zxi-yxi*zgamma) / jacobian
+ etay = (xxi*zgamma-xgamma*zxi) / jacobian
+ etaz = (xgamma*yxi-xxi*ygamma) / jacobian
+ gammax = (yxi*zeta-yeta*zxi) / jacobian
+ gammay = (xeta*zxi-xxi*zeta) / jacobian
+ gammaz = (xxi*yeta-xeta*yxi) / jacobian
+
+
+ ! resave the derivatives and the jacobian
+ ! distinguish between single and double precision for reals
+ if (ACTUALLY_STORE_ARRAYS) then
+ if(CUSTOM_REAL == SIZE_REAL) then
+ xixstore(i,j,k,ispec) = sngl(xix)
+ xiystore(i,j,k,ispec) = sngl(xiy)
+ xizstore(i,j,k,ispec) = sngl(xiz)
+ etaxstore(i,j,k,ispec) = sngl(etax)
+ etaystore(i,j,k,ispec) = sngl(etay)
+ etazstore(i,j,k,ispec) = sngl(etaz)
+ gammaxstore(i,j,k,ispec) = sngl(gammax)
+ gammaystore(i,j,k,ispec) = sngl(gammay)
+ gammazstore(i,j,k,ispec) = sngl(gammaz)
+ else
+ xixstore(i,j,k,ispec) = xix
+ xiystore(i,j,k,ispec) = xiy
+ xizstore(i,j,k,ispec) = xiz
+ etaxstore(i,j,k,ispec) = etax
+ etaystore(i,j,k,ispec) = etay
+ etazstore(i,j,k,ispec) = etaz
+ gammaxstore(i,j,k,ispec) = gammax
+ gammaystore(i,j,k,ispec) = gammay
+ gammazstore(i,j,k,ispec) = gammaz
+ endif
+ end if
+ enddo
+ enddo
+ enddo
+
+ end subroutine recalc_jacobian_gll3D
+
+
+
+ ! Hejun Zhu used to recalculate 2D jacobian according to gll points rather
+ ! than control nodes
+ ! Hejun Zhu JAN08, 2010
+
+ ! input parameters: myrank,ispecb,
+ ! xelm2D,yelm2D,zelm2D,
+ ! xigll,yigll,NSPEC2DMAX_AB,NGLLA,NGLLB
+
+ ! output results: jacobian2D,normal
+ subroutine recalc_jacobian_gll2D(myrank,ispecb, &
+ xelm2D,yelm2D,zelm2D,xigll,yigll,&
+ jacobian2D,normal,NGLLA,NGLLB,NSPEC2DMAX_AB)
+
+ implicit none
+ include "constants.h"
+ ! input parameters
+ integer::myrank,ispecb,NSPEC2DMAX_AB,NGLLA,NGLLB
+ double precision,dimension(NGLLA,NGLLB)::xelm2D,yelm2D,zelm2D
+ double precision,dimension(NGLLA)::xigll
+ double precision,dimension(NGLLB)::yigll
+
+ ! output results
+ real(kind=CUSTOM_REAL),dimension(NGLLA,NGLLB,NSPEC2DMAX_AB)::jacobian2D
+ real(kind=CUSTOM_REAL),dimension(3,NGLLA,NGLLB,NSPEC2DMAX_AB)::normal
+
+
+ ! local parameters in this subroutine
+ integer::i,j,i1,j1
+ double precision::xxi,xeta,yxi,yeta,zxi,zeta,&
+ xi,eta,xmesh,ymesh,zmesh,hlagrange,hlagrange_xi,hlagrange_eta,&
+ sumshape,sumdershapexi,sumdershapeeta,unx,uny,unz,jacobian
+ double precision,dimension(NGLLA)::hxir,hpxir
+ double precision,dimension(NGLLB)::hetar,hpetar
+
+ do j = 1,NGLLB
+ do i = 1,NGLLA
+ xxi = 0.0
+ xeta = 0.0
+ yxi = 0.0
+ yeta = 0.0
+ zxi = 0.0
+ zeta = 0.0
+
+ xi=xigll(i)
+ eta = yigll(j)
+
+ call lagrange_any(xi,NGLLA,xigll,hxir,hpxir)
+ call lagrange_any(eta,NGLLB,yigll,hetar,hpetar)
+
+
+ xmesh = 0.0
+ ymesh = 0.0
+ zmesh = 0.0
+ sumshape = 0.0
+ sumdershapexi = 0.0
+ sumdershapeeta = 0.0
+ do j1 = 1,NGLLB
+ do i1 = 1,NGLLA
+ hlagrange = hxir(i1)*hetar(j1)
+ hlagrange_xi = hpxir(i1)*hetar(j1)
+ hlagrange_eta = hxir(i1)*hpetar(j1)
+
+ xxi = xxi + xelm2D(i1,j1)*hlagrange_xi
+ xeta = xeta + xelm2D(i1,j1)*hlagrange_eta
+
+ yxi = yxi + yelm2D(i1,j1)*hlagrange_xi
+ yeta = yeta + yelm2D(i1,j1)*hlagrange_eta
+
+ zxi = zxi + zelm2D(i1,j1)*hlagrange_xi
+ zeta = zeta + zelm2D(i1,j1)*hlagrange_eta
+
+ xmesh = xmesh + xelm2D(i1,j1)*hlagrange
+ ymesh = ymesh + yelm2D(i1,j1)*hlagrange
+ zmesh = zmesh + zelm2D(i1,j1)*hlagrange
+ sumshape = sumshape + hlagrange
+ sumdershapexi = sumdershapexi + hlagrange_xi
+ sumdershapeeta = sumdershapeeta + hlagrange_eta
+ end do
+ end do
+
+
+ ! Check the lagrange polynomial
+ if (xmesh/=xelm2D(i,j).or.ymesh/=yelm2D(i,j).or.zmesh/=zelm2D(i,j)) then
+ call exit_MPI(myrank,'new boundary mesh is wrong in recalc_jacobian_gll2D')
+ end if
+
+ if (abs(sumshape-one) > TINYVAL) then
+ call exit_MPI(myrank,'error shape functions in recalc_jacobian_gll2D')
+ end if
+ if (abs(sumdershapexi) > TINYVAL) then
+ call exit_MPI(myrank,'error derivative xi in recalc_jacobian_gll2D')
+ end if
+ if (abs(sumdershapeeta) > TINYVAL) then
+ call exit_MPI(myrank,'error derivative eta in recalc_jacobian_gll2D')
+ end if
+
+ unx = yxi*zeta - yeta*zxi
+ uny = zxi*xeta - zeta*xxi
+ unz = xxi*yeta - xeta*yxi
+ jacobian = dsqrt(unx**2+uny**2+unz**2)
+ if (jacobian == ZERO) call exit_MPI(myrank,'2D Jacobian undefined in recalc_jacobian_gll2D')
+
+ if (CUSTOM_REAL == SIZE_REAL) then
+ jacobian2D(i,j,ispecb)=sngl(jacobian)
+ normal(1,i,j,ispecb)=sngl(unx/jacobian)
+ normal(2,i,j,ispecb)=sngl(uny/jacobian)
+ normal(3,i,j,ispecb)=sngl(unz/jacobian)
+ else
+ jacobian2D(i,j,ispecb)=jacobian
+ normal(1,i,j,ispecb)=unx/jacobian
+ normal(2,i,j,ispecb)=uny/jacobian
+ normal(3,i,j,ispecb)=unz/jacobian
+ endif
+ end do
+ end do
+
+ end subroutine recalc_jacobian_gll2D
+!< Hejun
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/compute_element_properties.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/compute_element_properties.f90 2010-01-21 16:01:30 UTC (rev 16155)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/compute_element_properties.f90 2010-01-21 16:14:52 UTC (rev 16156)
@@ -33,7 +33,7 @@
HETEROGEN_3D_MANTLE,CRUSTAL,ONE_CRUST, &
myrank,ibathy_topo,ATTENUATION,ATTENUATION_3D, &
ABSORBING_CONDITIONS,REFERENCE_1D_MODEL,THREE_D_MODEL, &
- RICB,RCMB,R670,RMOHO,RTOPDDOUBLEPRIME,R600,R220,R771,R400,R120,R80,RMIDDLE_CRUST,ROCEAN, &
+ RICB,RCMB,R670,RMOHO,RMOHO_FICTITIOUS_IN_MESHER,RTOPDDOUBLEPRIME,R600,R220,R771,R400,R120,R80,RMIDDLE_CRUST,ROCEAN, &
xelm,yelm,zelm,shape3D,dershape3D,rmin,rmax,rhostore,dvpstore, &
kappavstore,kappahstore,muvstore,muhstore,eta_anisostore, &
xixstore,xiystore,xizstore,etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore, &
@@ -45,7 +45,8 @@
numker,numhpa,numcof,ihpa,lmax,nylm, &
lmxhpa,itypehpa,ihpakern,numcoe,ivarkern, &
nconpt,iver,iconpt,conpt,xlaspl,xlospl,radspl, &
- coe,vercof,vercofd,ylmcof,wk1,wk2,wk3,kerstr,varstr,ACTUALLY_STORE_ARRAYS)
+ coe,vercof,vercofd,ylmcof,wk1,wk2,wk3,kerstr,varstr,ACTUALLY_STORE_ARRAYS,&
+ xigll,yigll,zigll)
implicit none
@@ -321,6 +322,9 @@
double precision RICB,RCMB,R670,RMOHO, &
RTOPDDOUBLEPRIME,R600,R220,R771,R400,R120,R80,RMIDDLE_CRUST,ROCEAN
+!> Hejun
+ integer:: RMOHO_FICTITIOUS_IN_MESHER
+!< Hejun
! use integer array to store values
integer, dimension(NX_BATHY,NY_BATHY) :: ibathy_topo
@@ -409,12 +413,33 @@
character(len=80) kerstr
character(len=40) varstr(maxker)
+!> Hejun
+! Parameters used to calculate Jacobian based upon 125 GLL points
+ double precision:: xigll(NGLLX)
+ double precision:: yigll(NGLLY)
+ double precision:: zigll(NGLLZ)
+! Parameter used to decide whether this element in crust or not
+ logical:: elem_in_crust
+ elem_in_crust=.false.
+!< Hejun
+
! **************
! add topography on the Moho *before* adding the 3D crustal model so that the streched
! mesh gets assigned the right model values
+! if(THREE_D_MODEL/=0 .and. (idoubling(ispec)==IFLAG_CRUST .or. idoubling(ispec)==IFLAG_220_80 &
+! .or. idoubling(ispec)==IFLAG_80_MOHO)) call moho_stretching(myrank,xelm,yelm,zelm,RMOHO,R220)
+
+
+!> Hejun
+! Stretch mesh to hornor smoothed moho thickness from crust2.0
if(THREE_D_MODEL/=0 .and. (idoubling(ispec)==IFLAG_CRUST .or. idoubling(ispec)==IFLAG_220_80 &
- .or. idoubling(ispec)==IFLAG_80_MOHO)) call moho_stretching(myrank,xelm,yelm,zelm,RMOHO,R220)
+ .or. idoubling(ispec)==IFLAG_80_MOHO)) then
+ call hornor_moho_crust2(myrank,xelm,yelm,zelm,RMOHO_FICTITIOUS_IN_MESHER,&
+ R220,RMIDDLE_CRUST,CM_V,elem_in_crust)
+ end if
+!< Hejun
+
! compute values for the Earth model
call get_model(myrank,iregion_code,nspec, &
kappavstore,kappahstore,muvstore,muhstore,eta_anisostore,rhostore,dvpstore,nspec_ani, &
@@ -433,39 +458,91 @@
numker,numhpa,numcof,ihpa,lmax,nylm, &
lmxhpa,itypehpa,ihpakern,numcoe,ivarkern, &
nconpt,iver,iconpt,conpt,xlaspl,xlospl,radspl, &
- coe,vercof,vercofd,ylmcof,wk1,wk2,wk3,kerstr,varstr)
+ coe,vercof,vercofd,ylmcof,wk1,wk2,wk3,kerstr,varstr,elem_in_crust)
-! add topography without the crustal model
- if(TOPOGRAPHY .and. (idoubling(ispec)==IFLAG_CRUST .or. idoubling(ispec)==IFLAG_220_80 &
- .or. idoubling(ispec)==IFLAG_80_MOHO)) call add_topography(myrank,xelm,yelm,zelm,ibathy_topo,R220)
+!> Hejun
+! Use GLL points or anchor points to capture TOPOGRAPHY and ELLIPCITY
+ if (.not.USE_GLL) then
+ ! add topography without the crustal model
+ if(TOPOGRAPHY .and. (idoubling(ispec)==IFLAG_CRUST .or. idoubling(ispec)==IFLAG_220_80 &
+ .or. idoubling(ispec)==IFLAG_80_MOHO)) call add_topography(myrank,xelm,yelm,zelm,ibathy_topo,R220)
-! add 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) &
- call add_topography_410_650(myrank,xelm,yelm,zelm,R220,R400,R670,R771, &
+ ! add 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) &
+ call add_topography_410_650(myrank,xelm,yelm,zelm,R220,R400,R670,R771, &
numker,numhpa,numcof,ihpa,lmax,nylm, &
lmxhpa,itypehpa,ihpakern,numcoe,ivarkern, &
nconpt,iver,iconpt,conpt,xlaspl,xlospl,radspl, &
coe,ylmcof,wk1,wk2,wk3,varstr)
-! 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)
+ ! 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)
+ ! 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) call get_ellipticity(xelm,yelm,zelm,nspl,rspl,espl,espl2)
+ ! make the Earth elliptical
+ if(ELLIPTICITY) call get_ellipticity(xelm,yelm,zelm,nspl,rspl,espl,espl2)
-! recompute coordinates and jacobian for real 3-D model
- call calc_jacobian(myrank,xixstore,xiystore,xizstore, &
- etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore, &
- xstore,ystore,zstore,xelm,yelm,zelm,shape3D,dershape3D,ispec,nspec,ACTUALLY_STORE_ARRAYS)
+ ! recompute coordinates and jacobian for real 3-D model
+ call calc_jacobian(myrank,xixstore,xiystore,xizstore, &
+ etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore, &
+ xstore,ystore,zstore,xelm,yelm,zelm,shape3D,dershape3D,ispec,nspec,ACTUALLY_STORE_ARRAYS)
+ else
+ ! first calculate coordinate for GLL points
+ call calc_jacobian(myrank,xixstore,xiystore,xizstore, &
+ etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore, &
+ xstore,ystore,zstore,xelm,yelm,zelm,shape3D,dershape3D,ispec,nspec,ACTUALLY_STORE_ARRAYS)
+
+ ! add topography without the crustal model, by using ALL 125 GLL points
+ if(TOPOGRAPHY .and. (idoubling(ispec)==IFLAG_CRUST .or. idoubling(ispec)==IFLAG_220_80 &
+ .or. idoubling(ispec)==IFLAG_80_MOHO))&
+ call add_topography_gll(myrank,xstore,ystore,zstore,ispec,nspec,ibathy_topo,R220)
+
+
+ ! add topography on 410 km and 650 km discontinuity in model S362ANI,
+ ! use 125 GLL points
+ 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) &
+ call add_topography_410_650_gll(myrank,xstore,ystore,zstore,ispec,nspec,R220,R400,R670,R771, &
+ numker,numhpa,numcof,ihpa,lmax,nylm, &
+ lmxhpa,itypehpa,ihpakern,numcoe,ivarkern, &
+ nconpt,iver,iconpt,conpt,xlaspl,xlospl,radspl, &
+ coe,ylmcof,wk1,wk2,wk3,varstr)
+
+
+ ! 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 ellipticity, use GLL points
+ if(ELLIPTICITY) call get_ellipticity_gll(xstore,ystore,zstore,ispec,nspec,nspl,rspl,espl,espl2)
+
+
+ ! recalculate jacobian according new coordinates
+ call recalc_jacobian_gll3D(myrank,xstore,ystore,zstore,xigll,yigll,zigll,&
+ ispec,nspec,ACTUALLY_STORE_ARRAYS,&
+ xixstore,xiystore,xizstore,&
+ etaxstore,etaystore,etazstore,&
+ gammaxstore,gammaystore,gammazstore)
+
+
+ end if
+ !< Hejun
end subroutine compute_element_properties
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_crust_mantle.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_crust_mantle.f90 2010-01-21 16:01:30 UTC (rev 16155)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_crust_mantle.f90 2010-01-21 16:14:52 UTC (rev 16156)
@@ -36,7 +36,7 @@
c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
ibool,idoubling,R_memory,epsilondev,epsilon_trace_over_3,one_minus_sum_beta, &
- alphaval,betaval,gammaval,factor_common,vx,vy,vz,vnspec,COMPUTE_AND_STORE_STRAIN, AM_V)
+ alphaval,betaval,gammaval,factor_common,vx,vy,vz,vnspec,COMPUTE_AND_STORE_STRAIN)
implicit none
@@ -65,8 +65,10 @@
integer :: Qn ! Number of points
end type attenuation_model_variables
- type (attenuation_model_variables) AM_V
+!> Hejun
+! type (attenuation_model_variables) AM_V
! attenuation_model_variables
+!< Hejun
! for forward or backward simulations
logical COMPUTE_AND_STORE_STRAIN
@@ -87,7 +89,9 @@
real(kind=CUSTOM_REAL) one_minus_sum_beta_use,minus_sum_beta
real(kind=CUSTOM_REAL), dimension(vx, vy, vz, vnspec) :: one_minus_sum_beta
- integer iregion_selected
+!>Hejun
+! integer iregion_selected
+!<Hejun
! for attenuation
real(kind=CUSTOM_REAL) R_xx_val,R_yy_val
@@ -172,7 +176,9 @@
! for gravity
integer int_radius
real(kind=CUSTOM_REAL) sigma_yx,sigma_zx,sigma_zy
- real(kind=CUSTOM_REAL) radius_cr
+!> Hejun
+! real(kind=CUSTOM_REAL) radius_cr
+!< Hejun
double precision radius,rho,minus_g,minus_dg
double precision minus_g_over_radius,minus_dg_plus_g_over_radius
double precision cos_theta,sin_theta,cos_phi,sin_phi
@@ -283,17 +289,19 @@
epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl
endif
+!> Hejun
! precompute terms for attenuation if needed
if(ATTENUATION_VAL) then
- if(ATTENUATION_3D_VAL) then
+! if(ATTENUATION_3D_VAL) then
one_minus_sum_beta_use = one_minus_sum_beta(i,j,k,ispec)
- else
- radius_cr = xstore(ibool(i,j,k,ispec))
- call get_attenuation_index(idoubling(ispec), dble(radius_cr), iregion_selected, .FALSE., AM_V)
- one_minus_sum_beta_use = one_minus_sum_beta(1,1,1,iregion_selected)
- endif
+! else
+! radius_cr = xstore(ibool(i,j,k,ispec))
+! call get_attenuation_index(idoubling(ispec), dble(radius_cr), iregion_selected, .FALSE., AM_V)
+! one_minus_sum_beta_use = one_minus_sum_beta(1,1,1,iregion_selected)
+! endif
minus_sum_beta = one_minus_sum_beta_use - 1.0
endif
+!< Hejun
!
! compute either isotropic or anisotropic elements
@@ -821,11 +829,13 @@
! IMPROVE we should probably use an average value instead
! reformatted R_memory to handle large factor_common and reduced [alpha,beta,gamma]val
- if(ATTENUATION_3D_VAL) then
+!> Hejun
+! if(ATTENUATION_3D_VAL) then
factor_common_c44_muv = factor_common(i_sls,:,:,:,ispec)
- else
- factor_common_c44_muv(:,:,:) = factor_common(i_sls,1,1,1,iregion_selected)
- endif
+! else
+! factor_common_c44_muv(:,:,:) = factor_common(i_sls,1,1,1,iregion_selected)
+! endif
+!< Hejun
if(ANISOTROPIC_3D_MANTLE_VAL) then
factor_common_c44_muv = factor_common_c44_muv * c44store(:,:,:,ispec)
else
@@ -867,7 +877,7 @@
c23store,c24store,c25store,c26store,c33store,c34store,c35store, &
c36store,c44store,c45store,c46store,c55store,c56store,c66store, &
ibool,idoubling,R_memory,epsilondev,epsilon_trace_over_3,one_minus_sum_beta, &
- alphaval,betaval,gammaval,factor_common,vx,vy,vz,vnspec,COMPUTE_AND_STORE_STRAIN, AM_V)
+ alphaval,betaval,gammaval,factor_common,vx,vy,vz,vnspec,COMPUTE_AND_STORE_STRAIN)
implicit none
@@ -896,8 +906,10 @@
integer :: Qn ! Number of points
end type attenuation_model_variables
- type (attenuation_model_variables) AM_V
+!> Hejun
+! type (attenuation_model_variables) AM_V
! attenuation_model_variables
+!< Hejun
! for forward or backward simulations
logical COMPUTE_AND_STORE_STRAIN
@@ -917,9 +929,9 @@
real(kind=CUSTOM_REAL) one_minus_sum_beta_use,minus_sum_beta
real(kind=CUSTOM_REAL), dimension(vx, vy, vz, vnspec) :: one_minus_sum_beta
-
- integer iregion_selected
-
+!> Hejun
+! integer iregion_selected
+!< Hejun
! for attenuation
real(kind=CUSTOM_REAL) R_xx_val,R_yy_val
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_memory
@@ -997,7 +1009,9 @@
! for gravity
integer int_radius
real(kind=CUSTOM_REAL) sigma_yx,sigma_zx,sigma_zy
- real(kind=CUSTOM_REAL) radius_cr
+!> Hejun
+! real(kind=CUSTOM_REAL) radius_cr
+!< Hejun
double precision radius,rho,minus_g,minus_dg
double precision minus_g_over_radius,minus_dg_plus_g_over_radius
double precision cos_theta,sin_theta,cos_phi,sin_phi
@@ -1188,19 +1202,19 @@
epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl
epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl
endif
-
+!> Hejun
! precompute terms for attenuation if needed
if(ATTENUATION_VAL) then
- if(ATTENUATION_3D_VAL) then
+! if(ATTENUATION_3D_VAL) then
one_minus_sum_beta_use = one_minus_sum_beta(i,j,k,ispec)
- else
- radius_cr = xstore(ibool(i,j,k,ispec))
- call get_attenuation_index(idoubling(ispec), dble(radius_cr), iregion_selected, .FALSE., AM_V)
- one_minus_sum_beta_use = one_minus_sum_beta(1,1,1,iregion_selected)
- endif
+! else
+! radius_cr = xstore(ibool(i,j,k,ispec))
+! call get_attenuation_index(idoubling(ispec), dble(radius_cr), iregion_selected, .FALSE., AM_V)
+! one_minus_sum_beta_use = one_minus_sum_beta(1,1,1,iregion_selected)
+! endif
minus_sum_beta = one_minus_sum_beta_use - 1.0
endif
-
+!< Hejun
!
! compute either isotropic or anisotropic elements
!
@@ -1753,13 +1767,14 @@
! get coefficients for that standard linear solid
! IMPROVE we use mu_v here even if there is some anisotropy
! IMPROVE we should probably use an average value instead
-
+!> Hejun
! reformatted R_memory to handle large factor_common and reduced [alpha,beta,gamma]val
- if(ATTENUATION_3D_VAL) then
+! if(ATTENUATION_3D_VAL) then
factor_common_c44_muv = factor_common(i_sls,:,:,:,ispec)
- else
- factor_common_c44_muv(:,:,:) = factor_common(i_sls,1,1,1,iregion_selected)
- endif
+! else
+! factor_common_c44_muv(:,:,:) = factor_common(i_sls,1,1,1,iregion_selected)
+! endif
+!< Hejun
if(ANISOTROPIC_3D_MANTLE_VAL) then
factor_common_c44_muv = factor_common_c44_muv * c44store(:,:,:,ispec)
else
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_inner_core.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_inner_core.f90 2010-01-21 16:01:30 UTC (rev 16155)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/compute_forces_inner_core.f90 2010-01-21 16:14:52 UTC (rev 16156)
@@ -34,7 +34,7 @@
kappavstore,muvstore,ibool,idoubling, &
c11store,c33store,c12store,c13store,c44store,R_memory,epsilondev,epsilon_trace_over_3,&
one_minus_sum_beta,alphaval,betaval,gammaval,factor_common, &
- vx,vy,vz,vnspec,COMPUTE_AND_STORE_STRAIN, AM_V)
+ vx,vy,vz,vnspec,COMPUTE_AND_STORE_STRAIN)
implicit none
@@ -63,8 +63,10 @@
integer :: Qn ! Number of points
end type attenuation_model_variables
- type (attenuation_model_variables) AM_V
+!> Hejun
+! type (attenuation_model_variables) AM_V
! attenuation_model_variables
+!< Hejun
! for forward or backward simulations
logical COMPUTE_AND_STORE_STRAIN
@@ -153,10 +155,11 @@
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: rho_s_H
double precision, dimension(NGLLX,NGLLY,NGLLZ) :: wgll_cube
real(kind=CUSTOM_REAL), dimension(NGLOB_INNER_CORE) :: xstore,ystore,zstore
- integer iregion_selected
+!> Hejun
+! integer iregion_selected
+! real(kind=CUSTOM_REAL) radius_cr
+!< Hejun
- real(kind=CUSTOM_REAL) radius_cr
-
! ****************************************************
! big loop over all spectral elements in the solid
! ****************************************************
@@ -261,15 +264,17 @@
epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl
endif
+!> Hejun
if(ATTENUATION_VAL) then
- if(ATTENUATION_3D_VAL) then
+! if(ATTENUATION_3D_VAL) then
minus_sum_beta = one_minus_sum_beta(i,j,k,ispec) - 1.0
- else
- radius_cr = xstore(ibool(i,j,k,ispec))
- call get_attenuation_index(idoubling(ispec), dble(radius_cr), iregion_selected, .TRUE., AM_V)
- minus_sum_beta = one_minus_sum_beta(1,1,1,iregion_selected) - 1.0
- endif ! ATTENUATION_3D_VAL
+! else
+! radius_cr = xstore(ibool(i,j,k,ispec))
+! call get_attenuation_index(idoubling(ispec), dble(radius_cr), iregion_selected, .TRUE., AM_V)
+! minus_sum_beta = one_minus_sum_beta(1,1,1,iregion_selected) - 1.0
+! endif ! ATTENUATION_3D_VAL
endif ! ATTENUATION_VAL
+!< Hejun
if(ANISOTROPIC_INNER_CORE_VAL) then
@@ -319,15 +324,16 @@
kappal = kappavstore(i,j,k,ispec)
mul = muvstore(i,j,k,ispec)
+!> Hejun
! use unrelaxed parameters if attenuation
if(ATTENUATION_VAL) then
- if(ATTENUATION_3D_VAL) then
+! if(ATTENUATION_3D_VAL) then
mul = mul * one_minus_sum_beta(i,j,k,ispec)
- else
- mul = mul * one_minus_sum_beta(1,1,1,iregion_selected)
- endif
+! else
+! mul = mul * one_minus_sum_beta(1,1,1,iregion_selected)
+! endif
endif
-
+!< Hejun
lambdalplus2mul = kappal + FOUR_THIRDS * mul
lambdal = lambdalplus2mul - 2.*mul
@@ -572,11 +578,13 @@
if(ATTENUATION_VAL) then
do i_sls = 1,N_SLS
- if(ATTENUATION_3D_VAL) then
+!> Hejun
+! if(ATTENUATION_3D_VAL) then
factor_common_use = factor_common(i_sls,:,:,:,ispec)
- else
- factor_common_use(:,:,:) = factor_common(i_sls,1,1,1,iregion_selected)
- endif
+! else
+! factor_common_use(:,:,:) = factor_common(i_sls,1,1,1,iregion_selected)
+! endif
+!< Hejun
do i_memory = 1,5
R_memory(i_memory,i_sls,:,:,:,ispec) = &
alphaval(i_sls) * &
@@ -613,7 +621,7 @@
kappavstore,muvstore,ibool,idoubling, &
c11store,c33store,c12store,c13store,c44store,R_memory,epsilondev,epsilon_trace_over_3,&
one_minus_sum_beta,alphaval,betaval,gammaval,factor_common, &
- vx,vy,vz,vnspec,COMPUTE_AND_STORE_STRAIN, AM_V)
+ vx,vy,vz,vnspec,COMPUTE_AND_STORE_STRAIN)
implicit none
@@ -642,8 +650,10 @@
integer :: Qn ! Number of points
end type attenuation_model_variables
- type (attenuation_model_variables) AM_V
+!> Hejun
+! type (attenuation_model_variables) AM_V
! attenuation_model_variables
+!< Hejun
! for forward or backward simulations
logical COMPUTE_AND_STORE_STRAIN
@@ -726,10 +736,11 @@
real(kind=CUSTOM_REAL), dimension(NDIM,NGLLX,NGLLY,NGLLZ) :: rho_s_H
double precision, dimension(NGLLX,NGLLY,NGLLZ) :: wgll_cube
real(kind=CUSTOM_REAL), dimension(NGLOB_INNER_CORE) :: xstore,ystore,zstore
- integer iregion_selected
+!> Hejun
+! integer iregion_selected
+! real(kind=CUSTOM_REAL) radius_cr
+!< Hejun
- real(kind=CUSTOM_REAL) radius_cr
-
! Deville
! manually inline the calls to the Deville et al. (2002) routines
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ) :: dummyx_loc,dummyy_loc,dummyz_loc, &
@@ -916,17 +927,17 @@
epsilondev_loc(4,i,j,k) = 0.5 * duzdxl_plus_duxdzl
epsilondev_loc(5,i,j,k) = 0.5 * duzdyl_plus_duydzl
endif
-
+!> Hejun
if(ATTENUATION_VAL) then
- if(ATTENUATION_3D_VAL) then
+! if(ATTENUATION_3D_VAL) then
minus_sum_beta = one_minus_sum_beta(i,j,k,ispec) - 1.0
- else
- radius_cr = xstore(ibool(i,j,k,ispec))
- call get_attenuation_index(idoubling(ispec), dble(radius_cr), iregion_selected, .TRUE., AM_V)
- minus_sum_beta = one_minus_sum_beta(1,1,1,iregion_selected) - 1.0
- endif ! ATTENUATION_3D_VAL
+! else
+! radius_cr = xstore(ibool(i,j,k,ispec))
+! call get_attenuation_index(idoubling(ispec), dble(radius_cr), iregion_selected, .TRUE., AM_V)
+! minus_sum_beta = one_minus_sum_beta(1,1,1,iregion_selected) - 1.0
+! endif ! ATTENUATION_3D_VAL
endif ! ATTENUATION_VAL
-
+!< Hejun
if(ANISOTROPIC_INNER_CORE_VAL) then
! elastic tensor for hexagonal symmetry in reduced notation:
!
@@ -972,16 +983,16 @@
! layer with no anisotropy, use kappav and muv for instance
kappal = kappavstore(i,j,k,ispec)
mul = muvstore(i,j,k,ispec)
-
+!> Hejun
! use unrelaxed parameters if attenuation
if(ATTENUATION_VAL) then
- if(ATTENUATION_3D_VAL) then
+! if(ATTENUATION_3D_VAL) then
mul = mul * one_minus_sum_beta(i,j,k,ispec)
- else
- mul = mul * one_minus_sum_beta(1,1,1,iregion_selected)
- endif
+! else
+! mul = mul * one_minus_sum_beta(1,1,1,iregion_selected)
+! endif
endif
-
+!< Hejun
lambdalplus2mul = kappal + FOUR_THIRDS * mul
lambdal = lambdalplus2mul - 2.*mul
@@ -1255,13 +1266,14 @@
! therefore Q_\alpha is not zero; for instance for V_p / V_s = sqrt(3)
! we get Q_\alpha = (9 / 4) * Q_\mu = 2.25 * Q_\mu
if(ATTENUATION_VAL) then
-
+!> Hejun
do i_sls = 1,N_SLS
- if(ATTENUATION_3D_VAL) then
+! if(ATTENUATION_3D_VAL) then
factor_common_use = factor_common(i_sls,:,:,:,ispec)
- else
- factor_common_use(:,:,:) = factor_common(i_sls,1,1,1,iregion_selected)
- endif
+! else
+! factor_common_use(:,:,:) = factor_common(i_sls,1,1,1,iregion_selected)
+! endif
+!< Hejun
do i_memory = 1,5
R_memory(i_memory,i_sls,:,:,:,ispec) = &
alphaval(i_sls) * &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in 2010-01-21 16:01:30 UTC (rev 16155)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/constants.h.in 2010-01-21 16:14:52 UTC (rev 16156)
@@ -102,6 +102,17 @@
!! pathname of the topography file
! character (len=*), parameter :: PATHNAME_TOPO_FILE = 'DATA/topo_bathy/topo_bathy_etopo2_smoothed_window7.dat'
+!!--- ETOPO1 1-minute model, implemented now, data must be created first
+!! size of topography and bathymetry file
+! integer, parameter :: NX_BATHY = 21600,NY_BATHY = 10800
+!! resolution of topography file in minutes
+! integer, parameter :: RESOLUTION_TOPO_FILE = 1
+!! pathname of the topography file
+! character (len=*), parameter :: PATHNAME_TOPO_FILE = 'DATA/topo_bathy/ETOPO1.xyz'
+
+! Use GLL points to capture TOPOGRAPHY and ELLIPTICITY
+ logical,parameter :: USE_GLL = .false.
+
! maximum depth of the oceans in trenches and height of topo in mountains
! to avoid taking into account spurious oscillations in global model ETOPO
logical, parameter :: USE_MAXIMUM_HEIGHT_TOPO = .false.
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90 2010-01-21 16:01:30 UTC (rev 16155)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90 2010-01-21 16:14:52 UTC (rev 16156)
@@ -41,7 +41,8 @@
rotation_matrix,ANGULAR_WIDTH_XI_RAD,ANGULAR_WIDTH_ETA_RAD,&
ATTENUATION,ATTENUATION_3D,SAVE_MESH_FILES, &
NCHUNKS,INCLUDE_CENTRAL_CUBE,ABSORBING_CONDITIONS,REFERENCE_1D_MODEL,THREE_D_MODEL, &
- R_CENTRAL_CUBE,RICB,RHO_OCEANS,RCMB,R670,RMOHO,RTOPDDOUBLEPRIME,R600,R220,R771,R400,R120,R80,RMIDDLE_CRUST,ROCEAN, &
+ R_CENTRAL_CUBE,RICB,RHO_OCEANS,RCMB,R670,RMOHO,RMOHO_FICTITIOUS_IN_MESHER,&
+ RTOPDDOUBLEPRIME,R600,R220,R771,R400,R120,R80,RMIDDLE_CRUST,ROCEAN, &
ner,ratio_sampling_array,doubling_index,r_bottom,r_top,this_region_has_a_doubling,CASE_3D, &
AMM_V,AM_V,QRFSI12_Q,M1066a_V,Mak135_V,Mref_V,SEA1DM_V,D3MM_V,HMM,JP3DM_V,SEA99M_V,CM_V, AM_S, AS_V, &
numker,numhpa,numcof,ihpa,lmax,nylm, &
@@ -353,6 +354,9 @@
double precision R_CENTRAL_CUBE,RICB,RHO_OCEANS,RCMB,R670,RMOHO, &
RTOPDDOUBLEPRIME,R600,R220,R771,R400,R120,R80,RMIDDLE_CRUST,ROCEAN
+!> Hejun
+ integer:: RMOHO_FICTITIOUS_IN_MESHER
+!< Hejun
character(len=150) LOCAL_PATH,errmsg
@@ -564,19 +568,30 @@
! create the name for the database of the current slide and region
call create_name_database(prname,myrank,iregion_code,LOCAL_PATH)
+!> Hejun
+! New Attenuation defination
! Attenuation
- if(ATTENUATION .and. ATTENUATION_3D) then
- T_c_source = AM_V%QT_c_source
- tau_s(:) = AM_V%Qtau_s(:)
- allocate(Qmu_store(NGLLX,NGLLY,NGLLZ,nspec))
- allocate(tau_e_store(N_SLS,NGLLX,NGLLY,NGLLZ,nspec))
- else
- allocate(Qmu_store(1,1,1,1))
- allocate(tau_e_store(N_SLS,1,1,1,1))
- Qmu_store(1,1,1,1) = 0.0d0
- tau_e_store(:,1,1,1,1) = 0.0d0
- endif
+! if(ATTENUATION .and. ATTENUATION_3D) then
+! T_c_source = AM_V%QT_c_source
+! tau_s(:) = AM_V%Qtau_s(:)
+! allocate(Qmu_store(NGLLX,NGLLY,NGLLZ,nspec))
+! allocate(tau_e_store(N_SLS,NGLLX,NGLLY,NGLLZ,nspec))
+! else
+! allocate(Qmu_store(1,1,1,1))
+! allocate(tau_e_store(N_SLS,1,1,1,1))
+! Qmu_store(1,1,1,1) = 0.0d0
+! tau_e_store(:,1,1,1,1) = 0.0d0
+! endif
+ if (ATTENUATION) then
+ T_c_source = AM_V%QT_c_source
+ tau_s(:) = AM_V%Qtau_s(:)
+ allocate(Qmu_store(NGLLX,NGLLY,NGLLZ,nspec))
+ allocate(tau_e_store(N_SLS,NGLLX,NGLLY,NGLLZ,nspec))
+ end if
+!< Hejun
+
+
! Gauss-Lobatto-Legendre points of integration
allocate(xigll(NGLLX))
allocate(yigll(NGLLY))
@@ -980,7 +995,7 @@
ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE,ISOTROPIC_3D_MANTLE,HETEROGEN_3D_MANTLE,CRUSTAL,ONE_CRUST, &
myrank,ibathy_topo,ATTENUATION,ATTENUATION_3D, &
ABSORBING_CONDITIONS,REFERENCE_1D_MODEL,THREE_D_MODEL, &
- RICB,RCMB,R670,RMOHO,RTOPDDOUBLEPRIME,R600,R220,R771,R400,R120,R80,RMIDDLE_CRUST,ROCEAN, &
+ RICB,RCMB,R670,RMOHO,RMOHO_FICTITIOUS_IN_MESHER,RTOPDDOUBLEPRIME,R600,R220,R771,R400,R120,R80,RMIDDLE_CRUST,ROCEAN, &
xelm,yelm,zelm,shape3D,dershape3D,rmin,rmax,rhostore,dvpstore,kappavstore,kappahstore,muvstore,muhstore,eta_anisostore, &
xixstore,xiystore,xizstore,etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore, &
c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
@@ -991,7 +1006,8 @@
numker,numhpa,numcof,ihpa,lmax,nylm, &
lmxhpa,itypehpa,ihpakern,numcoe,ivarkern, &
nconpt,iver,iconpt,conpt,xlaspl,xlospl,radspl, &
- coe,vercof,vercofd,ylmcof,wk1,wk2,wk3,kerstr,varstr,ACTUALLY_STORE_ARRAYS)
+ coe,vercof,vercofd,ylmcof,wk1,wk2,wk3,kerstr,varstr,ACTUALLY_STORE_ARRAYS,&
+ xigll,yigll,zigll)
! boundary mesh
if (ipass == 2 .and. SAVE_BOUNDARY_MESH .and. iregion_code == IREGION_CRUST_MANTLE) then
@@ -1182,7 +1198,7 @@
ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE,ISOTROPIC_3D_MANTLE,HETEROGEN_3D_MANTLE,CRUSTAL,ONE_CRUST, &
myrank,ibathy_topo,ATTENUATION,ATTENUATION_3D, &
ABSORBING_CONDITIONS,REFERENCE_1D_MODEL,THREE_D_MODEL, &
- RICB,RCMB,R670,RMOHO,RTOPDDOUBLEPRIME,R600,R220,R771,R400,R120,R80,RMIDDLE_CRUST,ROCEAN, &
+ RICB,RCMB,R670,RMOHO,RMOHO_FICTITIOUS_IN_MESHER,RTOPDDOUBLEPRIME,R600,R220,R771,R400,R120,R80,RMIDDLE_CRUST,ROCEAN, &
xelm,yelm,zelm,shape3D,dershape3D,rmin,rmax,rhostore,dvpstore,kappavstore,kappahstore,muvstore,muhstore,eta_anisostore, &
xixstore,xiystore,xizstore,etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore, &
c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
@@ -1193,7 +1209,8 @@
numker,numhpa,numcof,ihpa,lmax,nylm, &
lmxhpa,itypehpa,ihpakern,numcoe,ivarkern, &
nconpt,iver,iconpt,conpt,xlaspl,xlospl,radspl, &
- coe,vercof,vercofd,ylmcof,wk1,wk2,wk3,kerstr,varstr,ACTUALLY_STORE_ARRAYS)
+ coe,vercof,vercofd,ylmcof,wk1,wk2,wk3,kerstr,varstr,ACTUALLY_STORE_ARRAYS,&
+ xigll,yigll,zigll)
! boundary mesh
if (ipass == 2 .and. SAVE_BOUNDARY_MESH .and. iregion_code == IREGION_CRUST_MANTLE) then
@@ -1351,7 +1368,7 @@
ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE,ISOTROPIC_3D_MANTLE,HETEROGEN_3D_MANTLE,CRUSTAL,ONE_CRUST, &
myrank,ibathy_topo,ATTENUATION,ATTENUATION_3D, &
ABSORBING_CONDITIONS,REFERENCE_1D_MODEL,THREE_D_MODEL, &
- RICB,RCMB,R670,RMOHO,RTOPDDOUBLEPRIME,R600,R220,R771,R400,R120,R80,RMIDDLE_CRUST,ROCEAN, &
+ RICB,RCMB,R670,RMOHO,RMOHO_FICTITIOUS_IN_MESHER,RTOPDDOUBLEPRIME,R600,R220,R771,R400,R120,R80,RMIDDLE_CRUST,ROCEAN, &
xelm,yelm,zelm,shape3D,dershape3D,rmin,rmax,rhostore,dvpstore,kappavstore,kappahstore,muvstore,muhstore,eta_anisostore, &
xixstore,xiystore,xizstore,etaxstore,etaystore,etazstore,gammaxstore,gammaystore,gammazstore, &
c11store,c12store,c13store,c14store,c15store,c16store,c22store, &
@@ -1362,7 +1379,8 @@
numker,numhpa,numcof,ihpa,lmax,nylm, &
lmxhpa,itypehpa,ihpakern,numcoe,ivarkern, &
nconpt,iver,iconpt,conpt,xlaspl,xlospl,radspl, &
- coe,vercof,vercofd,ylmcof,wk1,wk2,wk3,kerstr,varstr,ACTUALLY_STORE_ARRAYS)
+ coe,vercof,vercofd,ylmcof,wk1,wk2,wk3,kerstr,varstr,ACTUALLY_STORE_ARRAYS,&
+ xigll,yigll,zigll)
enddo
enddo
enddo
@@ -1486,6 +1504,12 @@
RMIDDLE_CRUST,ROCEAN,M1066a_V,Mak135_V,Mref_V,SEA1DM_V)
call write_AVS_DX_surface_data(myrank,prname,nspec,iboun,ibool, &
idoubling,xstore,ystore,zstore,locval,ifseg,npointot)
+!> Hejun
+! Output material information for all GLL points
+! Can be use to check the mesh
+! call write_AVS_DX_global_data_gll(prname,nspec,xstore,ystore,zstore,&
+! rhostore,kappavstore,muvstore,Qmu_store)
+!< Hejun
endif
deallocate(locval,stat=ier); if(ier /= 0) stop 'error in deallocate'
@@ -1512,7 +1536,8 @@
normal_ymin,normal_ymax, &
normal_bottom,normal_top, &
NSPEC2D_BOTTOM,NSPEC2D_TOP, &
- NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX)
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,&
+ xigll,yigll,zigll)
! creating mass matrix in this slice (will be fully assembled in the solver)
allocate(rmass(nglob),stat=ier); if(ier /= 0) stop 'error in allocate'
@@ -1685,9 +1710,9 @@
jacobian2D_bottom,jacobian2D_top,nspec,nglob, &
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
TRANSVERSE_ISOTROPY,HETEROGEN_3D_MANTLE,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE,OCEANS, &
- tau_s,tau_e_store,Qmu_store,T_c_source,ATTENUATION,ATTENUATION_3D, &
+ tau_s,tau_e_store,Qmu_store,T_c_source,ATTENUATION, &
size(tau_e_store,2),size(tau_e_store,3),size(tau_e_store,4),size(tau_e_store,5),&
- NEX_PER_PROC_XI,NEX_PER_PROC_ETA,NEX_XI,ichunk,NCHUNKS,ABSORBING_CONDITIONS,AM_V)
+ NEX_PER_PROC_XI,NEX_PER_PROC_ETA,NEX_XI,ichunk,NCHUNKS,ABSORBING_CONDITIONS)
deallocate(rmass,stat=ier); if(ier /= 0) stop 'error in deallocate'
deallocate(rmass_ocean_load,stat=ier); if(ier /= 0) stop 'error in deallocate'
@@ -1831,8 +1856,12 @@
deallocate(nimin,nimax,njmin,njmax,nkmin_xi,nkmin_eta)
deallocate(rho_vp,rho_vs)
- deallocate(Qmu_store)
- deallocate(tau_e_store)
+!> Hejun
+ if (ATTENUATION) then
+ deallocate(Qmu_store)
+ deallocate(tau_e_store)
+ end if
+!< Hejun
end subroutine create_regions_mesh
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/crustal_model.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/crustal_model.f90 2010-01-21 16:01:30 UTC (rev 16155)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/crustal_model.f90 2010-01-21 16:14:52 UTC (rev 16156)
@@ -30,7 +30,7 @@
! based on software routines provided with the crust2.0 model by Bassin et al.
!
- subroutine crustal_model(lat,lon,x,vp,vs,rho,moho,found_crust,CM_V)
+ subroutine crustal_model(lat,lon,x,vp,vs,rho,moho,found_crust,CM_V,elem_in_crust)
implicit none
include "constants.h"
@@ -56,6 +56,9 @@
double precision x3,x4,x5,x6,x7,scaleval
double precision vps(NLAYERS_CRUST),vss(NLAYERS_CRUST),rhos(NLAYERS_CRUST),thicks(NLAYERS_CRUST)
+!> Hejun
+ logical::elem_in_crust
+!< Hejun
call crust(lat,lon,vps,vss,rhos,thicks,CM_V%abbreviation,CM_V%code,CM_V%thlr,CM_V%velocp,CM_V%velocs,CM_V%dens)
x3 = (R_EARTH-thicks(3)*1000.0d0)/R_EARTH
@@ -83,7 +86,7 @@
vp = vps(6)
vs = vss(6)
rho = rhos(6)
- else if(x > x7) then
+ else if(x > x7.or.elem_in_crust) then
vp = vps(7)
vs = vss(7)
rho = rhos(7)
@@ -97,9 +100,13 @@
vp = vp*1000.0d0/(R_EARTH*scaleval)
vs = vs*1000.0d0/(R_EARTH*scaleval)
rho = rho*1000.0d0/RHOAV
- moho = (h_uc+thicks(6)+thicks(7))*1000.0d0/R_EARTH
+! moho = (h_uc+thicks(6)+thicks(7))*1000.0d0/R_EARTH
endif
+ !> Hejun
+ ! No matter found_crust true or false, output moho thickness
+ moho = (h_uc+thicks(6)+thicks(7))*1000.0d0/R_EARTH
+ !< Hejun
end subroutine crustal_model
!---------------------------
@@ -171,11 +178,12 @@
implicit none
include "constants.h"
+!> Hejun
+! Change the CAP function to smooth crustal model
+ integer, parameter :: NTHETA = 4 !2
+ integer, parameter :: NPHI = 20 !10
+ double precision, parameter :: CAP = 1.0d0*PI/180.0d0 ! 2.0d0*PI/180.0d0
- integer, parameter :: NTHETA = 2
- integer, parameter :: NPHI = 10
- double precision, parameter :: CAP = 2.0d0*PI/180.0d0
-
! argument variables
double precision lat,lon
double precision rho(NLAYERS_CRUST),thick(NLAYERS_CRUST),velp(NLAYERS_CRUST),vels(NLAYERS_CRUST)
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/get_ellipticity.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/get_ellipticity.f90 2010-01-21 16:01:30 UTC (rev 16155)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/get_ellipticity.f90 2010-01-21 16:14:52 UTC (rev 16156)
@@ -63,3 +63,47 @@
end subroutine get_ellipticity
+ !> Hejun
+ ! get ellipticity according to GLL points
+ ! JAN08, 2010
+ subroutine get_ellipticity_gll(xstore,ystore,zstore,ispec,nspec,nspl,rspl,espl,espl2)
+
+ implicit none
+
+ include "constants.h"
+
+ integer nspl
+ integer::ispec,nspec
+ double precision,dimension(NGLLX,NGLLY,NGLLZ,nspec):: xstore,ystore,zstore
+ double precision rspl(NR),espl(NR),espl2(NR)
+
+ integer i,j,k
+
+ double precision ell
+ double precision r,theta,phi,factor
+ double precision cost,p20
+
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+
+ call xyz_2_rthetaphi_dble(xstore(i,j,k,ispec),ystore(i,j,k,ispec),zstore(i,j,k,ispec),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
+
+ xstore(i,j,k,ispec)=xstore(i,j,k,ispec)*factor
+ ystore(i,j,k,ispec)=ystore(i,j,k,ispec)*factor
+ zstore(i,j,k,ispec)=zstore(i,j,k,ispec)*factor
+
+ end do
+ end do
+ end do
+ end subroutine get_ellipticity_gll
+ !< Hejun
+
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/get_jacobian_boundaries.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/get_jacobian_boundaries.f90 2010-01-21 16:01:30 UTC (rev 16155)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/get_jacobian_boundaries.f90 2010-01-21 16:14:52 UTC (rev 16156)
@@ -36,7 +36,7 @@
normal_ymin,normal_ymax, &
normal_bottom,normal_top, &
NSPEC2D_BOTTOM,NSPEC2D_TOP, &
- NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX)
+ NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,xigll,yigll,zigll)
implicit none
@@ -83,6 +83,15 @@
double precision xelm(NGNOD2D),yelm(NGNOD2D),zelm(NGNOD2D)
+!> Hejun
+! Parameters used to calculate 2D Jacobian based upon 25 GLL points
+ integer:: i,j,k
+ double precision xelm2D(NGLLX,NGLLY),yelm2D(NGLLX,NGLLY),zelm2D(NGLLX,NGLLY)
+ double precision,dimension(NGLLX):: xigll
+ double precision,dimension(NGLLY):: yigll
+ double precision,dimension(NGLLZ):: zigll
+!< Hejun
+
! check that the parameter file is correct
if(NGNOD /= 27) call exit_MPI(myrank,'elements should have 27 control nodes')
if(NGNOD2D /= 9) call exit_MPI(myrank,'surface elements should have 9 control nodes')
@@ -94,6 +103,7 @@
ispecb5 = 0
ispecb6 = 0
+!> Hejun
do ispec=1,nspec
! determine if the element falls on a boundary
@@ -105,38 +115,52 @@
ispecb1=ispecb1+1
ibelm_xmin(ispecb1)=ispec
-! specify the 9 nodes for the 2-D boundary element
- xelm(1)=xstore(1,1,1,ispec)
- yelm(1)=ystore(1,1,1,ispec)
- zelm(1)=zstore(1,1,1,ispec)
- xelm(2)=xstore(1,NGLLY,1,ispec)
- yelm(2)=ystore(1,NGLLY,1,ispec)
- zelm(2)=zstore(1,NGLLY,1,ispec)
- xelm(3)=xstore(1,NGLLY,NGLLZ,ispec)
- yelm(3)=ystore(1,NGLLY,NGLLZ,ispec)
- zelm(3)=zstore(1,NGLLY,NGLLZ,ispec)
- xelm(4)=xstore(1,1,NGLLZ,ispec)
- yelm(4)=ystore(1,1,NGLLZ,ispec)
- zelm(4)=zstore(1,1,NGLLZ,ispec)
- xelm(5)=xstore(1,(NGLLY+1)/2,1,ispec)
- yelm(5)=ystore(1,(NGLLY+1)/2,1,ispec)
- zelm(5)=zstore(1,(NGLLY+1)/2,1,ispec)
- xelm(6)=xstore(1,NGLLY,(NGLLZ+1)/2,ispec)
- yelm(6)=ystore(1,NGLLY,(NGLLZ+1)/2,ispec)
- zelm(6)=zstore(1,NGLLY,(NGLLZ+1)/2,ispec)
- xelm(7)=xstore(1,(NGLLY+1)/2,NGLLZ,ispec)
- yelm(7)=ystore(1,(NGLLY+1)/2,NGLLZ,ispec)
- zelm(7)=zstore(1,(NGLLY+1)/2,NGLLZ,ispec)
- xelm(8)=xstore(1,1,(NGLLZ+1)/2,ispec)
- yelm(8)=ystore(1,1,(NGLLZ+1)/2,ispec)
- zelm(8)=zstore(1,1,(NGLLZ+1)/2,ispec)
- xelm(9)=xstore(1,(NGLLY+1)/2,(NGLLZ+1)/2,ispec)
- yelm(9)=ystore(1,(NGLLY+1)/2,(NGLLZ+1)/2,ispec)
- zelm(9)=zstore(1,(NGLLY+1)/2,(NGLLZ+1)/2,ispec)
+ if ( .not. USE_GLL) then
+ ! specify the 9 nodes for the 2-D boundary element
+ xelm(1)=xstore(1,1,1,ispec)
+ yelm(1)=ystore(1,1,1,ispec)
+ zelm(1)=zstore(1,1,1,ispec)
+ xelm(2)=xstore(1,NGLLY,1,ispec)
+ yelm(2)=ystore(1,NGLLY,1,ispec)
+ zelm(2)=zstore(1,NGLLY,1,ispec)
+ xelm(3)=xstore(1,NGLLY,NGLLZ,ispec)
+ yelm(3)=ystore(1,NGLLY,NGLLZ,ispec)
+ zelm(3)=zstore(1,NGLLY,NGLLZ,ispec)
+ xelm(4)=xstore(1,1,NGLLZ,ispec)
+ yelm(4)=ystore(1,1,NGLLZ,ispec)
+ zelm(4)=zstore(1,1,NGLLZ,ispec)
+ xelm(5)=xstore(1,(NGLLY+1)/2,1,ispec)
+ yelm(5)=ystore(1,(NGLLY+1)/2,1,ispec)
+ zelm(5)=zstore(1,(NGLLY+1)/2,1,ispec)
+ xelm(6)=xstore(1,NGLLY,(NGLLZ+1)/2,ispec)
+ yelm(6)=ystore(1,NGLLY,(NGLLZ+1)/2,ispec)
+ zelm(6)=zstore(1,NGLLY,(NGLLZ+1)/2,ispec)
+ xelm(7)=xstore(1,(NGLLY+1)/2,NGLLZ,ispec)
+ yelm(7)=ystore(1,(NGLLY+1)/2,NGLLZ,ispec)
+ zelm(7)=zstore(1,(NGLLY+1)/2,NGLLZ,ispec)
+ xelm(8)=xstore(1,1,(NGLLZ+1)/2,ispec)
+ yelm(8)=ystore(1,1,(NGLLZ+1)/2,ispec)
+ zelm(8)=zstore(1,1,(NGLLZ+1)/2,ispec)
+ xelm(9)=xstore(1,(NGLLY+1)/2,(NGLLZ+1)/2,ispec)
+ yelm(9)=ystore(1,(NGLLY+1)/2,(NGLLZ+1)/2,ispec)
+ zelm(9)=zstore(1,(NGLLY+1)/2,(NGLLZ+1)/2,ispec)
- call compute_jacobian_2D(myrank,ispecb1,xelm,yelm,zelm,dershape2D_x, &
+ call compute_jacobian_2D(myrank,ispecb1,xelm,yelm,zelm,dershape2D_x, &
jacobian2D_xmin,normal_xmin,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
-
+ else
+ ! get 25 GLL points for xmin
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ xelm2D(j,k) = xstore(1,j,k,ispec)
+ yelm2D(j,k) = ystore(1,j,k,ispec)
+ zelm2D(j,k) = zstore(1,j,k,ispec)
+ end do
+ end do
+ ! recalculate jacobian according to 2D GLL points
+ call recalc_jacobian_gll2D(myrank,ispecb1,xelm2D,yelm2D,zelm2D, &
+ yigll,zigll,jacobian2D_xmin,normal_xmin,&
+ NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
+ end if
endif
! on boundary: xmax
@@ -146,38 +170,53 @@
ispecb2=ispecb2+1
ibelm_xmax(ispecb2)=ispec
-! specify the 9 nodes for the 2-D boundary element
- xelm(1)=xstore(NGLLX,1,1,ispec)
- yelm(1)=ystore(NGLLX,1,1,ispec)
- zelm(1)=zstore(NGLLX,1,1,ispec)
- xelm(2)=xstore(NGLLX,NGLLY,1,ispec)
- yelm(2)=ystore(NGLLX,NGLLY,1,ispec)
- zelm(2)=zstore(NGLLX,NGLLY,1,ispec)
- xelm(3)=xstore(NGLLX,NGLLY,NGLLZ,ispec)
- yelm(3)=ystore(NGLLX,NGLLY,NGLLZ,ispec)
- zelm(3)=zstore(NGLLX,NGLLY,NGLLZ,ispec)
- xelm(4)=xstore(NGLLX,1,NGLLZ,ispec)
- yelm(4)=ystore(NGLLX,1,NGLLZ,ispec)
- zelm(4)=zstore(NGLLX,1,NGLLZ,ispec)
- xelm(5)=xstore(NGLLX,(NGLLY+1)/2,1,ispec)
- yelm(5)=ystore(NGLLX,(NGLLY+1)/2,1,ispec)
- zelm(5)=zstore(NGLLX,(NGLLY+1)/2,1,ispec)
- xelm(6)=xstore(NGLLX,NGLLY,(NGLLZ+1)/2,ispec)
- yelm(6)=ystore(NGLLX,NGLLY,(NGLLZ+1)/2,ispec)
- zelm(6)=zstore(NGLLX,NGLLY,(NGLLZ+1)/2,ispec)
- xelm(7)=xstore(NGLLX,(NGLLY+1)/2,NGLLZ,ispec)
- yelm(7)=ystore(NGLLX,(NGLLY+1)/2,NGLLZ,ispec)
- zelm(7)=zstore(NGLLX,(NGLLY+1)/2,NGLLZ,ispec)
- xelm(8)=xstore(NGLLX,1,(NGLLZ+1)/2,ispec)
- yelm(8)=ystore(NGLLX,1,(NGLLZ+1)/2,ispec)
- zelm(8)=zstore(NGLLX,1,(NGLLZ+1)/2,ispec)
- xelm(9)=xstore(NGLLX,(NGLLY+1)/2,(NGLLZ+1)/2,ispec)
- yelm(9)=ystore(NGLLX,(NGLLY+1)/2,(NGLLZ+1)/2,ispec)
- zelm(9)=zstore(NGLLX,(NGLLY+1)/2,(NGLLZ+1)/2,ispec)
+ if ( .not. USE_GLL) then
+ ! specify the 9 nodes for the 2-D boundary element
+ xelm(1)=xstore(NGLLX,1,1,ispec)
+ yelm(1)=ystore(NGLLX,1,1,ispec)
+ zelm(1)=zstore(NGLLX,1,1,ispec)
+ xelm(2)=xstore(NGLLX,NGLLY,1,ispec)
+ yelm(2)=ystore(NGLLX,NGLLY,1,ispec)
+ zelm(2)=zstore(NGLLX,NGLLY,1,ispec)
+ xelm(3)=xstore(NGLLX,NGLLY,NGLLZ,ispec)
+ yelm(3)=ystore(NGLLX,NGLLY,NGLLZ,ispec)
+ zelm(3)=zstore(NGLLX,NGLLY,NGLLZ,ispec)
+ xelm(4)=xstore(NGLLX,1,NGLLZ,ispec)
+ yelm(4)=ystore(NGLLX,1,NGLLZ,ispec)
+ zelm(4)=zstore(NGLLX,1,NGLLZ,ispec)
+ xelm(5)=xstore(NGLLX,(NGLLY+1)/2,1,ispec)
+ yelm(5)=ystore(NGLLX,(NGLLY+1)/2,1,ispec)
+ zelm(5)=zstore(NGLLX,(NGLLY+1)/2,1,ispec)
+ xelm(6)=xstore(NGLLX,NGLLY,(NGLLZ+1)/2,ispec)
+ yelm(6)=ystore(NGLLX,NGLLY,(NGLLZ+1)/2,ispec)
+ zelm(6)=zstore(NGLLX,NGLLY,(NGLLZ+1)/2,ispec)
+ xelm(7)=xstore(NGLLX,(NGLLY+1)/2,NGLLZ,ispec)
+ yelm(7)=ystore(NGLLX,(NGLLY+1)/2,NGLLZ,ispec)
+ zelm(7)=zstore(NGLLX,(NGLLY+1)/2,NGLLZ,ispec)
+ xelm(8)=xstore(NGLLX,1,(NGLLZ+1)/2,ispec)
+ yelm(8)=ystore(NGLLX,1,(NGLLZ+1)/2,ispec)
+ zelm(8)=zstore(NGLLX,1,(NGLLZ+1)/2,ispec)
+ xelm(9)=xstore(NGLLX,(NGLLY+1)/2,(NGLLZ+1)/2,ispec)
+ yelm(9)=ystore(NGLLX,(NGLLY+1)/2,(NGLLZ+1)/2,ispec)
+ zelm(9)=zstore(NGLLX,(NGLLY+1)/2,(NGLLZ+1)/2,ispec)
- call compute_jacobian_2D(myrank,ispecb2,xelm,yelm,zelm,dershape2D_x, &
+ call compute_jacobian_2D(myrank,ispecb2,xelm,yelm,zelm,dershape2D_x, &
jacobian2D_xmax,normal_xmax,NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
+ else
+ ! get 25 GLL points for xmax
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ xelm2D(j,k) = xstore(NGLLX,j,k,ispec)
+ yelm2D(j,k) = ystore(NGLLX,j,k,ispec)
+ zelm2D(j,k) = zstore(NGLLX,j,k,ispec)
+ end do
+ end do
+ ! recalculate jacobian according to 2D GLL points
+ call recalc_jacobian_gll2D(myrank,ispecb2,xelm2D,yelm2D,zelm2D,&
+ yigll,zigll,jacobian2D_xmax,normal_xmax,&
+ NGLLY,NGLLZ,NSPEC2DMAX_XMIN_XMAX)
+ end if
endif
! on boundary: ymin
@@ -187,38 +226,53 @@
ispecb3=ispecb3+1
ibelm_ymin(ispecb3)=ispec
-! specify the 9 nodes for the 2-D boundary element
- xelm(1)=xstore(1,1,1,ispec)
- yelm(1)=ystore(1,1,1,ispec)
- zelm(1)=zstore(1,1,1,ispec)
- xelm(2)=xstore(NGLLX,1,1,ispec)
- yelm(2)=ystore(NGLLX,1,1,ispec)
- zelm(2)=zstore(NGLLX,1,1,ispec)
- xelm(3)=xstore(NGLLX,1,NGLLZ,ispec)
- yelm(3)=ystore(NGLLX,1,NGLLZ,ispec)
- zelm(3)=zstore(NGLLX,1,NGLLZ,ispec)
- xelm(4)=xstore(1,1,NGLLZ,ispec)
- yelm(4)=ystore(1,1,NGLLZ,ispec)
- zelm(4)=zstore(1,1,NGLLZ,ispec)
- xelm(5)=xstore((NGLLX+1)/2,1,1,ispec)
- yelm(5)=ystore((NGLLX+1)/2,1,1,ispec)
- zelm(5)=zstore((NGLLX+1)/2,1,1,ispec)
- xelm(6)=xstore(NGLLX,1,(NGLLZ+1)/2,ispec)
- yelm(6)=ystore(NGLLX,1,(NGLLZ+1)/2,ispec)
- zelm(6)=zstore(NGLLX,1,(NGLLZ+1)/2,ispec)
- xelm(7)=xstore((NGLLX+1)/2,1,NGLLZ,ispec)
- yelm(7)=ystore((NGLLX+1)/2,1,NGLLZ,ispec)
- zelm(7)=zstore((NGLLX+1)/2,1,NGLLZ,ispec)
- xelm(8)=xstore(1,1,(NGLLZ+1)/2,ispec)
- yelm(8)=ystore(1,1,(NGLLZ+1)/2,ispec)
- zelm(8)=zstore(1,1,(NGLLZ+1)/2,ispec)
- xelm(9)=xstore((NGLLX+1)/2,1,(NGLLZ+1)/2,ispec)
- yelm(9)=ystore((NGLLX+1)/2,1,(NGLLZ+1)/2,ispec)
- zelm(9)=zstore((NGLLX+1)/2,1,(NGLLZ+1)/2,ispec)
+ if ( .not. USE_GLL) then
+ ! specify the 9 nodes for the 2-D boundary element
+ xelm(1)=xstore(1,1,1,ispec)
+ yelm(1)=ystore(1,1,1,ispec)
+ zelm(1)=zstore(1,1,1,ispec)
+ xelm(2)=xstore(NGLLX,1,1,ispec)
+ yelm(2)=ystore(NGLLX,1,1,ispec)
+ zelm(2)=zstore(NGLLX,1,1,ispec)
+ xelm(3)=xstore(NGLLX,1,NGLLZ,ispec)
+ yelm(3)=ystore(NGLLX,1,NGLLZ,ispec)
+ zelm(3)=zstore(NGLLX,1,NGLLZ,ispec)
+ xelm(4)=xstore(1,1,NGLLZ,ispec)
+ yelm(4)=ystore(1,1,NGLLZ,ispec)
+ zelm(4)=zstore(1,1,NGLLZ,ispec)
+ xelm(5)=xstore((NGLLX+1)/2,1,1,ispec)
+ yelm(5)=ystore((NGLLX+1)/2,1,1,ispec)
+ zelm(5)=zstore((NGLLX+1)/2,1,1,ispec)
+ xelm(6)=xstore(NGLLX,1,(NGLLZ+1)/2,ispec)
+ yelm(6)=ystore(NGLLX,1,(NGLLZ+1)/2,ispec)
+ zelm(6)=zstore(NGLLX,1,(NGLLZ+1)/2,ispec)
+ xelm(7)=xstore((NGLLX+1)/2,1,NGLLZ,ispec)
+ yelm(7)=ystore((NGLLX+1)/2,1,NGLLZ,ispec)
+ zelm(7)=zstore((NGLLX+1)/2,1,NGLLZ,ispec)
+ xelm(8)=xstore(1,1,(NGLLZ+1)/2,ispec)
+ yelm(8)=ystore(1,1,(NGLLZ+1)/2,ispec)
+ zelm(8)=zstore(1,1,(NGLLZ+1)/2,ispec)
+ xelm(9)=xstore((NGLLX+1)/2,1,(NGLLZ+1)/2,ispec)
+ yelm(9)=ystore((NGLLX+1)/2,1,(NGLLZ+1)/2,ispec)
+ zelm(9)=zstore((NGLLX+1)/2,1,(NGLLZ+1)/2,ispec)
- call compute_jacobian_2D(myrank,ispecb3,xelm,yelm,zelm,dershape2D_y, &
+ call compute_jacobian_2D(myrank,ispecb3,xelm,yelm,zelm,dershape2D_y, &
jacobian2D_ymin,normal_ymin,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
+ else
+ ! get 25 GLL points for ymin
+ do k =1 ,NGLLZ
+ do i = 1,NGLLX
+ xelm2D(i,k) = xstore(i,1,k,ispec)
+ yelm2D(i,k) = ystore(i,1,k,ispec)
+ zelm2D(i,k) = zstore(i,1,k,ispec)
+ end do
+ end do
+ ! recalcualte 2D jacobian according to GLL points
+ call recalc_jacobian_gll2D(myrank,ispecb3,xelm2D,yelm2D,zelm2D,&
+ xigll,zigll,jacobian2D_ymin,normal_ymin,&
+ NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
+ end if
endif
! on boundary: ymax
@@ -228,38 +282,53 @@
ispecb4=ispecb4+1
ibelm_ymax(ispecb4)=ispec
-! specify the 9 nodes for the 2-D boundary element
- xelm(1)=xstore(1,NGLLY,1,ispec)
- yelm(1)=ystore(1,NGLLY,1,ispec)
- zelm(1)=zstore(1,NGLLY,1,ispec)
- xelm(2)=xstore(NGLLX,NGLLY,1,ispec)
- yelm(2)=ystore(NGLLX,NGLLY,1,ispec)
- zelm(2)=zstore(NGLLX,NGLLY,1,ispec)
- xelm(3)=xstore(NGLLX,NGLLY,NGLLZ,ispec)
- yelm(3)=ystore(NGLLX,NGLLY,NGLLZ,ispec)
- zelm(3)=zstore(NGLLX,NGLLY,NGLLZ,ispec)
- xelm(4)=xstore(1,NGLLY,NGLLZ,ispec)
- yelm(4)=ystore(1,NGLLY,NGLLZ,ispec)
- zelm(4)=zstore(1,NGLLY,NGLLZ,ispec)
- xelm(5)=xstore((NGLLX+1)/2,NGLLY,1,ispec)
- yelm(5)=ystore((NGLLX+1)/2,NGLLY,1,ispec)
- zelm(5)=zstore((NGLLX+1)/2,NGLLY,1,ispec)
- xelm(6)=xstore(NGLLX,NGLLY,(NGLLZ+1)/2,ispec)
- yelm(6)=ystore(NGLLX,NGLLY,(NGLLZ+1)/2,ispec)
- zelm(6)=zstore(NGLLX,NGLLY,(NGLLZ+1)/2,ispec)
- xelm(7)=xstore((NGLLX+1)/2,NGLLY,NGLLZ,ispec)
- yelm(7)=ystore((NGLLX+1)/2,NGLLY,NGLLZ,ispec)
- zelm(7)=zstore((NGLLX+1)/2,NGLLY,NGLLZ,ispec)
- xelm(8)=xstore(1,NGLLY,(NGLLZ+1)/2,ispec)
- yelm(8)=ystore(1,NGLLY,(NGLLZ+1)/2,ispec)
- zelm(8)=zstore(1,NGLLY,(NGLLZ+1)/2,ispec)
- xelm(9)=xstore((NGLLX+1)/2,NGLLY,(NGLLZ+1)/2,ispec)
- yelm(9)=ystore((NGLLX+1)/2,NGLLY,(NGLLZ+1)/2,ispec)
- zelm(9)=zstore((NGLLX+1)/2,NGLLY,(NGLLZ+1)/2,ispec)
+ if ( .not. USE_GLL) then
+ ! specify the 9 nodes for the 2-D boundary element
+ xelm(1)=xstore(1,NGLLY,1,ispec)
+ yelm(1)=ystore(1,NGLLY,1,ispec)
+ zelm(1)=zstore(1,NGLLY,1,ispec)
+ xelm(2)=xstore(NGLLX,NGLLY,1,ispec)
+ yelm(2)=ystore(NGLLX,NGLLY,1,ispec)
+ zelm(2)=zstore(NGLLX,NGLLY,1,ispec)
+ xelm(3)=xstore(NGLLX,NGLLY,NGLLZ,ispec)
+ yelm(3)=ystore(NGLLX,NGLLY,NGLLZ,ispec)
+ zelm(3)=zstore(NGLLX,NGLLY,NGLLZ,ispec)
+ xelm(4)=xstore(1,NGLLY,NGLLZ,ispec)
+ yelm(4)=ystore(1,NGLLY,NGLLZ,ispec)
+ zelm(4)=zstore(1,NGLLY,NGLLZ,ispec)
+ xelm(5)=xstore((NGLLX+1)/2,NGLLY,1,ispec)
+ yelm(5)=ystore((NGLLX+1)/2,NGLLY,1,ispec)
+ zelm(5)=zstore((NGLLX+1)/2,NGLLY,1,ispec)
+ xelm(6)=xstore(NGLLX,NGLLY,(NGLLZ+1)/2,ispec)
+ yelm(6)=ystore(NGLLX,NGLLY,(NGLLZ+1)/2,ispec)
+ zelm(6)=zstore(NGLLX,NGLLY,(NGLLZ+1)/2,ispec)
+ xelm(7)=xstore((NGLLX+1)/2,NGLLY,NGLLZ,ispec)
+ yelm(7)=ystore((NGLLX+1)/2,NGLLY,NGLLZ,ispec)
+ zelm(7)=zstore((NGLLX+1)/2,NGLLY,NGLLZ,ispec)
+ xelm(8)=xstore(1,NGLLY,(NGLLZ+1)/2,ispec)
+ yelm(8)=ystore(1,NGLLY,(NGLLZ+1)/2,ispec)
+ zelm(8)=zstore(1,NGLLY,(NGLLZ+1)/2,ispec)
+ xelm(9)=xstore((NGLLX+1)/2,NGLLY,(NGLLZ+1)/2,ispec)
+ yelm(9)=ystore((NGLLX+1)/2,NGLLY,(NGLLZ+1)/2,ispec)
+ zelm(9)=zstore((NGLLX+1)/2,NGLLY,(NGLLZ+1)/2,ispec)
- call compute_jacobian_2D(myrank,ispecb4,xelm,yelm,zelm,dershape2D_y, &
+ call compute_jacobian_2D(myrank,ispecb4,xelm,yelm,zelm,dershape2D_y, &
jacobian2D_ymax,normal_ymax,NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
+ else
+ ! get 25 GLL points for ymax
+ do k =1,NGLLZ
+ do i = 1,NGLLX
+ xelm2D(i,k) = xstore(i,NGLLY,k,ispec)
+ yelm2D(i,k) = ystore(i,NGLLY,k,ispec)
+ zelm2D(i,k) = zstore(i,NGLLY,k,ispec)
+ end do
+ end do
+ ! recalculate jacobian for 2D GLL points
+ call recalc_jacobian_gll2D(myrank,ispecb4,xelm2D,yelm2D,zelm2D,&
+ xigll,zigll,jacobian2D_ymax,normal_ymax,&
+ NGLLX,NGLLZ,NSPEC2DMAX_YMIN_YMAX)
+ end if
endif
! on boundary: bottom
@@ -269,36 +338,52 @@
ispecb5=ispecb5+1
ibelm_bottom(ispecb5)=ispec
- xelm(1)=xstore(1,1,1,ispec)
- yelm(1)=ystore(1,1,1,ispec)
- zelm(1)=zstore(1,1,1,ispec)
- xelm(2)=xstore(NGLLX,1,1,ispec)
- yelm(2)=ystore(NGLLX,1,1,ispec)
- zelm(2)=zstore(NGLLX,1,1,ispec)
- xelm(3)=xstore(NGLLX,NGLLY,1,ispec)
- yelm(3)=ystore(NGLLX,NGLLY,1,ispec)
- zelm(3)=zstore(NGLLX,NGLLY,1,ispec)
- xelm(4)=xstore(1,NGLLY,1,ispec)
- yelm(4)=ystore(1,NGLLY,1,ispec)
- zelm(4)=zstore(1,NGLLY,1,ispec)
- xelm(5)=xstore((NGLLX+1)/2,1,1,ispec)
- yelm(5)=ystore((NGLLX+1)/2,1,1,ispec)
- zelm(5)=zstore((NGLLX+1)/2,1,1,ispec)
- xelm(6)=xstore(NGLLX,(NGLLY+1)/2,1,ispec)
- yelm(6)=ystore(NGLLX,(NGLLY+1)/2,1,ispec)
- zelm(6)=zstore(NGLLX,(NGLLY+1)/2,1,ispec)
- xelm(7)=xstore((NGLLX+1)/2,NGLLY,1,ispec)
- yelm(7)=ystore((NGLLX+1)/2,NGLLY,1,ispec)
- zelm(7)=zstore((NGLLX+1)/2,NGLLY,1,ispec)
- xelm(8)=xstore(1,(NGLLY+1)/2,1,ispec)
- yelm(8)=ystore(1,(NGLLY+1)/2,1,ispec)
- zelm(8)=zstore(1,(NGLLY+1)/2,1,ispec)
- xelm(9)=xstore((NGLLX+1)/2,(NGLLY+1)/2,1,ispec)
- yelm(9)=ystore((NGLLX+1)/2,(NGLLY+1)/2,1,ispec)
- zelm(9)=zstore((NGLLX+1)/2,(NGLLY+1)/2,1,ispec)
+ if ( .not. USE_GLL) then
+ xelm(1)=xstore(1,1,1,ispec)
+ yelm(1)=ystore(1,1,1,ispec)
+ zelm(1)=zstore(1,1,1,ispec)
+ xelm(2)=xstore(NGLLX,1,1,ispec)
+ yelm(2)=ystore(NGLLX,1,1,ispec)
+ zelm(2)=zstore(NGLLX,1,1,ispec)
+ xelm(3)=xstore(NGLLX,NGLLY,1,ispec)
+ yelm(3)=ystore(NGLLX,NGLLY,1,ispec)
+ zelm(3)=zstore(NGLLX,NGLLY,1,ispec)
+ xelm(4)=xstore(1,NGLLY,1,ispec)
+ yelm(4)=ystore(1,NGLLY,1,ispec)
+ zelm(4)=zstore(1,NGLLY,1,ispec)
+ xelm(5)=xstore((NGLLX+1)/2,1,1,ispec)
+ yelm(5)=ystore((NGLLX+1)/2,1,1,ispec)
+ zelm(5)=zstore((NGLLX+1)/2,1,1,ispec)
+ xelm(6)=xstore(NGLLX,(NGLLY+1)/2,1,ispec)
+ yelm(6)=ystore(NGLLX,(NGLLY+1)/2,1,ispec)
+ zelm(6)=zstore(NGLLX,(NGLLY+1)/2,1,ispec)
+ xelm(7)=xstore((NGLLX+1)/2,NGLLY,1,ispec)
+ yelm(7)=ystore((NGLLX+1)/2,NGLLY,1,ispec)
+ zelm(7)=zstore((NGLLX+1)/2,NGLLY,1,ispec)
+ xelm(8)=xstore(1,(NGLLY+1)/2,1,ispec)
+ yelm(8)=ystore(1,(NGLLY+1)/2,1,ispec)
+ zelm(8)=zstore(1,(NGLLY+1)/2,1,ispec)
+ xelm(9)=xstore((NGLLX+1)/2,(NGLLY+1)/2,1,ispec)
+ yelm(9)=ystore((NGLLX+1)/2,(NGLLY+1)/2,1,ispec)
+ zelm(9)=zstore((NGLLX+1)/2,(NGLLY+1)/2,1,ispec)
- call compute_jacobian_2D(myrank,ispecb5,xelm,yelm,zelm,dershape2D_bottom, &
+ call compute_jacobian_2D(myrank,ispecb5,xelm,yelm,zelm,dershape2D_bottom, &
jacobian2D_bottom,normal_bottom,NGLLX,NGLLY,NSPEC2D_BOTTOM)
+
+ else
+ ! get 25 GLL points for zmin
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+ xelm2D(i,j) = xstore(i,j,1,ispec)
+ yelm2D(i,j) = ystore(i,j,1,ispec)
+ zelm2D(i,j) = zstore(i,j,1,ispec)
+ end do
+ end do
+ ! recalcuate 2D jacobian according to GLL points
+ call recalc_jacobian_gll2D(myrank,ispecb5,xelm2D,yelm2D,zelm2D,&
+ xigll,yigll,jacobian2D_bottom,normal_bottom,&
+ NGLLX,NGLLY,NSPEC2D_BOTTOM)
+ end if
endif
@@ -309,41 +394,59 @@
ispecb6=ispecb6+1
ibelm_top(ispecb6)=ispec
- xelm(1)=xstore(1,1,NGLLZ,ispec)
- yelm(1)=ystore(1,1,NGLLZ,ispec)
- zelm(1)=zstore(1,1,NGLLZ,ispec)
- xelm(2)=xstore(NGLLX,1,NGLLZ,ispec)
- yelm(2)=ystore(NGLLX,1,NGLLZ,ispec)
- zelm(2)=zstore(NGLLX,1,NGLLZ,ispec)
- xelm(3)=xstore(NGLLX,NGLLY,NGLLZ,ispec)
- yelm(3)=ystore(NGLLX,NGLLY,NGLLZ,ispec)
- zelm(3)=zstore(NGLLX,NGLLY,NGLLZ,ispec)
- xelm(4)=xstore(1,NGLLY,NGLLZ,ispec)
- yelm(4)=ystore(1,NGLLY,NGLLZ,ispec)
- zelm(4)=zstore(1,NGLLY,NGLLZ,ispec)
- xelm(5)=xstore((NGLLX+1)/2,1,NGLLZ,ispec)
- yelm(5)=ystore((NGLLX+1)/2,1,NGLLZ,ispec)
- zelm(5)=zstore((NGLLX+1)/2,1,NGLLZ,ispec)
- xelm(6)=xstore(NGLLX,(NGLLY+1)/2,NGLLZ,ispec)
- yelm(6)=ystore(NGLLX,(NGLLY+1)/2,NGLLZ,ispec)
- zelm(6)=zstore(NGLLX,(NGLLY+1)/2,NGLLZ,ispec)
- xelm(7)=xstore((NGLLX+1)/2,NGLLY,NGLLZ,ispec)
- yelm(7)=ystore((NGLLX+1)/2,NGLLY,NGLLZ,ispec)
- zelm(7)=zstore((NGLLX+1)/2,NGLLY,NGLLZ,ispec)
- xelm(8)=xstore(1,(NGLLY+1)/2,NGLLZ,ispec)
- yelm(8)=ystore(1,(NGLLY+1)/2,NGLLZ,ispec)
- zelm(8)=zstore(1,(NGLLY+1)/2,NGLLZ,ispec)
- xelm(9)=xstore((NGLLX+1)/2,(NGLLY+1)/2,NGLLZ,ispec)
- yelm(9)=ystore((NGLLX+1)/2,(NGLLY+1)/2,NGLLZ,ispec)
- zelm(9)=zstore((NGLLX+1)/2,(NGLLY+1)/2,NGLLZ,ispec)
+ if ( .not. USE_GLL) then
+ xelm(1)=xstore(1,1,NGLLZ,ispec)
+ yelm(1)=ystore(1,1,NGLLZ,ispec)
+ zelm(1)=zstore(1,1,NGLLZ,ispec)
+ xelm(2)=xstore(NGLLX,1,NGLLZ,ispec)
+ yelm(2)=ystore(NGLLX,1,NGLLZ,ispec)
+ zelm(2)=zstore(NGLLX,1,NGLLZ,ispec)
+ xelm(3)=xstore(NGLLX,NGLLY,NGLLZ,ispec)
+ yelm(3)=ystore(NGLLX,NGLLY,NGLLZ,ispec)
+ zelm(3)=zstore(NGLLX,NGLLY,NGLLZ,ispec)
+ xelm(4)=xstore(1,NGLLY,NGLLZ,ispec)
+ yelm(4)=ystore(1,NGLLY,NGLLZ,ispec)
+ zelm(4)=zstore(1,NGLLY,NGLLZ,ispec)
+ xelm(5)=xstore((NGLLX+1)/2,1,NGLLZ,ispec)
+ yelm(5)=ystore((NGLLX+1)/2,1,NGLLZ,ispec)
+ zelm(5)=zstore((NGLLX+1)/2,1,NGLLZ,ispec)
+ xelm(6)=xstore(NGLLX,(NGLLY+1)/2,NGLLZ,ispec)
+ yelm(6)=ystore(NGLLX,(NGLLY+1)/2,NGLLZ,ispec)
+ zelm(6)=zstore(NGLLX,(NGLLY+1)/2,NGLLZ,ispec)
+ xelm(7)=xstore((NGLLX+1)/2,NGLLY,NGLLZ,ispec)
+ yelm(7)=ystore((NGLLX+1)/2,NGLLY,NGLLZ,ispec)
+ zelm(7)=zstore((NGLLX+1)/2,NGLLY,NGLLZ,ispec)
+ xelm(8)=xstore(1,(NGLLY+1)/2,NGLLZ,ispec)
+ yelm(8)=ystore(1,(NGLLY+1)/2,NGLLZ,ispec)
+ zelm(8)=zstore(1,(NGLLY+1)/2,NGLLZ,ispec)
+ xelm(9)=xstore((NGLLX+1)/2,(NGLLY+1)/2,NGLLZ,ispec)
+ yelm(9)=ystore((NGLLX+1)/2,(NGLLY+1)/2,NGLLZ,ispec)
+ zelm(9)=zstore((NGLLX+1)/2,(NGLLY+1)/2,NGLLZ,ispec)
- call compute_jacobian_2D(myrank,ispecb6,xelm,yelm,zelm,dershape2D_top, &
+ call compute_jacobian_2D(myrank,ispecb6,xelm,yelm,zelm,dershape2D_top, &
jacobian2D_top,normal_top,NGLLX,NGLLY,NSPEC2D_TOP)
+ else
+ ! get 25 GLL points for zmax
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+ xelm2D(i,j) = xstore(i,j,NGLLZ,ispec)
+ yelm2D(i,j) = ystore(i,j,NGLLZ,ispec)
+ zelm2D(i,j) = zstore(i,j,NGLLZ,ispec)
+ end do
+ end do
+ ! recalcuate jacobian according to 2D gll points
+ call recalc_jacobian_gll2D(myrank,ispecb6,xelm2D,yelm2D,zelm2D,&
+ xigll,yigll,jacobian2D_top,normal_top,&
+ NGLLX,NGLLY,NSPEC2D_TOP)
+ end if
+
endif
enddo
+!< Hejun
+
! check theoretical value of elements at the bottom
if(ispecb5 /= NSPEC2D_BOTTOM) then
call exit_MPI(myrank,'ispecb5 should equal NSPEC2D_BOTTOM')
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/get_model.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/get_model.f90 2010-01-21 16:01:30 UTC (rev 16155)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/get_model.f90 2010-01-21 16:14:52 UTC (rev 16156)
@@ -42,7 +42,7 @@
numker,numhpa,numcof,ihpa,lmax,nylm, &
lmxhpa,itypehpa,ihpakern,numcoe,ivarkern, &
nconpt,iver,iconpt,conpt,xlaspl,xlospl,radspl, &
- coe,vercof,vercofd,ylmcof,wk1,wk2,wk3,kerstr,varstr)
+ coe,vercof,vercofd,ylmcof,wk1,wk2,wk3,kerstr,varstr,elem_in_crust)
implicit none
@@ -400,6 +400,10 @@
character(len=80) kerstr
character(len=40) varstr(maxker)
+!> Hejun
+ logical::elem_in_crust
+!< Hejun
+
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
@@ -749,21 +753,23 @@
c66 = c44
endif
endif
-
+!> Hejun
+! Assign Attenuation after get 3-D crustal model
! This is here to identify how and where to include 3D attenuation
- if(ATTENUATION .and. ATTENUATION_3D) then
- call xyz_2_rthetaphi_dble(xmesh,ymesh,zmesh,r_dummy,theta,phi)
- call reduce(theta,phi)
- theta_degrees = theta / DEGREES_TO_RADIANS
- phi_degrees = phi / DEGREES_TO_RADIANS
- tau_e(:) = 0.0d0
- ! Get the value of Qmu (Attenuation) dependedent on
- ! the radius (r_prem) and idoubling flag
- !call attenuation_model_1D_PREM(r_prem, Qmu, idoubling)
- call attenuation_model_3D_QRFSI12(r_prem*R_EARTH_KM,theta_degrees,phi_degrees,Qmu,QRFSI12_Q,idoubling)
- ! Get tau_e from tau_s and Qmu
- call attenuation_conversion(Qmu, T_c_source, tau_s, tau_e, AM_V, AM_S, AS_V)
- endif
+! if(ATTENUATION .and. ATTENUATION_3D) then
+! call xyz_2_rthetaphi_dble(xmesh,ymesh,zmesh,r_dummy,theta,phi)
+! call reduce(theta,phi)
+! theta_degrees = theta / DEGREES_TO_RADIANS
+! phi_degrees = phi / DEGREES_TO_RADIANS
+! tau_e(:) = 0.0d0
+! ! Get the value of Qmu (Attenuation) dependedent on
+! ! the radius (r_prem) and idoubling flag
+! !call attenuation_model_1D_PREM(r_prem, Qmu, idoubling)
+! call attenuation_model_3D_QRFSI12(r_prem*R_EARTH_KM,theta_degrees,phi_degrees,Qmu,QRFSI12_Q,idoubling)
+! ! Get tau_e from tau_s and Qmu
+! call attenuation_conversion(Qmu, T_c_source, tau_s, tau_e, AM_V, AM_S, AS_V)
+! endif
+!<Hejun
! get the 3-D crustal model
if(CRUSTAL) then
@@ -811,7 +817,7 @@
lat=(PI/2.0d0-theta)*180.0d0/PI
lon=phi*180.0d0/PI
if(lon>180.0d0) lon=lon-360.0d0
- call crustal_model(lat,lon,r,vpc,vsc,rhoc,moho,found_crust,CM_V)
+ call crustal_model(lat,lon,r,vpc,vsc,rhoc,moho,found_crust,CM_V,elem_in_crust)
if (found_crust) then
vpv=vpc
vph=vpc
@@ -848,7 +854,7 @@
lat=(PI/2.0d0-theta)*180.0d0/PI
lon=phi*180.0d0/PI
if(lon>180.0d0) lon=lon-360.0d0
- call crustal_model(lat,lon,r,vpc,vsc,rhoc,moho,found_crust,CM_V)
+ call crustal_model(lat,lon,r,vpc,vsc,rhoc,moho,found_crust,CM_V,elem_in_crust)
if (found_crust) then
vpv=vpc
vph=vpc
@@ -884,6 +890,50 @@
endif
endif
+!> Hejun
+! New Attenuation assignment
+! Define 3D and 1D Attenuation after moho stretch
+! and before TOPOGRAPHY/ELLIPCITY
+ if (ATTENUATION) then
+ call xyz_2_rthetaphi_dble(xmesh,ymesh,zmesh,r_dummy,theta,phi)
+ call reduce(theta,phi)
+ lat=(PI/2.0d0-theta)*180.0d0/PI
+ lon=phi*180.0d0/PI
+ if(lon>180.0d0) lon=lon-360.0d0
+
+ theta_degrees = theta / DEGREES_TO_RADIANS
+ phi_degrees = phi / DEGREES_TO_RADIANS
+ tau_e(:) = 0.0d0
+ ! Get the value of Qmu (Attenuation) dependedent on
+ ! the radius (r_prem) and idoubling flag
+
+ if (ATTENUATION_3D) then
+ call attenuation_model_3D_QRFSI12(r_prem*R_EARTH_KM,theta_degrees,&
+ phi_degrees,Qmu,QRFSI12_Q,idoubling)
+ else
+ if (REFERENCE_1D_MODEL == REFERENCE_MODEL_PREM) then
+ call attenuation_model_1D_PREM(r_prem, Qmu)
+ else if (REFERENCE_1D_MODEL == REFERENCE_MODEL_REF) then
+ call attenuation_model_1D_REF(r_prem, Qmu)
+
+ if ( CRUSTAL ) then
+ ! extand the R80 to surface
+ if (r_prem > R80/R_EARTH) then
+ Qmu = 191.0d0
+ Qkappa = 943.0d0
+ end if
+ if ( r_prem > (ONE-moho) .or. elem_in_crust) then
+ Qmu=300.0d0
+ Qkappa=57822.5d0
+ end if
+ end if
+ end if
+ end if
+ ! Get tau_e from tau_s and Qmu
+ call attenuation_conversion(Qmu, T_c_source, tau_s, tau_e, AM_V, AM_S, AS_V)
+ endif
+!< Hejun
+
! define elastic parameters in the model
! distinguish between single and double precision for reals
@@ -1007,11 +1057,21 @@
endif
- if(ATTENUATION .and. ATTENUATION_3D) then
+ !> Hejun
+ ! No matter 1D or 3D attenuation, we save all gll point value
+ !if(ATTENUATION .and. ATTENUATION_3D) then
+ ! tau_e_store(:,i,j,k,ispec) = tau_e(:)
+ ! Qmu_store(i,j,k,ispec) = Qmu
+ !endif
+
+ if(ATTENUATION) then
tau_e_store(:,i,j,k,ispec) = tau_e(:)
Qmu_store(i,j,k,ispec) = Qmu
endif
+ !< Hejun
+
+
enddo
enddo
enddo
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D.f90 2010-01-21 16:01:30 UTC (rev 16155)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/meshfem3D.f90 2010-01-21 16:14:52 UTC (rev 16156)
@@ -1444,6 +1444,14 @@
call MPI_BCAST(AM_V%min_period, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
call MPI_BCAST(AM_V%max_period, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
+!> Hejun
+! New 1D Attenuation allocation
+ if(myrank /= 0) allocate(AM_V%Qtau_s(N_SLS))
+ call MPI_BCAST(AM_V%QT_c_source, 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
+ call MPI_BCAST(AM_V%Qtau_s(1), 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
+ call MPI_BCAST(AM_V%Qtau_s(2), 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
+ call MPI_BCAST(AM_V%Qtau_s(3), 1, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
+!< Hejun
call attenuation_model_setup(REFERENCE_1D_MODEL, RICB, RCMB, R670, R220, R80,AM_V,M1066a_V,Mak135_V,Mref_V,SEA1DM_V,AM_S,AS_V)
endif
@@ -1537,7 +1545,8 @@
rotation_matrix,ANGULAR_WIDTH_XI_RAD,ANGULAR_WIDTH_ETA_RAD, &
ATTENUATION,ATTENUATION_3D,SAVE_MESH_FILES, &
NCHUNKS,INCLUDE_CENTRAL_CUBE,ABSORBING_CONDITIONS,REFERENCE_1D_MODEL,THREE_D_MODEL, &
- R_CENTRAL_CUBE,RICB,RHO_OCEANS,RCMB,R670,RMOHO,RTOPDDOUBLEPRIME,R600,R220,R771,R400,R120,R80,RMIDDLE_CRUST,ROCEAN, &
+ R_CENTRAL_CUBE,RICB,RHO_OCEANS,RCMB,R670,RMOHO,RMOHO_FICTITIOUS_IN_MESHER,&
+ RTOPDDOUBLEPRIME,R600,R220,R771,R400,R120,R80,RMIDDLE_CRUST,ROCEAN, &
ner,ratio_sampling_array,doubling_index,r_bottom, r_top,this_region_has_a_doubling,CASE_3D, &
AMM_V,AM_V,QRFSI12_Q,M1066a_V,Mak135_V,Mref_V,SEA1DM_V,D3MM_V,HMM,JP3DM_V,SEA99M_V,CM_V,AM_S,AS_V, &
numker,numhpa,numcof,ihpa,lmax,nylm, &
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/moho_stretching.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/moho_stretching.f90 2010-01-21 16:01:30 UTC (rev 16155)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/moho_stretching.f90 2010-01-21 16:14:52 UTC (rev 16156)
@@ -299,3 +299,115 @@
end subroutine read_smooth_moho
+!> Hejun
+! stretching the moho according to the crust 2.0
+! input: myrank, xelm, yelm, zelm, RMOHO_FICTITIOUS_IN_MESHER R220,RMIDDLE_CRUST, CM_V
+! Dec, 30, 2009
+ subroutine hornor_moho_crust2(myrank,xelm,yelm,zelm,RMOHO_FICTITIOUS_IN_MESHER,R220,RMIDDLE_CRUST,CM_V,elem_in_crust)
+
+ implicit none
+
+ include "constants.h"
+
+ double precision xelm(NGNOD)
+ double precision yelm(NGNOD)
+ double precision zelm(NGNOD)
+ double precision R220,RMIDDLE_CRUST
+ integer::RMOHO_FICTITIOUS_IN_MESHER
+ integer::myrank
+ logical:: elem_in_crust
+
+! crustal_model_variables
+ type crustal_model_variables
+ sequence
+ double precision, dimension(NKEYS_CRUST,NLAYERS_CRUST) :: thlr
+ double precision, dimension(NKEYS_CRUST,NLAYERS_CRUST) :: velocp
+ double precision, dimension(NKEYS_CRUST,NLAYERS_CRUST) :: velocs
+ double precision, dimension(NKEYS_CRUST,NLAYERS_CRUST) :: dens
+ character(len=2) abbreviation(NCAP_CRUST/2,NCAP_CRUST)
+ character(len=2) code(NKEYS_CRUST)
+ end type crustal_model_variables
+
+ type(crustal_model_variables) CM_V
+
+
+ ! local parameters
+ integer:: ia,count
+ double precision:: r,theta,phi,lat,lon
+ double precision:: vpc,vsc,rhoc,moho,elevation,gamma
+ logical:: found_crust
+
+
+ count = 0
+ do ia= 1,NGNOD
+ call xyz_2_rthetaphi_dble(xelm(ia),yelm(ia),zelm(ia),r,theta,phi)
+ call reduce(theta,phi)
+
+ lat=(PI/2.0d0-theta)*180.0d0/PI
+ lon = phi*180.0d0/PI
+ if (lon>180.d0) lon = lon-360.0d0
+
+ call crustal_model(lat,lon,r,vpc,vsc,rhoc,moho,found_crust,CM_V,elem_in_crust)
+
+ moho = ONE-moho
+
+ elevation = 0.d0
+
+ if (moho<dble(RMOHO_FICTITIOUS_IN_MESHER)/R_EARTH) then
+
+ elevation = (moho-dble(RMOHO_FICTITIOUS_IN_MESHER)/R_EARTH)
+
+ if ( r >= dble(RMOHO_FICTITIOUS_IN_MESHER)/R_EARTH) then
+ gamma = (( R_UNIT_SPHERE -r )/(R_UNIT_SPHERE-dble(RMOHO_FICTITIOUS_IN_MESHER)/R_EARTH))
+ else if ( r < dble(RMOHO_FICTITIOUS_IN_MESHER)/R_EARTH) then
+ gamma = (( r-R220/R_EARTH)/(dble(RMOHO_FICTITIOUS_IN_MESHER)/R_EARTH-R220/R_EARTH))
+ ! since not all GLL points are exactlly at R220, use a small
+ ! tolerance for R220 detection, fix R220
+ if (abs(gamma) < SMALLVAL) then
+ gamma = 0.0
+ end if
+ end if
+ if(gamma < -0.0001 .or. gamma > 1.0001) call exit_MPI(myrank,'incorrect value of gamma for moho from crust 2.0')
+
+ xelm(ia)=xelm(ia)*(ONE+gamma*elevation/r)
+ yelm(ia)=yelm(ia)*(ONE+gamma*elevation/r)
+ zelm(ia)=zelm(ia)*(ONE+gamma*elevation/r)
+
+
+ else if ( moho> RMIDDLE_CRUST/R_EARTH) then
+ elevation = moho - RMIDDLE_CRUST/R_EARTH
+
+ if ( r > RMIDDLE_CRUST/R_EARTH) then
+ gamma= (R_UNIT_SPHERE-r)/(R_UNIT_SPHERE-RMIDDLE_CRUST/R_EARTH)
+ else if ( r <= RMIDDLE_CRUST/R_EARTH) then
+ gamma = (r-R220/R_EARTH)/(RMIDDLE_CRUST/R_EARTH-R220/R_EARTH)
+ ! since not all GLL points are exactlly at R220, use a small
+ ! tolerance for R220 detection, fix R220
+ if (abs(gamma) < SMALLVAL) then
+ gamma = 0.0
+ end if
+ end if
+ if(gamma < -0.0001 .or. gamma > 1.0001) call exit_MPI(myrank,'incorrect value of gamma for moho from crust 2.0')
+
+ xelm(ia)=xelm(ia)*(ONE+gamma*elevation/r)
+ yelm(ia)=yelm(ia)*(ONE+gamma*elevation/r)
+ zelm(ia)=zelm(ia)*(ONE+gamma*elevation/r)
+
+ end if
+
+
+ ! recalculate radius to decide whether this element is in the crust
+ call xyz_2_rthetaphi_dble(xelm(ia),yelm(ia),zelm(ia),r,theta,phi)
+
+ if ( r>= 0.9999d0*moho ) then
+ count = count + 1
+ end if
+
+ end do
+
+ if ( count == NGNOD) then
+ elem_in_crust = .true.
+ end if
+
+ end subroutine hornor_moho_crust2
+!< Hejun
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/read_compute_parameters.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/read_compute_parameters.f90 2010-01-21 16:01:30 UTC (rev 16155)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/read_compute_parameters.f90 2010-01-21 16:14:52 UTC (rev 16156)
@@ -142,6 +142,10 @@
integer :: tmp_sum_nglob2D_xi, tmp_sum_nglob2D_eta,divider,nglob_edges_h,nglob_edge_v,to_remove
+!> Hejun
+ double precision:: R80_FICTITIOUS_IN_MESHER
+!< Hejun
+
! get the base pathname for output files
call get_value_string(OUTPUT_FILES, 'OUTPUT_FILES', 'OUTPUT_FILES')
@@ -1046,8 +1050,14 @@
! from the d220 to the Earth surface
if(HONOR_1D_SPHERICAL_MOHO) then
RMOHO_FICTITIOUS_IN_MESHER = RMOHO
+ !> Hejun
+ R80_FICTITIOUS_IN_MESHER=R80
+ !< Hejun
else
RMOHO_FICTITIOUS_IN_MESHER = (R80 + R_EARTH) / 2
+ !> Hejun
+ R80_FICTITIOUS_IN_MESHER=R80-40000.d0
+ !< Hejun
endif
call read_value_double_precision(RECORD_LENGTH_IN_MINUTES, 'solver.RECORD_LENGTH_IN_MINUTES')
@@ -1340,15 +1350,15 @@
! the last region is in the inner core near the center of the Earth
r_top(1) = R_EARTH
- r_bottom(1) = R80
+ r_bottom(1) = R80_FICTITIOUS_IN_MESHER
r_top(2) = RMIDDLE_CRUST !!!! now fictitious
r_bottom(2) = RMOHO_FICTITIOUS_IN_MESHER !!!! now fictitious
r_top(3) = RMOHO_FICTITIOUS_IN_MESHER !!!! now fictitious
- r_bottom(3) = R80 !!!! now fictitious
+ r_bottom(3) = R80_FICTITIOUS_IN_MESHER !!!! now fictitious
- r_top(4) = R80
+ r_top(4) = R80_FICTITIOUS_IN_MESHER
r_bottom(4) = R220
r_top(5) = R220
@@ -1383,15 +1393,15 @@
! new definition of rmins & rmaxs
rmaxs(1) = ONE
- rmins(1) = R80 / R_EARTH
+ rmins(1) = R80_FICTITIOUS_IN_MESHER / R_EARTH
rmaxs(2) = RMIDDLE_CRUST / R_EARTH !!!! now fictitious
rmins(2) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH !!!! now fictitious
rmaxs(3) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH !!!! now fictitious
- rmins(3) = R80 / R_EARTH !!!! now fictitious
+ rmins(3) = R80_FICTITIOUS_IN_MESHER / R_EARTH !!!! now fictitious
- rmaxs(4) = R80 / R_EARTH
+ rmaxs(4) = R80_FICTITIOUS_IN_MESHER / R_EARTH
rmins(4) = R220 / R_EARTH
rmaxs(5) = R220 / R_EARTH
@@ -1474,9 +1484,9 @@
r_bottom(1) = RMOHO_FICTITIOUS_IN_MESHER
r_top(2) = RMOHO_FICTITIOUS_IN_MESHER
- r_bottom(2) = R80
+ r_bottom(2) = R80_FICTITIOUS_IN_MESHER
- r_top(3) = R80
+ r_top(3) = R80_FICTITIOUS_IN_MESHER
r_bottom(3) = R220
r_top(4) = R220
@@ -1514,9 +1524,9 @@
rmins(1) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
rmaxs(2) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
- rmins(2) = R80 / R_EARTH
+ rmins(2) = R80_FICTITIOUS_IN_MESHER / R_EARTH
- rmaxs(3) = R80 / R_EARTH
+ rmaxs(3) = R80_FICTITIOUS_IN_MESHER / R_EARTH
rmins(3) = R220 / R_EARTH
rmaxs(4) = R220 / R_EARTH
@@ -1601,9 +1611,9 @@
r_bottom(2) = RMOHO_FICTITIOUS_IN_MESHER
r_top(3) = RMOHO_FICTITIOUS_IN_MESHER
- r_bottom(3) = R80
+ r_bottom(3) = R80_FICTITIOUS_IN_MESHER
- r_top(4) = R80
+ r_top(4) = R80_FICTITIOUS_IN_MESHER
r_bottom(4) = R220
r_top(5) = R220
@@ -1644,9 +1654,9 @@
rmins(2) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
rmaxs(3) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
- rmins(3) = R80 / R_EARTH
+ rmins(3) = R80_FICTITIOUS_IN_MESHER / R_EARTH
- rmaxs(4) = R80 / R_EARTH
+ rmaxs(4) = R80_FICTITIOUS_IN_MESHER / R_EARTH
rmins(4) = R220 / R_EARTH
rmaxs(5) = R220 / R_EARTH
@@ -1729,15 +1739,15 @@
! the last region is in the inner core near the center of the Earth
r_top(1) = R_EARTH
- r_bottom(1) = R80
+ r_bottom(1) = R80_FICTITIOUS_IN_MESHER
r_top(2) = RMIDDLE_CRUST !!!! now fictitious
r_bottom(2) = RMOHO_FICTITIOUS_IN_MESHER !!!! now fictitious
r_top(3) = RMOHO_FICTITIOUS_IN_MESHER !!!! now fictitious
- r_bottom(3) = R80 !!!! now fictitious
+ r_bottom(3) = R80_FICTITIOUS_IN_MESHER !!!! now fictitious
- r_top(4) = R80
+ r_top(4) = R80_FICTITIOUS_IN_MESHER
r_bottom(4) = R220
r_top(5) = R220
@@ -1775,15 +1785,15 @@
! new definition of rmins & rmaxs
rmaxs(1) = ONE
- rmins(1) = R80 / R_EARTH
+ rmins(1) = R80_FICTITIOUS_IN_MESHER / R_EARTH
rmaxs(2) = RMIDDLE_CRUST / R_EARTH !!!! now fictitious
rmins(2) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH !!!! now fictitious
rmaxs(3) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH !!!! now fictitious
- rmins(3) = R80 / R_EARTH !!!! now fictitious
+ rmins(3) = R80_FICTITIOUS_IN_MESHER / R_EARTH !!!! now fictitious
- rmaxs(4) = R80 / R_EARTH
+ rmaxs(4) = R80_FICTITIOUS_IN_MESHER / R_EARTH
rmins(4) = R220 / R_EARTH
rmaxs(5) = R220 / R_EARTH
@@ -1869,9 +1879,9 @@
r_bottom(1) = RMOHO_FICTITIOUS_IN_MESHER
r_top(2) = RMOHO_FICTITIOUS_IN_MESHER
- r_bottom(2) = R80
+ r_bottom(2) = R80_FICTITIOUS_IN_MESHER
- r_top(3) = R80
+ r_top(3) = R80_FICTITIOUS_IN_MESHER
r_bottom(3) = R220
r_top(4) = R220
@@ -1912,9 +1922,9 @@
rmins(1) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
rmaxs(2) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
- rmins(2) = R80 / R_EARTH
+ rmins(2) = R80_FICTITIOUS_IN_MESHER / R_EARTH
- rmaxs(3) = R80 / R_EARTH
+ rmaxs(3) = R80_FICTITIOUS_IN_MESHER / R_EARTH
rmins(3) = R220 / R_EARTH
rmaxs(4) = R220 / R_EARTH
@@ -2001,9 +2011,9 @@
r_bottom(2) = RMOHO_FICTITIOUS_IN_MESHER
r_top(3) = RMOHO_FICTITIOUS_IN_MESHER
- r_bottom(3) = R80
+ r_bottom(3) = R80_FICTITIOUS_IN_MESHER
- r_top(4) = R80
+ r_top(4) = R80_FICTITIOUS_IN_MESHER
r_bottom(4) = R220
r_top(5) = R220
@@ -2047,9 +2057,9 @@
rmins(2) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
rmaxs(3) = RMOHO_FICTITIOUS_IN_MESHER / R_EARTH
- rmins(3) = R80 / R_EARTH
+ rmins(3) = R80_FICTITIOUS_IN_MESHER / R_EARTH
- rmaxs(4) = R80 / R_EARTH
+ rmaxs(4) = R80_FICTITIOUS_IN_MESHER / R_EARTH
rmins(4) = R220 / R_EARTH
rmaxs(5) = R220 / R_EARTH
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/save_arrays_solver.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/save_arrays_solver.f90 2010-01-21 16:01:30 UTC (rev 16155)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/save_arrays_solver.f90 2010-01-21 16:14:52 UTC (rev 16156)
@@ -46,8 +46,8 @@
NSPEC2DMAX_XMIN_XMAX,NSPEC2DMAX_YMIN_YMAX,NSPEC2D_BOTTOM,NSPEC2D_TOP, &
TRANSVERSE_ISOTROPY,HETEROGEN_3D_MANTLE,ANISOTROPIC_3D_MANTLE,ANISOTROPIC_INNER_CORE,OCEANS, &
tau_s,tau_e_store,Qmu_store,T_c_source, &
- ATTENUATION,ATTENUATION_3D,vx,vy,vz,vnspec, &
- NEX_PER_PROC_XI,NEX_PER_PROC_ETA,NEX_XI,ichunk,NCHUNKS,ABSORBING_CONDITIONS,AM_V)
+ ATTENUATION,vx,vy,vz,vnspec, &
+ NEX_PER_PROC_XI,NEX_PER_PROC_ETA,NEX_XI,ichunk,NCHUNKS,ABSORBING_CONDITIONS)
implicit none
@@ -72,10 +72,13 @@
integer :: Qn ! Number of points
end type attenuation_model_variables
- type (attenuation_model_variables) AM_V
+!> Hejun
+! type (attenuation_model_variables) AM_V
! attenuation_model_variables
+!< Hejun
- logical ATTENUATION,ATTENUATION_3D,ABSORBING_CONDITIONS
+ logical ATTENUATION,ABSORBING_CONDITIONS
+! logical ATTENUATION_3D
integer NEX_PER_PROC_XI,NEX_PER_PROC_ETA,NEX_XI,ichunk
@@ -364,16 +367,28 @@
close(27)
- if(ATTENUATION .and. ATTENUATION_3D) then
- open(unit=27, file=prname(1:len_trim(prname))//'attenuation3D.bin', status='unknown', form='unformatted',action='write')
+!> Hejun
+! No matter 1D or 3D Attenuation, we save value for gll points
+! if(ATTENUATION .and. ATTENUATION_3D) then
+! open(unit=27, file=prname(1:len_trim(prname))//'attenuation3D.bin', status='unknown', form='unformatted',action='write')
+! write(27) tau_s
+! write(27) tau_e_store
+! write(27) Qmu_store
+! write(27) T_c_source
+! close(27)
+! else if(ATTENUATION) then
+! call attenuation_save_arrays(prname, iregion_code, AM_V)
+! endif
+
+ if(ATTENUATION) then
+ open(unit=27, file=prname(1:len_trim(prname))//'attenuation.bin', status='unknown', form='unformatted',action='write')
write(27) tau_s
write(27) tau_e_store
write(27) Qmu_store
write(27) T_c_source
close(27)
- else if(ATTENUATION) then
- call attenuation_save_arrays(prname, iregion_code, AM_V)
endif
+!< Hejun
end subroutine save_arrays_solver
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/save_header_file.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/save_header_file.f90 2010-01-21 16:01:30 UTC (rev 16155)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/save_header_file.f90 2010-01-21 16:14:52 UTC (rev 16156)
@@ -415,20 +415,21 @@
write(IOUT,*) 'integer, parameter :: NUMMSGS_FACES_VAL = ',NPROC_XI*NUM_FACES*NUM_MSG_TYPES
write(IOUT,*) 'integer, parameter :: NCORNERSCHUNKS_VAL = ',NCORNERSCHUNKS
+!> Hejun
if(ATTENUATION) then
- if(ATTENUATION_3D) then
+! if(ATTENUATION_3D) then
att1 = NGLLX
att2 = NGLLY
att3 = NGLLZ
att4 = NSPEC(IREGION_CRUST_MANTLE)
att5 = NSPEC(IREGION_INNER_CORE)
- else
- att1 = 1
- att2 = 1
- att3 = 1
- att4 = NRAD_ATTENUATION
- att5 = NRAD_ATTENUATION
- endif
+! else
+! att1 = 1
+! att2 = 1
+! att3 = 1
+! att4 = NRAD_ATTENUATION
+! att5 = NRAD_ATTENUATION
+! endif
else
att1 = 1
att2 = 1
@@ -436,6 +437,7 @@
att4 = 1
att5 = 1
endif
+!< Hejun
write(IOUT,*) 'integer, parameter :: ATT1 = ',att1
write(IOUT,*) 'integer, parameter :: ATT2 = ',att2
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90 2010-01-21 16:01:30 UTC (rev 16155)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90 2010-01-21 16:14:52 UTC (rev 16156)
@@ -282,8 +282,10 @@
integer :: Qn ! Number of points
end type attenuation_model_variables
- type (attenuation_model_variables) AM_V
+!> Hejun
+! type (attenuation_model_variables) AM_V
! attenuation_model_variables
+!< Hejun
! memory variables and standard linear solids for attenuation
double precision, dimension(N_SLS) :: tau_sigma_dble
@@ -302,7 +304,9 @@
double precision, dimension(N_SLS,ATT1,ATT2,ATT3,ATT5) :: factor_common_inner_core_dble
double precision scale_factor,scale_factor_minus_one
- real(kind=CUSTOM_REAL) dist_cr
+!> Hejun
+! real(kind=CUSTOM_REAL) dist_cr
+!< Hejun
real(kind=CUSTOM_REAL), dimension(5,N_SLS,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ATTENUAT) :: R_memory_crust_mantle
real(kind=CUSTOM_REAL), dimension(5,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_STR_OR_ATT) :: epsilondev_crust_mantle
@@ -831,7 +835,9 @@
! if running on MareNostrum in Barcelona
character(len=400) system_command
- integer iregion_selected
+!> Hejun
+! integer iregion_selected
+!< Hejun
! computed in read_compute_parameters
integer, dimension(MAX_NUMBER_OF_MESH_LAYERS) :: ner,ratio_sampling_array
@@ -2376,27 +2382,41 @@
! get and store PREM attenuation model
+!> Hejun
+! Use the same reader for 1D and 3D attenuation model
! ATTENUATION_3D get values from mesher
- if(ATTENUATION_3D_VAL) then
- ! CRUST_MANTLE ATTENUATION
- call create_name_database(prname, myrank, IREGION_CRUST_MANTLE, LOCAL_PATH)
- call get_attenuation_model_3D(myrank, prname, omsb_crust_mantle_dble, &
+! if(ATTENUATION_3D_VAL) then
+! ! CRUST_MANTLE ATTENUATION
+! call create_name_database(prname, myrank, IREGION_CRUST_MANTLE, LOCAL_PATH)
+! call get_attenuation_model_3D(myrank, prname, omsb_crust_mantle_dble, &
+! factor_common_crust_mantle_dble, factor_scale_crust_mantle_dble, tau_sigma_dble, NSPEC_CRUST_MANTLE)
+! ! INNER_CORE ATTENUATION
+! call create_name_database(prname, myrank, IREGION_INNER_CORE, LOCAL_PATH)
+! call get_attenuation_model_3D(myrank, prname, omsb_inner_core_dble, &
+! factor_common_inner_core_dble, factor_scale_inner_core_dble, tau_sigma_dble, NSPEC_INNER_CORE)
+! else ! ATTENUATION = .true. .and. ATTENUATION_3D = .false.
+! call create_name_database(prname, myrank, IREGION_CRUST_MANTLE, LOCAL_PATH)
+! call get_attenuation_model_1D(myrank, prname, IREGION_CRUST_MANTLE, tau_sigma_dble, &
+! omsb_crust_mantle_dble, factor_common_crust_mantle_dble, &
+! factor_scale_crust_mantle_dble, NRAD_ATTENUATION,1,1,1, AM_V)
+! omsb_inner_core_dble(:,:,:,1:min(ATT4,ATT5)) = omsb_crust_mantle_dble(:,:,:,1:min(ATT4,ATT5))
+! factor_scale_inner_core_dble(:,:,:,1:min(ATT4,ATT5)) = factor_scale_crust_mantle_dble(:,:,:,1:min(ATT4,ATT5))
+! factor_common_inner_core_dble(:,:,:,:,1:min(ATT4,ATT5)) = factor_common_crust_mantle_dble(:,:,:,:,1:min(ATT4,ATT5))
+! ! Tell the Attenuation Code about the IDOUBLING regions within the Mesh
+! call set_attenuation_regions_1D(RICB, RCMB, R670, R220, R80, AM_V)
+! endif ! ATTENUATION_3D
+!< Hejun
+
+!> Hejun
+ ! CRUST_MANTLE ATTENUATION
+ call create_name_database(prname, myrank, IREGION_CRUST_MANTLE, LOCAL_PATH)
+ call get_attenuation_model_3D(myrank, prname, omsb_crust_mantle_dble, &
factor_common_crust_mantle_dble, factor_scale_crust_mantle_dble, tau_sigma_dble, NSPEC_CRUST_MANTLE)
- ! INNER_CORE ATTENUATION
- call create_name_database(prname, myrank, IREGION_INNER_CORE, LOCAL_PATH)
- call get_attenuation_model_3D(myrank, prname, omsb_inner_core_dble, &
+ ! INNER_CORE ATTENUATION
+ call create_name_database(prname, myrank, IREGION_INNER_CORE, LOCAL_PATH)
+ call get_attenuation_model_3D(myrank, prname, omsb_inner_core_dble, &
factor_common_inner_core_dble, factor_scale_inner_core_dble, tau_sigma_dble, NSPEC_INNER_CORE)
- else ! ATTENUATION = .true. .and. ATTENUATION_3D = .false.
- call create_name_database(prname, myrank, IREGION_CRUST_MANTLE, LOCAL_PATH)
- call get_attenuation_model_1D(myrank, prname, IREGION_CRUST_MANTLE, tau_sigma_dble, &
- omsb_crust_mantle_dble, factor_common_crust_mantle_dble, &
- factor_scale_crust_mantle_dble, NRAD_ATTENUATION,1,1,1, AM_V)
- omsb_inner_core_dble(:,:,:,1:min(ATT4,ATT5)) = omsb_crust_mantle_dble(:,:,:,1:min(ATT4,ATT5))
- factor_scale_inner_core_dble(:,:,:,1:min(ATT4,ATT5)) = factor_scale_crust_mantle_dble(:,:,:,1:min(ATT4,ATT5))
- factor_common_inner_core_dble(:,:,:,:,1:min(ATT4,ATT5)) = factor_common_crust_mantle_dble(:,:,:,:,1:min(ATT4,ATT5))
- ! Tell the Attenuation Code about the IDOUBLING regions within the Mesh
- call set_attenuation_regions_1D(RICB, RCMB, R670, R220, R80, AM_V)
- endif ! ATTENUATION_3D
+!< Hejun
if(CUSTOM_REAL == SIZE_REAL) then
factor_scale_crust_mantle = sngl(factor_scale_crust_mantle_dble)
@@ -2431,17 +2451,18 @@
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
-
+!> Hejun
! ATTENUATION_3D get scale_factor
- if(ATTENUATION_3D_VAL) then
+ !if(ATTENUATION_3D_VAL) then
! tau_mu and tau_sigma need to reference a point in the mesh
scale_factor = factor_scale_crust_mantle(i,j,k,ispec)
- else
- iglob = ibool_crust_mantle(i,j,k,ispec)
- dist_cr = xstore_crust_mantle(iglob)
- call get_attenuation_index(idoubling_crust_mantle(ispec), dble(dist_cr), iregion_selected, .FALSE., AM_V)
- scale_factor = factor_scale_crust_mantle(1,1,1,iregion_selected)
- endif ! ATTENUATION_3D
+ !else
+ ! iglob = ibool_crust_mantle(i,j,k,ispec)
+ ! dist_cr = xstore_crust_mantle(iglob)
+ ! call get_attenuation_index(idoubling_crust_mantle(ispec), dble(dist_cr), iregion_selected, .FALSE., AM_V)
+ ! scale_factor = factor_scale_crust_mantle(1,1,1,iregion_selected)
+ !endif ! ATTENUATION_3D
+!< Hejun
if(ANISOTROPIC_3D_MANTLE_VAL) then
scale_factor_minus_one = scale_factor - 1.
@@ -2486,15 +2507,16 @@
do k=1,NGLLZ
do j=1,NGLLY
do i=1,NGLLX
-
- if(ATTENUATION_3D_VAL) then
+!> Hejun
+! if(ATTENUATION_3D_VAL) then
scale_factor_minus_one = factor_scale_inner_core(i,j,k,ispec) - 1.0
- else
- iglob = ibool_inner_core(i,j,k,ispec)
- dist_cr = xstore_inner_core(iglob)
- call get_attenuation_index(idoubling_inner_core(ispec), dble(dist_cr), iregion_selected, .TRUE., AM_V)
- scale_factor_minus_one = factor_scale_inner_core(1,1,1,iregion_selected) - 1.
- endif
+! else
+! iglob = ibool_inner_core(i,j,k,ispec)
+! dist_cr = xstore_inner_core(iglob)
+! call get_attenuation_index(idoubling_inner_core(ispec), dble(dist_cr), iregion_selected, .TRUE., AM_V)
+! scale_factor_minus_one = factor_scale_inner_core(1,1,1,iregion_selected) - 1.
+! endif
+!< Hejun
if(ANISOTROPIC_INNER_CORE_VAL) then
mul = muvstore_inner_core(i,j,k,ispec)
@@ -2510,11 +2532,13 @@
+ scale_factor_minus_one * mul
endif
- if(ATTENUATION_3D_VAL) then
- muvstore_inner_core(i,j,k,ispec) = muvstore_inner_core(i,j,k,ispec) * factor_scale_inner_core(i,j,k,ispec)
- else
- muvstore_inner_core(i,j,k,ispec) = muvstore_inner_core(i,j,k,ispec) * factor_scale_inner_core(1,1,1,iregion_selected)
- endif
+!> Hejun
+ !if(ATTENUATION_3D_VAL) then
+ muvstore_inner_core(i,j,k,ispec) = muvstore_inner_core(i,j,k,ispec) * factor_scale_inner_core(i,j,k,ispec)
+ !else
+ ! muvstore_inner_core(i,j,k,ispec)=muvstore_inner_core(i,j,k,ispec)*factor_scale_inner_core(1,1,1,iregion_selected)
+ !endif
+!< Hejun
enddo
enddo
@@ -3667,7 +3691,7 @@
R_memory_crust_mantle,epsilondev_crust_mantle,eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
alphaval,betaval,gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),COMPUTE_AND_STORE_STRAIN,AM_V)
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),COMPUTE_AND_STORE_STRAIN)
else
call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
displ_crust_mantle,accel_crust_mantle, &
@@ -3691,7 +3715,7 @@
R_memory_crust_mantle,epsilondev_crust_mantle,eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
alphaval,betaval,gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),COMPUTE_AND_STORE_STRAIN,AM_V)
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),COMPUTE_AND_STORE_STRAIN)
endif
if (SIMULATION_TYPE == 3) then
@@ -3719,7 +3743,7 @@
b_R_memory_crust_mantle,b_epsilondev_crust_mantle,b_eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),COMPUTE_AND_STORE_STRAIN,AM_V)
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),COMPUTE_AND_STORE_STRAIN)
else
call compute_forces_crust_mantle(minus_gravity_table,density_table,minus_deriv_gravity_table, &
b_displ_crust_mantle,b_accel_crust_mantle, &
@@ -3743,7 +3767,7 @@
b_R_memory_crust_mantle,b_epsilondev_crust_mantle,b_eps_trace_over_3_crust_mantle,one_minus_sum_beta_crust_mantle, &
b_alphaval,b_betaval,b_gammaval,factor_common_crust_mantle, &
size(factor_common_crust_mantle,2), size(factor_common_crust_mantle,3), &
- size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),COMPUTE_AND_STORE_STRAIN,AM_V)
+ size(factor_common_crust_mantle,4), size(factor_common_crust_mantle,5),COMPUTE_AND_STORE_STRAIN)
endif
endif
@@ -3996,7 +4020,7 @@
alphaval,betaval,gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),COMPUTE_AND_STORE_STRAIN,AM_V)
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5),COMPUTE_AND_STORE_STRAIN)
else
call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
displ_inner_core,accel_inner_core, &
@@ -4013,7 +4037,7 @@
alphaval,betaval,gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),COMPUTE_AND_STORE_STRAIN,AM_V)
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5),COMPUTE_AND_STORE_STRAIN)
endif
if (SIMULATION_TYPE == 3) then
@@ -4033,7 +4057,7 @@
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),COMPUTE_AND_STORE_STRAIN,AM_V)
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5),COMPUTE_AND_STORE_STRAIN)
else
call compute_forces_inner_core(minus_gravity_table,density_table,minus_deriv_gravity_table, &
b_displ_inner_core,b_accel_inner_core, &
@@ -4050,7 +4074,7 @@
b_alphaval,b_betaval,b_gammaval, &
factor_common_inner_core, &
size(factor_common_inner_core,2), size(factor_common_inner_core,3), &
- size(factor_common_inner_core,4), size(factor_common_inner_core,5),COMPUTE_AND_STORE_STRAIN,AM_V)
+ size(factor_common_inner_core,4), size(factor_common_inner_core,5),COMPUTE_AND_STORE_STRAIN)
endif
endif
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/topo_bathy.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/topo_bathy.f90 2010-01-21 16:01:30 UTC (rev 16155)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/topo_bathy.f90 2010-01-21 16:14:52 UTC (rev 16156)
@@ -44,6 +44,10 @@
double precision samples_per_degree_topo
double precision xlo
+!> Hejun
+ double precision:: lon_corner,lat_corner,ratio_lon,ratio_lat
+!< Hejun
+
xlo = xlon
if(xlon < 0.d0) xlo = xlo + 360.d0
@@ -58,9 +62,30 @@
iel1 = int(xlo/samples_per_degree_topo)
if(iel1 <= 0 .or. iel1 > NX_BATHY) iel1 = NX_BATHY
+!> Hejun
+! Use bilinear interpolation rather nearest point interpolation
! convert integer value to double precision
- value = dble(ibathy_topo(iel1,iadd1))
+! value = dble(ibathy_topo(iel1,iadd1))
+
+ lon_corner=iel1*samples_per_degree_topo
+ lat_corner=90.d0-iadd1*samples_per_degree_topo
+
+ ratio_lon = (xlo-lon_corner)/samples_per_degree_topo
+ ratio_lat = (xlat-lat_corner)/samples_per_degree_topo
+
+ if(ratio_lon<0.0) ratio_lon=0.0
+ if(ratio_lon>1.0) ratio_lon=1.0
+ if(ratio_lat<0.0) ratio_lat=0.0
+ if(ratio_lat>1.0) ratio_lat=1.0
+
+! convert integer value to double precision
+ value = dble(ibathy_topo(iel1,iadd1))*(1-ratio_lon)*(1.-ratio_lat)+&
+ dble(ibathy_topo(iel1+1,iadd1))*ratio_lon*(1.-ratio_lat)+&
+ dble(ibathy_topo(iel1+1,iadd1+1))*ratio_lon*ratio_lat+&
+ dble(ibathy_topo(iel1,iadd1+1))*(1.-ratio_lon)*ratio_lat
+ !< Hejun
+
end subroutine get_topo_bathy
! -------------------------------------------
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/write_AVS_DX_global_data.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/write_AVS_DX_global_data.f90 2010-01-21 16:01:30 UTC (rev 16155)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/write_AVS_DX_global_data.f90 2010-01-21 16:14:52 UTC (rev 16156)
@@ -193,3 +193,166 @@
end subroutine write_AVS_DX_global_data
+!> Hejun
+! write material information for gll points
+ subroutine write_AVS_DX_global_data_gll(prname,nspec, &
+ xstore,ystore,zstore,rhostore,kappavstore,muvstore,Qmustore)
+
+ implicit none
+
+ include "constants.h"
+
+ integer nspec
+ character(len=150) prname
+
+ double precision xstore(NGLLX,NGLLY,NGLLZ,nspec)
+ double precision ystore(NGLLX,NGLLY,NGLLZ,nspec)
+ double precision zstore(NGLLX,NGLLY,NGLLZ,nspec)
+
+ real(kind=CUSTOM_REAL) kappavstore(NGLLX,NGLLY,NGLLZ,nspec)
+ real(kind=CUSTOM_REAL) muvstore(NGLLX,NGLLY,NGLLZ,nspec)
+ real(kind=CUSTOM_REAL) rhostore(NGLLX,NGLLY,NGLLZ,nspec)
+ double precision:: Qmustore(NGLLX,NGLLY,NGLLZ,nspec)
+
+ ! local parameters
+ double precision,dimension(8):: vp,vs,rho,Qmu
+ double precision:: vp_average,vs_average,rho_average,Qmu_average
+
+ integer flag(NGLLX,NGLLY,NGLLZ,nspec)
+
+ integer ispec,i,j,k
+ integer iglob1,iglob2,iglob3,iglob4,iglob5,iglob6,iglob7,iglob8
+ integer numpoin,nelem
+
+
+! writing points
+ open(unit=10,file=prname(1:len_trim(prname))//'AVS_DXpoints_gll.txt',status='unknown')
+
+! number of points in AVS or DX file
+ write(10,*) nspec*NGLLX*NGLLY*NGLLZ
+
+
+! output global AVS or DX points
+ numpoin = 0
+ do ispec=1,nspec
+ do k = 1,NGLLZ
+ do j = 1,NGLLY
+ do i = 1,NGLLX
+ numpoin = numpoin + 1
+ write(10,*) numpoin,sngl(xstore(i,j,k,ispec)),&
+ sngl(ystore(i,j,k,ispec)),sngl(zstore(i,j,k,ispec))
+ flag(i,j,k,ispec) = numpoin
+ end do
+ end do
+ end do
+ enddo
+
+ close(10)
+
+! writing elements
+ open(unit=10,file=prname(1:len_trim(prname))//'AVS_DXelements_gll.txt',status='unknown')
+
+
+! number of elements in AVS or DX file
+ write(10,*) nspec*(NGLLX-1)*(NGLLY-1)*(NGLLZ-1)
+
+ nelem = 0
+! output global AVS or DX elements
+ do ispec=1,nspec
+ do k = 1,NGLLZ-1
+ do j = 1,NGLLY-1
+ do i = 1,NGLLX-1
+ nelem = nelem + 1
+ iglob1=flag(i,j,k,ispec)
+ iglob2=flag(i+1,j,k,ispec)
+ iglob3=flag(i+1,j+1,k,ispec)
+ iglob4=flag(i,j+1,k,ispec)
+ iglob5=flag(i,j,k+1,ispec)
+ iglob6=flag(i+1,j,k+1,ispec)
+ iglob7=flag(i+1,j+1,k+1,ispec)
+ iglob8=flag(i,j+1,k+1,ispec)
+
+ write(10,*) nelem,iglob1, &
+ iglob2,iglob3,iglob4,&
+ iglob5,iglob6,iglob7,iglob8
+ end do
+ end do
+ end do
+ enddo
+
+ close(10)
+
+! writing elements properity
+ open(unit=1001,file=prname(1:len_trim(prname))//'AVS_DXmaterials_gll.txt',status='unknown')
+
+! number of elements in AVS or DX file
+ write(1001,*) nspec*(NGLLX-1)*(NGLLY-1)*(NGLLZ-1)
+
+ nelem = 0
+! output global AVS or DX elements
+ do ispec=1,nspec
+ do k = 1,NGLLZ-1
+ do j = 1,NGLLY-1
+ do i = 1,NGLLX-1
+ nelem = nelem + 1
+ rho(1)=dble(rhostore(i,j,k,ispec))
+ vs(1)=dble(sqrt(muvstore(i,j,k,ispec)/rhostore(i,j,k,ispec)))
+ vp(1)=dble(sqrt(kappavstore(i,j,k,ispec)/rhostore(i,j,k,ispec)+4.d0*vs(1)*vs(1)/3.d0))
+ Qmu(1)=dble(Qmustore(i,j,k,ispec))
+
+ rho(2)=dble(rhostore(i+1,j,k,ispec))
+ vs(2)=dble(sqrt(muvstore(i+1,j,k,ispec)/rhostore(i+1,j,k,ispec)))
+ vp(2)=dble(sqrt(kappavstore(i+1,j,k,ispec)/rhostore(i+1,j,k,ispec)+4.d0*vs(2)*vs(2)/3.d0))
+ Qmu(2)=dble(Qmustore(i+1,j,k,ispec))
+
+ rho(3)=dble(rhostore(i+1,j+1,k,ispec))
+ vs(3)=dble(sqrt(muvstore(i+1,j+1,k,ispec)/rhostore(i+1,j+1,k,ispec)))
+ vp(3)=dble(sqrt(kappavstore(i+1,j+1,k,ispec)/rhostore(i+1,j+1,k,ispec)+4.d0*vs(3)*vs(3)/3.d0))
+ Qmu(3)=dble(Qmustore(i+1,j+1,k,ispec))
+
+ rho(4)=dble(rhostore(i,j+1,k,ispec))
+ vs(4)=dble(sqrt(muvstore(i,j+1,k,ispec)/rhostore(i,j+1,k,ispec)))
+ vp(4)=dble(sqrt(kappavstore(i,j+1,k,ispec)/rhostore(i,j+1,k,ispec)+4.d0*vs(4)*vs(4)/3.d0))
+ Qmu(4)=dble(Qmustore(i,j+1,k,ispec))
+
+ rho(5)=dble(rhostore(i,j,k+1,ispec))
+ vs(5)=dble(sqrt(muvstore(i,j,k+1,ispec)/rhostore(i,j,k+1,ispec)))
+ vp(5)=dble(sqrt(kappavstore(i,j,k+1,ispec)/rhostore(i,j,k+1,ispec)+4.d0*vs(5)*vs(5)/3.d0))
+ Qmu(5)=dble(Qmustore(i,j,k+1,ispec))
+
+ rho(6)=dble(rhostore(i+1,j,k+1,ispec))
+ vs(6)=dble(sqrt(muvstore(i+1,j,k+1,ispec)/rhostore(i+1,j,k+1,ispec)))
+ vp(6)=dble(sqrt(kappavstore(i+1,j,k+1,ispec)/rhostore(i+1,j,k+1,ispec)+4.d0*vs(6)*vs(6)/3.d0))
+ Qmu(6)=dble(Qmustore(i+1,j,k+1,ispec))
+
+ rho(7)=dble(rhostore(i+1,j+1,k+1,ispec))
+ vs(7)=dble(sqrt(muvstore(i+1,j+1,k+1,ispec)/rhostore(i+1,j+1,k+1,ispec)))
+ vp(7)=dble(sqrt(kappavstore(i+1,j+1,k+1,ispec)/rhostore(i+1,j+1,k+1,ispec)+4.d0*vs(7)*vs(7)/3.d0))
+ Qmu(7)=dble(Qmustore(i+1,j+1,k+1,ispec))
+
+ rho(8)=dble(rhostore(i,j+1,k+1,ispec))
+ vs(8)=dble(sqrt(muvstore(i,j+1,k+1,ispec)/rhostore(i,j+1,k+1,ispec)))
+ vp(8)=dble(sqrt(kappavstore(i,j+1,k+1,ispec)/rhostore(i,j+1,k+1,ispec)+4.d0*vs(8)*vs(8)/3.d0))
+ Qmu(8)=dble(Qmustore(i,j+1,k+1,ispec))
+
+ !rho_average=sum(rho(1:4))/4.d0
+ !vp_average=sum(vp(1:4))/4.d0
+ !vs_average=sum(vs(1:4))/4.d0
+ !Qmu_average=sum(Qmu(1:4))/4.d0
+ rho_average=rho(1)
+ vp_average=vp(1)
+ vs_average=vs(1)
+ Qmu_average=Qmu(1)
+
+ write(1001,*) nelem,rho_average,vp_average,vs_average,Qmu_average
+ end do
+ end do
+ end do
+ enddo
+
+ close(1001)
+
+ end subroutine write_AVS_DX_global_data_gll
+!< Hejun
+
+
More information about the CIG-COMMITS
mailing list