[cig-commits] r17086 - in seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY: . UTILS/seis_process
yangl at geodynamics.org
yangl at geodynamics.org
Thu Aug 12 15:31:38 PDT 2010
Author: yangl
Date: 2010-08-12 15:31:38 -0700 (Thu, 12 Aug 2010)
New Revision: 17086
Added:
seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/UTILS/seis_process/asc2sac/
Modified:
seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/auto_ner.f90
seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/combine_AVS_DX.f90
seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/create_regions_mesh.f90
seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/specfem3D.f90
seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/todo_list_please_dont_remove.txt
seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/write_AVS_DX_global_faces_data.f90
seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/write_AVS_DX_surface_data.f90
Log:
NOISE: merging changes from trunk (rev:17049,17055,17056,17067) (Done)
Copied: seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/UTILS/seis_process/asc2sac (from rev 17049, seismo/3D/SPECFEM3D_GLOBE/trunk/UTILS/seis_process/asc2sac)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/auto_ner.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/auto_ner.f90 2010-08-12 22:24:52 UTC (rev 17085)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/auto_ner.f90 2010-08-12 22:31:38 UTC (rev 17086)
@@ -120,7 +120,7 @@
MIN_ATTENUATION_PERIOD = TMP
if(N_SLS < 2 .OR. N_SLS > 5) then
- call exit_MPI_without_rank('N_SLS must be greater than 1 or less than 6')
+ stop 'N_SLS must be greater than 1 or less than 6'
endif
! Compute Max Attenuation Period
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/combine_AVS_DX.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/combine_AVS_DX.f90 2010-08-12 22:24:52 UTC (rev 17085)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/combine_AVS_DX.f90 2010-08-12 22:31:38 UTC (rev 17086)
@@ -549,7 +549,9 @@
if(ivalue == 1) then
open(unit=10,file=prname(1:len_trim(prname))//'AVS_DXelementsfaces.txt',status='old',action='read')
open(unit=12,file=prname(1:len_trim(prname))//'AVS_DXpointsfaces.txt',status='old',action='read')
- else if(ivalue == 2) then
+ if(icolor == 5 .or. icolor == 6) &
+ open(unit=13,file=prname(1:len_trim(prname))//'AVS_DXelementsfaces_dvp_dvs.txt',status='old',action='read')
+ else if(ivalue == 2) then
open(unit=10,file=prname(1:len_trim(prname))//'AVS_DXelementschunks.txt',status='old',action='read')
if(icolor == 5 .or. icolor == 6) &
open(unit=13,file=prname(1:len_trim(prname))//'AVS_DXelementschunks_dvp_dvs.txt',status='old',action='read')
@@ -557,8 +559,8 @@
else if(ivalue == 3) then
open(unit=10,file=prname(1:len_trim(prname))//'AVS_DXelementssurface.txt',status='old',action='read')
open(unit=12,file=prname(1:len_trim(prname))//'AVS_DXpointssurface.txt',status='old',action='read')
- !if(icolor == 5 .or. icolor == 6) &
- ! open(unit=13,file=prname(1:len_trim(prname))//'AVS_DXelementschunks_dvp_dvs.txt',status='old',action='read')
+ if(icolor == 5 .or. icolor == 6) &
+ open(unit=13,file=prname(1:len_trim(prname))//'AVS_DXelementssurface_dvp_dvs.txt',status='old',action='read')
endif
@@ -809,8 +811,10 @@
endif
- if(ISOTROPIC_3D_MANTLE) then
+ if(icolor == 5 .or. icolor == 6) then
+ if(ISOTROPIC_3D_MANTLE) then
+
! compute absolute maximum for dvp
rnorm_factor = maxval(dabs(dvp(:)))
@@ -843,6 +847,7 @@
endif
enddo
+ endif
endif
! ************* generate element data values ******************
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/create_regions_mesh.f90 2010-08-12 22:24:52 UTC (rev 17085)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/create_regions_mesh.f90 2010-08-12 22:31:38 UTC (rev 17086)
@@ -660,7 +660,11 @@
call write_AVS_DX_global_data(myrank,prname,nspec,ibool,idoubling,xstore,ystore,zstore,locval,ifseg,npointot)
call write_AVS_DX_global_faces_data(myrank,prname,nspec,iMPIcut_xi,iMPIcut_eta,ibool, &
- idoubling,xstore,ystore,zstore,locval,ifseg,npointot)
+ idoubling,xstore,ystore,zstore,locval,ifseg,npointot, &
+ rhostore,kappavstore,muvstore,nspl,rspl,espl,espl2, &
+ ELLIPTICITY,ISOTROPIC_3D_MANTLE, &
+ RICB,RCMB,RTOPDDOUBLEPRIME,R600,R670,R220,R771,R400,R120,R80,RMOHO, &
+ RMIDDLE_CRUST,ROCEAN,iregion_code)
call write_AVS_DX_global_chunks_data(myrank,prname,nspec,iboun,ibool, &
idoubling,xstore,ystore,zstore,locval,ifseg,npointot, &
@@ -670,7 +674,11 @@
RMIDDLE_CRUST,ROCEAN,iregion_code)
call write_AVS_DX_surface_data(myrank,prname,nspec,iboun,ibool, &
- idoubling,xstore,ystore,zstore,locval,ifseg,npointot)
+ idoubling,xstore,ystore,zstore,locval,ifseg,npointot, &
+ rhostore,kappavstore,muvstore,nspl,rspl,espl,espl2, &
+ ELLIPTICITY,ISOTROPIC_3D_MANTLE, &
+ RICB,RCMB,RTOPDDOUBLEPRIME,R600,R670,R220,R771,R400,R120,R80,RMOHO, &
+ RMIDDLE_CRUST,ROCEAN,iregion_code)
!> Hejun
! Output material information for all GLL points
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/specfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/specfem3D.f90 2010-08-12 22:24:52 UTC (rev 17085)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/specfem3D.f90 2010-08-12 22:31:38 UTC (rev 17086)
@@ -1389,7 +1389,7 @@
if(MOVIE_COARSE .and. NOISE_TOMOGRAPHY ==0) then ! only output corners !for noise tomography, must NOT be coarse
nmovie_points = 2 * 2 * NSPEC2D_TOP(IREGION_CRUST_MANTLE)
if(NGLLX /= NGLLY) &
- call MPI_exit(myrank,'MOVIE_COARSE together with MOVIE_SURFACE requires NGLLX=NGLLY')
+ call exit_MPI(myrank,'MOVIE_COARSE together with MOVIE_SURFACE requires NGLLX=NGLLY')
NIT = NGLLX - 1
else
nmovie_points = NGLLX * NGLLY * NSPEC2D_TOP(IREGION_CRUST_MANTLE)
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/todo_list_please_dont_remove.txt
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/todo_list_please_dont_remove.txt 2010-08-12 22:24:52 UTC (rev 17085)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/todo_list_please_dont_remove.txt 2010-08-12 22:31:38 UTC (rev 17086)
@@ -4,6 +4,46 @@
Things that could be done in a future version:
+Subject: Re: SEM attenuation
+Date: Fri, 23 Jul 2010 12:16:40 +0200
+From: Tarje Nissen-Meyer <tarjen at ethz.ch>
+Organization: ETH Zurich
+To: Dimitri Komatitsch <dimitri.komatitsch at univ-pau.fr>
+CC: <yingz at vt.edu>, Jeroen Tromp <jtromp at princeton.edu>, Min Chen <mchen at gps.caltech.edu>, Vala Hjorleifsdottir <vala at ldeo.columbia.edu>, "Brian Savage" <savage at uri.edu>, Shiann-Jong Lee <sjlee at earth.sinica.edu.tw>, "Roland Martin" <roland.martin at univ-pau.fr>, Bernhard Schuberth <mail at bernhard-schuberth.de>, Carl Tape <carltape at fas.harvard.edu>, "Anne Sieminski" <anne.sieminski at obs.ujf-grenoble.fr>, Paul Friberg <p.friberg at isti.com>, Kasper van Wijk <kasper at cgiss.boisestate.edu>, "Dylan Mikesell" <dmikesell at cgiss.boisestate.edu>, Federica Magnoni <federica.magnoni at ingv.it>, Ebru Bozdag <bozdag at princeton.edu>, Hejun Zhu <hejunzhu at princeton.edu>, Pieyre Le Loher <pieyre.le_loher at inria.fr>, Christina Morency <cmorency at Princeton.EDU>, Emanuele Casarotti <emanuele.casarotti at gmail.com>, Piero Basini <basini at bo.ingv.it>, "Emmanuel Chaljub" <Emmanuel.Chaljub at obs.ujf-grenoble.fr>, Qinya Liu <liuqy at physics.utoronto.ca>, Yang Luo <yangl at Princeton.EDU>, Daniel Peter <dpeter at Princeton.EDU>
+
+Hi all,
+
+Just to add that symplectic schemes are based on conservative systems
+(they approximate the Hamiltonian) and we don't know yet how they behave
+for dissipative media. But still of course something to be added and
+looked at. Should also give credit to Jean-Paul Ampuero who really first
+suggested trying these schemes for elastodynamics with the SEM.
+
+Best regards,
+Tarje
+
+On 23/07/10 02:21, Dimitri Komatitsch wrote:
+>
+> Hi Ying,
+>
+> What I know for sure if that back in 1999 when I developed the time
+> scheme for attenuation I used a trick that made implementation much
+> easier but that also makes the 4th order Runge Kutta (RK4) time scheme
+> become second order only, i.e. RK2 instead of RK4. This means that for
+> very long simulations (for instance multi-orbit surface waves)
+> attenuation can become inaccurate because of a lack of accuracy in the
+> time scheme (in which case reducing Delta_t purposely solves the
+> problem, but of course makes the simulation more expensive).
+>
+> A way of solving this problem could be to switch to better time
+> schemes such as the symplectic time scheme that Tarje introduced in a
+> GJI paper a few years ago. Tarje and I should probably implement that
+> in the code one day... (therefore I cc him; we talked about this at
+> the EGU meeting in May)
+>
+> Thank you,
+> Dimitri.
+
- symplectic time scheme (will be done by Tarje Nissen-Meyer and/or Jean-Paul Ampuero) (would be useful in the 2D version of the code as well):
Hi Jeroen, Perfect. I think talking to Jean-Paul Ampuero would be useful
as well because in Utrecht last year he had told us that
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/write_AVS_DX_global_faces_data.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/write_AVS_DX_global_faces_data.f90 2010-08-12 22:24:52 UTC (rev 17085)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/write_AVS_DX_global_faces_data.f90 2010-08-12 22:31:38 UTC (rev 17086)
@@ -30,7 +30,10 @@
subroutine write_AVS_DX_global_faces_data(myrank,prname,nspec,iMPIcut_xi,iMPIcut_eta, &
ibool,idoubling,xstore,ystore,zstore,num_ibool_AVS_DX,mask_ibool, &
- npointot)
+ npointot,rhostore,kappavstore,muvstore,nspl,rspl,espl,espl2, &
+ ELLIPTICITY,ISOTROPIC_3D_MANTLE, &
+ RICB,RCMB,RTOPDDOUBLEPRIME,R600,R670,R220,R771,R400,R120,R80,RMOHO, &
+ RMIDDLE_CRUST,ROCEAN,iregion_code)
implicit none
@@ -41,13 +44,22 @@
integer idoubling(nspec)
+ logical ELLIPTICITY,ISOTROPIC_3D_MANTLE
+
logical iMPIcut_xi(2,nspec)
logical iMPIcut_eta(2,nspec)
+ double precision RICB,RCMB,RTOPDDOUBLEPRIME,R600,R670,R220,R771, &
+ R400,R120,R80,RMOHO,RMIDDLE_CRUST,ROCEAN
+
double precision xstore(NGLLX,NGLLY,NGLLZ,nspec)
double precision ystore(NGLLX,NGLLY,NGLLZ,nspec)
double precision zstore(NGLLX,NGLLY,NGLLZ,nspec)
+ real(kind=CUSTOM_REAL) kappavstore(NGLLX,NGLLY,NGLLZ,nspec)
+ real(kind=CUSTOM_REAL) muvstore(NGLLX,NGLLY,NGLLZ,nspec)
+ real(kind=CUSTOM_REAL) rhostore(NGLLX,NGLLY,NGLLZ,nspec)
+
! logical mask used to output global points only once
integer npointot
logical mask_ibool(npointot)
@@ -56,12 +68,24 @@
integer num_ibool_AVS_DX(npointot)
integer ispec
+ integer i,j,k,np
integer iglob1,iglob2,iglob3,iglob4,iglob5,iglob6,iglob7,iglob8
integer npoin,numpoin,nspecface,ispecface
+ double precision r,rho,vp,vs,Qkappa,Qmu
+ double precision vpv,vph,vsv,vsh,eta_aniso
+ double precision x,y,z,theta,phi_dummy,cost,p20,ell,factor
+ real(kind=CUSTOM_REAL) dvp,dvs
+
+! for ellipticity
+ integer nspl
+ double precision rspl(NR),espl(NR),espl2(NR)
+
! processor identification
character(len=150) prname
+ integer iregion_code
+
! writing points
open(unit=10,file=prname(1:len_trim(prname))//'AVS_DXpointsfaces.txt',status='unknown')
@@ -288,6 +312,8 @@
! writing elements
open(unit=10,file=prname(1:len_trim(prname))//'AVS_DXelementsfaces.txt',status='unknown')
+ if(ISOTROPIC_3D_MANTLE) &
+ open(unit=11,file=prname(1:len_trim(prname))//'AVS_DXelementsfaces_dvp_dvs.txt',status='unknown')
! number of elements in AVS or DX file
write(10,*) nspecface
@@ -306,12 +332,82 @@
iglob7=ibool(NGLLX,NGLLY,NGLLZ,ispec)
iglob8=ibool(1,NGLLY,NGLLZ,ispec)
+! include lateral variations if needed
+
+ if(ISOTROPIC_3D_MANTLE) then
+ ! pick a point within the element and get its radius
+ r=dsqrt(xstore(2,2,2,ispec)**2+ystore(2,2,2,ispec)**2+zstore(2,2,2,ispec)**2)
+
+ if(r > RCMB/R_EARTH .and. r < R_UNIT_SPHERE) then
+ ! average over the element
+ dvp = 0.0
+ dvs = 0.0
+ np =0
+ do k=2,NGLLZ-1
+ do j=2,NGLLY-1
+ do i=2,NGLLX-1
+ np=np+1
+ x=xstore(i,j,k,ispec)
+ y=ystore(i,j,k,ispec)
+ z=zstore(i,j,k,ispec)
+ r=dsqrt(x*x+y*y+z*z)
+ ! take out ellipticity
+ if(ELLIPTICITY) then
+ call xyz_2_rthetaphi_dble(x,y,z,r,theta,phi_dummy)
+ cost=dcos(theta)
+ p20=0.5d0*(3.0d0*cost*cost-1.0d0)
+ call spline_evaluation(rspl,espl,espl2,nspl,r,ell)
+ factor=ONE-(TWO/3.0d0)*ell*p20
+ r=r/factor
+ endif
+
+
+ ! gets reference model values: rho,vpv,vph,vsv,vsh and eta_aniso
+ call meshfem3D_models_get1D_val(myrank,iregion_code,idoubling(ispec), &
+ r,rho,vpv,vph,vsv,vsh,eta_aniso, &
+ Qkappa,Qmu,RICB,RCMB, &
+ RTOPDDOUBLEPRIME,R80,R120,R220,R400,R600,R670,R771, &
+ RMOHO,RMIDDLE_CRUST,ROCEAN)
+
+ ! calculates isotropic values
+ vp = sqrt(((8.d0+4.d0*eta_aniso)*vph*vph + 3.d0*vpv*vpv &
+ + (8.d0 - 8.d0*eta_aniso)*vsv*vsv)/15.d0)
+ vs = sqrt(((1.d0-2.d0*eta_aniso)*vph*vph + vpv*vpv &
+ + 5.d0*vsh*vsh + (6.d0+4.d0*eta_aniso)*vsv*vsv)/15.d0)
+
+ if( abs(rhostore(i,j,k,ispec))< 1.e-20 ) then
+ print*,' attention: rhostore close to zero',rhostore(i,j,k,ispec),r,i,j,k,ispec
+ dvp = 0.0
+ dvs = 0.0
+ else if( abs(sngl(vp))< 1.e-20 ) then
+ print*,' attention: vp close to zero',sngl(vp),r,i,j,k,ispec
+ dvp = 0.0
+ else if( abs(sngl(vs))< 1.e-20 ) then
+ print*,' attention: vs close to zero',sngl(vs),r,i,j,k,ispec
+ dvs = 0.0
+ else
+ dvp = dvp + (sqrt((kappavstore(i,j,k,ispec)+4.*muvstore(i,j,k,ispec)/3.)/rhostore(i,j,k,ispec)) - sngl(vp))/sngl(vp)
+ dvs = dvs + (sqrt(muvstore(i,j,k,ispec)/rhostore(i,j,k,ispec)) - sngl(vs))/sngl(vs)
+ endif
+
+ enddo
+ enddo
+ enddo
+ dvp = dvp / np
+ dvs = dvs / np
+ else
+ dvp = 0.0
+ dvs = 0.0
+ endif
+ endif
+
! face xi = xi_min
if(iMPIcut_xi(1,ispec)) then
ispecface = ispecface + 1
write(10,*) ispecface,idoubling(ispec),num_ibool_AVS_DX(iglob1), &
num_ibool_AVS_DX(iglob4),num_ibool_AVS_DX(iglob8), &
num_ibool_AVS_DX(iglob5)
+ if(ISOTROPIC_3D_MANTLE) write(11,*) ispecface,dvp,dvs
endif
! face xi = xi_max
@@ -320,6 +416,7 @@
write(10,*) ispecface,idoubling(ispec),num_ibool_AVS_DX(iglob2), &
num_ibool_AVS_DX(iglob3),num_ibool_AVS_DX(iglob7), &
num_ibool_AVS_DX(iglob6)
+ if(ISOTROPIC_3D_MANTLE) write(11,*) ispecface,dvp,dvs
endif
! face eta = eta_min
@@ -328,6 +425,7 @@
write(10,*) ispecface,idoubling(ispec),num_ibool_AVS_DX(iglob1), &
num_ibool_AVS_DX(iglob2),num_ibool_AVS_DX(iglob6), &
num_ibool_AVS_DX(iglob5)
+ if(ISOTROPIC_3D_MANTLE) write(11,*) ispecface,dvp,dvs
endif
! face eta = eta_max
@@ -336,6 +434,7 @@
write(10,*) ispecface,idoubling(ispec),num_ibool_AVS_DX(iglob4), &
num_ibool_AVS_DX(iglob3),num_ibool_AVS_DX(iglob7), &
num_ibool_AVS_DX(iglob8)
+ if(ISOTROPIC_3D_MANTLE) write(11,*) ispecface,dvp,dvs
endif
endif
@@ -346,6 +445,7 @@
call exit_MPI(myrank,'incorrect number of surface elements in AVS or DX file creation')
close(10)
+ if(ISOTROPIC_3D_MANTLE) close(11)
end subroutine write_AVS_DX_global_faces_data
Modified: seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/write_AVS_DX_surface_data.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/write_AVS_DX_surface_data.f90 2010-08-12 22:24:52 UTC (rev 17085)
+++ seismo/3D/SPECFEM3D_GLOBE/branches/NOISE_TOMOGRAPHY/write_AVS_DX_surface_data.f90 2010-08-12 22:31:38 UTC (rev 17086)
@@ -28,7 +28,11 @@
! create AVS or DX 2D data for the surface of the model
! to be recombined in postprocessing
subroutine write_AVS_DX_surface_data(myrank,prname,nspec,iboun, &
- ibool,idoubling,xstore,ystore,zstore,num_ibool_AVS_DX,mask_ibool,npointot)
+ ibool,idoubling,xstore,ystore,zstore,num_ibool_AVS_DX,mask_ibool,npointot,&
+ rhostore,kappavstore,muvstore,nspl,rspl,espl,espl2, &
+ ELLIPTICITY,ISOTROPIC_3D_MANTLE, &
+ RICB,RCMB,RTOPDDOUBLEPRIME,R600,R670,R220,R771,R400,R120,R80,RMOHO, &
+ RMIDDLE_CRUST,ROCEAN,iregion_code)
implicit none
@@ -40,11 +44,24 @@
integer idoubling(nspec)
logical iboun(6,nspec)
+ logical ELLIPTICITY,ISOTROPIC_3D_MANTLE
+ double precision RICB,RCMB,RTOPDDOUBLEPRIME,R600,R670,R220,R771, &
+ R400,R120,R80,RMOHO,RMIDDLE_CRUST,ROCEAN
+
+ double precision r,rho,vp,vs,Qkappa,Qmu
+ double precision vpv,vph,vsv,vsh,eta_aniso
+ double precision x,y,z,theta,phi_dummy,cost,p20,ell,factor
+ real(kind=CUSTOM_REAL) dvp,dvs
+
double precision xstore(NGLLX,NGLLY,NGLLZ,nspec)
double precision ystore(NGLLX,NGLLY,NGLLZ,nspec)
double precision zstore(NGLLX,NGLLY,NGLLZ,nspec)
+ real(kind=CUSTOM_REAL) kappavstore(NGLLX,NGLLY,NGLLZ,nspec)
+ real(kind=CUSTOM_REAL) muvstore(NGLLX,NGLLY,NGLLZ,nspec)
+ real(kind=CUSTOM_REAL) rhostore(NGLLX,NGLLY,NGLLZ,nspec)
+
! logical mask used to output global points only once
integer npointot
logical mask_ibool(npointot)
@@ -53,12 +70,19 @@
integer num_ibool_AVS_DX(npointot)
integer ispec
+ integer i,j,k,np
integer, dimension(8) :: iglobval
integer npoin,numpoin,nspecface,ispecface
+! for ellipticity
+ integer nspl
+ double precision rspl(NR),espl(NR),espl2(NR)
+
! processor identification
character(len=150) prname
+ integer iregion_code
+
! writing points
open(unit=10,file=prname(1:len_trim(prname))//'AVS_DXpointssurface.txt',status='unknown')
@@ -158,6 +182,8 @@
! writing elements
open(unit=10,file=prname(1:len_trim(prname))//'AVS_DXelementssurface.txt',status='unknown')
+ if(ISOTROPIC_3D_MANTLE) &
+ open(unit=11,file=prname(1:len_trim(prname))//'AVS_DXelementssurface_dvp_dvs.txt',status='unknown')
! number of elements in AVS or DX file
write(10,*) nspecface
@@ -165,20 +191,89 @@
ispecface = 0
do ispec=1,nspec
! only if at the surface
- if(iboun(6,ispec)) then
+ if(iboun(6,ispec)) then
- iglobval(5)=ibool(1,1,NGLLZ,ispec)
- iglobval(6)=ibool(NGLLX,1,NGLLZ,ispec)
- iglobval(7)=ibool(NGLLX,NGLLY,NGLLZ,ispec)
- iglobval(8)=ibool(1,NGLLY,NGLLZ,ispec)
+ iglobval(5)=ibool(1,1,NGLLZ,ispec)
+ iglobval(6)=ibool(NGLLX,1,NGLLZ,ispec)
+ iglobval(7)=ibool(NGLLX,NGLLY,NGLLZ,ispec)
+ iglobval(8)=ibool(1,NGLLY,NGLLZ,ispec)
-! top face
- ispecface = ispecface + 1
- write(10,*) ispecface,idoubling(ispec),num_ibool_AVS_DX(iglobval(5)), &
- num_ibool_AVS_DX(iglobval(6)),num_ibool_AVS_DX(iglobval(7)), &
- num_ibool_AVS_DX(iglobval(8))
+ if(ISOTROPIC_3D_MANTLE) then
+ ! pick a point within the element and get its radius
+ r=dsqrt(xstore(2,2,2,ispec)**2+ystore(2,2,2,ispec)**2+zstore(2,2,2,ispec)**2)
- endif
+ if(r > RCMB/R_EARTH .and. r < R_UNIT_SPHERE) then
+ ! average over the element
+ dvp = 0.0
+ dvs = 0.0
+ np =0
+ do k=2,NGLLZ-1
+ do j=2,NGLLY-1
+ do i=2,NGLLX-1
+ np=np+1
+ x=xstore(i,j,k,ispec)
+ y=ystore(i,j,k,ispec)
+ z=zstore(i,j,k,ispec)
+ r=dsqrt(x*x+y*y+z*z)
+ ! take out ellipticity
+ if(ELLIPTICITY) then
+ call xyz_2_rthetaphi_dble(x,y,z,r,theta,phi_dummy)
+ cost=dcos(theta)
+ p20=0.5d0*(3.0d0*cost*cost-1.0d0)
+ call spline_evaluation(rspl,espl,espl2,nspl,r,ell)
+ factor=ONE-(TWO/3.0d0)*ell*p20
+ r=r/factor
+ endif
+
+
+ ! gets reference model values: rho,vpv,vph,vsv,vsh and eta_aniso
+ call meshfem3D_models_get1D_val(myrank,iregion_code,idoubling(ispec), &
+ r,rho,vpv,vph,vsv,vsh,eta_aniso, &
+ Qkappa,Qmu,RICB,RCMB, &
+ RTOPDDOUBLEPRIME,R80,R120,R220,R400,R600,R670,R771, &
+ RMOHO,RMIDDLE_CRUST,ROCEAN)
+
+ ! calculates isotropic values
+ vp = sqrt(((8.d0+4.d0*eta_aniso)*vph*vph + 3.d0*vpv*vpv &
+ + (8.d0 - 8.d0*eta_aniso)*vsv*vsv)/15.d0)
+ vs = sqrt(((1.d0-2.d0*eta_aniso)*vph*vph + vpv*vpv &
+ + 5.d0*vsh*vsh + (6.d0+4.d0*eta_aniso)*vsv*vsv)/15.d0)
+
+ if( abs(rhostore(i,j,k,ispec))< 1.e-20 ) then
+ print*,' attention: rhostore close to zero',rhostore(i,j,k,ispec),r,i,j,k,ispec
+ dvp = 0.0
+ dvs = 0.0
+ else if( abs(sngl(vp))< 1.e-20 ) then
+ print*,' attention: vp close to zero',sngl(vp),r,i,j,k,ispec
+ dvp = 0.0
+ else if( abs(sngl(vs))< 1.e-20 ) then
+ print*,' attention: vs close to zero',sngl(vs),r,i,j,k,ispec
+ dvs = 0.0
+ else
+ dvp = dvp + (sqrt((kappavstore(i,j,k,ispec)+4.*muvstore(i,j,k,ispec)/3.) &
+ /rhostore(i,j,k,ispec)) - sngl(vp))/sngl(vp)
+ dvs = dvs + (sqrt(muvstore(i,j,k,ispec)/rhostore(i,j,k,ispec)) - sngl(vs))/sngl(vs)
+ endif
+
+ enddo
+ enddo
+ enddo
+ dvp = dvp / np
+ dvs = dvs / np
+ else
+ dvp = 0.0
+ dvs = 0.0
+ endif
+ endif
+
+ ! top face
+ ispecface = ispecface + 1
+ write(10,*) ispecface,idoubling(ispec),num_ibool_AVS_DX(iglobval(5)), &
+ num_ibool_AVS_DX(iglobval(6)),num_ibool_AVS_DX(iglobval(7)), &
+ num_ibool_AVS_DX(iglobval(8))
+ if(ISOTROPIC_3D_MANTLE) write(11,*) ispecface,dvp,dvs
+
+ endif
enddo
! check that number of surface elements output is okay
@@ -186,6 +281,7 @@
call exit_MPI(myrank,'incorrect number of surface elements in AVS or DX file creation')
close(10)
+ if(ISOTROPIC_3D_MANTLE) close(11)
end subroutine write_AVS_DX_surface_data
More information about the CIG-COMMITS
mailing list