[cig-commits] commit: First pass at a helper package for custom Earth models.
Mercurial
hg at geodynamics.org
Sun Jul 3 20:06:52 PDT 2011
changeset: 0:9b613f27327a
user: Walter Landry <wlandry at caltech.edu>
date: Sat Jul 02 17:09:21 2011 -0700
files: DATA/PPM/readme.txt README create_tgz src/meshfem3D/model_aniso_mantle.f90 src/meshfem3D/model_atten3D_QRFSI12.f90 src/meshfem3D/model_crust.f90 src/meshfem3D/model_s20rts.f90 test_compilation
description:
First pass at a helper package for custom Earth models.
diff -r 000000000000 -r 9b613f27327a DATA/PPM/readme.txt
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/DATA/PPM/readme.txt Sat Jul 02 17:09:21 2011 -0700
@@ -0,0 +1,31 @@
+----------------------------------------------------------------------
+PPM - point-profile models
+----------------------------------------------------------------------
+
+ for generic Vs models given as depth profiles at lon/lat using a text-file format like:
+
+ #lon(deg), lat(deg), depth(km), Vs-perturbation wrt PREM(%), Vs-PREM (km/s)
+ -10.00000 31.00000 40.00000 -1.775005 4.400000
+ -10.00000 32.00000 40.00000 -1.056823 4.400000
+ ...
+
+ (the first line in the file is ignored)
+
+
+1. in order to use it, create a symbolic link 'model.txt'
+ in this directory DATA/PPM/:
+ > cd DATA/PPM/
+ > ln -s my_ppm_model model.txt
+
+
+2. change in DATA/Par_file the corresponding line to:
+ #..
+ MODEL = PPM
+ ..
+
+3. any specifics can be added to the routine implementations in
+ the subroutine file: 'model_ppm.f90'
+
+ note: by default, the routine assumes that data given are Vs perturbations
+ and that the data have a regular spacing, i.e. constant increments,
+ for each lon,lat and depth column
\ No newline at end of file
diff -r 000000000000 -r 9b613f27327a README
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/README Sat Jul 02 17:09:21 2011 -0700
@@ -0,0 +1,35 @@
+These scripts assume a Unix-like environment: Linux, Mac OS X, or
+Windows with Cygwin.
+
+To use a custom Earth model, you must either
+
+1) Modify one of the existing Earth models. This involves modifying
+ one of the f90 files in src/meshfem3D. The exact details of what
+ must be done are spelled in Chapter 10 of the manual.
+
+ https://geodynamics.org/portals/seismo/doc/manual_SPECFEM3D_GLOBE.pdf
+
+ Do not forget to modify src/meshfem3D/Makefile if you add any extra
+ source files.
+
+2) Specify your model as a deviation from PREM. See the file
+ DATA/PPM/readme.txt for instructions.
+
+Once you have completed your modifications, you should do a simple
+test to make sure the code compiles. To do this, download a copy
+of Specfem 3D Globe from
+
+ http://geodynamics.org/cig/software/specfem3d-globe/SPECFEM3D_GLOBE_V5.1.1.tar.gz
+
+Configure and build Specfem 3D Globe. Then run
+
+ ./test_compilation /path/to/SPECFEM3D_GLOBE
+
+You should also make sure your modified version gives correct answers.
+
+To create a file suitable for uploading, write a description of the
+Earth model in the file "description.txt" or "description.html".
+Then, run
+
+ ./create_tgz Backend_Name
+
diff -r 000000000000 -r 9b613f27327a create_tgz
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/create_tgz Sat Jul 02 17:09:21 2011 -0700
@@ -0,0 +1,21 @@
+#!/bin/bash
+
+if [ $# -ne 1 ]; then
+ echo "Usage: create_tgz Backend_Name"
+ exit 1
+fi
+
+if [ ! -s description.txt ]; then
+ if [ ! -s description.html ]; then
+ echo "You must put something into description.txt or description.html file"
+ exit 1
+ fi
+fi
+
+echo "Creating the tarfile ${1}.tgz"
+echo ""
+
+tar -zcvf "${1}.tgz" DATA/ src/ description.*[a-z]
+
+echo ""
+echo "Now upload the file ${1}.tgz to the portal."
diff -r 000000000000 -r 9b613f27327a src/meshfem3D/model_aniso_mantle.f90
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/meshfem3D/model_aniso_mantle.f90 Sat Jul 02 17:09:21 2011 -0700
@@ -0,0 +1,907 @@
+!=====================================================================
+!
+! 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 University of Pau / CNRS / INRIA, France
+! (c) Princeton University / California Institute of Technology and University of Pau / CNRS / INRIA
+! 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.
+!
+!=====================================================================
+
+!--------------------------------------------------------------------------------------------------
+!
+! Jean-Paul Montagner, January 2002
+! modified by Min Chen, Caltech, February 2002
+!
+! input is (r, theta, phi), output is the matrix cij(6x6)
+! 0 <= r <= 1, 0 <= theta <= pi, 0 <= phi <= 2 pi
+!
+! returns non-dimensionalized cij
+!
+! creates parameters p(i=1,14,r,theta,phi)
+! from model glob-prem3sm01 globpreman3sm01 (Montagner, 2002)
+!
+!--------------------------------------------------------------------------------------------------
+
+
+ subroutine model_aniso_mantle_broadcast(myrank,AMM_V)
+
+! standard routine to setup model
+
+ implicit none
+
+ include "constants.h"
+ ! standard include of the MPI library
+ include 'mpif.h'
+
+ ! model_aniso_mantle_variables
+ type model_aniso_mantle_variables
+ sequence
+ double precision beta(14,34,37,73)
+ double precision pro(47)
+ integer npar1
+ integer dummy_pad ! padding 4 bytes to align the structure
+ end type model_aniso_mantle_variables
+
+ type (model_aniso_mantle_variables) AMM_V
+ ! model_aniso_mantle_variables
+
+ integer :: myrank
+ integer :: ier
+
+ ! the variables read are declared and stored in structure AMM_V
+ if(myrank == 0) call read_aniso_mantle_model(AMM_V)
+
+ ! broadcast the information read on the master to the nodes
+ call MPI_BCAST(AMM_V%npar1,1,MPI_INTEGER,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST(AMM_V%beta,14*34*37*73,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST(AMM_V%pro,47,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+
+
+ end subroutine model_aniso_mantle_broadcast
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+
+ subroutine model_aniso_mantle(r,theta,phi,rho, &
+ c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66,&
+ AMM_V)
+
+ implicit none
+
+ include "constants.h"
+
+! model_aniso_mantle_variables
+ type model_aniso_mantle_variables
+ sequence
+ double precision beta(14,34,37,73)
+ double precision pro(47)
+ integer npar1
+ integer dummy_pad ! padding 4 bytes to align the structure
+ end type model_aniso_mantle_variables
+
+ type (model_aniso_mantle_variables) AMM_V
+! model_aniso_mantle_variables
+
+ double precision r,theta,phi
+ double precision rho
+ double precision c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26, &
+ c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
+
+ double precision d11,d12,d13,d14,d15,d16,d22,d23,d24,d25,d26, &
+ d33,d34,d35,d36,d44,d45,d46,d55,d56,d66
+
+ double precision colat,lon
+
+ lon = phi / DEGREES_TO_RADIANS
+ colat = theta / DEGREES_TO_RADIANS
+
+! uncomment this line to suppress the anisotropic mantle model
+! call exit_MPI_without_rank('please provide an anisotropic mantle model for subroutine aniso_mantle_model')
+
+! assign the local (d_ij) or global (c_ij) anisotropic parameters.
+! The c_ij are the coefficients in the global
+! reference frame used in SPECFEM3D.
+ call build_cij(AMM_V%pro,AMM_V%npar1,rho,AMM_V%beta,r,colat,lon,&
+ d11,d12,d13,d14,d15,d16,d22,d23,d24,d25,d26,d33,d34,d35,d36,&
+ d44,d45,d46,d55,d56,d66)
+
+ call rotate_aniso_tensor(theta,phi,d11,d12,d13,d14,d15,d16,d22,d23,d24,d25,d26,&
+ d33,d34,d35,d36,d44,d45,d46,d55,d56,d66,&
+ c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,&
+ c33,c34,c35,c36,c44,c45,c46,c55,c56,c66)
+
+ end subroutine model_aniso_mantle
+
+!--------------------------------------------------------------------
+
+ subroutine build_cij(pro,npar1,rho,beta,r,theta,phi,&
+ d11,d12,d13,d14,d15,d16,d22,d23,d24,d25,d26,d33,d34,d35,d36,&
+ d44,d45,d46,d55,d56,d66)
+
+ implicit none
+
+ include "constants.h"
+
+ integer npar1,ndepth,idep,ipar,itheta,ilon,icz0,nx0,ny0,nz0,&
+ ict0,ict1,icp0,icp1,icz1
+
+ double precision d11,d12,d13,d14,d15,d16,d22,d23,d24,d25,d26, &
+ d33,d34,d35,d36,d44,d45,d46,d55,d56,d66
+ double precision r,theta,phi,rho,depth,tei,tet,ph,fi,x0,y0,pxy0
+ double precision d1,d2,d3,d4,sd,thickness,dprof1,dprof2,eps,pc1,pc2,pc3,pc4,&
+ dpr1,dpr2,param,scale_GPa,scaleval
+ double precision A,C,F,AL,AN,BC,BS,GC,GS,HC,HS,EC,ES,C1p,C1sv,C1sh,C3,S1p,S1sv,S1sh,S3
+ double precision beta(14,34,37,73),pro(47)
+ double precision anispara(14,2,4),elpar(14)
+
+ ndepth = npar1
+ pxy0 = 5.
+ x0 = 0.
+ y0 = 0.
+ nx0 = 37
+ ny0 = 73
+ nz0 = 34
+
+! avoid edge effects
+ if(theta==0.0d0) theta=0.000001d0
+ if(theta==180.d0) theta=0.999999d0*theta
+ if(phi==0.0d0) phi=0.000001d0
+ if(phi==360.d0) phi=0.999999d0*phi
+
+! dimensionalize
+ depth = R_EARTH_KM*(R_UNIT_SPHERE - r)
+ if(depth <= pro(nz0) .or. depth >= pro(1)) call exit_MPI_without_rank('r out of range in build_cij')
+ itheta = int(theta + pxy0)/pxy0
+ ilon = int(phi + pxy0)/pxy0
+ tet = theta
+ ph = phi
+
+ icz0 = 0
+ do idep = 1,ndepth
+ if(pro(idep) > depth) icz0 = icz0 + 1
+ enddo
+
+!
+! Interpolation for depth between dep1(iz0) and dep2(iz1)
+!
+! 1 (ict0,icp0) 2 (ict0,icp1)
+! 3 (ict1,icp0) 4 (ict1,icp1)
+!
+
+ ict0 = itheta
+ ict1 = ict0 + 1
+ icp0 = ilon
+ icp1 = icp0 + 1
+ icz1 = icz0 + 1
+
+! check that parameters make sense
+ if(ict0 < 1 .or. ict0 > nx0) call exit_MPI_without_rank('ict0 out of range')
+ if(ict1 < 1 .or. ict1 > nx0) call exit_MPI_without_rank('ict1 out of range')
+ if(icp0 < 1 .or. icp0 > ny0) call exit_MPI_without_rank('icp0 out of range')
+ if(icp1 < 1 .or. icp1 > ny0) call exit_MPI_without_rank('icp1 out of range')
+ if(icz0 < 1 .or. icz0 > nz0) call exit_MPI_without_rank('icz0 out of range')
+ if(icz1 < 1 .or. icz1 > nz0) call exit_MPI_without_rank('icz1 out of range')
+
+ do ipar = 1,14
+ anispara(ipar,1,1) = beta(ipar,icz0,ict0,icp0)
+ anispara(ipar,2,1) = beta(ipar,icz1,ict0,icp0)
+ anispara(ipar,1,2) = beta(ipar,icz0,ict0,icp1)
+ anispara(ipar,2,2) = beta(ipar,icz1,ict0,icp1)
+ anispara(ipar,1,3) = beta(ipar,icz0,ict1,icp0)
+ anispara(ipar,2,3) = beta(ipar,icz1,ict1,icp0)
+ anispara(ipar,1,4) = beta(ipar,icz0,ict1,icp1)
+ anispara(ipar,2,4) = beta(ipar,icz1,ict1,icp1)
+ enddo
+
+!
+! calculation of distances between the selected point and grid points
+!
+ tei = pxy0*ict0 + x0 - pxy0
+ fi = pxy0*icp0 + y0 - pxy0
+
+!*** d1=de(tet,ph,tei,fi)
+
+ d1 = dsqrt(((tei - tet)**2) + ((fi - ph)**2)*(dsin((tet + tei)*DEGREES_TO_RADIANS/2.)**2))
+
+!*** d2=de(tet,ph,tei+pxy0,fi)
+
+ d2 = dsqrt(((tei - tet + pxy0)**2) + ((fi - ph)**2)*(dsin((tet + tei + pxy0)*DEGREES_TO_RADIANS/2.)**2))
+
+!*** d3=de(tet,ph,tei,fi+pxy0)
+
+ d3 = dsqrt(((tei - tet)**2) + ((fi - ph + pxy0)**2)*(dsin((tet + tei)*DEGREES_TO_RADIANS/2.)**2))
+
+!*** d4=de(tet,ph,tei+pxy0,fi+pxy0)
+
+ d4 = dsqrt(((tei - tet + pxy0)**2) + ((fi - ph + pxy0)**2)*(dsin((tet + tei + pxy0)*DEGREES_TO_RADIANS/2.)**2))
+
+ sd = d2*d3*d4 + d1*d2*d4 + d1*d3*d4 + d1*d2*d3
+ thickness = pro(icz0) - pro(icz1)
+ dprof1 = pro(icz0) - depth
+ dprof2 = depth - pro(icz1)
+ eps = 0.01
+
+ do ipar = 1,14
+ if(thickness < eps)then
+ pc1 = anispara(ipar,1,1)
+ pc2 = anispara(ipar,1,2)
+ pc3 = anispara(ipar,1,3)
+ pc4 = anispara(ipar,1,4)
+ else
+ dpr1 = dprof1/thickness
+ dpr2 = dprof2/thickness
+ pc1 = anispara(ipar,1,1)*dpr2+anispara(ipar,2,1)*dpr1
+ pc2 = anispara(ipar,1,2)*dpr2+anispara(ipar,2,2)*dpr1
+ pc3 = anispara(ipar,1,3)*dpr2+anispara(ipar,2,3)*dpr1
+ pc4 = anispara(ipar,1,4)*dpr2+anispara(ipar,2,4)*dpr1
+ endif
+ param = pc1*d2*d3*d4 + pc2*d1*d3*d4 + pc3*d1*d2*d4 + pc4*d1*d2*d3
+ param = param/sd
+ elpar(ipar) = param
+ enddo
+
+ d11 = ZERO
+ d12 = ZERO
+ d13 = ZERO
+ d14 = ZERO
+ d15 = ZERO
+ d16 = ZERO
+ d22 = ZERO
+ d23 = ZERO
+ d24 = ZERO
+ d25 = ZERO
+ d26 = ZERO
+ d33 = ZERO
+ d34 = ZERO
+ d35 = ZERO
+ d36 = ZERO
+ d44 = ZERO
+ d45 = ZERO
+ d46 = ZERO
+ d55 = ZERO
+ d56 = ZERO
+ d66 = ZERO
+!
+! create dij
+!
+ rho = elpar(1)
+ A = elpar(2)
+ C = elpar(3)
+ F = elpar(4)
+ AL = elpar(5)
+ AN = elpar(6)
+ BC = elpar(7)
+ BS = elpar(8)
+ GC = elpar(9)
+ GS = elpar(10)
+ HC = elpar(11)
+ HS = elpar(12)
+ EC = elpar(13)
+ ES = elpar(14)
+ C1p = 0.0d0
+ S1p = 0.0d0
+ C1sv = 0.0d0
+ S1sv = 0.0d0
+ C1sh = 0.0d0
+ S1sh = 0.0d0
+ C3 = 0.0d0
+ S3 = 0.0d0
+
+ d11 = A + EC + BC
+ d12 = A - 2.*AN - EC
+ d13 = F + HC
+ d14 = S3 + 2.*S1sh + 2.*S1p
+ d15 = 2.*C1p + C3
+ d16 = -BS/2. - ES
+ d22 = A + EC - BC
+ d23 = F - HC
+ d24 = 2.*S1p - S3
+ d25 = 2.*C1p - 2.*C1sh - C3
+ d26 = -BS/2. + ES
+ d33 = C
+ d34 = 2.*(S1p - S1sv)
+ d35 = 2.*(C1p - C1sv)
+ d36 = -HS
+ d44 = AL - GC
+ d45 = -GS
+ d46 = C1sh - C3
+ d55 = AL + GC
+ d56 = S3 - S1sh
+ d66 = AN - EC
+
+! non-dimensionalize the elastic coefficients using
+! the scale of GPa--[g/cm^3][(km/s)^2]
+ scaleval = dsqrt(PI*GRAV*RHOAV)
+ scale_GPa =(RHOAV/1000.d0)*((R_EARTH*scaleval/1000.d0)**2)
+ d11 = d11/scale_GPa
+ d12 = d12/scale_GPa
+ d13 = d13/scale_GPa
+ d14 = d14/scale_GPa
+ d15 = d15/scale_GPa
+ d16 = d16/scale_GPa
+ d22 = d22/scale_GPa
+ d23 = d23/scale_GPa
+ d24 = d24/scale_GPa
+ d25 = d25/scale_GPa
+ d26 = d26/scale_GPa
+ d33 = d33/scale_GPa
+ d34 = d34/scale_GPa
+ d35 = d35/scale_GPa
+ d36 = d36/scale_GPa
+ d44 = d44/scale_GPa
+ d45 = d45/scale_GPa
+ d46 = d46/scale_GPa
+ d55 = d55/scale_GPa
+ d56 = d56/scale_GPa
+ d66 = d66/scale_GPa
+
+! non-dimensionalize
+ rho = rho*1000.d0/RHOAV
+
+ end subroutine build_cij
+
+!--------------------------------------------------------------
+
+ subroutine read_aniso_mantle_model(AMM_V)
+
+ implicit none
+
+ include "constants.h"
+
+! model_aniso_mantle_variables
+ type model_aniso_mantle_variables
+ sequence
+ double precision beta(14,34,37,73)
+ double precision pro(47)
+ integer npar1
+ integer dummy_pad ! padding 4 bytes to align the structure
+ end type model_aniso_mantle_variables
+
+ type (model_aniso_mantle_variables) AMM_V
+! model_aniso_mantle_variables
+
+ integer nx,ny,np1,np2,ipar,ipa1,ipa,ilat,ilon,il,idep,nfin,nfi0,nf,nri
+ double precision xinf,yinf,pxy,ppp,angle,A,A2L,AL,af
+ double precision ra(47),pari(14,47)
+ double precision bet2(14,34,37,73)
+ double precision alph(73,37),ph(73,37)
+ character(len=150) glob_prem3sm01, globpreman3sm01
+
+ np1 = 1
+ np2 = 34
+ AMM_V%npar1 = (np2 - np1 + 1)
+
+!
+! glob-prem3sm01: model with rho,A,L,xi-1,1-phi,eta
+!
+ call get_value_string(glob_prem3sm01, 'model.glob_prem3sm01', 'DATA/Montagner_model/glob-prem3sm01')
+ open(19,file=glob_prem3sm01,status='old',action='read')
+
+!
+! read the models
+!
+! reference model: PREM or ACY400
+!
+ call lecmod(nri,pari,ra)
+!
+! read tomographic model (equivalent T.I. model)
+!
+ ipa = 0
+ nfi0 = 6
+ nfin = 14
+ do nf = 1,nfi0
+ ipa = ipa + 1
+ do idep = 1,AMM_V%npar1
+ il = idep + np1 - 1
+ read(19,"(2f4.0,2i3,f4.0)",end = 88) xinf,yinf,nx,ny,pxy
+
+ ppp = 1.
+ read(19,"(f5.0,f8.4)",end = 88) AMM_V%pro(idep),ppp
+
+ if(nf == 1) pari(nf,il) = ppp
+ if(nf == 2) pari(nf,il) = ppp
+ if(nf == 3) pari(nf,il) = ppp
+ if(nf == 4) ppp = pari(nf,il)
+ if(nf == 5) ppp = pari(nf,il)
+ do ilat = 1,nx
+ read(19,"(17f7.2)",end = 88) (AMM_V%beta(ipa,idep,ilat,ilon),ilon = 1,ny)
+!
+! calculation of A,C,F,L,N
+!
+! bet2(1,...)=rho, bet2(2,...)=A,bet2(3,...)=L,bet2(4,...)=xi
+! bet2(5,...)=phi=C/A, bet2(6,...)=eta=F/A-2L
+! bet2(7,...)=Bc, bet2(8,...)=Bs,bet2(9,...)=Gc,bet2(10,...)=Gs
+! bet2(11,...)=Hc, bet2(12,...)=Hs,bet2(13,...)=Ec,bet2(14,...)=Es
+!
+ do ilon = 1,ny
+ if(nf <= 3 .or. nf >= 6)then
+ bet2(ipa,idep,ilat,ilon) = AMM_V%beta(ipa,idep,ilat,ilon)*0.01*ppp + ppp
+ else
+ if(nf == 4)bet2(ipa,idep,ilat,ilon) = AMM_V%beta(ipa,idep,ilat,ilon)*0.01 + 1.
+ if(nf == 5)bet2(ipa,idep,ilat,ilon) = - AMM_V%beta(ipa,idep,ilat,ilon)*0.01 + 1.
+ endif
+ enddo
+
+ enddo
+ enddo
+ enddo
+88 close(19)
+
+!
+! read anisotropic azimuthal parameters
+!
+
+!
+! beta(ipa,idep,ilat,ilon) are sorted in (amplitude, phase)
+! normalized, in percents: 100 G/L
+!
+ call get_value_string(globpreman3sm01, 'model.globpreman3sm01', 'DATA/Montagner_model/globpreman3sm01')
+ open(unit=15,file=globpreman3sm01,status='old',action='read')
+
+ do nf = 7,nfin,2
+ ipa = nf
+ ipa1 = ipa + 1
+ do idep = 1,AMM_V%npar1
+ il = idep + np1 - 1
+ read(15,"(2f4.0,2i3,f4.0)",end = 888) xinf,yinf,nx,ny,pxy
+ read(15,"(f5.0,f8.4)",end = 888) AMM_V%pro(idep),ppp
+ if(nf == 7) ppp = pari(2,il)
+ if(nf == 9) ppp = pari(3,il)
+ af = pari(6,il)*(pari(2,il) - 2.*pari(3,il))
+ if(nf == 11) ppp = af
+ if(nf == 13) ppp = (pari(4,il) + 1.)*pari(3,il)
+
+ do ilat = 1,nx
+ read(15,"(17f7.2)",end = 888) (alph(ilon,ilat),ilon = 1,ny)
+ enddo
+
+ do ilat=1,nx
+ read(15,"(17f7.2)",end = 888) (ph(ilon,ilat),ilon = 1,ny)
+ enddo
+
+ do ilat = 1,nx
+ do ilon = 1,ny
+ angle = 2.*DEGREES_TO_RADIANS*ph(ilon,ilat)
+ AMM_V%beta(ipa,idep,ilat,ilon) = alph(ilon,ilat)*ppp*0.01d0
+ AMM_V%beta(ipa1,idep,ilat,ilon) = ph(ilon,ilat)
+ bet2(ipa,idep,ilat,ilon) = alph(ilon,ilat)*dcos(angle)*ppp*0.01d0
+ bet2(ipa1,idep,ilat,ilon) = alph(ilon,ilat)*dsin(angle)*ppp*0.01d0
+ enddo
+ enddo
+
+ enddo
+ enddo
+
+888 close(15)
+
+ do idep = 1,AMM_V%npar1
+ do ilat = 1,nx
+ do ilon = 1,ny
+
+! rho
+ AMM_V%beta(1,idep,ilat,ilon) = bet2(1,idep,ilat,ilon)
+
+! A
+ AMM_V%beta(2,idep,ilat,ilon) = bet2(2,idep,ilat,ilon)
+ A=bet2(2,idep,ilat,ilon)
+
+! C
+ AMM_V%beta(3,idep,ilat,ilon) = bet2(5,idep,ilat,ilon)*A
+
+! F
+ A2L = A - 2.*bet2(3,idep,ilat,ilon)
+ AMM_V%beta(4,idep,ilat,ilon) = bet2(6,idep,ilat,ilon)*A2L
+
+! L
+ AMM_V%beta(5,idep,ilat,ilon) = bet2(3,idep,ilat,ilon)
+ AL = bet2(3,idep,ilat,ilon)
+
+! N
+ AMM_V%beta(6,idep,ilat,ilon) = bet2(4,idep,ilat,ilon)*AL
+
+! azimuthal terms
+ do ipar = 7,14
+ AMM_V%beta(ipar,idep,ilat,ilon) = bet2(ipar,idep,ilat,ilon)
+ enddo
+
+ enddo
+ enddo
+ enddo
+
+ end subroutine read_aniso_mantle_model
+
+!--------------------------------------------------------------------
+
+ subroutine lecmod(nri,pari,ra)
+
+ implicit none
+
+! read the reference Earth model: rho, Vph, Vsv, XI, PHI, ETA
+! array par(i,nlayer)
+! output: array pari(ipar, nlayer): rho, A, L, xi-1, phi-1, eta-1
+
+ integer i,j,k,ip,ifanis,idum1,idum2,idum3,nlayer,nout,neff,&
+ nband,nri,minlay,moho,kiti
+ double precision pari(14,47),qkappa(47),qshear(47),par(6,47)
+ double precision epa(14,47),ra(47),dcori(47),ri(47)
+ double precision corpar(21,47)
+ double precision aa,an,al,af,ac,vpv,vph,vsv,vsh,rho,red,a2l
+ character(len=80) null
+ character(len=150) Adrem119
+
+ ifanis = 1
+ nri = 47
+
+ call get_value_string(Adrem119, 'model.Adrem119', 'DATA/Montagner_model/Adrem119')
+ open(unit=13,file=Adrem119,status='old',action='read')
+ read(13,*,end = 77) nlayer,minlay,moho,nout,neff,nband,kiti,null
+
+ if(kiti == 0) read(13,"(20a4)",end = 77) idum1
+ read(13,"(20a4)",end = 77) idum2
+ read(13,"(20a4)",end = 77) idum3
+
+ do i = 1,nlayer
+ read(13,"(4x,f11.1,8d12.5)",end = 77) ra(i),(par(k,i),k = 1,6),qshear(i),qkappa(i)
+ enddo
+
+ do i = 1,nlayer
+ ri(i) = 0.001*ra(i)
+ enddo
+
+ do i = 1,nlayer
+ rho = par(1,i)
+ pari(1,i) = rho
+! A : pari(2,i)
+ pari(2,i) = rho*(par(2,i)**2)
+ aa = pari(2,i)
+! L : pari(3,i)
+ pari(3,i) = rho*(par(3,i)**2)
+ al = pari(3,i)
+! Xi : pari(4,i)= (N-L)/L
+ an = al*par(4,i)
+ pari(4,i) = 0.
+ pari(4,i) = par(4,i) - 1.
+! Phi : pari(5,i)=(a-c)/a
+ pari(5,i) = - par(5,i) + 1.
+ ac = par(5,i)*aa
+! f : pari(4,i)
+ af = par(6,i)*(aa - 2.*al)
+ pari(6,i) = par(6,i)
+ do ip = 7,14
+ pari(ip,i) = 0.
+ enddo
+ vsv = 0.
+ vsh = 0.
+ if(al < 0.0001 .or. an < 0.0001) goto 12
+ vsv = dsqrt(al/rho)
+ vsh = dsqrt(an/rho)
+ 12 vpv = dsqrt(ac/rho)
+ vph = dsqrt(aa/rho)
+ enddo
+
+ red = 1.
+ do i = 1,nlayer
+ read(13,"(15x,6e12.5,f11.1)",end = 77) (epa(j,i),j = 1,6),dcori(i)
+ epa(7,i) = epa(2,i)
+ epa(8,i) = epa(2,i)
+ epa(9,i) = epa(3,i)
+ epa(10,i) = epa(3,i)
+
+ a2l = pari(2,i) - 2.*pari(3,i)
+ epa(11,i) = epa(6,i)*a2l
+ epa(12,i) = epa(6,i)*a2l
+ epa(13,i) = epa(3,i)
+ epa(14,i) = epa(3,i)
+
+ do j = 1,14
+ epa(j,i) = red*epa(j,i)
+ enddo
+
+ read(13,"(21f7.3)",end = 77) (corpar(j,i),j = 1,21)
+
+ enddo
+
+77 close(13)
+
+ end subroutine lecmod
+
+!--------------------------------------------------------------------
+
+ subroutine rotate_aniso_tensor(theta,phi,d11,d12,d13,d14,d15,d16,&
+ d22,d23,d24,d25,d26,&
+ d33,d34,d35,d36,d44,d45,d46,d55,d56,d66,&
+ c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,&
+ c33,c34,c35,c36,c44,c45,c46,c55,c56,c66)
+
+ implicit none
+
+ include "constants.h"
+
+ double precision theta,phi
+ double precision c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26, &
+ c33,c34,c35,c36,c44,c45,c46,c55,c56,c66
+ double precision d11,d12,d13,d14,d15,d16,d22,d23,d24,d25,d26, &
+ d33,d34,d35,d36,d44,d45,d46,d55,d56,d66
+ double precision costheta,sintheta,cosphi,sinphi
+ double precision costhetasq,sinthetasq,cosphisq,sinphisq
+ double precision costwotheta,sintwotheta,costwophi,sintwophi
+ double precision cosfourtheta,sinfourtheta
+ double precision costhetafour,sinthetafour,cosphifour,sinphifour
+ double precision sintwophisq,sintwothetasq
+
+ costheta = dcos(theta)
+ sintheta = dsin(theta)
+ cosphi = dcos(phi)
+ sinphi = dsin(phi)
+
+ costhetasq = costheta * costheta
+ sinthetasq = sintheta * sintheta
+ cosphisq = cosphi * cosphi
+ sinphisq = sinphi * sinphi
+
+ costhetafour = costhetasq * costhetasq
+ sinthetafour = sinthetasq * sinthetasq
+ cosphifour = cosphisq * cosphisq
+ sinphifour = sinphisq * sinphisq
+
+ costwotheta = dcos(2.d0*theta)
+ sintwotheta = dsin(2.d0*theta)
+ costwophi = dcos(2.d0*phi)
+ sintwophi = dsin(2.d0*phi)
+
+ cosfourtheta = dcos(4.d0*theta)
+ sinfourtheta = dsin(4.d0*theta)
+ sintwothetasq = sintwotheta * sintwotheta
+ sintwophisq = sintwophi * sintwophi
+
+! recompute 21 anisotropic coefficients for full anisotropoc model using Mathematica
+
+c11 = d22*sinphifour - 2.*sintwophi*sinphisq*(d26*costheta + d24*sintheta) - &
+ 2.*cosphisq*sintwophi*(d16*costhetasq*costheta + &
+ (d14 + 2*d56)*costhetasq*sintheta + &
+ (d36 + 2*d45)*costheta*sinthetasq + d34*sintheta*sinthetasq) + &
+ cosphifour*(d11*costhetafour + 2.*d15*costhetasq*sintwotheta + &
+ (d13 + 2.*d55)*sintwothetasq/2. + &
+ 2.*d35*sintwotheta*sinthetasq + d33*sinthetafour) + &
+ (sintwophisq/4.)*(d12 + d23 + 2.*(d44 + d66) + &
+ (d12 - d23 - 2.*d44 + 2.*d66)*costwotheta + &
+ 2.*(d25 + 2.*d46)*sintwotheta)
+
+c12 = -((sintwophi/2.)*sinphisq*((3.*d16 - 4.*d26 + d36 + 2.*d45)*costheta + &
+ (d16 - d36 - 2.*d45)*(4.*costhetasq*costheta - 3.*costheta) + &
+ 2.*(d14 - 2.*d24 + d34 + 2.*d56 + &
+ (d14 - d34 + 2.*d56)*costwotheta)*sintheta))/2. + &
+ cosphisq*sintwophi*(d16*costhetasq*costheta - d24*sintheta + &
+ (d14 + 2.*d56)*costhetasq*sintheta + d34*sintheta*sinthetasq + &
+ costheta*(-d26 + (d36 + 2.*d45)*sinthetasq)) + &
+ (sintwophisq/4.)*(d22 + d11*costhetafour + &
+ 2.*d15*costhetasq*sintwotheta - 4.*d44*sinthetasq + &
+ d33*sinthetafour + costhetasq*(-4.*d66 + &
+ 2.*(d13 + 2.*d55)*sinthetasq) + &
+ costheta*(-8.*d46*sintheta + 4.*d35*sintheta*sinthetasq)) + &
+ (cosphifour + sinphifour)*(d12*costhetasq + &
+ d23*sinthetasq + d25*sintwotheta)
+
+c13 = sinphisq*(d23*costhetasq - d25*sintwotheta + d12*sinthetasq) - &
+ sintwophi*(d36*costhetasq*costheta + &
+ (d34 - 2.*d56)*costhetasq*sintheta + &
+ (d16 - 2.*d45)*costheta*sinthetasq + d14*sintheta*sinthetasq) + &
+ (cosphisq*(d11 + 6.*d13 + d33 - 4.*d55 - &
+ (d11 - 2.*d13 + d33 - 4.*d55)*cosfourtheta + &
+ 4.*(-d15 + d35)*sinfourtheta))/8.
+
+c14 = (-4.*cosphi*sinphisq*((-d14 - 2.*d24 + d34 + 2.*d56)*costheta + &
+ (d14 - d34 + 2.*d56)*(4.*costhetasq*costheta - 3.*costheta) + &
+ 2.*(-d16 + d26 + d36 + (-d16 + d36 + 2.*d45)*costwotheta)*sintheta) + &
+ 8.*cosphisq*cosphi*(d14*costhetasq*costheta - &
+ (d16 - 2.*d45)*costhetasq*sintheta + &
+ (d34 - 2.*d56)*costheta*sinthetasq - d36*sintheta*sinthetasq) + &
+ 4.*sinphi*sinphisq*(2.*d25*costwotheta + (-d12 + d23)*sintwotheta) + &
+ cosphisq*sinphi*(4.*(d15 + d35 - 4*d46)*costwotheta + &
+ 4.*(d15 - d35)*cosfourtheta - &
+ 2.*(d11 - d33 + 4.*d44 - 4.*d66 + &
+ (d11 - 2.*d13 + d33 - 4.*d55)*costwotheta)*sintwotheta))/8.
+
+c15 = (8.*sinphi*sinphisq*(-(d24*costheta) + d26*sintheta) + &
+ 4.*cosphi*sinphisq*(2.*(d25 + 2.*d46)*costwotheta + &
+ (-d12 + d23 + 2.*d44 - 2.*d66)*sintwotheta) + &
+ cosphisq*cosphi*(4.*(d15 + d35)*costwotheta + &
+ 4.*(d15 - d35)*cosfourtheta - 2.*(d11 - d33 + &
+ (d11 - 2.*d13 + d33 - 4.*d55)*costwotheta)*sintwotheta) - &
+ 2.*cosphisq*sinphi*((d14 + 3.*d34 + 2.*d56)*costheta + &
+ 3.*(d14 - d34 + 2.*d56)*(4.*costhetasq*costheta - 3.*costheta) - &
+ (3.*d16 + d36 + 2.*d45)*sintheta + &
+ 3.*(-d16 + d36 + 2.*d45)*(-4.*sinthetasq*sintheta + 3.*sintheta)))/8.
+
+c16 = -(sinphifour*(d26*costheta + d24*sintheta)) - &
+ (3.*(sintwophisq/4.)*((3.*d16 - 4.*d26 + d36 + 2.*d45)*costheta + &
+ (d16 - d36 - 2.*d45)*(4.*costhetasq*costheta - 3.*costheta) + &
+ 2.*(d14 - 2.*d24 + d34 + 2.*d56 + &
+ (d14 - d34 + 2.*d56)*costwotheta)*sintheta))/4. + &
+ cosphifour*(d16*costhetasq*costheta + &
+ (d14 + 2.*d56)*costhetasq*sintheta + &
+ (d36 + 2.*d45)*costheta*sinthetasq + d34*sintheta*sinthetasq) + &
+ (sintwophi/2.)*sinphisq*(-d22 + (d12 + 2.*d66)*costhetasq + &
+ 2.*d46*sintwotheta + (d23 + 2.*d44)*sinthetasq + d25*sintwotheta) + &
+ cosphisq*(sintwophi/2.)*(d11*costhetafour + &
+ 2.*d15*costhetasq*sintwotheta - (d23 + 2.*d44)*sinthetasq + &
+ d33*sinthetafour - costhetasq*(d12 + &
+ 2.*d66 - 2.*(d13 + 2.*d55)*sinthetasq) - &
+ (d25 - d35 + 2.*d46 + d35*costwotheta)*sintwotheta)
+
+c22 = d22*cosphifour + 2.*cosphisq*sintwophi*(d26*costheta + d24*sintheta) + &
+ 2.*sintwophi*sinphisq*(d16*costhetasq*costheta + &
+ (d14 + 2.*d56)*costhetasq*sintheta + &
+ (d36 + 2.*d45)*costheta*sinthetasq + d34*sintheta*sinthetasq) + &
+ sinphifour*(d11*costhetafour + 2.*d15*costhetasq*sintwotheta + &
+ (d13 + 2.*d55)*sintwothetasq/2. + &
+ 2.*d35*sintwotheta*sinthetasq + d33*sinthetafour) + &
+ (sintwophisq/4.)*(d12 + d23 + 2.*(d44 + d66) + &
+ (d12 - d23 - 2.*d44 + 2.*d66)*costwotheta + &
+ 2.*(d25 + 2.*d46)*sintwotheta)
+
+c23 = d13*costhetafour*sinphisq + &
+ sintheta*sinthetasq*(d14*sintwophi + d13*sinphisq*sintheta) + &
+ costheta*sinthetasq*((d16 - 2.*d45)*sintwophi + &
+ 2.*(d15 - d35)*sinphisq*sintheta) + &
+ costhetasq*costheta*(d36*sintwophi + &
+ 2.*(-d15 + d35)*sinphisq*sintheta) + &
+ costhetasq*sintheta*((d34 - 2.*d56)*sintwophi + &
+ (d11 + d33 - 4.*d55)*sinphisq*sintheta) + &
+ cosphisq*(d23*costhetasq - d25*sintwotheta + d12*sinthetasq)
+
+c24 = (8.*cosphisq*cosphi*(d24*costheta - d26*sintheta) + &
+ 4.*cosphisq*sinphi*(2.*(d25 + 2.*d46)*costwotheta + &
+ (-d12 + d23 + 2.*d44 - 2.*d66)*sintwotheta) + &
+ sinphi*sinphisq*(4.*(d15 + d35)*costwotheta + &
+ 4.*(d15 - d35)*cosfourtheta - &
+ 2.*(d11 - d33 + (d11 - 2.*d13 + &
+ d33 - 4.*d55)*costwotheta)*sintwotheta) + &
+ 2.*cosphi*sinphisq*((d14 + 3.*d34 + 2.*d56)*costheta + &
+ 3.*(d14 - d34 + 2.*d56)*(4.*costhetasq*costheta - 3.*costheta) - &
+ (3.*d16 + d36 + 2.*d45)*sintheta + &
+ 3.*(-d16 + d36 + 2.*d45)*(-4.*sinthetasq*sintheta + 3.*sintheta)))/8.
+
+c25 = (4.*cosphisq*sinphi*((-d14 - 2.*d24 + d34 + 2.*d56)*costheta + &
+ (d14 - d34 + 2.*d56)*(4.*costhetasq*costheta - 3.*costheta) + &
+ 2.*(-d16 + d26 + d36 + (-d16 + d36 + 2.*d45)*costwotheta)*sintheta) - &
+ 8.*sinphi*sinphisq*(d14*costhetasq*costheta - &
+ (d16 - 2.*d45)*costhetasq*sintheta + &
+ (d34 - 2.*d56)*costheta*sinthetasq - d36*sintheta*sinthetasq) + &
+ 4.*cosphisq*cosphi*(2.*d25*costwotheta + (-d12 + d23)*sintwotheta) + &
+ cosphi*sinphisq*(4.*(d15 + d35 - 4.*d46)*costwotheta + &
+ 4.*(d15 - d35)*cosfourtheta - 2.*(d11 - d33 + 4.*d44 - 4.*d66 + &
+ (d11 - 2.*d13 + d33 - 4.*d55)*costwotheta)*sintwotheta))/8.
+
+c26 = cosphifour*(d26*costheta + d24*sintheta) + &
+ (3.*(sintwophisq/4.)*((3.*d16 - 4.*d26 + d36 + 2.*d45)*costheta + &
+ (d16 - d36 - 2.*d45)*(4.*costhetasq*costheta - 3.*costheta) + &
+ 2.*(d14 - 2.*d24 + d34 + 2.*d56 + &
+ (d14 - d34 + 2.*d56)*costwotheta)*sintheta))/4. - &
+ sinphifour*(d16*costhetasq*costheta + &
+ (d14 + 2.*d56)*costhetasq*sintheta + &
+ (d36 + 2.*d45)*costheta*sinthetasq + d34*sintheta*sinthetasq) + &
+ cosphisq*(sintwophi/2.)*(-d22 + (d12 + 2.*d66)*costhetasq + &
+ 2.*d46*sintwotheta + (d23 + 2.*d44)*sinthetasq + &
+ d25*sintwotheta) + (sintwophi/2.)*sinphisq*(d11*costhetafour + &
+ 2.*d15*costhetasq*sintwotheta - (d23 + 2.*d44)*sinthetasq + &
+ d33*sinthetafour - costhetasq*(d12 + &
+ 2.*d66 - 2.*(d13 + 2.*d55)*sinthetasq) - &
+ (d25 - d35 + 2.*d46 + d35*costwotheta)*sintwotheta)
+
+c33 = d33*costhetafour - 2.*d35*costhetasq*sintwotheta + &
+ (d13 + 2.*d55)*sintwothetasq/2. - &
+ 2.*d15*sintwotheta*sinthetasq + d11*sinthetafour
+
+c34 = cosphi*(d34*costhetasq*costheta - (d36 + 2.*d45)*costhetasq*sintheta + &
+ (d14 + 2.*d56)*costheta*sinthetasq - d16*sintheta*sinthetasq) + &
+ (sinphi*(4.*(d15 + d35)*costwotheta + 4.*(-d15 + d35)*cosfourtheta + &
+ 2.*(-d11 + d33)*sintwotheta + &
+ (d11 - 2.*d13 + d33 - 4.*d55)*sinfourtheta))/8.
+
+c35 = sinphi*(-(d34*costhetasq*costheta) + &
+ (d36 + 2.*d45)*costhetasq*sintheta - &
+ (d14 + 2.*d56)*costheta*sinthetasq + d16*sintheta*sinthetasq) + &
+ (cosphi*(4.*(d15 + d35)*costwotheta + 4.*(-d15 + d35)*cosfourtheta + &
+ 2.*(-d11 + d33)*sintwotheta + &
+ (d11 - 2.*d13 + d33 - 4.*d55)*sinfourtheta))/8.
+
+c36 = (4.*costwophi*((d16 + 3.*d36 - 2.*d45)*costheta + &
+ (-d16 + d36 + 2.*d45)*(4.*costhetasq*costheta - 3.*costheta) + &
+ (3.*d14 + d34 - 2.*d56)*sintheta + &
+ (-d14 + d34 - 2.*d56)*(-4.*sinthetasq*sintheta + 3.*sintheta)) + &
+ sintwophi*(d11 - 4.*d12 + 6.*d13 - 4.*d23 + d33 - 4.*d55 + &
+ 4.*(d12 - d23)*costwotheta - &
+ (d11 - 2.*d13 + d33 - 4.*d55)*cosfourtheta + &
+ 8.*d25*sintwotheta + 4.*(-d15 + d35)*sinfourtheta))/16.
+
+c44 = (d11 - 2.*d13 + d33 + 4.*(d44 + d55 + d66) - &
+ (d11 - 2.*d13 + d33 - 4.*(d44 - d55 + d66))*costwophi + &
+ 4.*sintwophi*((d16 - d36 + 2.*d45)*costheta + &
+ (-d16 + d36 + 2.*d45)*(4.*costhetasq*costheta - 3.*costheta) - &
+ 2.*(d14 - d34 + (d14 - d34 + 2.*d56)*costwotheta)*sintheta) + &
+ 8.*cosphisq*((d44 - d66)*costwotheta - 2.*d46*sintwotheta) + &
+ 2.*sinphisq*(-((d11 - 2.*d13 + d33 - 4.*d55)*cosfourtheta) + &
+ 4.*(-d15 + d35)*sinfourtheta))/16.
+
+c45 = (4.*costwophi*((d16 - d36 + 2.*d45)*costheta + &
+ (-d16 + d36 + 2.*d45)*(4.*costhetasq*costheta - 3.*costheta) - &
+ 2.*(d14 - d34 + (d14 - d34 + 2.*d56)*costwotheta)*sintheta) + &
+ sintwophi*(d11 - 2.*d13 + d33 - 4.*(d44 - d55 + d66) + &
+ 4.*(-d44 + d66)*costwotheta - &
+ (d11 - 2.*d13 + d33 - 4.*d55)*cosfourtheta + 8.*d46*sintwotheta + &
+ 4.*(-d15 + d35)*sinfourtheta))/16.
+
+c46 = (-2.*sinphi*sinphisq*((-d14 + d34 + 2.*d56)*costheta + &
+ (d14 - d34 + 2.*d56)*(4.*costhetasq*costheta - 3.*costheta) + &
+ 2.*(-d16 + d36 + (-d16 + d36 + 2.*d45)*costwotheta)*sintheta) + &
+ 4.*cosphisq*cosphi*(2.*d46*costwotheta + (d44 - d66)*sintwotheta) + &
+ cosphi*sinphisq*(4.*(d15 - 2.*d25 + d35 - 2.*d46)*costwotheta + &
+ 4.*(d15 - d35)*cosfourtheta - &
+ 2.*(d11 - 2.*d12 + 2.*d23 - d33 + 2.*d44 - 2.*d66 + &
+ (d11 - 2.*d13 + d33 - 4.*d55)*costwotheta)*sintwotheta) + &
+ 4.*cosphisq*sinphi*((d14 - 2.*d24 + d34)*costheta + &
+ (d14 - d34 + 2.*d56)*(4.*costhetasq*costheta - 3.*costheta) - &
+ (d16 - 2.*d26 + d36)*sintheta + &
+ (-d16 + d36 + 2.*d45)*(-4.*sinthetasq*sintheta + 3.*sintheta)))/8.
+
+c55 = d66*sinphisq*sinthetasq + (sintwotheta/2.)*(-2.*d46*sinphisq + &
+ (d36 + d45)*sintwophi*sintheta) + &
+ costhetasq*(d44*sinphisq + (d14 + d56)*sintwophi*sintheta) - &
+ sintwophi*(d45*costhetasq*costheta + d34*costhetasq*sintheta + &
+ d16*costheta*sinthetasq + d56*sintheta*sinthetasq) + &
+ (cosphisq*(d11 - 2.*d13 + d33 + 4.*d55 - &
+ (d11 - 2.*d13 + d33 - 4.*d55)*cosfourtheta + &
+ 4.*(-d15 + d35)*sinfourtheta))/8.
+
+c56 = (8.*cosphisq*cosphi*(d56*costhetasq*costheta - &
+ (d16 - d36 - d45)*costhetasq*sintheta - &
+ (d14 - d34 + d56)*costheta*sinthetasq - d45*sintheta*sinthetasq) + &
+ 4.*sinphi*sinphisq*(2.*d46*costwotheta + (d44 - d66)*sintwotheta) + &
+ cosphisq*sinphi*(4.*(d15 - 2.*d25 + d35 - 2.*d46)*costwotheta + &
+ 4.*(d15 - d35)*cosfourtheta - &
+ 2.*(d11 - 2.*d12 + 2.*d23 - d33 + 2.*d44 - 2.*d66 + &
+ (d11 - 2.*d13 + d33 - 4.*d55)*costwotheta)*sintwotheta) - &
+ 4.*cosphi*sinphisq*((d14 - 2.*d24 + d34)*costheta + &
+ (d14 - d34 + 2.*d56)*(4.*costhetasq*costheta - 3.*costheta) - &
+ (d16 - 2.*d26 + d36)*sintheta + &
+ (-d16 + d36 + 2.*d45)*(-4.*sinthetasq*sintheta + 3.*sintheta)))/8.
+
+c66 = -((sintwophi/2.)*sinphisq*((3.*d16 - 4.*d26 + d36 + 2.*d45)*costheta + &
+ (d16 - d36 - 2.*d45)*(4.*costhetasq*costheta - 3.*costheta) + &
+ 2.*(d14 - 2.*d24 + d34 + 2.*d56 + &
+ (d14 - d34 + 2.*d56)*costwotheta)*sintheta))/2. + &
+ cosphisq*sintwophi*(d16*costhetasq*costheta - d24*sintheta + &
+ (d14 + 2.*d56)*costhetasq*sintheta + d34*sintheta*sinthetasq + &
+ costheta*(-d26 + (d36 + 2.*d45)*sinthetasq)) + &
+ (sintwophisq/4.)*(d22 + d11*costhetafour + &
+ 2.*d15*costhetasq*sintwotheta - 2.*(d23 + d44)*sinthetasq + &
+ d33*sinthetafour - 2.*sintwotheta*(d25 + d46 - d35*sinthetasq) - &
+ 2.*costhetasq*(d12 + d66 - (d13 + 2.*d55)*sinthetasq)) + &
+ (cosphifour + sinphifour)*(d66*costhetasq + &
+ d44*sinthetasq + d46*sintwotheta)
+
+
+end subroutine rotate_aniso_tensor
+!--------------------------------------------------------------------
+
diff -r 000000000000 -r 9b613f27327a src/meshfem3D/model_atten3D_QRFSI12.f90
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/meshfem3D/model_atten3D_QRFSI12.f90 Sat Jul 02 17:09:21 2011 -0700
@@ -0,0 +1,736 @@
+!=====================================================================
+!
+! 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 University of Pau / CNRS / INRIA, France
+! (c) Princeton University / California Institute of Technology and University of Pau / CNRS / INRIA
+! 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.
+!
+!=====================================================================
+
+!--------------------------------------------------------------------------------------------------
+!
+! This file contains subroutines to read in and get values for
+! 3-D attenuation model QRFSI12 (Dalton, Ekstrom, & Dziewonski, 2008)
+!
+! C.A. Dalton, G. Ekstr\"om and A.M. Dziewonski, 2008.
+! The global attenuation structure of the upper mantle,
+! J. Geophys. Res., 113, B05317,10.1029/2006JB004394
+!
+! Last edit: Colleen Dalton, March 25, 2008
+!
+! Q1: what are theta and phi?
+! Q2: units for radius?
+! Q3: what to do about core?
+!--------------------------------------------------------------------------------------------------
+
+
+ subroutine model_atten3D_QRFSI12_broadcast(myrank,QRFSI12_Q)
+
+! standard routine to setup model
+
+ implicit none
+
+ include "constants.h"
+ ! standard include of the MPI library
+ include 'mpif.h'
+
+ ! model_atten3D_QRFSI12_variables
+ type model_atten3D_QRFSI12_variables
+ sequence
+ double precision dqmu(NKQ,NSQ)
+ double precision spknt(NKQ)
+ double precision refdepth(NDEPTHS_REFQ)
+ double precision refqmu(NDEPTHS_REFQ)
+ end type model_atten3D_QRFSI12_variables
+
+ type (model_atten3D_QRFSI12_variables) QRFSI12_Q
+ ! model_atten3D_QRFSI12_variables
+
+ integer :: myrank
+ integer :: ier
+
+ if(myrank == 0) call read_atten_model_3D_QRFSI12(QRFSI12_Q)
+
+ call MPI_BCAST(QRFSI12_Q%dqmu, NKQ*NSQ,MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
+ call MPI_BCAST(QRFSI12_Q%spknt, NKQ,MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
+ call MPI_BCAST(QRFSI12_Q%refdepth, NDEPTHS_REFQ,MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
+ call MPI_BCAST(QRFSI12_Q%refqmu, NDEPTHS_REFQ,MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ier)
+
+ if(myrank == 0) write(IMAIN,*) 'read 3D attenuation model'
+
+
+ end subroutine
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine read_atten_model_3D_QRFSI12(QRFSI12_Q)
+
+ implicit none
+
+ include "constants.h"
+
+! three_d_model_atten3D_QRFSI12_variables
+ type model_atten3D_QRFSI12_variables
+ sequence
+ double precision dqmu(NKQ,NSQ)
+ double precision spknt(NKQ)
+ double precision refdepth(NDEPTHS_REFQ)
+ double precision refqmu(NDEPTHS_REFQ)
+ end type model_atten3D_QRFSI12_variables
+
+ type (model_atten3D_QRFSI12_variables) QRFSI12_Q
+! three_d_model_atten3D_QRFSI12_variables
+
+ integer j,k,l,m
+ integer index,ll,mm
+ double precision v1,v2
+
+ character(len=150) QRFSI12,QRFSI12_ref
+
+! read in QRFSI12
+! hard-wire for now
+ QRFSI12='DATA/QRFSI12/QRFSI12.dat'
+ QRFSI12_ref='DATA/QRFSI12/ref_QRFSI12'
+
+! get the dq model coefficients
+ open(unit=10,file=QRFSI12,status='old',action='read')
+ do k=1,NKQ
+ read(10,*)index
+ j=0
+ do l=0,MAXL_Q
+ do m=0,l
+ if(m.eq.0)then
+ j=j+1
+ read(10,*)ll,mm,v1
+ QRFSI12_Q%dqmu(k,j)=v1
+ else
+ j=j+2
+ read(10,*)ll,mm,v1,v2
+ ! write(*,*) 'k,l,m,ll,mm:',k,l,m,ll,mm,v1
+ QRFSI12_Q%dqmu(k,j-1)=2.*v1
+ QRFSI12_Q%dqmu(k,j)=-2.*v2
+ endif
+ enddo
+ enddo
+ enddo
+ close(10)
+
+! get the depths (km) of the spline knots
+ QRFSI12_Q%spknt(1) = 24.4
+ QRFSI12_Q%spknt(2) = 75.0
+ QRFSI12_Q%spknt(3) = 150.0
+ QRFSI12_Q%spknt(4) = 225.0
+ QRFSI12_Q%spknt(5) = 300.0
+ QRFSI12_Q%spknt(6) = 410.0
+ QRFSI12_Q%spknt(7) = 530.0
+ QRFSI12_Q%spknt(8) = 650.0
+
+! get the depths and 1/Q values of the reference model
+ open(11,file=QRFSI12_ref,status='old',action='read')
+ do j=1,NDEPTHS_REFQ
+ read(11,*)QRFSI12_Q%refdepth(j),QRFSI12_Q%refqmu(j)
+ enddo
+ close(11)
+
+
+ end subroutine read_atten_model_3D_QRFSI12
+
+!----------------------------------
+!----------------------------------
+
+ subroutine model_atten3D_QRFSI12(radius,theta,phi,Qmu,QRFSI12_Q,idoubling)
+
+ implicit none
+
+ include "constants.h"
+
+! model_atten3D_QRFSI12_variables
+ type model_atten3D_QRFSI12_variables
+ sequence
+ double precision dqmu(NKQ,NSQ)
+ double precision spknt(NKQ)
+ double precision refdepth(NDEPTHS_REFQ)
+ double precision refqmu(NDEPTHS_REFQ)
+ end type model_atten3D_QRFSI12_variables
+
+ type (model_atten3D_QRFSI12_variables) QRFSI12_Q
+! model_atten3D_QRFSI12_variables
+
+ integer i,j,k,n,idoubling
+ integer ifnd
+ double precision radius,theta,phi,Qmu,smallq,dqmu,smallq_ref
+ real(kind=4) splpts(NKQ),splcon(NKQ),splcond(NKQ)
+ real(kind=4) depth,ylat,xlon
+ real(kind=4) shdep(NSQ)
+ real(kind=4) xlmvec(NSQ),wk1(NSQ),wk2(NSQ),wk3(NSQ)
+ double precision, parameter :: rmoho_prem = 6371.0-24.4
+ double precision, parameter :: rcmb = 3480.0
+
+ !in Colleen's original code theta refers to the latitude. Here we have redefined theta to be colatitude
+ !to agree with the rest of specfem
+! print *,'entering QRFSI12 subroutine'
+
+ ylat=90.0d0-theta
+ xlon=phi
+
+! only checks radius for crust, idoubling is missleading for oceanic crust when we want to expand mantle up to surface...
+! !if(idoubling == IFLAG_CRUST .or. radius >= rmoho) then
+ if( radius >= rmoho_prem ) then
+ ! print *,'QRFSI12: we are in the crust'
+ Qmu = 600.0d0
+ else if(idoubling == IFLAG_INNER_CORE_NORMAL .or. idoubling == IFLAG_MIDDLE_CENTRAL_CUBE .or. &
+ idoubling == IFLAG_BOTTOM_CENTRAL_CUBE .or. idoubling == IFLAG_TOP_CENTRAL_CUBE .or. &
+ idoubling == IFLAG_IN_FICTITIOUS_CUBE) then
+ ! print *,'QRFSI12: we are in the inner core'
+ Qmu = 84.6d0
+ else if(idoubling == IFLAG_OUTER_CORE_NORMAL) then
+ ! print *,'QRFSI12: we are in the outer core'
+ Qmu = 0.0d0
+ else !we are in the mantle
+ depth = 6371.-radius
+! print *,'QRFSI12: we are in the mantle at depth',depth
+ ifnd=0
+ do i=2,NDEPTHS_REFQ
+ if(depth >= QRFSI12_Q%refdepth(i-1) .and. depth < QRFSI12_Q%refdepth(i))then
+ ifnd=i
+ endif
+ enddo
+ if(ifnd == 0)then
+ write(6,"('problem finding reference Q value at depth: ',f8.3)") depth
+ stop
+ endif
+ smallq_ref=QRFSI12_Q%refqmu(ifnd)
+ smallq = smallq_ref
+
+ if(depth < 650.d0) then !Colleen's model is only defined between depths of 24.4 and 650km
+ do j=1,NSQ
+ shdep(j)=0.
+ enddo
+ do n=1,NKQ
+ splpts(n)=QRFSI12_Q%spknt(n)
+ enddo
+ call vbspl(depth,NKQ,splpts,splcon,splcond)
+ do n=1,NKQ
+ do j=1,NSQ
+ shdep(j)=shdep(j)+(splcon(n)*QRFSI12_Q%dqmu(n,j))
+ enddo
+ enddo
+ call ylm(ylat,xlon,MAXL_Q,xlmvec,wk1,wk2,wk3)
+ dqmu=0.
+ do k=1,NSQ
+ dqmu=dqmu+xlmvec(k)*shdep(k)
+ enddo
+ smallq = smallq_ref + dqmu
+ endif
+ ! if smallq is small and negative (due to numerical error), Qmu is very large:
+ if(smallq < 0.0d0) smallq = 1.0d0/ATTENUATION_COMP_MAXIMUM
+ Qmu = 1/smallq
+ ! Qmu is larger than MAX_ATTENUATION_VALUE, set it to ATTENUATION_COMP_MAXIMUM. This assumes that this
+ ! value is high enough that at this point there is almost no attenuation at all.
+ if(Qmu >= ATTENUATION_COMP_MAXIMUM) Qmu = 0.99d0*ATTENUATION_COMP_MAXIMUM
+
+ endif
+
+ end subroutine model_atten3D_QRFSI12
+
+!----------------------------------
+!----------------------------------
+
+!!$ subroutine vbspl(x,np,xarr,splcon,splcond)
+!!$!
+!!$!---- this subroutine returns the spline contributions at a particular value of x
+!!$!
+!!$ implicit none
+!!$
+!!$ integer :: np
+!!$
+!!$ real(kind=4) :: xarr(np),x
+!!$ real(kind=4) :: splcon(np)
+!!$ real(kind=4) :: splcond(np)
+!!$
+!!$ real(kind=4) :: r1,r2,r3,r4,r5,r6,r7,r8,r9,r10,r11,r12,r13
+!!$ real(kind=4) :: r1d,r2d,r3d,r4d,r5d,r6d,r7d,r8d,r9d,r10d,r11d,r12d,r13d,val,vald
+!!$
+!!$ real(kind=4) :: rr1,rr2,rr3,rr4,rr5,rr6,rr7,rr8,rr9,rr10,rr11,rr12
+!!$ real(kind=4) :: rr1d,rr2d,rr3d,rr4d,rr5d,rr6d,rr7d,rr8d,rr9d,rr10d,rr11d,rr12d
+!!$
+!!$ integer :: iflag,interval,ik,ib
+!!$
+!!$!
+!!$!---- iflag=1 ==>> second derivative is 0 at end points
+!!$!---- iflag=0 ==>> first derivative is 0 at end points
+!!$!
+!!$ iflag=1
+!!$!
+!!$!---- first, find out within which interval x falls
+!!$!
+!!$ interval=0
+!!$ ik=1
+!!$ do while(interval == 0.and.ik < np)
+!!$ ik=ik+1
+!!$ if(x >= xarr(ik-1).and.x <= xarr(ik)) interval=ik-1
+!!$ enddo
+!!$ if(x > xarr(np)) then
+!!$ interval=np
+!!$ endif
+!!$
+!!$ if(interval == 0) then
+!!$! write(6,"('low value:',2f10.3)") x,xarr(1)
+!!$ else if(interval > 0.and.interval < np) then
+!!$! write(6,"('bracket:',i5,3f10.3)") interval,xarr(interval),x,xarr(interval+1)
+!!$ else
+!!$! write(6,"('high value:',2f10.3)") xarr(np),x
+!!$ endif
+!!$
+!!$ do ib=1,np
+!!$ val=0.
+!!$ vald=0.
+!!$ if(ib == 1) then
+!!$
+!!$ r1=(x-xarr(1))/(xarr(2)-xarr(1))
+!!$ r2=(xarr(3)-x)/(xarr(3)-xarr(1))
+!!$ r4=(xarr(2)-x)/(xarr(2)-xarr(1))
+!!$ r5=(x-xarr(1))/(xarr(2)-xarr(1))
+!!$ r6=(xarr(3)-x)/(xarr(3)-xarr(1))
+!!$ r10=(xarr(2)-x)/(xarr(2)-xarr(1))
+!!$ r11=(x-xarr(1)) /(xarr(2)-xarr(1))
+!!$ r12=(xarr(3)-x)/(xarr(3)-xarr(2))
+!!$ r13=(xarr(2)-x)/(xarr(2)-xarr(1))
+!!$
+!!$ r1d=1./(xarr(2)-xarr(1))
+!!$ r2d=-1./(xarr(3)-xarr(1))
+!!$ r4d=-1./(xarr(2)-xarr(1))
+!!$ r5d=1./(xarr(2)-xarr(1))
+!!$ r6d=-1./(xarr(3)-xarr(1))
+!!$ r10d=-1./(xarr(2)-xarr(1))
+!!$ r11d=1./(xarr(2)-xarr(1))
+!!$ r12d=-1./(xarr(3)-xarr(2))
+!!$ r13d=-1./(xarr(2)-xarr(1))
+!!$
+!!$ if(interval == ib.or.interval == 0) then
+!!$ if(iflag == 0) then
+!!$ val=r1*r4*r10 + r2*r5*r10 + r2*r6*r11 +r13**3
+!!$ vald=r1d*r4*r10+r1*r4d*r10+r1*r4*r10d
+!!$ vald=vald+r2d*r5*r10+r2*r5d*r10+r2*r5*r10d
+!!$ vald=vald+r2d*r6*r11+r2*r6d*r11+r2*r6*r11d
+!!$ vald=vald+3.*r13d*r13**2
+!!$ else if(iflag == 1) then
+!!$ val=0.6667*(r1*r4*r10 + r2*r5*r10 + r2*r6*r11 &
+!!$ + 1.5*r13**3)
+!!$ vald=r1d*r4*r10+r1*r4d*r10+r1*r4*r10d
+!!$ vald=vald+r2d*r5*r10+r2*r5d*r10+r2*r5*r10d
+!!$ vald=vald+r2d*r6*r11+r2*r6d*r11+r2*r6*r11d
+!!$ vald=vald+4.5*r13d*r13**2
+!!$ vald=0.6667*vald
+!!$ endif
+!!$ else if(interval == ib+1) then
+!!$ if(iflag == 0) then
+!!$ val=r2*r6*r12
+!!$ vald=r2d*r6*r12+r2*r6d*r12+r2*r6*r12d
+!!$ else if(iflag == 1) then
+!!$ val=0.6667*r2*r6*r12
+!!$ vald=0.6667*(r2d*r6*r12+r2*r6d*r12+r2*r6*r12d)
+!!$ endif
+!!$ else
+!!$ val=0.
+!!$ endif
+!!$
+!!$ else if(ib == 2) then
+!!$
+!!$ rr1=(x-xarr(1))/(xarr(2)-xarr(1))
+!!$ rr2=(xarr(3)-x)/(xarr(3)-xarr(1))
+!!$ rr4=(xarr(2)-x)/(xarr(2)-xarr(1))
+!!$ rr5=(x-xarr(1))/(xarr(2)-xarr(1))
+!!$ rr6=(xarr(3)-x)/(xarr(3)-xarr(1))
+!!$ rr10=(xarr(2)-x)/(xarr(2)-xarr(1))
+!!$ rr11=(x-xarr(1)) /(xarr(2)-xarr(1))
+!!$ rr12=(xarr(3)-x)/(xarr(3)-xarr(2))
+!!$
+!!$ rr1d=1./(xarr(2)-xarr(1))
+!!$ rr2d=-1./(xarr(3)-xarr(1))
+!!$ rr4d=-1./(xarr(2)-xarr(1))
+!!$ rr5d=1./(xarr(2)-xarr(1))
+!!$ rr6d=-1./(xarr(3)-xarr(1))
+!!$ rr10d=-1./(xarr(2)-xarr(1))
+!!$ rr11d=1./(xarr(2)-xarr(1))
+!!$ rr12d=-1./(xarr(3)-xarr(2))
+!!$
+!!$ r1=(x-xarr(ib-1))/(xarr(ib+1)-xarr(ib-1))
+!!$ r2=(xarr(ib+2)-x)/(xarr(ib+2)-xarr(ib-1))
+!!$ r3=(x-xarr(ib-1))/(xarr(ib)-xarr(ib-1))
+!!$ r4=(xarr(ib+1)-x)/(xarr(ib+1)-xarr(ib-1))
+!!$ r5=(x-xarr(ib-1))/(xarr(ib+1)-xarr(ib-1))
+!!$ r6=(xarr(ib+2)-x)/(xarr(ib+2)-xarr(ib))
+!!$ r8=(xarr(ib)-x)/ (xarr(ib)-xarr(ib-1))
+!!$ r9=(x-xarr(ib-1))/(xarr(ib)-xarr(ib-1))
+!!$ r10=(xarr(ib+1)-x)/(xarr(ib+1)-xarr(ib))
+!!$ r11=(x-xarr(ib)) /(xarr(ib+1)-xarr(ib))
+!!$ r12=(xarr(ib+2)-x)/(xarr(ib+2)-xarr(ib+1))
+!!$
+!!$ r1d=1./(xarr(ib+1)-xarr(ib-1))
+!!$ r2d=-1./(xarr(ib+2)-xarr(ib-1))
+!!$ r3d=1./(xarr(ib)-xarr(ib-1))
+!!$ r4d=-1./(xarr(ib+1)-xarr(ib-1))
+!!$ r5d=1./(xarr(ib+1)-xarr(ib-1))
+!!$ r6d=-1./(xarr(ib+2)-xarr(ib))
+!!$ r8d=-1./ (xarr(ib)-xarr(ib-1))
+!!$ r9d=1./(xarr(ib)-xarr(ib-1))
+!!$ r10d=-1./(xarr(ib+1)-xarr(ib))
+!!$ r11d=1./(xarr(ib+1)-xarr(ib))
+!!$ r12d=-1./(xarr(ib+2)-xarr(ib+1))
+!!$
+!!$ if(interval == ib-1.or.interval == 0) then
+!!$ val=r1*r3*r8 + r1*r4*r9 + r2*r5*r9
+!!$ vald=r1d*r3*r8+r1*r3d*r8+r1*r3*r8d
+!!$ vald=vald+r1d*r4*r9+r1*r4d*r9+r1*r4*r9d
+!!$ vald=vald+r2d*r5*r9+r2*r5d*r9+r2*r5*r9d
+!!$ if(iflag == 1) then
+!!$ val=val+0.3333*(rr1*rr4*rr10 + rr2*rr5*rr10 + &
+!!$ rr2*rr6*rr11)
+!!$ vald=vald+0.3333*(rr1d*rr4*rr10+rr1*rr4d*rr10+ &
+!!$ rr1*rr4*rr10d)
+!!$ vald=vald+0.3333*(rr2d*rr5*rr10+rr2*rr5d*rr10+ &
+!!$ rr2*rr5*rr10d)
+!!$ vald=vald+0.3333*(rr2d*rr6*rr11+rr2*rr6d*rr11+ &
+!!$ rr2*rr6*rr11d)
+!!$ endif
+!!$ else if(interval == ib) then
+!!$ val=r1*r4*r10 + r2*r5*r10 + r2*r6*r11
+!!$ vald=r1d*r4*r10+r1*r4d*r10+r1*r4*r10d
+!!$ vald=vald+r2d*r5*r10+r2*r5d*r10+r2*r5*r10d
+!!$ vald=vald+r2d*r6*r11+r2*r6d*r11+r2*r6*r11d
+!!$ if(iflag == 1) then
+!!$ val=val+0.3333*rr2*rr6*rr12
+!!$ vald=vald+0.3333*(rr2d*rr6*rr12+rr2*rr6d*rr12+ &
+!!$ rr2*rr6*rr12d)
+!!$ endif
+!!$ else if(interval == ib+1) then
+!!$ val=r2*r6*r12
+!!$ vald=r2d*r6*r12+r2*r6d*r12+r2*r6*r12d
+!!$ else
+!!$ val=0.
+!!$ endif
+!!$ else if(ib == np-1) then
+!!$
+!!$ rr1=(x-xarr(np-2))/(xarr(np)-xarr(np-2))
+!!$ rr2=(xarr(np)-x)/(xarr(np)-xarr(np-1))
+!!$ rr3=(x-xarr(np-2))/(xarr(np)-xarr(np-2))
+!!$ rr4=(xarr(np)-x)/(xarr(np)-xarr(np-1))
+!!$ rr5=(x-xarr(np-1))/(xarr(np)-xarr(np-1))
+!!$ rr7=(x-xarr(np-2))/(xarr(np-1)-xarr(np-2))
+!!$ rr8=(xarr(np)-x)/ (xarr(np)-xarr(np-1))
+!!$ rr9=(x-xarr(np-1))/(xarr(np)-xarr(np-1))
+!!$
+!!$ rr1d=1./(xarr(np)-xarr(np-2))
+!!$ rr2d=-1./(xarr(np)-xarr(np-1))
+!!$ rr3d=1./(xarr(np)-xarr(np-2))
+!!$ rr4d=-1./(xarr(np)-xarr(np-1))
+!!$ rr5d=1./(xarr(np)-xarr(np-1))
+!!$ rr7d=1./(xarr(np-1)-xarr(np-2))
+!!$ rr8d=-1./ (xarr(np)-xarr(np-1))
+!!$ rr9d=1./(xarr(np)-xarr(np-1))
+!!$
+!!$ r1=(x-xarr(ib-2))/(xarr(ib+1)-xarr(ib-2))
+!!$ r2=(xarr(ib+1)-x)/(xarr(ib+1)-xarr(ib-1))
+!!$ r3=(x-xarr(ib-2))/(xarr(ib)-xarr(ib-2))
+!!$ r4=(xarr(ib+1)-x)/(xarr(ib+1)-xarr(ib-1))
+!!$ r5=(x-xarr(ib-1))/(xarr(ib+1)-xarr(ib-1))
+!!$ r6=(xarr(ib+1)-x)/(xarr(ib+1)-xarr(ib))
+!!$ r7=(x-xarr(ib-2))/(xarr(ib-1)-xarr(ib-2))
+!!$ r8=(xarr(ib)-x)/ (xarr(ib)-xarr(ib-1))
+!!$ r9=(x-xarr(ib-1))/(xarr(ib)-xarr(ib-1))
+!!$ r10=(xarr(ib+1)-x)/(xarr(ib+1)-xarr(ib))
+!!$ r11=(x-xarr(ib)) /(xarr(ib+1)-xarr(ib))
+!!$
+!!$ r1d=1./(xarr(ib+1)-xarr(ib-2))
+!!$ r2d=-1./(xarr(ib+1)-xarr(ib-1))
+!!$ r3d=1./(xarr(ib)-xarr(ib-2))
+!!$ r4d=-1./(xarr(ib+1)-xarr(ib-1))
+!!$ r5d=1./(xarr(ib+1)-xarr(ib-1))
+!!$ r6d=-1./(xarr(ib+1)-xarr(ib))
+!!$ r7d=1./(xarr(ib-1)-xarr(ib-2))
+!!$ r8d=-1./(xarr(ib)-xarr(ib-1))
+!!$ r9d=1./(xarr(ib)-xarr(ib-1))
+!!$ r10d=-1./(xarr(ib+1)-xarr(ib))
+!!$ r11d=1./(xarr(ib+1)-xarr(ib))
+!!$
+!!$ if(interval == ib-2) then
+!!$ val=r1*r3*r7
+!!$ vald=r1d*r3*r7+r1*r3d*r7+r1*r3*r7d
+!!$ else if(interval == ib-1) then
+!!$ val=r1*r3*r8 + r1*r4*r9 + r2*r5*r9
+!!$ vald=r1d*r3*r8+r1*r3d*r8+r1*r3*r8d
+!!$ vald=vald+r1d*r4*r9+r1*r4d*r9+r1*r4*r9d
+!!$ vald=vald+r2d*r5*r9+r2*r5d*r9+r2*r5*r9d
+!!$ if(iflag == 1) then
+!!$ val=val+0.3333*rr1*rr3*rr7
+!!$ vald=vald+0.3333*(rr1d*rr3*rr7+rr1*rr3d*rr7+ &
+!!$ rr1*rr3*rr7d)
+!!$ endif
+!!$ else if(interval == ib.or.interval == np) then
+!!$ val=r1*r4*r10 + r2*r5*r10 + r2*r6*r11
+!!$ vald=r1d*r4*r10+r1*r4d*r10+r1*r4*r10d
+!!$ vald=vald+r2d*r5*r10+r2*r5d*r10+r2*r5*r10d
+!!$ vald=vald+r2d*r6*r11+r2*r6d*r11+r2*r6*r11d
+!!$ if(iflag == 1) then
+!!$ val=val+0.3333*(rr1*rr3*rr8 + rr1*rr4*rr9 + &
+!!$ rr2*rr5*rr9)
+!!$ vald=vald+0.3333*(rr1d*rr3*rr8+rr1*rr3d*rr8+ &
+!!$ rr1*rr3*rr8d)
+!!$ vald=vald+0.3333*(rr1d*rr4*rr9+rr1*rr4d*rr9+ &
+!!$ rr1*rr4*rr9d)
+!!$ vald=vald+0.3333*(rr2d*rr5*rr9+rr2*rr5d*rr9+ &
+!!$ rr2*rr5*rr9d)
+!!$ endif
+!!$ else
+!!$ val=0.
+!!$ endif
+!!$ else if(ib == np) then
+!!$
+!!$ r1=(x-xarr(np-2))/(xarr(np)-xarr(np-2))
+!!$ r2=(xarr(np)-x)/(xarr(np)-xarr(np-1))
+!!$ r3=(x-xarr(np-2))/(xarr(np)-xarr(np-2))
+!!$ r4=(xarr(np)-x)/(xarr(np)-xarr(np-1))
+!!$ r5=(x-xarr(np-1))/(xarr(np)-xarr(np-1))
+!!$ r7=(x-xarr(np-2))/(xarr(np-1)-xarr(np-2))
+!!$ r8=(xarr(np)-x)/ (xarr(np)-xarr(np-1))
+!!$ r9=(x-xarr(np-1))/(xarr(np)-xarr(np-1))
+!!$ r13=(x-xarr(np-1))/(xarr(np)-xarr(np-1))
+!!$
+!!$ r1d=1./(xarr(np)-xarr(np-2))
+!!$ r2d=-1./(xarr(np)-xarr(np-1))
+!!$ r3d=1./(xarr(np)-xarr(np-2))
+!!$ r4d=-1./(xarr(np)-xarr(np-1))
+!!$ r5d=1./(xarr(np)-xarr(np-1))
+!!$ r7d=1./(xarr(np-1)-xarr(np-2))
+!!$ r8d=-1./ (xarr(np)-xarr(np-1))
+!!$ r9d=1./(xarr(np)-xarr(np-1))
+!!$ r13d=1./(xarr(np)-xarr(np-1))
+!!$
+!!$ if(interval == np-2) then
+!!$ if(iflag == 0) then
+!!$ val=r1*r3*r7
+!!$ vald=r1d*r3*r7+r1*r3d*r7+r1*r3*r7d
+!!$ else if(iflag == 1) then
+!!$ val=0.6667*r1*r3*r7
+!!$ vald=0.6667*(r1d*r3*r7+r1*r3d*r7+r1*r3*r7d)
+!!$ endif
+!!$ else if(interval == np-1.or.interval == np) then
+!!$ if(iflag == 0) then
+!!$ val=r1*r3*r8 + r1*r4*r9 + r2*r5*r9 + r13**3
+!!$ vald=r1d*r3*r8+r1*r3d*r8+r1*r3*r8d
+!!$ vald=vald+r1d*r4*r9+r1*r4d*r9+r1*r4*r9d
+!!$ vald=vald+r2d*r5*r9+r2*r5d*r9+r2*r5*r9d
+!!$ vald=vald+3.*r13d*r13**2
+!!$ else if(iflag == 1) then
+!!$ val=0.6667*(r1*r3*r8 + r1*r4*r9 + r2*r5*r9 + &
+!!$ 1.5*r13**3)
+!!$ vald=r1d*r3*r8+r1*r3d*r8+r1*r3*r8d
+!!$ vald=vald+r1d*r4*r9+r1*r4d*r9+r1*r4*r9d
+!!$ vald=vald+r2d*r5*r9+r2*r5d*r9+r2*r5*r9d
+!!$ vald=vald+4.5*r13d*r13**2
+!!$ vald=0.6667*vald
+!!$ endif
+!!$ else
+!!$ val=0.
+!!$ endif
+!!$ else
+!!$
+!!$ r1=(x-xarr(ib-2))/(xarr(ib+1)-xarr(ib-2))
+!!$ r2=(xarr(ib+2)-x)/(xarr(ib+2)-xarr(ib-1))
+!!$ r3=(x-xarr(ib-2))/(xarr(ib)-xarr(ib-2))
+!!$ r4=(xarr(ib+1)-x)/(xarr(ib+1)-xarr(ib-1))
+!!$ r5=(x-xarr(ib-1))/(xarr(ib+1)-xarr(ib-1))
+!!$ r6=(xarr(ib+2)-x)/(xarr(ib+2)-xarr(ib))
+!!$ r7=(x-xarr(ib-2))/(xarr(ib-1)-xarr(ib-2))
+!!$ r8=(xarr(ib)-x)/ (xarr(ib)-xarr(ib-1))
+!!$ r9=(x-xarr(ib-1))/(xarr(ib)-xarr(ib-1))
+!!$ r10=(xarr(ib+1)-x)/(xarr(ib+1)-xarr(ib))
+!!$ r11=(x-xarr(ib)) /(xarr(ib+1)-xarr(ib))
+!!$ r12=(xarr(ib+2)-x)/(xarr(ib+2)-xarr(ib+1))
+!!$
+!!$ r1d=1./(xarr(ib+1)-xarr(ib-2))
+!!$ r2d=-1./(xarr(ib+2)-xarr(ib-1))
+!!$ r3d=1./(xarr(ib)-xarr(ib-2))
+!!$ r4d=-1./(xarr(ib+1)-xarr(ib-1))
+!!$ r5d=1./(xarr(ib+1)-xarr(ib-1))
+!!$ r6d=-1./(xarr(ib+2)-xarr(ib))
+!!$ r7d=1./(xarr(ib-1)-xarr(ib-2))
+!!$ r8d=-1./ (xarr(ib)-xarr(ib-1))
+!!$ r9d=1./(xarr(ib)-xarr(ib-1))
+!!$ r10d=-1./(xarr(ib+1)-xarr(ib))
+!!$ r11d=1./(xarr(ib+1)-xarr(ib))
+!!$ r12d=-1./(xarr(ib+2)-xarr(ib+1))
+!!$
+!!$ if(interval == ib-2) then
+!!$ val=r1*r3*r7
+!!$ vald=r1d*r3*r7+r1*r3d*r7+r1*r3*r7d
+!!$ else if(interval == ib-1) then
+!!$ val=r1*r3*r8 + r1*r4*r9 + r2*r5*r9
+!!$ vald=r1d*r3*r8+r1*r3d*r8+r1*r3*r8d
+!!$ vald=vald+r1d*r4*r9+r1*r4d*r9+r1*r4*r9d
+!!$ vald=vald+r2d*r5*r9+r2*r5d*r9+r2*r5*r9d
+!!$ else if(interval == ib) then
+!!$ val=r1*r4*r10 + r2*r5*r10 + r2*r6*r11
+!!$ vald=r1d*r4*r10+r1*r4d*r10+r1*r4*r10d
+!!$ vald=vald+r2d*r5*r10+r2*r5d*r10+r2*r5*r10d
+!!$ vald=vald+r2d*r6*r11+r2*r6d*r11+r2*r6*r11d
+!!$ else if(interval == ib+1) then
+!!$ val=r2*r6*r12
+!!$ vald=r2d*r6*r12+r2*r6d*r12+r2*r6*r12d
+!!$ else
+!!$ val=0.
+!!$ endif
+!!$ endif
+!!$ splcon(ib)=val
+!!$ splcond(ib)=vald
+!!$ enddo
+!!$
+!!$ end subroutine vbspl
+
+!----------------------------------
+!----------------------------------
+
+!!$ subroutine ylm(XLAT,XLON,LMAX,Y,WK1,WK2,WK3)
+!!$
+!!$ implicit none
+!!$
+!!$ complex TEMP,FAC,DFAC
+!!$
+!!$ real(kind=4) WK1(1),WK2(1),WK3(1),Y(1),XLAT,XLON
+!!$
+!!$ integer :: LMAX
+!!$
+!!$!
+!!$! WK1,WK2,WK3 SHOULD BE DIMENSIONED AT LEAST (LMAX+1)*4
+!!$!
+!!$ real(kind=4), parameter :: RADIAN = 57.2957795
+!!$
+!!$ integer :: IM,IL1,IND,LM1,L
+!!$
+!!$ real(kind=4) :: THETA,PHI
+!!$
+!!$ THETA=(90.-XLAT)/RADIAN
+!!$ PHI=XLON/RADIAN
+!!$
+!!$ IND=0
+!!$ LM1=LMAX+1
+!!$
+!!$ DO IL1=1,LM1
+!!$
+!!$ L=IL1-1
+!!$ CALL legndr(THETA,L,L,WK1,WK2,WK3)
+!!$
+!!$ FAC=(1.,0.)
+!!$ DFAC=CEXP(CMPLX(0.,PHI))
+!!$
+!!$ do IM=1,IL1
+!!$ TEMP=FAC*CMPLX(WK1(IM),0.)
+!!$ IND=IND+1
+!!$ Y(IND)=REAL(TEMP)
+!!$ IF(IM == 1) GOTO 20
+!!$ IND=IND+1
+!!$ Y(IND)=AIMAG(TEMP)
+!!$ 20 FAC=FAC*DFAC
+!!$ enddo
+!!$
+!!$ enddo
+!!$
+!!$ end subroutine ylm
+
+!!$ subroutine legndr(THETA,L,M,X,XP,XCOSEC)
+!!$ implicit none
+!!$
+!!$ integer :: L,M,i,k,LP1,MP1
+!!$ real(kind=4) :: THETA,X,XP,XCOSEC,SFL3
+!!$
+!!$ DIMENSION X(2),XP(2),XCOSEC(2)
+!!$ DOUBLE PRECISION SMALL,SUM,COMPAR,CT,ST,FCT,COT,FPI,X1,X2,X3,F1,F2,XM,TH,DSFL3,COSEC
+!!$ DATA FPI/12.56637062D0/
+!!$! DFLOAT(I)=FLOAT(I)
+!!$ SUM=0.D0
+!!$ LP1=L+1
+!!$ TH=THETA
+!!$ CT=DCOS(TH)
+!!$ ST=DSIN(TH)
+!!$ MP1=M+1
+!!$ FCT=DSQRT(dble(FLOAT(2*L+1))/FPI)
+!!$ SFL3=SQRT(FLOAT(L*(L+1)))
+!!$ COMPAR=dble(FLOAT(2*L+1))/FPI
+!!$ DSFL3=SFL3
+!!$ SMALL=1.D-16*COMPAR
+!!$ do I=1,MP1
+!!$ X(I)=0.
+!!$ XCOSEC(I)=0.
+!!$ XP(I)=0.
+!!$ enddo
+!!$ IF(L.GT.1.AND.ABS(THETA).GT.1.E-5) GO TO 3
+!!$ X(1)=FCT
+!!$ IF(L.EQ.0) RETURN
+!!$ X(1)=CT*FCT
+!!$ X(2)=-ST*FCT/DSFL3
+!!$ XP(1)=-ST*FCT
+!!$ XP(2)=-.5D0*CT*FCT*DSFL3
+!!$ IF(ABS(THETA).LT.1.E-5) XCOSEC(2)=XP(2)
+!!$ IF(ABS(THETA).GE.1.E-5) XCOSEC(2)=X(2)/ST
+!!$ RETURN
+!!$ 3 X1=1.D0
+!!$ X2=CT
+!!$ DO I=2,L
+!!$ X3=(dble(FLOAT(2*I-1))*CT*X2-dble(FLOAT(I-1))*X1)/dble(FLOAT(I))
+!!$ X1=X2
+!!$ X2=X3
+!!$ enddo
+!!$ COT=CT/ST
+!!$ COSEC=1./ST
+!!$ X3=X2*FCT
+!!$ X2=dble(FLOAT(L))*(X1-CT*X2)*FCT/ST
+!!$ X(1)=X3
+!!$ X(2)=X2
+!!$ SUM=X3*X3
+!!$ XP(1)=-X2
+!!$ XP(2)=dble(FLOAT(L*(L+1)))*X3-COT*X2
+!!$ X(2)=-X(2)/SFL3
+!!$ XCOSEC(2)=X(2)*COSEC
+!!$ XP(2)=-XP(2)/SFL3
+!!$ SUM=SUM+2.D0*X(2)*X(2)
+!!$ IF(SUM-COMPAR.GT.SMALL) RETURN
+!!$ X1=X3
+!!$ X2=-X2/DSQRT(dble(FLOAT(L*(L+1))))
+!!$ DO I=3,MP1
+!!$ K=I-1
+!!$ F1=DSQRT(dble(FLOAT((L+I-1)*(L-I+2))))
+!!$ F2=DSQRT(dble(FLOAT((L+I-2)*(L-I+3))))
+!!$ XM=K
+!!$ X3=-(2.D0*COT*(XM-1.D0)*X2+F2*X1)/F1
+!!$ SUM=SUM+2.D0*X3*X3
+!!$ IF(SUM-COMPAR.GT.SMALL.AND.I.NE.LP1) RETURN
+!!$ X(I)=X3
+!!$ XCOSEC(I)=X(I)*COSEC
+!!$ X1=X2
+!!$ XP(I)=-(F1*X2+XM*COT*X3)
+!!$ X2=X3
+!!$ enddo
+!!$ RETURN
+!!$ end subroutine legndr
+
diff -r 000000000000 -r 9b613f27327a src/meshfem3D/model_crust.f90
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/meshfem3D/model_crust.f90 Sat Jul 02 17:09:21 2011 -0700
@@ -0,0 +1,742 @@
+!=====================================================================
+!
+! 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 University of Pau / CNRS / INRIA, France
+! (c) Princeton University / California Institute of Technology and University of Pau / CNRS / INRIA
+! 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.
+!
+!=====================================================================
+
+!--------------------------------------------------------------------------------------------------
+! CRUST 2.0 model by Bassin et al. (2000)
+!
+! C. Bassin, G. Laske, and G. Masters.
+! The current limits of resolution for surface wave tomography in North America.
+! EOS, 81: F897, 2000.
+!
+! reads and smooths crust2.0 model
+!--------------------------------------------------------------------------------------------------
+
+
+ subroutine model_crust_broadcast(myrank,CM_V)
+
+! standard routine to setup model
+
+ implicit none
+
+ include "constants.h"
+ ! standard include of the MPI library
+ include 'mpif.h'
+
+ ! model_crust_variables
+ type model_crust_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)
+ character(len=2) dummy_pad ! padding 2 bytes to align the structure
+ end type model_crust_variables
+
+ type (model_crust_variables) CM_V
+ ! model_crust_variables
+
+ integer :: myrank
+ integer :: ier
+
+ ! the variables read are declared and stored in structure CM_V
+ if(myrank == 0) call read_crust_model(CM_V)
+
+ ! broadcast the information read on the master to the nodes
+ call MPI_BCAST(CM_V%thlr,NKEYS_CRUST*NLAYERS_CRUST,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST(CM_V%velocp,NKEYS_CRUST*NLAYERS_CRUST,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST(CM_V%velocs,NKEYS_CRUST*NLAYERS_CRUST,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST(CM_V%dens,NKEYS_CRUST*NLAYERS_CRUST,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST(CM_V%abbreviation,NCAP_CRUST*NCAP_CRUST,MPI_CHARACTER,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST(CM_V%code,2*NKEYS_CRUST,MPI_CHARACTER,0,MPI_COMM_WORLD,ier)
+
+
+ end subroutine model_crust_broadcast
+
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+
+ subroutine model_crust(lat,lon,x,vp,vs,rho,moho,found_crust,CM_V,elem_in_crust)
+
+ implicit none
+ include "constants.h"
+
+! model_crust_variables
+ type model_crust_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)
+ character(len=2) dummy_pad ! padding 2 bytes to align the structure
+ end type model_crust_variables
+
+ type (model_crust_variables) CM_V
+! model_crust_variables
+
+ double precision lat,lon,x,vp,vs,rho,moho
+ logical found_crust,elem_in_crust
+
+ ! local parameters
+ double precision h_sed,h_uc
+ double precision x3,x4,x5,x6,x7,scaleval
+ double precision vps(NLAYERS_CRUST),vss(NLAYERS_CRUST),rhos(NLAYERS_CRUST),thicks(NLAYERS_CRUST)
+
+ ! initializes
+ vp = 0.d0
+ vs = 0.d0
+ rho = 0.d0
+
+ ! gets smoothed crust2.0 structure
+ call crust_CAPsmoothed(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
+ h_sed = thicks(3) + thicks(4)
+ x4 = (R_EARTH-h_sed*1000.0d0)/R_EARTH
+ h_uc = h_sed + thicks(5)
+ x5 = (R_EARTH-h_uc*1000.0d0)/R_EARTH
+ x6 = (R_EARTH-(h_uc+thicks(6))*1000.0d0)/R_EARTH
+ x7 = (R_EARTH-(h_uc+thicks(6)+thicks(7))*1000.0d0)/R_EARTH
+
+ found_crust = .true.
+
+! if(x > x3 .and. INCLUDE_SEDIMENTS_CRUST &
+! .and. h_sed >= MINIMUM_SEDIMENT_THICKNESS) then
+ if(x > x3 .and. INCLUDE_SEDIMENTS_CRUST ) then
+ vp = vps(3)
+ vs = vss(3)
+ rho = rhos(3)
+! else if(x > x4 .and. INCLUDE_SEDIMENTS_CRUST &
+! .and. h_sed >= MINIMUM_SEDIMENT_THICKNESS) then
+ else if(x > x4 .and. INCLUDE_SEDIMENTS_CRUST ) then
+ vp = vps(4)
+ vs = vss(4)
+ rho = rhos(4)
+ else if(x > x5) then
+ vp = vps(5)
+ vs = vss(5)
+ rho = rhos(5)
+ else if(x > x6) then
+ vp = vps(6)
+ vs = vss(6)
+ rho = rhos(6)
+ else if(x > x7 .or. elem_in_crust) then
+ ! takes lower crustal values only if x is slightly above moho depth or
+ ! if elem_in_crust is set
+ !
+ ! note: it looks like this does distinguish between GLL points at the exact moho boundary,
+ ! where the point is on the interface between both,
+ ! oceanic elements and mantle elements below
+ vp = vps(7)
+ vs = vss(7)
+ rho = rhos(7)
+ else
+ ! note: if x is exactly the moho depth this will return false
+ found_crust = .false.
+ endif
+
+ ! non-dimensionalize
+ if (found_crust) then
+ scaleval = dsqrt(PI*GRAV*RHOAV)
+ vp = vp*1000.0d0/(R_EARTH*scaleval)
+ vs = vs*1000.0d0/(R_EARTH*scaleval)
+ rho = rho*1000.0d0/RHOAV
+ endif
+
+ ! checks moho value
+ !moho = h_uc + thicks(6) + thicks(7)
+ !if( moho /= thicks(NLAYERS_CRUST) ) then
+ ! print*,'moho:',moho,thicks(NLAYERS_CRUST)
+ ! print*,' lat/lon/x:',lat,lon,x
+ !endif
+
+ ! No matter found_crust true or false, output moho thickness
+ moho = (h_uc+thicks(6)+thicks(7))*1000.0d0/R_EARTH
+
+ end subroutine model_crust
+
+!---------------------------
+
+ subroutine read_crust_model(CM_V)
+
+ implicit none
+ include "constants.h"
+
+! model_crust_variables
+ type model_crust_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)
+ character(len=2) dummy_pad ! padding 2 bytes to align the structure
+ end type model_crust_variables
+
+ type (model_crust_variables) CM_V
+! model_crust_variables
+
+! local variables
+ integer i
+ integer ila,icolat
+ integer ikey
+
+ double precision h_moho_min,h_moho_max
+
+ character(len=150) CNtype2, CNtype2_key_modif
+
+ call get_value_string(CNtype2, 'model.CNtype2', 'DATA/crust2.0/CNtype2.txt')
+ call get_value_string(CNtype2_key_modif, 'model.CNtype2_key_modif', 'DATA/crust2.0/CNtype2_key_modif.txt')
+
+ open(unit=1,file=CNtype2,status='old',action='read')
+ do ila=1,NCAP_CRUST/2
+ read(1,*) icolat,(CM_V%abbreviation(ila,i),i=1,NCAP_CRUST)
+ enddo
+ close(1)
+
+ open(unit=1,file=CNtype2_key_modif,status='old',action='read')
+ h_moho_min=HUGEVAL
+ h_moho_max=-HUGEVAL
+ do ikey=1,NKEYS_CRUST
+ read (1,"(a2)") CM_V%code(ikey)
+ read (1,*) (CM_V%velocp(ikey,i),i=1,NLAYERS_CRUST)
+ read (1,*) (CM_V%velocs(ikey,i),i=1,NLAYERS_CRUST)
+ read (1,*) (CM_V%dens(ikey,i),i=1,NLAYERS_CRUST)
+ read (1,*) (CM_V%thlr(ikey,i),i=1,NLAYERS_CRUST-1),CM_V%thlr(ikey,NLAYERS_CRUST)
+ if(CM_V%thlr(ikey,NLAYERS_CRUST) > h_moho_max) h_moho_max=CM_V%thlr(ikey,NLAYERS_CRUST)
+ if(CM_V%thlr(ikey,NLAYERS_CRUST) < h_moho_min) h_moho_min=CM_V%thlr(ikey,NLAYERS_CRUST)
+ enddo
+ close(1)
+
+ if(h_moho_min == HUGEVAL .or. h_moho_max == -HUGEVAL) &
+ stop 'incorrect moho depths in read_crust_model'
+
+ end subroutine read_crust_model
+
+!---------------------------
+
+ subroutine crust_CAPsmoothed(lat,lon,velp,vels,rho,thick,abbreviation,&
+ code,thlr,velocp,velocs,dens)
+
+! crustal vp and vs in km/s, layer thickness in km
+!
+! crust2.0 is smoothed with a cap of size CAP using NTHETA points
+! in the theta direction and NPHI in the phi direction.
+! The cap is rotated to the North Pole.
+
+ implicit none
+ include "constants.h"
+
+ ! sampling rate for CAP points
+ integer, parameter :: NTHETA = 4
+ integer, parameter :: NPHI = 20
+
+ ! argument variables
+ double precision lat,lon
+ double precision rho(NLAYERS_CRUST),thick(NLAYERS_CRUST),velp(NLAYERS_CRUST),vels(NLAYERS_CRUST)
+ double precision thlr(NKEYS_CRUST,NLAYERS_CRUST),velocp(NKEYS_CRUST,NLAYERS_CRUST)
+ double precision velocs(NKEYS_CRUST,NLAYERS_CRUST),dens(NKEYS_CRUST,NLAYERS_CRUST)
+ character(len=2) code(NKEYS_CRUST),abbreviation(NCAP_CRUST/2,NCAP_CRUST)
+
+ !-------------------------------
+ ! work-around to avoid jacobian problems when stretching mesh elements;
+ ! one could also try to slightly change the shape of the doulbing element bricks (which cause the problem)...
+ !
+ ! defines a "critical" region around the andes to have at least a 2-degree smoothing;
+ ! critical region can lead to negative jacobians for mesh stretching when CAP smoothing is too small
+ double precision,parameter :: LAT_CRITICAL_ANDES = -20.0d0
+ double precision,parameter :: LON_CRITICAL_ANDES = -70.0d0
+ double precision,parameter :: CRITICAL_RANGE = 70.0d0
+ !-------------------------------
+
+ ! local variables
+ double precision xlon(NTHETA*NPHI),xlat(NTHETA*NPHI),weight(NTHETA*NPHI)
+ double precision rhol(NLAYERS_CRUST),thickl(NLAYERS_CRUST),velpl(NLAYERS_CRUST),velsl(NLAYERS_CRUST)
+ double precision weightl,cap_degree,dist
+ double precision h_sed
+ integer i,icolat,ilon,ierr
+ character(len=2) crustaltype
+
+ ! checks latitude/longitude
+ if(lat > 90.0d0 .or. lat < -90.0d0 .or. lon > 180.0d0 .or. lon < -180.0d0) &
+ stop 'error in latitude/longitude range in crust'
+
+ ! makes sure lat/lon are within crust2.0 range
+ if(lat==90.0d0) lat=89.9999d0
+ if(lat==-90.0d0) lat=-89.9999d0
+ if(lon==180.0d0) lon=179.9999d0
+ if(lon==-180.0d0) lon=-179.9999d0
+
+ ! sets up smoothing points
+ ! by default uses CAP smoothing with 1 degree
+ cap_degree = 1.0d0
+
+ ! checks if inside/outside of critical region for mesh stretching
+ if( SMOOTH_CRUST ) then
+ dist = dsqrt( (lon-LON_CRITICAL_ANDES)**2 + (lat-LAT_CRITICAL_ANDES )**2 )
+ if( dist < CRITICAL_RANGE ) then
+ ! increases cap smoothing degree
+ ! scales between -1 at center and 0 at border
+ dist = dist / CRITICAL_RANGE - 1.0d0
+ ! shifts value to 1 at center and 0 to the border with exponential decay
+ dist = 1.0d0 - exp( - dist*dist*10.0d0 )
+ ! increases smoothing degree inside of critical region to 2 degree
+ cap_degree = cap_degree + dist
+ endif
+ endif
+
+ ! gets smoothing points and weights
+ call CAP_vardegree(lon,lat,xlon,xlat,weight,cap_degree,NTHETA,NPHI)
+
+ ! initializes
+ velp(:) = 0.0d0
+ vels(:) = 0.0d0
+ rho(:) = 0.0d0
+ thick(:) = 0.0d0
+
+ ! loops over weight points
+ do i=1,NTHETA*NPHI
+ ! gets crust values
+ call icolat_ilon(xlat(i),xlon(i),icolat,ilon)
+ crustaltype = abbreviation(icolat,ilon)
+ call get_crust_structure(crustaltype,velpl,velsl,rhol,thickl, &
+ code,thlr,velocp,velocs,dens,ierr)
+ if(ierr /= 0) stop 'error in routine get_crust_structure'
+
+ ! sediment thickness
+ h_sed = thickl(3) + thickl(4)
+
+ ! takes upper crust value if sediment too thin
+ if( h_sed < MINIMUM_SEDIMENT_THICKNESS ) then
+ velpl(3) = velpl(5)
+ velpl(4) = velpl(5)
+ velsl(3) = velsl(5)
+ velsl(4) = velsl(5)
+ rhol(3) = rhol(5)
+ rhol(4) = rhol(5)
+ endif
+
+ ! weighting value
+ weightl = weight(i)
+
+ ! total, smoothed values
+ rho(:) = rho(:) + weightl*rhol(:)
+ thick(:) = thick(:) + weightl*thickl(:)
+ velp(:) = velp(:) + weightl*velpl(:)
+ vels(:) = vels(:) + weightl*velsl(:)
+ enddo
+
+ end subroutine crust_CAPsmoothed
+
+
+!------------------------------------------------------
+
+ subroutine icolat_ilon(xlat,xlon,icolat,ilon)
+
+ implicit none
+
+
+! argument variables
+ double precision xlat,xlon
+ integer icolat,ilon
+
+ if(xlat > 90.0d0 .or. xlat < -90.0d0 .or. xlon > 180.0d0 .or. xlon < -180.0d0) &
+ stop 'error in latitude/longitude range in icolat_ilon'
+ icolat=int(1+((90.d0-xlat)/2.d0))
+ if(icolat == 91) icolat=90
+ ilon=int(1+((180.d0+xlon)/2.d0))
+ if(ilon == 181) ilon=1
+
+ if(icolat>90 .or. icolat<1) stop 'error in routine icolat_ilon'
+ if(ilon<1 .or. ilon>180) stop 'error in routine icolat_ilon'
+
+ end subroutine icolat_ilon
+
+!---------------------------------------------------------------------
+
+ subroutine get_crust_structure(type,vptyp,vstyp,rhtyp,thtp, &
+ code,thlr,velocp,velocs,dens,ierr)
+
+ implicit none
+ include "constants.h"
+
+
+! argument variables
+ integer ierr
+ double precision rhtyp(NLAYERS_CRUST),thtp(NLAYERS_CRUST)
+ double precision vptyp(NLAYERS_CRUST),vstyp(NLAYERS_CRUST)
+ character(len=2) type,code(NKEYS_CRUST)
+ double precision thlr(NKEYS_CRUST,NLAYERS_CRUST),velocp(NKEYS_CRUST,NLAYERS_CRUST)
+ double precision velocs(NKEYS_CRUST,NLAYERS_CRUST),dens(NKEYS_CRUST,NLAYERS_CRUST)
+
+! local variables
+ integer i,ikey
+
+ ierr=1
+ do ikey=1,NKEYS_CRUST
+ if (code(ikey) == type) then
+ do i=1,NLAYERS_CRUST
+ vptyp(i)=velocp(ikey,i)
+ vstyp(i)=velocs(ikey,i)
+ rhtyp(i)=dens(ikey,i)
+ enddo
+ do i=1,NLAYERS_CRUST-1
+ thtp(i)=thlr(ikey,i)
+ enddo
+ ! get distance to Moho from the bottom of the ocean or the ice
+ thtp(NLAYERS_CRUST)=thlr(ikey,NLAYERS_CRUST)-thtp(1)-thtp(2)
+ ierr=0
+ endif
+ enddo
+
+ end subroutine get_crust_structure
+
+
+!---------------------------
+
+ subroutine CAP_vardegree(lon,lat,xlon,xlat,weight,CAP_DEGREE,NTHETA,NPHI)
+
+! calculates weighting points to smooth around lon/lat location with
+! a smoothing range of CAP_DEGREE
+!
+! The cap is rotated to the North Pole.
+!
+! returns: xlon,xlat,weight
+
+ implicit none
+ include "constants.h"
+
+ ! sampling rate
+ integer :: NTHETA
+ integer :: NPHI
+ ! smoothing size (in degrees)
+ double precision :: CAP_DEGREE
+
+ ! argument variables
+ double precision lat,lon
+ double precision xlon(NTHETA*NPHI),xlat(NTHETA*NPHI),weight(NTHETA*NPHI)
+
+ ! local variables
+ double precision CAP
+ double precision theta,phi,sint,cost,sinp,cosp,wght,total
+ double precision r_rot,theta_rot,phi_rot
+ double precision rotation_matrix(3,3),x(3),xc(3)
+ double precision dtheta,dphi,cap_area,dweight,pi_over_nphi
+ integer i,j,k
+ integer itheta,iphi
+
+ double precision, parameter :: RADIANS_TO_DEGREES = 180.d0 / PI
+ double precision, parameter :: PI_OVER_TWO = PI / 2.0d0
+
+ ! initializes
+ xlon(:) = 0.d0
+ xlat(:) = 0.d0
+ weight(:) = 0.d0
+
+ ! checks cap degree size
+ if( CAP_DEGREE < TINYVAL ) then
+ ! no cap smoothing
+ print*,'error cap:',CAP_DEGREE
+ print*,' lat/lon:',lat,lon
+ stop 'error cap_degree too small'
+ endif
+
+ ! pre-compute parameters
+ CAP = CAP_DEGREE * PI/180.0d0
+ dtheta = 0.5d0 * CAP / dble(NTHETA)
+ dphi = TWO_PI / dble(NPHI)
+ cap_area = TWO_PI * (1.0d0 - dcos(CAP))
+ dweight = CAP / dble(NTHETA) * dphi / cap_area
+ pi_over_nphi = PI/dble(NPHI)
+
+ ! colatitude/longitude in radian
+ theta = (90.0d0 - lat ) * DEGREES_TO_RADIANS
+ phi = lon * DEGREES_TO_RADIANS
+
+ sint = dsin(theta)
+ cost = dcos(theta)
+ sinp = dsin(phi)
+ cosp = dcos(phi)
+
+ ! set up rotation matrix to go from cap at North pole
+ ! to cap around point of interest
+ rotation_matrix(1,1) = cosp*cost
+ rotation_matrix(1,2) = -sinp
+ rotation_matrix(1,3) = cosp*sint
+ rotation_matrix(2,1) = sinp*cost
+ rotation_matrix(2,2) = cosp
+ rotation_matrix(2,3) = sinp*sint
+ rotation_matrix(3,1) = -sint
+ rotation_matrix(3,2) = 0.0d0
+ rotation_matrix(3,3) = cost
+
+ ! calculates points over a cap at the North pole and rotates them to specified lat/lon point
+ i = 0
+ total = 0.0d0
+ do itheta = 1,NTHETA
+
+ theta = dble(2*itheta-1)*dtheta
+ cost = dcos(theta)
+ sint = dsin(theta)
+ wght = sint*dweight
+
+ do iphi = 1,NPHI
+
+ i = i+1
+
+ ! get the weight associated with this integration point (same for all phi)
+ weight(i) = wght
+
+ total = total + weight(i)
+ phi = dble(2*iphi-1)*pi_over_nphi
+ cosp = dcos(phi)
+ sinp = dsin(phi)
+
+ ! x,y,z coordinates of integration point in cap at North pole
+ xc(1) = sint*cosp
+ xc(2) = sint*sinp
+ xc(3) = cost
+
+ ! get x,y,z coordinates in cap around point of interest
+ do j=1,3
+ x(j) = 0.0d0
+ do k=1,3
+ x(j) = x(j)+rotation_matrix(j,k)*xc(k)
+ enddo
+ enddo
+
+ ! get latitude and longitude (degrees) of integration point
+ call xyz_2_rthetaphi_dble(x(1),x(2),x(3),r_rot,theta_rot,phi_rot)
+ call reduce(theta_rot,phi_rot)
+ xlat(i) = (PI_OVER_TWO - theta_rot) * RADIANS_TO_DEGREES
+ xlon(i) = phi_rot * RADIANS_TO_DEGREES
+ if(xlon(i) > 180.0d0) xlon(i) = xlon(i) - 360.0d0
+
+ enddo
+
+ enddo
+ if(abs(total-1.0d0) > 0.001d0) then
+ print*,'error cap:',total,CAP_DEGREE
+ stop 'error in cap integration for variable degree'
+ endif
+
+ end subroutine
+
+
+!---------------------------
+! unused routines...
+!
+! subroutine crust_singlevalue(lat,lon,velp,vels,rho,thick,abbreviation,&
+! code,thlr,velocp,velocs,dens)
+!
+!! crustal vp and vs in km/s, layer thickness in km
+!
+!! uses crust2.0 as is, without smoothing
+!
+! implicit none
+! include "constants.h"
+!
+!! argument variables
+! double precision lat,lon
+! double precision rho(NLAYERS_CRUST),thick(NLAYERS_CRUST),velp(NLAYERS_CRUST),vels(NLAYERS_CRUST)
+! double precision thlr(NKEYS_CRUST,NLAYERS_CRUST),velocp(NKEYS_CRUST,NLAYERS_CRUST)
+! double precision velocs(NKEYS_CRUST,NLAYERS_CRUST),dens(NKEYS_CRUST,NLAYERS_CRUST)
+! character(len=2) code(NKEYS_CRUST),abbreviation(NCAP_CRUST/2,NCAP_CRUST)
+!
+!! local variables
+! integer icolat,ilon,ierr
+! character(len=2) crustaltype
+!
+!
+!! get integer colatitude and longitude of crustal cap
+!! -90<lat<90 -180<lon<180
+! if(lat > 90.0d0 .or. lat < -90.0d0 .or. lon > 180.0d0 .or. lon < -180.0d0) &
+! stop 'error in latitude/longitude range in crust'
+! if(lat==90.0d0) lat=89.9999d0
+! if(lat==-90.0d0) lat=-89.9999d0
+! if(lon==180.0d0) lon=179.9999d0
+! if(lon==-180.0d0) lon=-179.9999d0
+!
+! call icolat_ilon(lat,lon,icolat,ilon)
+! crustaltype = abbreviation(icolat,ilon)
+! call get_crust_structure(crustaltype,velp,vels,rho,thick, &
+! code,thlr,velocp,velocs,dens,ierr)
+! if( ierr /= 0 ) stop 'error in routine get_crust_structure'
+!
+! end subroutine crust_singlevalue
+!
+!---------------------------
+!
+!
+! subroutine crust_org(lat,lon,velp,vels,rho,thick,abbreviation,code,thlr,velocp,velocs,dens)
+!
+!! crustal vp and vs in km/s, layer thickness in km
+!! crust2.0 is smoothed with a cap of size CAP using NTHETA points
+!! in the theta direction and NPHI in the phi direction.
+!! The cap is rotated to the North Pole.
+!
+! implicit none
+! include "constants.h"
+!! 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
+!
+!! argument variables
+! double precision lat,lon
+! double precision rho(NLAYERS_CRUST),thick(NLAYERS_CRUST),velp(NLAYERS_CRUST),vels(NLAYERS_CRUST)
+! double precision thlr(NKEYS_CRUST,NLAYERS_CRUST),velocp(NKEYS_CRUST,NLAYERS_CRUST)
+! double precision velocs(NKEYS_CRUST,NLAYERS_CRUST),dens(NKEYS_CRUST,NLAYERS_CRUST)
+! character(len=2) code(NKEYS_CRUST),abbreviation(NCAP_CRUST/2,NCAP_CRUST)
+!
+!! local variables
+! integer i,j,k,icolat,ilon,ierr
+! integer itheta,iphi,npoints
+! double precision theta,phi,sint,cost,sinp,cosp,dtheta,dphi,cap_area,wght,total
+! double precision r_rot,theta_rot,phi_rot
+! double precision rotation_matrix(3,3),x(3),xc(3)
+! double precision xlon(NTHETA*NPHI),xlat(NTHETA*NPHI),weight(NTHETA*NPHI)
+! double precision rhol(NLAYERS_CRUST),thickl(NLAYERS_CRUST),velpl(NLAYERS_CRUST),velsl(NLAYERS_CRUST)
+! character(len=2) crustaltype
+!
+!! get integer colatitude and longitude of crustal cap
+!! -90<lat<90 -180<lon<180
+! if(lat > 90.0d0 .or. lat < -90.0d0 .or. lon > 180.0d0 .or. lon < -180.0d0) &
+! stop 'error in latitude/longitude range in crust'
+! if(lat==90.0d0) lat=89.9999d0
+! if(lat==-90.0d0) lat=-89.9999d0
+! if(lon==180.0d0) lon=179.9999d0
+! if(lon==-180.0d0) lon=-179.9999d0
+!
+! call icolat_ilon(lat,lon,icolat,ilon)
+! crustaltype=abbreviation(icolat,ilon)
+! call get_crust_structure(crustaltype,velp,vels,rho,thick, &
+! code,thlr,velocp,velocs,dens,ierr)
+!
+!! uncomment the following line to use crust2.0 as is, without smoothing
+!!
+!! return
+!
+! theta = (90.0-lat)*PI/180.0
+! phi = lon*PI/180.0
+!
+! sint = sin(theta)
+! cost = cos(theta)
+! sinp = sin(phi)
+! cosp = cos(phi)
+!
+!! set up rotation matrix to go from cap at North pole
+!! to cap around point of interest
+! rotation_matrix(1,1) = cosp*cost
+! rotation_matrix(1,2) = -sinp
+! rotation_matrix(1,3) = cosp*sint
+! rotation_matrix(2,1) = sinp*cost
+! rotation_matrix(2,2) = cosp
+! rotation_matrix(2,3) = sinp*sint
+! rotation_matrix(3,1) = -sint
+! rotation_matrix(3,2) = 0.0
+! rotation_matrix(3,3) = cost
+!
+! dtheta = CAP/dble(NTHETA)
+! dphi = 2.0*PI/dble(NPHI)
+! cap_area = 2.0*PI*(1.0-cos(CAP))
+!
+!! integrate over a cap at the North pole
+! i = 0
+! total = 0.0
+! do itheta = 1,NTHETA
+!
+! theta = 0.5*dble(2*itheta-1)*CAP/dble(NTHETA)
+! cost = cos(theta)
+! sint = sin(theta)
+! wght = sint*dtheta*dphi/cap_area
+!
+! do iphi = 1,NPHI
+!
+! i = i+1
+!! get the weight associated with this integration point (same for all phi)
+! weight(i) = wght
+! total = total + weight(i)
+! phi = dble(2*iphi-1)*PI/dble(NPHI)
+! cosp = cos(phi)
+! sinp = sin(phi)
+!! x,y,z coordinates of integration point in cap at North pole
+! xc(1) = sint*cosp
+! xc(2) = sint*sinp
+! xc(3) = cost
+!! get x,y,z coordinates in cap around point of interest
+! do j=1,3
+! x(j) = 0.0
+! do k=1,3
+! x(j) = x(j)+rotation_matrix(j,k)*xc(k)
+! enddo
+! enddo
+!! get latitude and longitude (degrees) of integration point
+! call xyz_2_rthetaphi_dble(x(1),x(2),x(3),r_rot,theta_rot,phi_rot)
+! call reduce(theta_rot,phi_rot)
+! xlat(i) = (PI/2.0-theta_rot)*180.0/PI
+! xlon(i) = phi_rot*180.0/PI
+! if(xlon(i) > 180.0) xlon(i) = xlon(i)-360.0
+!
+! enddo
+!
+! enddo
+!
+! if(abs(total-1.0) > 0.001) stop 'error in cap integration for crust2.0'
+!
+! npoints = i
+!
+! do j=1,NLAYERS_CRUST
+! rho(j)=0.0d0
+! thick(j)=0.0d0
+! velp(j)=0.0d0
+! vels(j)=0.0d0
+! enddo
+!
+! do i=1,npoints
+! call icolat_ilon(xlat(i),xlon(i),icolat,ilon)
+! crustaltype=abbreviation(icolat,ilon)
+! call get_crust_structure(crustaltype,velpl,velsl,rhol,thickl, &
+! code,thlr,velocp,velocs,dens,ierr)
+! if(ierr /= 0) stop 'error in routine get_crust_structure'
+! do j=1,NLAYERS_CRUST
+! rho(j)=rho(j)+weight(i)*rhol(j)
+! thick(j)=thick(j)+weight(i)*thickl(j)
+! velp(j)=velp(j)+weight(i)*velpl(j)
+! vels(j)=vels(j)+weight(i)*velsl(j)
+! enddo
+! enddo
+!
+! end subroutine crust_org
+
diff -r 000000000000 -r 9b613f27327a src/meshfem3D/model_s20rts.f90
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/src/meshfem3D/model_s20rts.f90 Sat Jul 02 17:09:21 2011 -0700
@@ -0,0 +1,515 @@
+!=====================================================================
+!
+! 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 University of Pau / CNRS / INRIA, France
+! (c) Princeton University / California Institute of Technology and University of Pau / CNRS / INRIA
+! 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.
+!
+!=====================================================================
+
+!--------------------------------------------------------------------------------------------------
+! S20rts
+!
+! 3D mantle model S20RTS [Ritsema et al., 1999]
+!
+! Note that S20RTS uses transversely isotropic PREM as a background
+! model, and that we use the PREM radial attenuation model when ATTENUATION is incorporated.
+!--------------------------------------------------------------------------------------------------
+
+
+ subroutine model_s20rts_broadcast(myrank,S20RTS_V)
+
+! standard routine to setup model
+
+ implicit none
+
+ include "constants.h"
+ ! standard include of the MPI library
+ include 'mpif.h'
+
+! model_s20rts_variables s20rts
+ type model_s20rts_variables
+ sequence
+ double precision dvs_a(0:NK_20,0:NS_20,0:NS_20)
+ double precision dvs_b(0:NK_20,0:NS_20,0:NS_20)
+ double precision dvp_a(0:NK_20,0:NS_20,0:NS_20)
+ double precision dvp_b(0:NK_20,0:NS_20,0:NS_20)
+ double precision spknt(NK_20+1)
+ double precision qq0(NK_20+1,NK_20+1)
+ double precision qq(3,NK_20+1,NK_20+1)
+ end type model_s20rts_variables
+
+ type (model_s20rts_variables) S20RTS_V
+! model_s20rts_variables
+
+ integer :: myrank
+ integer :: ier
+
+ ! the variables read are declared and stored in structure S20RTS_V
+ if(myrank == 0) call read_model_s20rts(S20RTS_V)
+
+ ! broadcast the information read on the master to the nodes
+ call MPI_BCAST(S20RTS_V%dvs_a,(NK_20+1)*(NS_20+1)*(NS_20+1),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST(S20RTS_V%dvs_b,(NK_20+1)*(NS_20+1)*(NS_20+1),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST(S20RTS_V%dvp_a,(NK_20+1)*(NS_20+1)*(NS_20+1),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST(S20RTS_V%dvp_b,(NK_20+1)*(NS_20+1)*(NS_20+1),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST(S20RTS_V%spknt,NK_20+1,MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST(S20RTS_V%qq0,(NK_20+1)*(NK_20+1),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+ call MPI_BCAST(S20RTS_V%qq,3*(NK_20+1)*(NK_20+1),MPI_DOUBLE_PRECISION,0,MPI_COMM_WORLD,ier)
+
+ end subroutine model_s20rts_broadcast
+!
+!-------------------------------------------------------------------------------------------------
+!
+
+ subroutine read_model_s20rts(S20RTS_V)
+
+ implicit none
+
+ include "constants.h"
+
+! model_s20rts_variables
+ type model_s20rts_variables
+ sequence
+ double precision dvs_a(0:NK_20,0:NS_20,0:NS_20)
+ double precision dvs_b(0:NK_20,0:NS_20,0:NS_20)
+ double precision dvp_a(0:NK_20,0:NS_20,0:NS_20)
+ double precision dvp_b(0:NK_20,0:NS_20,0:NS_20)
+ double precision spknt(NK_20+1)
+ double precision qq0(NK_20+1,NK_20+1)
+ double precision qq(3,NK_20+1,NK_20+1)
+ end type model_s20rts_variables
+
+ type (model_s20rts_variables) S20RTS_V
+! model_s20rts_variables
+
+ integer k,l,m
+
+ character(len=150) S20RTS, P12
+
+ call get_value_string(S20RTS, 'model.S20RTS', 'DATA/s20rts/S20RTS.dat')
+ call get_value_string(P12, 'model.P12', 'DATA/s20rts/P12.dat')
+
+! S20RTS degree 20 S model from Ritsema
+ open(unit=10,file=S20RTS,status='old',action='read')
+ do k=0,NK_20
+ do l=0,NS_20
+ read(10,*) S20RTS_V%dvs_a(k,l,0),(S20RTS_V%dvs_a(k,l,m),S20RTS_V%dvs_b(k,l,m),m=1,l)
+ enddo
+ enddo
+ close(10)
+
+! P12 degree 12 P model from Ritsema
+ open(unit=10,file=P12,status='old',action='read')
+ do k=0,NK_20
+ do l=0,12
+ read(10,*) S20RTS_V%dvp_a(k,l,0),(S20RTS_V%dvp_a(k,l,m),S20RTS_V%dvp_b(k,l,m),m=1,l)
+ enddo
+ do l=13,NS_20
+ S20RTS_V%dvp_a(k,l,0) = 0.0d0
+ do m=1,l
+ S20RTS_V%dvp_a(k,l,m) = 0.0d0
+ S20RTS_V%dvp_b(k,l,m) = 0.0d0
+ enddo
+ enddo
+ enddo
+ close(10)
+
+! set up the splines used as radial basis functions by Ritsema
+ call s20rts_splhsetup(S20RTS_V)
+
+ end subroutine read_model_s20rts
+
+!---------------------------
+
+ subroutine mantle_s20rts(radius,theta,phi,dvs,dvp,drho,S20RTS_V)
+
+ implicit none
+
+ include "constants.h"
+
+! model_s20rts_variables
+ type model_s20rts_variables
+ sequence
+ double precision dvs_a(0:NK_20,0:NS_20,0:NS_20)
+ double precision dvs_b(0:NK_20,0:NS_20,0:NS_20)
+ double precision dvp_a(0:NK_20,0:NS_20,0:NS_20)
+ double precision dvp_b(0:NK_20,0:NS_20,0:NS_20)
+ double precision spknt(NK_20+1)
+ double precision qq0(NK_20+1,NK_20+1)
+ double precision qq(3,NK_20+1,NK_20+1)
+ end type model_s20rts_variables
+
+ type (model_s20rts_variables) S20RTS_V
+! model_s20rts_variables
+
+! factor to convert perturbations in shear speed to perturbations in density
+ double precision, parameter :: SCALE_RHO = 0.40d0
+
+ double precision radius,theta,phi,dvs,dvp,drho
+
+ double precision, parameter :: RMOHO_ = 6346600.d0
+ double precision, parameter :: RCMB_ = 3480000.d0
+ double precision, parameter :: R_EARTH_ = 6371000.d0
+ double precision, parameter :: ZERO_ = 0.d0
+
+ integer l,m,k
+ double precision r_moho,r_cmb,xr
+ double precision dvs_alm,dvs_blm
+ double precision dvp_alm,dvp_blm
+ double precision s20rts_rsple,radial_basis(0:NK_20)
+ double precision sint,cost,x(2*NS_20+1),dx(2*NS_20+1)
+
+ dvs = ZERO_
+ dvp = ZERO_
+ drho = ZERO_
+
+ r_moho = RMOHO_ / R_EARTH_
+ r_cmb = RCMB_ / R_EARTH_
+ if(radius>=r_moho .or. radius <= r_cmb) return
+
+ xr=-1.0d0+2.0d0*(radius-r_cmb)/(r_moho-r_cmb)
+ do k=0,NK_20
+ radial_basis(k)=s20rts_rsple(1,NK_20+1,S20RTS_V%spknt(1),S20RTS_V%qq0(1,NK_20+1-k),S20RTS_V%qq(1,1,NK_20+1-k),xr)
+ enddo
+
+ do l=0,NS_20
+ sint=dsin(theta)
+ cost=dcos(theta)
+ call lgndr(l,cost,sint,x,dx)
+
+ dvs_alm=0.0d0
+ dvp_alm=0.0d0
+ do k=0,NK_20
+ dvs_alm=dvs_alm+radial_basis(k)*S20RTS_V%dvs_a(k,l,0)
+ dvp_alm=dvp_alm+radial_basis(k)*S20RTS_V%dvp_a(k,l,0)
+ enddo
+ dvs=dvs+dvs_alm*x(1)
+ dvp=dvp+dvp_alm*x(1)
+
+ do m=1,l
+ dvs_alm=0.0d0
+ dvp_alm=0.0d0
+ dvs_blm=0.0d0
+ dvp_blm=0.0d0
+ do k=0,NK_20
+ dvs_alm=dvs_alm+radial_basis(k)*S20RTS_V%dvs_a(k,l,m)
+ dvp_alm=dvp_alm+radial_basis(k)*S20RTS_V%dvp_a(k,l,m)
+ dvs_blm=dvs_blm+radial_basis(k)*S20RTS_V%dvs_b(k,l,m)
+ dvp_blm=dvp_blm+radial_basis(k)*S20RTS_V%dvp_b(k,l,m)
+ enddo
+ dvs=dvs+(dvs_alm*dcos(dble(m)*phi)+dvs_blm*dsin(dble(m)*phi))*x(m+1)
+ dvp=dvp+(dvp_alm*dcos(dble(m)*phi)+dvp_blm*dsin(dble(m)*phi))*x(m+1)
+ enddo
+
+ enddo
+
+ drho = SCALE_RHO*dvs
+
+ end subroutine mantle_s20rts
+
+!----------------------------------
+
+ subroutine s20rts_splhsetup(S20RTS_V)!!!!!!!!!!!!!!(spknt,qq0,qq)
+
+ implicit none
+ include "constants.h"
+
+!!!!!!!!!!!!!!!!!!! double precision spknt(NK_20+1),qq0(NK_20+1,NK_20+1),qq(3,NK_20+1,NK_20+1)
+
+! model_s20rts_variables
+ type model_s20rts_variables
+ sequence
+ double precision dvs_a(0:NK_20,0:NS_20,0:NS_20)
+ double precision dvs_b(0:NK_20,0:NS_20,0:NS_20)
+ double precision dvp_a(0:NK_20,0:NS_20,0:NS_20)
+ double precision dvp_b(0:NK_20,0:NS_20,0:NS_20)
+ double precision spknt(NK_20+1)
+ double precision qq0(NK_20+1,NK_20+1)
+ double precision qq(3,NK_20+1,NK_20+1)
+ end type model_s20rts_variables
+
+ type (model_s20rts_variables) S20RTS_V
+! model_s20rts_variables
+
+
+ integer i,j
+ double precision qqwk(3,NK_20+1)
+
+ S20RTS_V%spknt(1) = -1.00000d0
+ S20RTS_V%spknt(2) = -0.78631d0
+ S20RTS_V%spknt(3) = -0.59207d0
+ S20RTS_V%spknt(4) = -0.41550d0
+ S20RTS_V%spknt(5) = -0.25499d0
+ S20RTS_V%spknt(6) = -0.10909d0
+ S20RTS_V%spknt(7) = 0.02353d0
+ S20RTS_V%spknt(8) = 0.14409d0
+ S20RTS_V%spknt(9) = 0.25367d0
+ S20RTS_V%spknt(10) = 0.35329d0
+ S20RTS_V%spknt(11) = 0.44384d0
+ S20RTS_V%spknt(12) = 0.52615d0
+ S20RTS_V%spknt(13) = 0.60097d0
+ S20RTS_V%spknt(14) = 0.66899d0
+ S20RTS_V%spknt(15) = 0.73081d0
+ S20RTS_V%spknt(16) = 0.78701d0
+ S20RTS_V%spknt(17) = 0.83810d0
+ S20RTS_V%spknt(18) = 0.88454d0
+ S20RTS_V%spknt(19) = 0.92675d0
+ S20RTS_V%spknt(20) = 0.96512d0
+ S20RTS_V%spknt(21) = 1.00000d0
+
+ do i=1,NK_20+1
+ do j=1,NK_20+1
+ if(i == j) then
+ S20RTS_V%qq0(j,i)=1.0d0
+ else
+ S20RTS_V%qq0(j,i)=0.0d0
+ endif
+ enddo
+ enddo
+ do i=1,NK_20+1
+ call s20rts_rspln(1,NK_20+1,S20RTS_V%spknt(1),S20RTS_V%qq0(1,i),S20RTS_V%qq(1,1,i),qqwk(1,1))
+ enddo
+
+ end subroutine s20rts_splhsetup
+
+!----------------------------------
+
+! changed the obsolecent f77 features in the two routines below
+! now still awful Fortran, but at least conforms to f90 standard
+
+ double precision function s20rts_rsple(I1,I2,X,Y,Q,S)
+
+ implicit none
+
+! rsple returns the value of the function y(x) evaluated at point S
+! using the cubic spline coefficients computed by rspln and saved in Q.
+! If S is outside the interval (x(i1),x(i2)) rsple extrapolates
+! using the first or last interpolation polynomial. The arrays must
+! be dimensioned at least - x(i2), y(i2), and q(3,i2).
+
+ integer i1,i2
+ double precision X(*),Y(*),Q(3,*),s
+
+ integer i,ii
+ double precision h
+
+ i = 1
+ II=I2-1
+
+! GUARANTEE I WITHIN BOUNDS.
+ I=MAX0(I,I1)
+ I=MIN0(I,II)
+
+! SEE IF X IS INCREASING OR DECREASING.
+ IF(X(I2)-X(I1) < 0) goto 1
+ IF(X(I2)-X(I1) >= 0) goto 2
+
+! X IS DECREASING. CHANGE I AS NECESSARY.
+ 1 IF(S-X(I) <= 0) goto 3
+ IF(S-X(I) > 0) goto 4
+
+ 4 I=I-1
+
+ IF(I-I1 < 0) goto 11
+ IF(I-I1 == 0) goto 6
+ IF(I-I1 > 0) goto 1
+
+ 3 IF(S-X(I+1) < 0) goto 5
+ IF(S-X(I+1) >= 0) goto 6
+
+ 5 I=I+1
+
+ IF(I-II < 0) goto 3
+ IF(I-II == 0) goto 6
+ IF(I-II > 0) goto 7
+
+! X IS INCREASING. CHANGE I AS NECESSARY.
+ 2 IF(S-X(I+1) <= 0) goto 8
+ IF(S-X(I+1) > 0) goto 9
+
+ 9 I=I+1
+
+ IF(I-II < 0) goto 2
+ IF(I-II == 0) goto 6
+ IF(I-II > 0) goto 7
+
+ 8 IF(S-X(I) < 0) goto 10
+ IF(S-X(I) >= 0) goto 6
+
+ 10 I=I-1
+ IF(I-I1 < 0) goto 11
+ IF(I-I1 == 0) goto 6
+ IF(I-I1 > 0) goto 8
+
+ 7 I=II
+ GOTO 6
+ 11 I=I1
+
+! CALCULATE RSPLE USING SPLINE COEFFICIENTS IN Y AND Q.
+ 6 H=S-X(I)
+ S20RTS_RSPLE=Y(I)+H*(Q(1,I)+H*(Q(2,I)+H*Q(3,I)))
+
+ end function s20rts_rsple
+
+!----------------------------------
+
+ subroutine s20rts_rspln(I1,I2,X,Y,Q,F)
+
+ implicit none
+
+! Subroutine rspln computes cubic spline interpolation coefficients
+! for y(x) between grid points i1 and i2 saving them in q.The
+! interpolation is continuous with continuous first and second
+! derivatives. It agrees exactly with y at grid points and with the
+! three point first derivatives at both end points (i1 and i2).
+! X must be monotonic but if two successive values of x are equal
+! a discontinuity is assumed and separate interpolation is done on
+! each strictly monotonic segment. The arrays must be dimensioned at
+! least - x(i2), y(i2), q(3,i2), and f(3,i2).
+! F is working storage for rspln.
+
+ integer i1,i2
+ double precision X(*),Y(*),Q(3,*),F(3,*)
+
+ integer i,j,k,j1,j2
+ double precision y0,a0,b0,b1,h,h2,ha,h2a,h3a,h2b
+ double precision YY(3),small
+ equivalence (YY(1),Y0)
+ data SMALL/1.0d-08/,YY/0.0d0,0.0d0,0.0d0/
+
+ J1=I1+1
+ Y0=0.0d0
+
+! BAIL OUT IF THERE ARE LESS THAN TWO POINTS TOTAL
+ IF(I2-I1 < 0) return
+ IF(I2-I1 == 0) goto 17
+ IF(I2-I1 > 0) goto 8
+
+ 8 A0=X(J1-1)
+! SEARCH FOR DISCONTINUITIES.
+ DO 3 I=J1,I2
+ B0=A0
+ A0=X(I)
+ IF(DABS((A0-B0)/DMAX1(A0,B0)) < SMALL) GOTO 4
+ 3 CONTINUE
+ 17 J1=J1-1
+ J2=I2-2
+ GOTO 5
+ 4 J1=J1-1
+ J2=I-3
+! SEE IF THERE ARE ENOUGH POINTS TO INTERPOLATE (AT LEAST THREE).
+ 5 IF(J2+1-J1 < 0) goto 9
+ IF(J2+1-J1 == 0) goto 10
+ IF(J2+1-J1 > 0) goto 11
+
+! ONLY TWO POINTS. USE LINEAR INTERPOLATION.
+ 10 J2=J2+2
+ Y0=(Y(J2)-Y(J1))/(X(J2)-X(J1))
+ DO J=1,3
+ Q(J,J1)=YY(J)
+ Q(J,J2)=YY(J)
+ enddo
+ GOTO 12
+
+! MORE THAN TWO POINTS. DO SPLINE INTERPOLATION.
+ 11 A0=0.
+ H=X(J1+1)-X(J1)
+ H2=X(J1+2)-X(J1)
+ Y0=H*H2*(H2-H)
+ H=H*H
+ H2=H2*H2
+! CALCULATE DERIVITIVE AT NEAR END.
+ B0=(Y(J1)*(H-H2)+Y(J1+1)*H2-Y(J1+2)*H)/Y0
+ B1=B0
+
+! EXPLICITLY REDUCE BANDED MATRIX TO AN UPPER BANDED MATRIX.
+ DO I=J1,J2
+ H=X(I+1)-X(I)
+ Y0=Y(I+1)-Y(I)
+ H2=H*H
+ HA=H-A0
+ H2A=H-2.0d0*A0
+ H3A=2.0d0*H-3.0d0*A0
+ H2B=H2*B0
+ Q(1,I)=H2/HA
+ Q(2,I)=-HA/(H2A*H2)
+ Q(3,I)=-H*H2A/H3A
+ F(1,I)=(Y0-H*B0)/(H*HA)
+ F(2,I)=(H2B-Y0*(2.0d0*H-A0))/(H*H2*H2A)
+ F(3,I)=-(H2B-3.0d0*Y0*HA)/(H*H3A)
+ A0=Q(3,I)
+ B0=F(3,I)
+ enddo
+
+! TAKE CARE OF LAST TWO ROWS.
+ I=J2+1
+ H=X(I+1)-X(I)
+ Y0=Y(I+1)-Y(I)
+ H2=H*H
+ HA=H-A0
+ H2A=H*HA
+ H2B=H2*B0-Y0*(2.0d0*H-A0)
+ Q(1,I)=H2/HA
+ F(1,I)=(Y0-H*B0)/H2A
+ HA=X(J2)-X(I+1)
+ Y0=-H*HA*(HA+H)
+ HA=HA*HA
+
+! CALCULATE DERIVATIVE AT FAR END.
+ Y0=(Y(I+1)*(H2-HA)+Y(I)*HA-Y(J2)*H2)/Y0
+ Q(3,I)=(Y0*H2A+H2B)/(H*H2*(H-2.0d0*A0))
+ Q(2,I)=F(1,I)-Q(1,I)*Q(3,I)
+
+! SOLVE UPPER BANDED MATRIX BY REVERSE ITERATION.
+ DO J=J1,J2
+ K=I-1
+ Q(1,I)=F(3,K)-Q(3,K)*Q(2,I)
+ Q(3,K)=F(2,K)-Q(2,K)*Q(1,I)
+ Q(2,K)=F(1,K)-Q(1,K)*Q(3,K)
+ I=K
+ enddo
+ Q(1,I)=B1
+! FILL IN THE LAST POINT WITH A LINEAR EXTRAPOLATION.
+ 9 J2=J2+2
+ DO J=1,3
+ Q(J,J2)=YY(J)
+ enddo
+
+! SEE IF THIS DISCONTINUITY IS THE LAST.
+ 12 IF(J2-I2 < 0) then
+ goto 6
+ else
+ return
+ endif
+
+! NO. GO BACK FOR MORE.
+ 6 J1=J2+2
+ IF(J1-I2 <= 0) goto 8
+ IF(J1-I2 > 0) goto 7
+
+! THERE IS ONLY ONE POINT LEFT AFTER THE LATEST DISCONTINUITY.
+ 7 DO J=1,3
+ Q(J,I2)=YY(J)
+ enddo
+
+ end subroutine s20rts_rspln
+
diff -r 000000000000 -r 9b613f27327a test_compilation
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/test_compilation Sat Jul 02 17:09:21 2011 -0700
@@ -0,0 +1,20 @@
+#!/bin/bash
+
+if [ $# -ne 1 -o ! -d "$1" ]; then
+ echo "Usage: test_compilation /path/to/SPECFEM3D_GLOBE"
+ echo ""
+ echo "The Specfem 3D Globe directory must be already configured"
+ exit 1
+fi
+
+if [ ! -f "$1/Makefile" ]; then
+ echo "ERROR: The Specfem 3D Globe directory must be already configured."
+ echo "The directory"
+ echo " $1"
+ echo "is missing the file 'Makefile'. Please run 'configure' in that directory."
+ exit 1
+fi
+
+cp -a DATA/ src/ "$1"
+cd "$1"
+make
More information about the CIG-COMMITS
mailing list