[cig-commits] r17067 - seismo/3D/SPECFEM3D_GLOBE/trunk
pieyre at geodynamics.org
pieyre at geodynamics.org
Thu Aug 5 10:17:20 PDT 2010
Author: pieyre
Date: 2010-08-05 10:17:20 -0700 (Thu, 05 Aug 2010)
New Revision: 17067
Modified:
seismo/3D/SPECFEM3D_GLOBE/trunk/auto_ner.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/combine_AVS_DX.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/write_AVS_DX_global_faces_data.f90
seismo/3D/SPECFEM3D_GLOBE/trunk/write_AVS_DX_surface_data.f90
Log:
modifications, bug fixes related to the creation of OpenDX and AVS files and fixes of a few wrong references found in the code
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/auto_ner.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/auto_ner.f90 2010-08-04 14:42:45 UTC (rev 17066)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/auto_ner.f90 2010-08-05 17:17:20 UTC (rev 17067)
@@ -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/trunk/combine_AVS_DX.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/combine_AVS_DX.f90 2010-08-04 14:42:45 UTC (rev 17066)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/combine_AVS_DX.f90 2010-08-05 17:17:20 UTC (rev 17067)
@@ -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/trunk/create_regions_mesh.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90 2010-08-04 14:42:45 UTC (rev 17066)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/create_regions_mesh.f90 2010-08-05 17:17:20 UTC (rev 17067)
@@ -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/trunk/specfem3D.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90 2010-08-04 14:42:45 UTC (rev 17066)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/specfem3D.f90 2010-08-05 17:17:20 UTC (rev 17067)
@@ -1382,7 +1382,7 @@
if(MOVIE_COARSE) then !only output corners
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
print *,' type: surface - coarse'
print *,' points: ',nmovie_points
Modified: seismo/3D/SPECFEM3D_GLOBE/trunk/write_AVS_DX_global_faces_data.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/write_AVS_DX_global_faces_data.f90 2010-08-04 14:42:45 UTC (rev 17066)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/write_AVS_DX_global_faces_data.f90 2010-08-05 17:17:20 UTC (rev 17067)
@@ -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/trunk/write_AVS_DX_surface_data.f90
===================================================================
--- seismo/3D/SPECFEM3D_GLOBE/trunk/write_AVS_DX_surface_data.f90 2010-08-04 14:42:45 UTC (rev 17066)
+++ seismo/3D/SPECFEM3D_GLOBE/trunk/write_AVS_DX_surface_data.f90 2010-08-05 17:17:20 UTC (rev 17067)
@@ -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