[cig-commits] r21020 - in seismo/3D: SPECFEM3D/trunk SPECFEM3D_GLOBE/trunk/setup SPECFEM3D_GLOBE/trunk/src/specfem3D
elliott.sales.de.andrade at geodynamics.org
elliott.sales.de.andrade at geodynamics.org
Mon Nov 12 03:51:41 PST 2012
Author: elliott.sales.de.andrade
Date: 2012-11-12 03:51:41 -0800 (Mon, 12 Nov 2012)
New Revision: 21020
Added:
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_regular_points.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/save_regular_kernels.f90
Modified:
seismo/3D/SPECFEM3D/trunk/todo_list_please_dont_remove.txt
seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/Makefile.in
seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90
Log:
Remove an old TODO item that's finished a while ago.
Modified: seismo/3D/SPECFEM3D/trunk/todo_list_please_dont_remove.txt
===================================================================
--- seismo/3D/SPECFEM3D/trunk/todo_list_please_dont_remove.txt 2012-11-12 08:57:29 UTC (rev 21019)
+++ seismo/3D/SPECFEM3D/trunk/todo_list_please_dont_remove.txt 2012-11-12 11:51:41 UTC (rev 21020)
@@ -442,18 +442,6 @@
------------------------------------------------
-sensitivity kernel output:
-------------------------------------------------
-
-- suggestion 31:
-----------------
-
-Elliott Sales de Andrade and Qinya will implement a routine to interpolate and output the sensitiviy kernels on a topologically-regular grid
-instead of the non-structured SEM grid (this can be useful for many processing packages that require a topologically-regular grid
-of points as input).
-
-
-------------------------------------------------
databases in SPECFEM2D:
------------------------------------------------
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in 2012-11-12 08:57:29 UTC (rev 21019)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/setup/constants.h.in 2012-11-12 11:51:41 UTC (rev 21020)
@@ -279,6 +279,18 @@
! (default is to use transverse isotropy only between MOHO and 220 )
logical, parameter :: USE_FULL_TISO_MANTLE = .false.
+! output kernels on a regular grid instead of on the mesh points
+ logical, parameter :: SAVE_REGULAR_KL = .false.
+
+! regular kernel parameters
+ character(len=*), parameter :: PATHNAME_KL_REG = 'DATA/kl_reg_grid.txt'
+ integer, parameter :: NM_KL_REG_LAYER = 100
+ integer, parameter :: NM_KL_REG_PTS = 300000
+ real, parameter :: KL_REG_MIN_LON = 0.0
+ real, parameter :: KL_REG_MAX_LON = 360.0
+ real, parameter :: KL_REG_MIN_LAT = -90.0
+ real, parameter :: KL_REG_MAX_LAT = +90.0
+
!!-----------------------------------------------------------
!!
!! time stamp information
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/Makefile.in
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/Makefile.in 2012-11-12 08:57:29 UTC (rev 21019)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/Makefile.in 2012-11-12 11:51:41 UTC (rev 21020)
@@ -103,6 +103,7 @@
$O/intgrl.o \
$O/lagrange_poly.o \
$O/locate_receivers.o \
+ $O/locate_regular_points.o \
$O/locate_sources.o \
$O/make_ellipticity.o \
$O/make_gravity.o \
@@ -119,6 +120,7 @@
$O/recompute_jacobian.o \
$O/reduce.o \
$O/rthetaphi_xyz.o \
+ $O/save_regular_kernels.o \
$O/write_c_binary.o \
$O/write_seismograms.o \
$O/write_output_ASCII.o \
@@ -379,6 +381,9 @@
$O/save_kernels.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/save_kernels.f90
${FCCOMPILE_NO_CHECK} -c -o $O/save_kernels.o ${FCFLAGS_f90} $S/save_kernels.f90
+$O/save_regular_kernels.o: ${SETUP}/constants.h ${OUTPUT}/values_from_mesher.h $S/save_regular_kernels.f90
+ ${FCCOMPILE_NO_CHECK} -c -o $O/save_regular_kernels.o ${FCFLAGS_f90} $S/save_regular_kernels.f90
+
###
### specfem3D - regular compilation options here
###
@@ -485,6 +490,9 @@
$O/locate_receivers.o: ${SETUP}/constants.h $S/locate_receivers.f90
${MPIFCCOMPILE_CHECK} -c -o $O/locate_receivers.o ${FCFLAGS_f90} $S/locate_receivers.f90
+$O/locate_regular_points.o: ${SETUP}/constants.h $S/locate_regular_points.f90
+ ${MPIFCCOMPILE_CHECK} -c -o $O/locate_regular_points.o ${FCFLAGS_f90} $S/locate_regular_points.f90
+
$O/locate_sources.o: ${SETUP}/constants.h $S/locate_sources.f90
${MPIFCCOMPILE_CHECK} -c -o $O/locate_sources.o ${FCFLAGS_f90} $S/locate_sources.f90
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_regular_points.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_regular_points.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/locate_regular_points.f90 2012-11-12 11:51:41 UTC (rev 21020)
@@ -0,0 +1,491 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 5 . 1
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Princeton University, USA
+! and CNRS / INRIA / University of Pau, France
+! (c) Princeton University and CNRS / INRIA / University of Pau
+! April 2011
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+subroutine read_kl_regular_grid(GRID)
+
+ implicit none
+ include 'constants.h'
+
+ type kl_reg_grid_variables
+ sequence
+ real dlat
+ real dlon
+ integer nlayer
+ real rlayer(NM_KL_REG_LAYER)
+ integer ndoubling(NM_KL_REG_LAYER)
+ integer nlat(NM_KL_REG_LAYER)
+ integer nlon(NM_KL_REG_LAYER)
+ integer npts_total
+ integer npts_before_layer(NM_KL_REG_LAYER+1)
+ end type kl_reg_grid_variables
+
+ type (kl_reg_grid_variables), intent(inout) :: GRID
+
+ integer :: ios,nlayer,i,nlat,nlon,npts_this_layer
+
+ ! improvements to make: read-in by master and broadcast to all slaves
+ open(10,file=PATHNAME_KL_REG,iostat=ios,status='old',action='read')
+
+ read(10,*) GRID%dlat, GRID%dlon
+
+ nlayer = 1
+ do while (nlayer <= NM_KL_REG_LAYER)
+ read(10,*,iostat=ios) GRID%rlayer(nlayer), GRID%ndoubling(nlayer)
+ if (ios.ne.0) exit
+ nlayer = nlayer + 1
+ enddo
+ close(10)
+
+ if (nlayer > NM_KL_REG_LAYER) then
+ call exit_MPI('Increase NM_KL_REG_LAYER limit')
+ endif
+
+ GRID%nlayer = nlayer
+
+ GRID%npts_total = 0
+ GRID%npts_before_layer = 0
+ do i = 1, nlayer
+ nlon = floor((KL_REG_MAX_LON-KL_REG_MIN_LON)/(GRID%dlon*GRID%ndoubling(i)))+1
+ GRID%nlon(i) = nlon
+ nlat = floor((KL_REG_MAX_LAT-KL_REG_MIN_LAT)/(GRID%dlat*GRID%ndoubling(i)))+1
+ GRID%nlat(i) = nlat
+ npts_this_layer = nlon * nlat
+ GRID%npts_total = GRID%npts_total + npts_this_layer
+ GRID%npts_before_layer(i+1) = GRID%npts_before_layer(i) + npts_this_layer
+ enddo
+ if (GRID%npts_total <= 0) then
+ call exit_MPI('No Model points read in')
+ endif
+
+end subroutine read_kl_regular_grid
+
+!==============================================================
+
+subroutine find_regular_grid_slice_number(slice_number, GRID, &
+ NCHUNKS, NPROC_XI, NPROC_ETA)
+
+ implicit none
+ include 'constants.h'
+
+ integer, intent(out) :: slice_number(*)
+
+ type kl_reg_grid_variables
+ sequence
+ real dlat
+ real dlon
+ integer nlayer
+ real rlayer(NM_KL_REG_LAYER)
+ integer ndoubling(NM_KL_REG_LAYER)
+ integer nlat(NM_KL_REG_LAYER)
+ integer nlon(NM_KL_REG_LAYER)
+ integer npts_total
+ integer npts_before_layer(NM_KL_REG_LAYER+1)
+ end type kl_reg_grid_variables
+ type (kl_reg_grid_variables), intent(in) :: GRID
+
+ integer, intent(in) :: NCHUNKS,NPROC_XI,NPROC_ETA
+
+ real(kind=CUSTOM_REAL) :: xi_width, eta_width
+ integer :: nproc, ilayer, isp, ilat, ilon, k, chunk_isp
+ integer :: iproc_xi, iproc_eta
+ real :: lat,lon,th,ph,x,y,z,xik,etak,xi_isp,eta_isp,xi1,eta1
+
+ ! assuming 6 chunks full global simulations right now
+ if (NCHUNKS /= 6 .or. NPROC_XI /= NPROC_ETA) then
+ call exit_MPI('Only deal with 6 chunks at this moment')
+ endif
+
+ xi_width=PI/2; eta_width=PI/2; nproc=NPROC_XI
+ ilayer=0
+
+ do isp = 1,GRID%npts_total
+ if (isp == GRID%npts_before_layer(ilayer+1)+1) ilayer=ilayer+1
+ ilat = (isp - GRID%npts_before_layer(ilayer) - 1) / GRID%nlat(ilayer)
+ ilon = (isp - GRID%npts_before_layer(ilayer) - 1) - ilat * GRID%nlat(ilayer)
+
+ ! (lat,lon,radius) for isp point
+ lat = KL_REG_MIN_LAT + ilat * GRID%dlat * GRID%ndoubling(ilayer)
+ th = (90 - lat) * DEGREES_TO_RADIANS
+ lon = KL_REG_MIN_LON + ilon * GRID%dlon * GRID%ndoubling(ilayer)
+ ph = lon * DEGREES_TO_RADIANS
+ x = sin(th) * cos(ph); y = sin(th) * sin(ph); z = cos(th)
+
+ ! figure out slice number
+ chunk_isp = 1; xi_isp = 0; eta_isp = 0
+ do k = 1, NCHUNKS
+ call chunk_map(k, x, y, z, xik, etak)
+ if (abs(xik) <= PI/4 .and. abs(etak) <= PI/4) then
+ chunk_isp = k; xi_isp = xik; eta_isp = etak; exit
+ endif
+ enddo
+ xi1 = xi_isp / xi_width * 2; eta1 = eta_isp / eta_width * 2
+ iproc_xi = floor((xi1+1)/2 * nproc)
+ iproc_eta = floor((eta1+1)/2 * nproc)
+ slice_number(isp) = nproc * nproc * (chunk_isp-1) + nproc * iproc_eta + iproc_xi
+ enddo
+
+end subroutine find_regular_grid_slice_number
+
+!==============================================================
+
+! how about using single precision for the iterations?
+subroutine locate_reg_points(myrank,npoints_slice,points_slice,GRID, &
+ NEX_XI,nspec,xstore,ystore,zstore,ibool, &
+ xigll,yigll,zigll,ispec_reg, &
+ hxir_reg,hetar_reg,hgammar_reg)
+
+ implicit none
+ include 'constants.h'
+
+ ! declarations of regular grid model
+ integer, intent(in) :: myrank, npoints_slice
+ integer, dimension(NM_KL_REG_PTS), intent(in) :: points_slice
+
+ type kl_reg_grid_variables
+ sequence
+ real dlat
+ real dlon
+ integer nlayer
+ real rlayer(NM_KL_REG_LAYER)
+ integer ndoubling(NM_KL_REG_LAYER)
+ integer nlat(NM_KL_REG_LAYER)
+ integer nlon(NM_KL_REG_LAYER)
+ integer npts_total
+ integer npts_before_layer(NM_KL_REG_LAYER+1)
+ end type kl_reg_grid_variables
+ type (kl_reg_grid_variables), intent(in) :: GRID
+
+ ! simulation geometry
+ integer, intent(in) :: NEX_XI, nspec
+ real(kind=CUSTOM_REAL), dimension(*), intent(in) :: xstore,ystore,zstore
+ integer, dimension(NGLLX,NGLLY,NGLLZ,*), intent(in) :: ibool
+
+ ! Gauss-Lobatto-Legendre points of integration and weights
+ double precision, dimension(NGLLX), intent(in) :: xigll
+ double precision, dimension(NGLLY), intent(in) :: yigll
+ double precision, dimension(NGLLZ), intent(in) :: zigll
+
+ ! output
+ integer, dimension(NM_KL_REG_PTS), intent(out) :: ispec_reg
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NM_KL_REG_PTS), intent(out) :: hxir_reg
+ real(kind=CUSTOM_REAL), dimension(NGLLY,NM_KL_REG_PTS), intent(out) :: hetar_reg
+ real(kind=CUSTOM_REAL), dimension(NGLLZ,NM_KL_REG_PTS), intent(out) :: hgammar_reg
+
+ ! GLL number of anchors
+ integer, dimension(NGNOD) :: iaddx, iaddy, iaddr
+
+ integer :: i, j, k, isp, ilayer, ilat, ilon, iglob, ix_in, iy_in, iz_in
+ integer :: ispec_in, ispec, iter_loop, ia, ipoint
+ double precision :: lat, lon, radius, th, ph, x,y,z
+ double precision :: x_target, y_target, z_target
+ double precision :: distmin,dist,typical_size
+ double precision :: xi,eta,gamma,dx,dy,dz,dxi,deta,dgamma
+ double precision :: xix,xiy,xiz
+ double precision :: etax,etay,etaz
+ double precision :: gammax,gammay,gammaz
+
+ logical locate_target
+ double precision, dimension(NGNOD) :: xelm, yelm, zelm
+
+ double precision, dimension(NGLLX) :: hxir
+ double precision, dimension(NGLLY) :: hetar
+ double precision, dimension(NGLLZ) :: hgammar
+
+ ! DEBUG
+ !real(kind=CUSTOM_REAL), dimension(npoints_slice) :: dist_final
+
+ !---------------------------
+
+ call hex_nodes2(iaddx,iaddy,iaddr)
+
+ ! compute typical size of elements at the surface
+ typical_size = TWO_PI * R_UNIT_SPHERE / (4.*NEX_XI)
+
+ ! use 10 times the distance as a criterion for source detection
+ typical_size = 10. * typical_size
+
+ ! DEBUG
+ !dist_final=HUGEVAL
+
+ do ipoint = 1, npoints_slice
+ isp = points_slice(ipoint)
+ do ilayer = 1, GRID%nlayer
+ if (isp <= GRID%npts_before_layer(ilayer+1)) exit
+ enddo
+
+ ilat = (isp - GRID%npts_before_layer(ilayer) - 1) / GRID%nlat(ilayer)
+ ilon = (isp - GRID%npts_before_layer(ilayer) - 1) - ilat * GRID%nlat(ilayer)
+
+ ! (lat,lon,radius) for isp point
+ lat = KL_REG_MIN_LAT + ilat * GRID%dlat * GRID%ndoubling(ilayer)
+ lon = KL_REG_MIN_LON + ilon * GRID%dlon * GRID%ndoubling(ilayer)
+ ! convert radius to meters and then scale
+ radius = GRID%rlayer(ilayer) * 1000.0 / R_EARTH
+ ! (x,y,z) for isp point
+ th = (90 - lat) * DEGREES_TO_RADIANS; ph = lon * DEGREES_TO_RADIANS
+ x_target = radius * sin(th) * cos(ph)
+ y_target = radius * sin(th) * sin(ph)
+ z_target = radius * cos(th)
+
+ ! first exclude elements too far away
+ locate_target = .false.; distmin = HUGEVAL
+ do ispec = 1,nspec
+ iglob = ibool(1,1,1,ispec)
+ dist = dsqrt((x_target - xstore(iglob))**2 &
+ + (y_target - ystore(iglob))**2 &
+ + (z_target - zstore(iglob))**2)
+ if (dist > typical_size) cycle
+
+ locate_target = .true.
+ ! loop only on points inside the element
+ ! exclude edges to ensure this point is not
+ ! shared with other elements
+ ! can be improved if we have a better algorithm of determining if a point
+ ! exists inside a 3x3x3 specfem element ???
+
+ do k = 2, NGLLZ-1
+ do j = 2, NGLLY-1
+ do i = 2, NGLLX-1
+ iglob = ibool(i,j,k,ispec)
+ dist = dsqrt((x_target - xstore(iglob))**2 &
+ +(y_target - ystore(iglob))**2 &
+ +(z_target - zstore(iglob))**2)
+ if (dist < distmin) then
+ ix_in=i; iy_in=j; iz_in=k; ispec_in=ispec; distmin=dist
+ endif
+ enddo
+ enddo
+ enddo
+
+ enddo
+ if (.not. locate_target) stop 'error in point_source() array'
+
+ xi = xigll(ix_in)
+ eta = yigll(iy_in)
+ gamma = zigll(iz_in)
+ ispec_reg(ipoint) = ispec_in
+
+ ! anchors
+ do ia = 1, NGNOD
+ iglob = ibool(iaddx(ia), iaddy(ia), iaddr(ia), ispec_in)
+ xelm(ia) = dble(xstore(iglob))
+ yelm(ia) = dble(ystore(iglob))
+ zelm(ia) = dble(zstore(iglob))
+ enddo
+
+ ! iterate to solve the nonlinear system
+ do iter_loop = 1,NUM_ITER
+
+ ! recompute jacobian for the new point
+ call recompute_jacobian(xelm,yelm,zelm, xi,eta,gamma, x,y,z, &
+ xix,xiy,xiz, etax,etay,etaz, gammax,gammay,gammaz)
+
+ ! compute distance to target location
+ dx = - (x - x_target)
+ dy = - (y - y_target)
+ dz = - (z - z_target)
+
+ ! compute increments
+ dxi = xix*dx + xiy*dy + xiz*dz
+ deta = etax*dx + etay*dy + etaz*dz
+ dgamma = gammax*dx + gammay*dy + gammaz*dz
+
+ ! update values
+ xi = xi + dxi
+ eta = eta + deta
+ gamma = gamma + dgamma
+
+ ! Debugging
+ !if (abs(xi) > 1.d0+TINYVAL .or. abs(eta) > 1.d0+TINYVAL &
+ ! .or. abs(gamma) > 1.0d0+TINYVAL) then
+ ! print *, 'Outside the element ', myrank, ipoint,' : ', &
+ ! iter_loop,xi,eta,gamma
+ !endif
+
+ ! impose that we stay in that element
+ ! (useful if user gives a source outside the mesh for instance)
+ if (xi > 1.d0) xi = 1.d0
+ if (xi < -1.d0) xi = -1.d0
+ if (eta > 1.d0) eta = 1.d0
+ if (eta < -1.d0) eta = -1.d0
+ if (gamma > 1.d0) gamma = 1.d0
+ if (gamma < -1.d0) gamma = -1.d0
+
+ enddo
+
+ ! DEBUG: recompute jacobian for the new point (can be commented after debug)
+ !call recompute_jacobian(xelm,yelm,zelm,xi,eta,gamma,x,y,z,xix,xiy,xiz,etax,etay,etaz,gammax,gammay,gammaz)
+ !dist_final(ipoint)=dsqrt((x_target-x)**2+(y_target-y)**2+(z_target-z)**2)
+
+ ! store l(xi),l(eta),l(gamma)
+ call lagrange_any2(xi, NGLLX, xigll, hxir)
+ call lagrange_any2(eta, NGLLY, yigll, hetar)
+ call lagrange_any2(gamma, NGLLZ, zigll, hgammar)
+ hxir_reg(:,ipoint) = hxir
+ hetar_reg(:,ipoint) = hetar
+ hgammar_reg(:,ipoint) = hgammar
+
+ enddo ! ipoint
+
+! DEBUG
+! print *, 'Maximum distance discrepancy ', maxval(dist_final(1:npoints_slice))
+
+end subroutine locate_reg_points
+
+!==============================================================
+
+subroutine hex_nodes2(iaddx,iaddy,iaddz)
+
+ implicit none
+ include 'constants.h'
+
+ integer, dimension(NGNOD), intent(out) :: iaddx,iaddy,iaddz
+ integer :: ia
+
+ ! define topology of the control element
+ call hex_nodes(iaddx,iaddy,iaddz)
+
+ ! define coordinates of the control points of the element
+ do ia=1,NGNOD
+
+ if (iaddx(ia) == 0) then
+ iaddx(ia) = 1
+ else if (iaddx(ia) == 1) then
+ iaddx(ia) = (NGLLX+1)/2
+ else if (iaddx(ia) == 2) then
+ iaddx(ia) = NGLLX
+ else
+ stop 'incorrect value of iaddx'
+ endif
+
+ if (iaddy(ia) == 0) then
+ iaddy(ia) = 1
+ else if (iaddy(ia) == 1) then
+ iaddy(ia) = (NGLLY+1)/2
+ else if (iaddy(ia) == 2) then
+ iaddy(ia) = NGLLY
+ else
+ stop 'incorrect value of iaddy'
+ endif
+
+ if (iaddz(ia) == 0) then
+ iaddz(ia) = 1
+ else if (iaddz(ia) == 1) then
+ iaddz(ia) = (NGLLZ+1)/2
+ else if (iaddz(ia) == 2) then
+ iaddz(ia) = NGLLZ
+ else
+ stop 'incorrect value of iaddz'
+ endif
+
+ enddo
+
+end subroutine hex_nodes2
+
+!==============================================================
+
+subroutine lagrange_any2(xi,NGLL,xigll,h)
+
+! subroutine to compute the Lagrange interpolants based upon the GLL points
+! and their first derivatives at any point xi in [-1,1]
+
+ implicit none
+
+ double precision, intent(in) :: xi
+ integer, intent(in) :: NGLL
+ double precision, dimension(NGLL), intent(in) :: xigll
+ double precision, dimension(NGLL), intent(out) :: h
+
+ integer :: dgr,i
+ double precision :: prod1,prod2
+
+ do dgr=1,NGLL
+ prod1 = 1.0d0
+ prod2 = 1.0d0
+
+ do i=1,NGLL
+ if (i /= dgr) then
+ prod1 = prod1 * (xi - xigll(i))
+ prod2 = prod2 * (xigll(dgr) - xigll(i))
+ endif
+ enddo
+
+ h(dgr) = prod1 / prod2
+ enddo
+
+end subroutine lagrange_any2
+
+!==============================================================
+
+subroutine chunk_map(k,xx,yy,zz,xi,eta)
+
+ ! this program get the xi,eta for (xx,yy,zz)
+ ! point under the k'th chunk coordinate
+ ! transformation
+
+ implicit none
+ include 'constants.h'
+
+ integer, intent(in) :: k
+ real, intent(in) :: xx, yy, zz
+ real, intent(out) :: xi, eta
+
+ real :: x, y, z
+ real, parameter :: EPS=1e-6
+
+ x = xx; y = yy; z = zz
+ if (0 <= x .and. x < EPS) x = EPS
+ if (-EPS < x .and. x < 0) x = -EPS
+ if (0 <= y .and. y < EPS) y = EPS
+ if (-EPS < y .and. y < 0) y = -EPS
+ if (0 <= z .and. z < EPS) z = EPS
+ if (-EPS < z .and. z < 0) z = -EPS
+
+ if (k == CHUNK_AB) then
+ xi = atan(y/z); eta = atan(-x/z)
+ if (z < 0) xi = 10
+ else if (k == CHUNK_AC) then
+ xi = atan(-z/y); eta = atan(x/y)
+ if (y > 0) xi = 10
+ else if (k == CHUNK_BC) then
+ xi = atan(-z/x); eta = atan(-y/x)
+ if (x > 0) xi = 10
+ else if (k == CHUNK_AC_ANTIPODE) then
+ xi = atan(-z/y); eta = atan(-x/y)
+ if (y < 0) xi = 10
+ else if (k == CHUNK_BC_ANTIPODE) then
+ xi = atan(z/x); eta = atan(-y/x)
+ if (x < 0) xi = 10
+ else if (k == CHUNK_AB_ANTIPODE) then
+ xi = atan(y/z); eta = atan(x/z)
+ if (z > 0) xi = 10
+ else
+ stop 'chunk number k < 6'
+ end if
+
+end subroutine chunk_map
+
Added: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/save_regular_kernels.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/save_regular_kernels.f90 (rev 0)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/save_regular_kernels.f90 2012-11-12 11:51:41 UTC (rev 21020)
@@ -0,0 +1,513 @@
+!=====================================================================
+!
+! S p e c f e m 3 D G l o b e V e r s i o n 5 . 1
+! --------------------------------------------------
+!
+! Main authors: Dimitri Komatitsch and Jeroen Tromp
+! Princeton University, USA
+! and CNRS / INRIA / University of Pau, France
+! (c) Princeton University and CNRS / INRIA / University of Pau
+! April 2011
+!
+! This program is free software; you can redistribute it and/or modify
+! it under the terms of the GNU General Public License as published by
+! the Free Software Foundation; either version 2 of the License, or
+! (at your option) any later version.
+!
+! This program is distributed in the hope that it will be useful,
+! but WITHOUT ANY WARRANTY; without even the implied warranty of
+! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+! GNU General Public License for more details.
+!
+! You should have received a copy of the GNU General Public License along
+! with this program; if not, write to the Free Software Foundation, Inc.,
+! 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+!
+!=====================================================================
+
+ subroutine save_regular_kernels_crust_mantle(myrank, &
+ GRID, npoints_slice, &
+ hxir_reg, hetar_reg, hgammar_reg, ispec_reg, &
+ scale_t,scale_displ, &
+ cijkl_kl_crust_mantle,rho_kl_crust_mantle, &
+ alpha_kl_crust_mantle,beta_kl_crust_mantle, &
+ ystore_crust_mantle,zstore_crust_mantle, &
+ rhostore_crust_mantle,muvstore_crust_mantle, &
+ kappavstore_crust_mantle,ibool_crust_mantle, &
+ kappahstore_crust_mantle,muhstore_crust_mantle, &
+ eta_anisostore_crust_mantle,ispec_is_tiso_crust_mantle, &
+ ! --idoubling_crust_mantle, &
+ LOCAL_PATH)
+
+ implicit none
+
+ include "constants.h"
+ include "OUTPUT_FILES/values_from_mesher.h"
+
+ integer myrank
+
+ type kl_reg_grid_variables
+ sequence
+ real dlat
+ real dlon
+ integer nlayer
+ real rlayer(NM_KL_REG_LAYER)
+ integer ndoubling(NM_KL_REG_LAYER)
+ integer nlat(NM_KL_REG_LAYER)
+ integer nlon(NM_KL_REG_LAYER)
+ integer npts_total
+ integer npts_before_layer(NM_KL_REG_LAYER+1)
+ end type kl_reg_grid_variables
+ type (kl_reg_grid_variables), intent(in) :: GRID
+
+ integer, intent(in) :: npoints_slice
+ real, dimension(NGLLX, NM_KL_REG_PTS), intent(in) :: hxir_reg
+ real, dimension(NGLLY, NM_KL_REG_PTS), intent(in) :: hetar_reg
+ real, dimension(NGLLZ, NM_KL_REG_PTS), intent(in) :: hgammar_reg
+ integer, dimension(NM_KL_REG_PTS), intent(in) :: ispec_reg
+
+ double precision :: scale_t,scale_displ
+
+ real(kind=CUSTOM_REAL), dimension(21,NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT) :: &
+ cijkl_kl_crust_mantle
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE_ADJOINT) :: &
+ rho_kl_crust_mantle, beta_kl_crust_mantle, alpha_kl_crust_mantle
+
+ real(kind=CUSTOM_REAL), dimension(NGLOB_CRUST_MANTLE) :: &
+ ystore_crust_mantle,zstore_crust_mantle
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_ISO_MANTLE) :: &
+ rhostore_crust_mantle,kappavstore_crust_mantle,muvstore_crust_mantle
+
+ real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPECMAX_TISO_MANTLE) :: &
+ kappahstore_crust_mantle,muhstore_crust_mantle,eta_anisostore_crust_mantle
+
+! integer, dimension(NSPEC_CRUST_MANTLE) :: idoubling_crust_mantle
+ logical, dimension(NSPEC_CRUST_MANTLE) :: ispec_is_tiso_crust_mantle
+
+ integer, dimension(NGLLX,NGLLY,NGLLZ,NSPEC_CRUST_MANTLE) :: ibool_crust_mantle
+
+ character(len=150) LOCAL_PATH
+
+ ! local parameters
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: &
+ cijkl_kl_crust_mantle_reg
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
+ rho_kl_crust_mantle_reg, beta_kl_crust_mantle_reg, alpha_kl_crust_mantle_reg
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:), allocatable :: &
+ mu_kl_crust_mantle, kappa_kl_crust_mantle, rhonotprime_kl_crust_mantle
+ real(kind=CUSTOM_REAL),dimension(21) :: cijkl_kl_local
+ real(kind=CUSTOM_REAL) :: scale_kl,scale_kl_ani,scale_kl_rho
+ real(kind=CUSTOM_REAL) :: rhol,mul,kappal,rho_kl,alpha_kl,beta_kl
+ integer :: ispec,i,j,k,iglob
+ character(len=150) prname
+ double precision :: hlagrange
+ integer :: ipoint
+
+ ! transverse isotropic parameters
+ real(kind=CUSTOM_REAL), dimension(21) :: an_kl
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:),allocatable :: &
+ alphav_kl_crust_mantle,alphah_kl_crust_mantle, &
+ betav_kl_crust_mantle,betah_kl_crust_mantle, &
+ eta_kl_crust_mantle
+
+ ! bulk parameterization
+ real(kind=CUSTOM_REAL), dimension(:,:,:,:),allocatable :: &
+ bulk_c_kl_crust_mantle,bulk_beta_kl_crust_mantle, &
+ bulk_betav_kl_crust_mantle,bulk_betah_kl_crust_mantle
+ real(kind=CUSTOM_REAL) :: A,C,F,L,N,eta
+ real(kind=CUSTOM_REAL) :: muvl,kappavl,muhl,kappahl
+ real(kind=CUSTOM_REAL) :: alphav_sq,alphah_sq,betav_sq,betah_sq,bulk_sq
+
+ ! scaling factors
+ scale_kl = scale_t/scale_displ * 1.d9
+ ! For anisotropic kernels
+ ! final unit : [s km^(-3) GPa^(-1)]
+ scale_kl_ani = scale_t**3 / (RHOAV*R_EARTH**3) * 1.d18
+ ! final unit : [s km^(-3) (kg/m^3)^(-1)]
+ scale_kl_rho = scale_t / scale_displ / RHOAV * 1.d9
+
+ ! allocates temporary arrays
+ allocate(cijkl_kl_crust_mantle_reg(21,NGLLX,NGLLY,NGLLZ,npoints_slice), &
+ rho_kl_crust_mantle_reg(NGLLX,NGLLY,NGLLZ,npoints_slice), &
+ beta_kl_crust_mantle_reg(NGLLX,NGLLY,NGLLZ,npoints_slice), &
+ alpha_kl_crust_mantle_reg(NGLLX,NGLLY,NGLLZ,npoints_slice))
+
+ if( SAVE_TRANSVERSE_KL ) then
+ ! transverse isotropic kernel arrays for file output
+ allocate(alphav_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,npoints_slice), &
+ alphah_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,npoints_slice), &
+ betav_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,npoints_slice), &
+ betah_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,npoints_slice), &
+ eta_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,npoints_slice))
+
+ ! isotropic kernel arrays for file output
+ allocate(bulk_c_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,npoints_slice), &
+ bulk_betav_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,npoints_slice), &
+ bulk_betah_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,npoints_slice), &
+ bulk_beta_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,npoints_slice))
+ endif
+
+ if( .not. ANISOTROPIC_KL ) then
+ ! allocates temporary isotropic kernel arrays for file output
+ allocate(bulk_c_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,npoints_slice), &
+ bulk_beta_kl_crust_mantle(NGLLX,NGLLY,NGLLZ,npoints_slice))
+ endif
+
+ ! crust_mantle
+ do ipoint = 1, npoints_slice
+ ispec = ispec_reg(ipoint)
+ do k = 1, NGLLZ
+ do j = 1, NGLLY
+ do i = 1, NGLLX
+
+ hlagrange = hxir_reg(i,ipoint)*hetar_reg(j,ipoint)*hgammar_reg(k,ipoint)
+
+ if (ANISOTROPIC_KL) then
+
+ ! For anisotropic kernels
+ iglob = ibool_crust_mantle(i,j,k,ispec)
+
+ ! The cartesian global cijkl_kl are rotated into the spherical local cijkl_kl
+ ! ystore and zstore are thetaval and phival (line 2252) -- dangerous
+ call rotate_kernels_dble(cijkl_kl_crust_mantle(:,i,j,k,ispec),cijkl_kl_local, &
+ ystore_crust_mantle(iglob),zstore_crust_mantle(iglob))
+
+ cijkl_kl_crust_mantle_reg(:,i,j,k,ipoint) = cijkl_kl_local * scale_kl_ani * hlagrange
+ rho_kl_crust_mantle_reg(i,j,k,ipoint) = rho_kl_crust_mantle(i,j,k,ispec) * scale_kl_rho * hlagrange
+
+ ! transverse isotropic kernel calculations
+ if( SAVE_TRANSVERSE_KL ) then
+ ! note: transverse isotropic kernels are calculated for all elements
+ !
+ ! however, the factors A,C,L,N,F are based only on transverse elements
+ ! in between Moho and 220 km layer, otherwise they are taken from isotropic values
+
+ rhol = rhostore_crust_mantle(i,j,k,ispec)
+
+ ! transverse isotropic parameters from compute_force_crust_mantle.f90
+ ! C=rhovpvsq A=rhovphsq L=rhovsvsq N=rhovshsq eta=F/(A - 2 L)
+
+ ! Get A,C,F,L,N,eta from kappa,mu
+ ! element can have transverse isotropy if between d220 and Moho
+ !if( .not. (TRANSVERSE_ISOTROPY_VAL .and. &
+ ! (idoubling_crust_mantle(ispec) == IFLAG_80_MOHO .or. &
+ ! idoubling_crust_mantle(ispec) == IFLAG_220_80))) then
+ if( .not. ispec_is_tiso_crust_mantle(ispec) ) then
+
+ ! layer with no transverse isotropy
+ ! A,C,L,N,F from isotropic model
+
+ mul = muvstore_crust_mantle(i,j,k,ispec)
+ kappal = kappavstore_crust_mantle(i,j,k,ispec)
+ muvl = mul
+ muhl = mul
+ kappavl = kappal
+ kappahl = kappal
+
+ A = kappal + FOUR_THIRDS * mul
+ C = A
+ L = mul
+ N = mul
+ F = kappal - 2._CUSTOM_REAL/3._CUSTOM_REAL * mul
+ eta = 1._CUSTOM_REAL
+
+ else
+
+ ! A,C,L,N,F from transverse isotropic model
+ kappavl = kappavstore_crust_mantle(i,j,k,ispec)
+ kappahl = kappahstore_crust_mantle(i,j,k,ispec)
+ muvl = muvstore_crust_mantle(i,j,k,ispec)
+ muhl = muhstore_crust_mantle(i,j,k,ispec)
+ kappal = kappavl
+
+ A = kappahl + FOUR_THIRDS * muhl
+ C = kappavl + FOUR_THIRDS * muvl
+ L = muvl
+ N = muhl
+ eta = eta_anisostore_crust_mantle(i,j,k,ispec) ! that is F / (A - 2 L)
+ F = eta * ( A - 2._CUSTOM_REAL * L )
+
+ endif
+
+ ! note: cijkl_kl_local() is fully anisotropic C_ij kernel components (non-dimensionalized)
+ ! for GLL point at (i,j,k,ispec)
+
+ ! Purpose : compute the kernels for the An coeffs (an_kl)
+ ! from the kernels for Cij (cijkl_kl_local)
+ ! At r,theta,phi fixed
+ ! kernel def : dx = kij * dcij + krho * drho
+ ! = kAn * dAn + krho * drho
+
+ ! Definition of the input array cij_kl :
+ ! cij_kl(1) = C11 ; cij_kl(2) = C12 ; cij_kl(3) = C13
+ ! cij_kl(4) = C14 ; cij_kl(5) = C15 ; cij_kl(6) = C16
+ ! cij_kl(7) = C22 ; cij_kl(8) = C23 ; cij_kl(9) = C24
+ ! cij_kl(10) = C25 ; cij_kl(11) = C26 ; cij_kl(12) = C33
+ ! cij_kl(13) = C34 ; cij_kl(14) = C35 ; cij_kl(15) = C36
+ ! cij_kl(16) = C44 ; cij_kl(17) = C45 ; cij_kl(18) = C46
+ ! cij_kl(19) = C55 ; cij_kl(20) = C56 ; cij_kl(21) = C66
+ ! where the Cij (Voigt's notation) are defined as function of
+ ! the components of the elastic tensor in spherical coordinates
+ ! by eq. (A.1) of Chen & Tromp, GJI 168 (2007)
+
+ ! From the relations giving Cij in function of An
+ ! Checked with Min Chen's results (routine build_cij)
+
+ an_kl(1) = cijkl_kl_local(1)+cijkl_kl_local(2)+cijkl_kl_local(7) !A
+ an_kl(2) = cijkl_kl_local(12) !C
+ an_kl(3) = -2*cijkl_kl_local(2)+cijkl_kl_local(21) !N
+ an_kl(4) = cijkl_kl_local(16)+cijkl_kl_local(19) !L
+ an_kl(5) = cijkl_kl_local(3)+cijkl_kl_local(8) !F
+
+ ! not used yet
+ !an_kl(6)=2*cijkl_kl_local(5)+2*cijkl_kl_local(10)+2*cijkl_kl_local(14) !Jc
+ !an_kl(7)=2*cijkl_kl_local(4)+2*cijkl_kl_local(9)+2*cijkl_kl_local(13) !Js
+ !an_kl(8)=-2*cijkl_kl_local(14) !Kc
+ !an_kl(9)=-2*cijkl_kl_local(13) !Ks
+ !an_kl(10)=-2*cijkl_kl_local(10)+cijkl_kl_local(18) !Mc
+ !an_kl(11)=2*cijkl_kl_local(4)-cijkl_kl_local(20) !Ms
+ !an_kl(12)=cijkl_kl_local(1)-cijkl_kl_local(7) !Bc
+ !an_kl(13)=-1./2.*(cijkl_kl_local(6)+cijkl_kl_local(11)) !Bs
+ !an_kl(14)=cijkl_kl_local(3)-cijkl_kl_local(8) !Hc
+ !an_kl(15)=-cijkl_kl_local(15) !Hs
+ !an_kl(16)=-cijkl_kl_local(16)+cijkl_kl_local(19) !Gc
+ !an_kl(17)=-cijkl_kl_local(17) !Gs
+ !an_kl(18)=cijkl_kl_local(5)-cijkl_kl_local(10)-cijkl_kl_local(18) !Dc
+ !an_kl(19)=cijkl_kl_local(4)-cijkl_kl_local(9)+cijkl_kl_local(20) !Ds
+ !an_kl(20)=cijkl_kl_local(1)-cijkl_kl_local(2)+cijkl_kl_local(7)-cijkl_kl_local(21) !Ec
+ !an_kl(21)=-cijkl_kl_local(6)+cijkl_kl_local(11) !Es
+
+ ! K_rho (primary kernel, for a parameterization (A,C,L,N,F,rho) )
+ rhonotprime_kl_crust_mantle(i,j,k,ipoint) = rhol * rho_kl_crust_mantle_reg(i,j,k,ipoint) / scale_kl_rho
+
+ ! note: transverse isotropic kernels are calculated for ALL elements,
+ ! and not just transverse elements
+ !
+ ! note: the kernels are for relative perturbations (delta ln (m_i) = (m_i - m_0)/m_i )
+ !
+ ! Gets transverse isotropic kernels
+ ! (see Appendix B of Sieminski et al., GJI 171, 2007)
+
+ ! for parameterization: ( alpha_v, alpha_h, beta_v, beta_h, eta, rho )
+ ! K_alpha_v
+ alphav_kl_crust_mantle(i,j,k,ipoint) = 2*C*an_kl(2) * hlagrange
+ ! K_alpha_h
+ alphah_kl_crust_mantle(i,j,k,ipoint) = (2*A*an_kl(1) + 2*A*eta*an_kl(5)) * hlagrange
+ ! K_beta_v
+ betav_kl_crust_mantle(i,j,k,ipoint) = (2*L*an_kl(4) - 4*L*eta*an_kl(5)) * hlagrange
+ ! K_beta_h
+ betah_kl_crust_mantle(i,j,k,ipoint) = 2*N*an_kl(3) * hlagrange
+ ! K_eta
+ eta_kl_crust_mantle(i,j,k,ipoint) = F*an_kl(5) * hlagrange
+ ! K_rhoprime (for a parameterization (alpha_v, ..., rho) )
+ rho_kl_crust_mantle_reg(i,j,k,ipoint) = (A*an_kl(1) + C*an_kl(2) &
+ + N*an_kl(3) + L*an_kl(4) + F*an_kl(5)) * hlagrange &
+ + rhonotprime_kl_crust_mantle(i,j,k,ipoint)
+
+ ! write the kernel in physical units (01/05/2006)
+ rhonotprime_kl_crust_mantle(i,j,k,ipoint) = - rhonotprime_kl_crust_mantle(i,j,k,ispec) * scale_kl * hlagrange
+
+ alphav_kl_crust_mantle(i,j,k,ipoint) = - alphav_kl_crust_mantle(i,j,k,ipoint) * scale_kl
+ alphah_kl_crust_mantle(i,j,k,ipoint) = - alphah_kl_crust_mantle(i,j,k,ipoint) * scale_kl
+ betav_kl_crust_mantle(i,j,k,ipoint) = - betav_kl_crust_mantle(i,j,k,ipoint) * scale_kl
+ betah_kl_crust_mantle(i,j,k,ipoint) = - betah_kl_crust_mantle(i,j,k,ipoint) * scale_kl
+ eta_kl_crust_mantle(i,j,k,ipoint) = - eta_kl_crust_mantle(i,j,k,ipoint) * scale_kl
+ rho_kl_crust_mantle_reg(i,j,k,ipoint) = - rho_kl_crust_mantle_reg(i,j,k,ipoint) * scale_kl
+
+ ! for parameterization: ( bulk, beta_v, beta_h, eta, rho )
+ ! where kappa_v = kappa_h = kappa and bulk c = sqrt( kappa / rho )
+ betav_sq = muvl / rhol
+ betah_sq = muhl / rhol
+ alphav_sq = ( kappal + FOUR_THIRDS * muvl ) / rhol
+ alphah_sq = ( kappal + FOUR_THIRDS * muhl ) / rhol
+ bulk_sq = kappal / rhol
+
+ bulk_c_kl_crust_mantle(i,j,k,ipoint) = &
+ bulk_sq / alphav_sq * alphav_kl_crust_mantle(i,j,k,ipoint) + &
+ bulk_sq / alphah_sq * alphah_kl_crust_mantle(i,j,k,ipoint)
+
+ bulk_betah_kl_crust_mantle(i,j,k,ipoint) = &
+ betah_kl_crust_mantle(i,j,k,ipoint) + &
+ FOUR_THIRDS * betah_sq / alphah_sq * alphah_kl_crust_mantle(i,j,k,ipoint)
+
+ bulk_betav_kl_crust_mantle(i,j,k,ipoint) = &
+ betav_kl_crust_mantle(i,j,k,ipoint) + &
+ FOUR_THIRDS * betav_sq / alphav_sq * alphav_kl_crust_mantle(i,j,k,ipoint)
+ ! the rest, K_eta and K_rho are the same as above
+
+ ! to check: isotropic kernels from transverse isotropic ones
+ alpha_kl_crust_mantle_reg(i,j,k,ipoint) = alphav_kl_crust_mantle(i,j,k,ipoint) &
+ + alphah_kl_crust_mantle(i,j,k,ipoint)
+ beta_kl_crust_mantle_reg(i,j,k,ipoint) = betav_kl_crust_mantle(i,j,k,ipoint) &
+ + betah_kl_crust_mantle(i,j,k,ipoint)
+ !rho_kl_crust_mantle_reg(i,j,k,ipoint) = rhonotprime_kl_crust_mantle(i,j,k,ipoint) &
+ ! + alpha_kl_crust_mantle_reg(i,j,k,ipoint) &
+ ! + beta_kl_crust_mantle_reg(i,j,k,ipoint)
+ bulk_beta_kl_crust_mantle(i,j,k,ipoint) = bulk_betah_kl_crust_mantle(i,j,k,ipoint) &
+ + bulk_betav_kl_crust_mantle(i,j,k,ipoint)
+
+ endif ! SAVE_TRANSVERSE_KL
+
+ else
+
+ ! isotropic kernels
+
+ rhol = rhostore_crust_mantle(i,j,k,ispec)
+ mul = muvstore_crust_mantle(i,j,k,ispec)
+ kappal = kappavstore_crust_mantle(i,j,k,ispec)
+
+ ! kernel values for rho, kappa, mu (primary kernel values)
+ rho_kl = - rhol * rho_kl_crust_mantle(i,j,k,ispec)
+ alpha_kl = - kappal * alpha_kl_crust_mantle(i,j,k,ispec) ! note: alpha_kl corresponds to K_kappa
+ beta_kl = - 2 * mul * beta_kl_crust_mantle(i,j,k,ispec) ! note: beta_kl corresponds to K_mu
+
+ ! for a parameterization: (rho,mu,kappa) "primary" kernels
+ rhonotprime_kl_crust_mantle(i,j,k,ipoint) = rho_kl * scale_kl * hlagrange
+ mu_kl_crust_mantle(i,j,k,ipoint) = beta_kl * scale_kl * hlagrange
+ kappa_kl_crust_mantle(i,j,k,ipoint) = alpha_kl * scale_kl * hlagrange
+
+ ! for a parameterization: (rho,alpha,beta)
+ ! kernels rho^prime, beta, alpha
+ rho_kl_crust_mantle_reg(i,j,k,ipoint) = (rho_kl + alpha_kl + beta_kl) * scale_kl * hlagrange
+ beta_kl_crust_mantle_reg(i,j,k,ipoint) = &
+ 2._CUSTOM_REAL * (beta_kl - FOUR_THIRDS * mul * alpha_kl / kappal) * scale_kl * hlagrange
+ alpha_kl_crust_mantle_reg(i,j,k,ipoint) = &
+ 2._CUSTOM_REAL * (1 + FOUR_THIRDS * mul / kappal) * alpha_kl * scale_kl * hlagrange
+
+ ! for a parameterization: (rho,bulk, beta)
+ ! where bulk wave speed is c = sqrt( kappa / rho)
+ ! note: rhoprime is the same as for (rho,alpha,beta) parameterization
+ bulk_c_kl_crust_mantle(i,j,k,ipoint) = 2._CUSTOM_REAL * alpha_kl * scale_kl * hlagrange
+ bulk_beta_kl_crust_mantle(i,j,k,ipoint) = 2._CUSTOM_REAL * beta_kl * scale_kl * hlagrange
+
+ endif
+
+ enddo
+ enddo
+ enddo
+ enddo
+
+ call create_name_database(prname,myrank,IREGION_CRUST_MANTLE,LOCAL_PATH)
+
+ ! For anisotropic kernels
+ if (ANISOTROPIC_KL) then
+
+ ! outputs transverse isotropic kernels only
+ if( SAVE_TRANSVERSE_KL ) then
+ ! transverse isotropic kernels
+ ! (alpha_v, alpha_h, beta_v, beta_h, eta, rho ) parameterization
+ open(unit=27,file=trim(prname)//'alphav_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) alphav_kl_crust_mantle
+ close(27)
+ open(unit=27,file=trim(prname)//'alphah_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) alphah_kl_crust_mantle
+ close(27)
+ open(unit=27,file=trim(prname)//'betav_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) betav_kl_crust_mantle
+ close(27)
+ open(unit=27,file=trim(prname)//'betah_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) betah_kl_crust_mantle
+ close(27)
+ open(unit=27,file=trim(prname)//'eta_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) eta_kl_crust_mantle
+ close(27)
+ open(unit=27,file=trim(prname)//'rho_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) rho_kl_crust_mantle_reg
+ close(27)
+
+ ! in case one is interested in primary kernel K_rho
+ !open(unit=27,file=trim(prname)//'rhonotprime_kernel.bin',status='unknown',form='unformatted',action='write')
+ !write(27) rhonotprime_kl_crust_mantle
+ !close(27)
+
+ ! (bulk, beta_v, beta_h, eta, rho ) parameterization: K_eta and K_rho same as above
+ open(unit=27,file=trim(prname)//'bulk_c_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) bulk_c_kl_crust_mantle
+ close(27)
+ open(unit=27,file=trim(prname)//'bulk_betav_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) bulk_betav_kl_crust_mantle
+ close(27)
+ open(unit=27,file=trim(prname)//'bulk_betah_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) bulk_betah_kl_crust_mantle
+ close(27)
+
+ ! to check: isotropic kernels
+ open(unit=27,file=trim(prname)//'alpha_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) alpha_kl_crust_mantle_reg
+ close(27)
+ open(unit=27,file=trim(prname)//'beta_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) beta_kl_crust_mantle_reg
+ close(27)
+ open(unit=27,file=trim(prname)//'bulk_beta_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) bulk_beta_kl_crust_mantle
+ close(27)
+
+ else
+
+ ! fully anisotropic kernels
+ ! note: the C_ij and density kernels are not for relative perturbations (delta ln( m_i) = delta m_i / m_i),
+ ! but absolute perturbations (delta m_i = m_i - m_0)
+ open(unit=27,file=trim(prname)//'rho_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) - rho_kl_crust_mantle_reg
+ close(27)
+ open(unit=27,file=trim(prname)//'cijkl_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) - cijkl_kl_crust_mantle_reg
+ close(27)
+
+ endif
+
+ else
+ ! primary kernels: (rho,kappa,mu) parameterization
+ open(unit=27,file=trim(prname)//'rhonotprime_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) rhonotprime_kl_crust_mantle
+ close(27)
+ open(unit=27,file=trim(prname)//'kappa_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) kappa_kl_crust_mantle
+ close(27)
+ open(unit=27,file=trim(prname)//'mu_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) mu_kl_crust_mantle
+ close(27)
+
+ ! (rho, alpha, beta ) parameterization
+ open(unit=27,file=trim(prname)//'rho_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) rho_kl_crust_mantle_reg
+ close(27)
+ open(unit=27,file=trim(prname)//'alpha_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) alpha_kl_crust_mantle_reg
+ close(27)
+ open(unit=27,file=trim(prname)//'beta_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) beta_kl_crust_mantle_reg
+ close(27)
+
+ ! (rho, bulk, beta ) parameterization, K_rho same as above
+ open(unit=27,file=trim(prname)//'bulk_c_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) bulk_c_kl_crust_mantle
+ close(27)
+ open(unit=27,file=trim(prname)//'bulk_beta_kernel.bin',status='unknown',form='unformatted',action='write')
+ write(27) bulk_beta_kl_crust_mantle
+ close(27)
+
+
+ endif
+
+ ! cleans up temporary kernel arrays
+ if( SAVE_TRANSVERSE_KL ) then
+ deallocate(alphav_kl_crust_mantle,alphah_kl_crust_mantle, &
+ betav_kl_crust_mantle,betah_kl_crust_mantle, &
+ eta_kl_crust_mantle)
+ deallocate(bulk_c_kl_crust_mantle,bulk_betah_kl_crust_mantle, &
+ bulk_betav_kl_crust_mantle,bulk_beta_kl_crust_mantle)
+ endif
+ if( .not. ANISOTROPIC_KL ) then
+ deallocate(bulk_c_kl_crust_mantle,bulk_beta_kl_crust_mantle)
+ endif
+ deallocate(cijkl_kl_crust_mantle_reg, &
+ rho_kl_crust_mantle_reg, &
+ beta_kl_crust_mantle_reg, &
+ alpha_kl_crust_mantle_reg)
+
+ end subroutine save_regular_kernels_crust_mantle
+
+
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90 2012-11-12 08:57:29 UTC (rev 21019)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/src/specfem3D/specfem3D.F90 2012-11-12 11:51:41 UTC (rev 21020)
@@ -671,6 +671,29 @@
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLY,NGLLZ,NSPEC_INNER_CORE_ADJOINT) :: rho_kl_inner_core, &
beta_kl_inner_core, alpha_kl_inner_core
+! For saving kernels on a regular grid
+ type kl_reg_grid_variables
+ sequence
+ real dlat
+ real dlon
+ integer nlayer
+ real rlayer(NM_KL_REG_LAYER)
+ integer ndoubling(NM_KL_REG_LAYER)
+ integer nlat(NM_KL_REG_LAYER)
+ integer nlon(NM_KL_REG_LAYER)
+ integer npts_total
+ integer npts_before_layer(NM_KL_REG_LAYER+1)
+ end type kl_reg_grid_variables
+ type (kl_reg_grid_variables) KL_REG_GRID
+
+ integer isp, npoints_slice
+ integer, dimension(:), allocatable :: slice_number
+ integer, dimension(NM_KL_REG_PTS) :: points_slice
+ integer, dimension(NM_KL_REG_PTS) :: ispec_reg
+ real, dimension(NGLLX, NM_KL_REG_PTS) :: hxir_reg
+ real, dimension(NGLLY, NM_KL_REG_PTS) :: hetar_reg
+ real, dimension(NGLLZ, NM_KL_REG_PTS) :: hgammar_reg
+
real(kind=CUSTOM_REAL), dimension(:,:,:,:,:), allocatable :: absorb_xmin_crust_mantle5, &
absorb_xmax_crust_mantle5, absorb_ymin_crust_mantle5, absorb_ymax_crust_mantle5
@@ -943,9 +966,10 @@
#ifdef USE_SERIAL_CASCADE_FOR_IOs
logical :: you_can_start_doing_IOs
+#endif
integer msg_status(MPI_STATUS_SIZE)
-#endif
+! *************************************************
! ************** PROGRAM STARTS HERE **************
!
!-------------------------------------------------------------------------------------------------
@@ -1317,6 +1341,63 @@
endif
+ if (SAVE_REGULAR_KL) then
+ call read_kl_regular_grid(KL_REG_GRID)
+
+ if (myrank==0) then
+ allocate(slice_number(KL_REG_GRID%npts_total))
+
+! print *, 'slice npts =', KL_REG_GRID%npts_total
+ call find_regular_grid_slice_number(slice_number, KL_REG_GRID, NCHUNKS_VAL, &
+ NPROC_XI_VAL, NPROC_ETA_VAL)
+ do i = NPROCTOT_VAL-1,0,-1
+ npoints_slice = 0
+ do isp = 1,KL_REG_GRID%npts_total
+ if (slice_number(isp) == i) then
+ npoints_slice = npoints_slice + 1
+ if (npoints_slice > NM_KL_REG_PTS) stop 'Exceeding NM_KL_REG_PTS limit'
+ points_slice(npoints_slice) = isp
+ endif
+ enddo
+
+ if (i /= 0) then
+ call MPI_Send(npoints_slice,1,MPI_INTEGER,i,i,MPI_COMM_WORLD,ier)
+ if (npoints_slice > 0) then
+ call MPI_Send(points_slice,npoints_slice,MPI_INTEGER,i,2*i,MPI_COMM_WORLD,ier)
+ endif
+ endif
+ enddo
+
+ open(unit=IOUT,file=trim(OUTPUT_FILES)//'/kl_grid_slice.txt',status='unknown',action='write')
+ write(IOUT,*) slice_number
+ close(IOUT)
+
+ deallocate(slice_number)
+ else
+ call MPI_Recv(npoints_slice,1,MPI_INTEGER,0,myrank,MPI_COMM_WORLD,msg_status,ier)
+ if (npoints_slice > 0) then
+ call MPI_Recv(points_slice,npoints_slice,MPI_INTEGER,0,2*myrank,MPI_COMM_WORLD,msg_status,ier)
+ endif
+ endif
+
+ ! this is the core part that takes up most of the computation time,
+ ! and presumably the more processors involved the faster.
+ if (npoints_slice > 0) then
+ call locate_reg_points(myrank, npoints_slice, points_slice, KL_REG_GRID, &
+ NEX_XI, NSPEC_CRUST_MANTLE, &
+ xstore_crust_mantle, ystore_crust_mantle, zstore_crust_mantle, &
+ ibool_crust_mantle, &
+ xigll, yigll, zigll, &
+ ispec_reg, hxir_reg, hetar_reg, hgammar_reg)
+ endif
+
+ if (myrank==0) then
+ write(IMAIN,*) ' '
+ write(IMAIN,*) 'Finish locating kernel output regular grid'
+ write(IMAIN,*) ' '
+ endif
+ endif
+
#ifdef USE_SERIAL_CASCADE_FOR_IOs
you_can_start_doing_IOs = .true.
if (myrank < NPROC_XI_VAL*NPROC_ETA_VAL-1) &
@@ -4671,7 +4752,22 @@
! dump kernel arrays
if (SIMULATION_TYPE == 3) then
+
! crust mantle
+ if (SAVE_REGULAR_KL) then
+ call save_regular_kernels_crust_mantle(myrank, &
+ KL_REG_GRID, npoints_slice, hxir_reg, hetar_reg, hgammar_reg, &
+ scale_t,scale_displ, &
+ cijkl_kl_crust_mantle,rho_kl_crust_mantle, &
+ alpha_kl_crust_mantle,beta_kl_crust_mantle, &
+ ystore_crust_mantle,zstore_crust_mantle, &
+ rhostore_crust_mantle,muvstore_crust_mantle, &
+ kappavstore_crust_mantle,ibool_crust_mantle, &
+ kappahstore_crust_mantle,muhstore_crust_mantle, &
+ eta_anisostore_crust_mantle,ispec_is_tiso_crust_mantle, &
+ ! --idoubling_crust_mantle, &
+ LOCAL_PATH)
+ else
call save_kernels_crust_mantle(myrank,scale_t,scale_displ, &
cijkl_kl_crust_mantle,rho_kl_crust_mantle, &
alpha_kl_crust_mantle,beta_kl_crust_mantle, &
@@ -4682,6 +4778,7 @@
eta_anisostore_crust_mantle,ispec_is_tiso_crust_mantle, &
! --idoubling_crust_mantle, &
LOCAL_PATH)
+ endif
!<YANGL
! noise strength kernel
More information about the CIG-COMMITS
mailing list