[cig-commits] r11913 - seismo/2D/SPECFEM2D/trunk
dkomati1 at geodynamics.org
dkomati1 at geodynamics.org
Sat May 3 17:41:04 PDT 2008
Author: dkomati1
Date: 2008-05-03 17:41:04 -0700 (Sat, 03 May 2008)
New Revision: 11913
Modified:
seismo/2D/SPECFEM2D/trunk/compute_curl_one_element.f90
seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90
seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90
seismo/2D/SPECFEM2D/trunk/paco_beyond_critical.f90
seismo/2D/SPECFEM2D/trunk/paco_convolve_fft.f90
seismo/2D/SPECFEM2D/trunk/specfem2D.F90
Log:
- fixed a small inconsistency in the definition of the amplitude of a pressure source in an acoustic element
- converted Paco's subroutines to Fortran95: double complex and CDEXP not allowed anymore
- removed several unused variables
Modified: seismo/2D/SPECFEM2D/trunk/compute_curl_one_element.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_curl_one_element.f90 2008-05-03 17:24:53 UTC (rev 11912)
+++ seismo/2D/SPECFEM2D/trunk/compute_curl_one_element.f90 2008-05-04 00:41:04 UTC (rev 11913)
@@ -66,7 +66,7 @@
real(kind=CUSTOM_REAL), dimension(NGLLZ,NGLLZ) :: hprime_zz
! local variables
- integer :: i,j,k,iglob
+ integer :: i,j,k
! jacobian
real(kind=CUSTOM_REAL) :: xixl,xizl,gammaxl,gammazl
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90 2008-05-03 17:24:53 UTC (rev 11912)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_acoustic.f90 2008-05-04 00:41:04 UTC (rev 11913)
@@ -41,16 +41,14 @@
!========================================================================
subroutine compute_forces_acoustic(npoin,nspec,nelemabs,numat, &
- iglob_source,ispec_selected_source,is_proc_source,source_type,it,NSTEP,anyabs, &
- assign_external_model,initialfield,ibool,kmato,numabs, &
+ anyabs,assign_external_model,ibool,kmato,numabs, &
elastic,codeabs,potential_dot_dot_acoustic,potential_dot_acoustic, &
potential_acoustic,density,elastcoef,xix,xiz,gammax,gammaz,jacobian, &
- vpext,rhoext,source_time_function,hprime_xx,hprimewgll_xx, &
+ vpext,rhoext,hprime_xx,hprimewgll_xx, &
hprime_zz,hprimewgll_zz,wxgll,wzgll, &
ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
jbegin_left,jend_left,jbegin_right,jend_right, &
- nspec_inner_outer, ispec_inner_outer_to_glob, num_phase_inner_outer &
- )
+ nspec_inner_outer, ispec_inner_outer_to_glob, num_phase_inner_outer)
! compute forces for the acoustic elements
@@ -58,9 +56,9 @@
include "constants.h"
- integer :: npoin,nspec,nelemabs,numat,iglob_source,ispec_selected_source,is_proc_source,source_type,it,NSTEP
+ integer :: npoin,nspec,nelemabs,numat
- logical :: anyabs,assign_external_model,initialfield
+ logical :: anyabs,assign_external_model
integer, dimension(NGLLX,NGLLZ,nspec) :: ibool
integer, dimension(nspec) :: kmato
@@ -75,7 +73,6 @@
double precision, dimension(4,numat) :: elastcoef
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLZ,nspec) :: xix,xiz,gammax,gammaz,jacobian
double precision, dimension(NGLLX,NGLLZ,nspec) :: vpext,rhoext
- real(kind=CUSTOM_REAL), dimension(NSTEP) :: source_time_function
! derivatives of Lagrange polynomials
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: hprime_xx,hprimewgll_xx
@@ -338,31 +335,6 @@
endif ! end of absorbing boundaries
-
-! --- add the source
- if(.not. initialfield) then
-
- if (is_proc_source == 1 ) then
-! collocated force
-! beware, for acoustic medium, source is a pressure source
-! the sign is negative because pressure p = - Chi_dot_dot therefore we need
-! to add minus the source to Chi_dot_dot to get plus the source in pressure
- if(source_type == 1) then
- if(.not. elastic(ispec_selected_source)) then
- potential_dot_dot_acoustic(iglob_source) = potential_dot_dot_acoustic(iglob_source) - source_time_function(it)
- endif
-
-! moment tensor
- else if(source_type == 2) then
- if(.not. elastic(ispec_selected_source)) then
- call exit_MPI('cannot have moment tensor source in acoustic element')
- endif
- endif
- endif
- else
- call exit_MPI('wrong source type')
- endif
-
endif ! end of computation that needs to be done only once, during the first call to compute_forces_acoustic
end subroutine compute_forces_acoustic
Modified: seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90 2008-05-03 17:24:53 UTC (rev 11912)
+++ seismo/2D/SPECFEM2D/trunk/compute_forces_elastic.f90 2008-05-04 00:41:04 UTC (rev 11913)
@@ -40,7 +40,7 @@
!
!========================================================================
- subroutine compute_forces_elastic(npoin,nspec,nelemabs,numat,iglob_source, &
+ subroutine compute_forces_elastic(npoin,nspec,nelemabs,numat, &
ispec_selected_source,is_proc_source,source_type,it,NSTEP,anyabs,assign_external_model, &
initialfield,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,angleforce,deltatcube, &
deltatfourth,twelvedeltat,fourdeltatsquare,ibool,kmato,numabs,elastic,codeabs, &
@@ -60,7 +60,7 @@
include "constants.h"
- integer :: npoin,nspec,nelemabs,numat,iglob_source,ispec_selected_source,is_proc_source,source_type,it,NSTEP
+ integer :: npoin,nspec,nelemabs,numat,ispec_selected_source,is_proc_source,source_type,it,NSTEP
logical :: anyabs,assign_external_model,initialfield,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,add_Bielak_conditions
@@ -323,7 +323,6 @@
!--- left absorbing boundary
if(codeabs(ILEFT,ispecabs)) then
-!!$ if(.false.) then
i = 1
@@ -392,7 +391,6 @@
!--- right absorbing boundary
if(codeabs(IRIGHT,ispecabs)) then
-!!$ if(.false.) then
i = NGLLX
Modified: seismo/2D/SPECFEM2D/trunk/paco_beyond_critical.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/paco_beyond_critical.f90 2008-05-03 17:24:53 UTC (rev 11912)
+++ seismo/2D/SPECFEM2D/trunk/paco_beyond_critical.f90 2008-05-04 00:41:04 UTC (rev 11913)
@@ -52,12 +52,12 @@
double precision :: X, Z, xmin, xmax, zmin, zmax
integer :: inode
- double complex CAKA,CAQA,UI,UR
- double complex UX,UZ,SX,SZ,SXZ,A2,B2,AL,AK,AM
+ complex(selected_real_kind(15,300)) :: CAKA,CAQA,UI,UR
+ complex(selected_real_kind(15,300)) :: UX,UZ,SX,SZ,SXZ,A2,B2,AL,AK,AM
- double complex :: TX,TZ
+ complex(selected_real_kind(15,300)) :: TX,TZ
- double complex, dimension(:),allocatable::Field_Ux,Field_Uz,Field_Tx,Field_Tz
+ complex(selected_real_kind(15,300)), dimension(:),allocatable::Field_Ux,Field_Uz,Field_Tx,Field_Tz
double precision :: TS
@@ -227,7 +227,7 @@
TOTO=0.0d0
- CALL DESFXY(TOTO,TOTO,source_type,UX,UZ,SX,SZ,SXZ,A1,B1,A2,B2,AL,AK,AM,RLM,ANU)
+ CALL DESFXY(TOTO,TOTO,source_type,UX,UZ,SX,SZ,SXZ,A1,B1,A2,B2,AL,AK,AM,RLM)
!! DK DK write the frequency seismograms
TX = SX *VNX+SXZ*VNZ
@@ -267,7 +267,7 @@
IF (source_type==2) CALL ONDASS(GAMR,AKA,AQA,A1,B1,A2,B2,AL,AK,AM,ANU,BEALF)
IF (source_type==3) CALL ONDASR(AQA,A1,B1,A2,B2,AL,AK,AM,ANU,BEALF)
- CALL DESFXY(X,Z,source_type,UX,UZ,SX,SZ,SXZ,A1,B1,A2,B2,AL,AK,AM,RLM,ANU)
+ CALL DESFXY(X,Z,source_type,UX,UZ,SX,SZ,SXZ,A1,B1,A2,B2,AL,AK,AM,RLM)
!! DK DK write the frequency seismograms
TX = SX *VNX+SXZ*VNZ
@@ -351,37 +351,36 @@
!---
-SUBROUTINE DESFXY(X,Z,ICAS,UX,UZ,SX,SZ,SXZ,A1,B1,A2,B2,AL,AK,AM,RLM,ANU)
+SUBROUTINE DESFXY(X,Z,ICAS,UX,UZ,SX,SZ,SXZ,A1,B1,A2,B2,AL,AK,AM,RLM)
implicit none
- double precision A1,B1,RLM,ANU,X,Z
+ double precision A1,B1,RLM,X,Z
integer ICAS
- double complex UX,UZ,SX,SZ,SXZ,A2,B2,AL,AK,AM
- double complex UI,FAC
- double complex AUX1,AUX2,FI1,FI2,PS1,PS2
+ complex(selected_real_kind(15,300)) :: UX,UZ,SX,SZ,SXZ,A2,B2,AL,AK,AM
+ complex(selected_real_kind(15,300)) :: UI,FAC
+ complex(selected_real_kind(15,300)) :: AUX1,AUX2,FI1,FI2,PS1,PS2
-
UI=(0.0d0,1.0d0)
if (A1/=0.0d0) then
- AUX1=A1*CDEXP(UI*(AM*Z-AL*X)) !campo P incidente
+ AUX1=A1*EXP(UI*(AM*Z-AL*X)) !campo P incidente
else
AUX1=CMPLX(0.0d0)
end if
if (A2/=0.0d0) then
- AUX2=A2*CDEXP(-UI*(AM*Z+AL*X)) *1.0d0 !campo P reflejado
+ AUX2=A2*EXP(-UI*(AM*Z+AL*X)) *1.0d0 !campo P reflejado
else
AUX2=CMPLX(0.0d0)
end if
FI1=AUX1+AUX2
FI2=AUX1-AUX2
if (B1/=0.0d0) then
- AUX1=B1*CDEXP(UI*(AK*Z-AL*X)) !campo S incidente
+ AUX1=B1*EXP(UI*(AK*Z-AL*X)) !campo S incidente
else
AUX1=CMPLX(0.0d0)
end if
if (B2/=0.0d0) then
- AUX2=B2*CDEXP(-UI*(AK*Z+AL*X)) *1.0d0 !campo S reflejado
+ AUX2=B2*EXP(-UI*(AK*Z+AL*X)) *1.0d0 !campo S reflejado
else
AUX2=CMPLX(0.0d0)
end if
@@ -415,7 +414,7 @@
implicit none
double precision CA,CB,A,B
- double complex FA,FB,ZER,UI
+ complex(selected_real_kind(15,300)) :: FA,FB,ZER,UI
ZER=(0.0d0,0.0d0)
UI=(0.0d0,1.0d0)
@@ -442,7 +441,7 @@
implicit none
- double complex FA,FB,A2,B2,DEN,AUX
+ complex(selected_real_kind(15,300)) :: FA,FB,A2,B2,DEN,AUX
AUX=FB*FB-1.0d0
DEN=4.0d0*FA*FB+AUX*AUX
@@ -457,7 +456,7 @@
implicit none
double precision A1,B1,ANU,CA,CB,GP,AQB,BEALF
- double complex A2,B2,FA,FB,ZER,AL,AK,AM
+ complex(selected_real_kind(15,300)) :: A2,B2,FA,FB,ZER,AL,AK,AM
ZER=(0.0d0,0.0d0)
BEALF=SQRT((1.0d0-2.0d0*ANU)/2.0d0/(1.0d0-ANU))
@@ -493,7 +492,7 @@
implicit none
double precision A1,B1,ANU,CA,CB,GS,AQB,BEALF,AKB
- double complex A2,B2,FA,FB,ZER,AL,AK,AM
+ complex(selected_real_kind(15,300)) :: A2,B2,FA,FB,ZER,AL,AK,AM
ZER=(0.0d0,0.0d0)
BEALF=SQRT((1.0d0-2.0d0*ANU)/2.0d0/(1.0d0-ANU))
@@ -548,7 +547,7 @@
implicit none
double precision A1,B1,ANU,CA,CB,AQB,BEALF,ba2
- double complex A2,B2,FA,FB,ZER,AL,AK,AM
+ complex(selected_real_kind(15,300)) :: A2,B2,FA,FB,ZER,AL,AK,AM
double precision, external :: crb
Modified: seismo/2D/SPECFEM2D/trunk/paco_convolve_fft.f90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/paco_convolve_fft.f90 2008-05-03 17:24:53 UTC (rev 11912)
+++ seismo/2D/SPECFEM2D/trunk/paco_convolve_fft.f90 2008-05-04 00:41:04 UTC (rev 11913)
@@ -15,22 +15,18 @@
integer :: NFREC,N,NSTEP
- double complex, dimension(NFREC+1) :: Field
+ complex(selected_real_kind(15,300)), dimension(NFREC+1) :: Field
- double complex CR(2*NFREC)
+ complex(selected_real_kind(15,300)) :: CR(2*NFREC)
double precision, dimension(NSTEP) :: output_field
- integer :: J,i,label
+ integer :: J,label
- integer :: beg_UV, end_UV
+ double precision :: AN,FUN,RAIZ,dt,tp,ts
- double precision :: AN,FUN,RAIZ,X,Z,dt,tp,ts
-
double precision, external :: RIC, deRIC, de2RIC
- integer inode
-
N=2*NFREC
AN = N
@@ -67,22 +63,21 @@
implicit none
- integer NSTEP, j,jn,nar,N,jstep,jbegin,jend,label,nfrec,mult,delay
+ integer NSTEP, j,jn,N,label,nfrec,mult,delay
double precision :: RAIZ
- double complex VC
+ complex(selected_real_kind(15,300)) :: VC
double precision VT(2*NFREC)
- double precision :: filt,eta,dt
+ double precision :: filt,dt
-
double precision, dimension(NSTEP) :: output_field
- double complex, dimension(NFREC+1) :: V
+ complex(selected_real_kind(15,300)), dimension(NFREC+1) :: V
- double complex CY(2*NFREC),CR(2*NFREC)
+ complex(selected_real_kind(15,300)) :: CY(2*NFREC),CR(2*NFREC)
N=2*NFREC
@@ -152,7 +147,7 @@
include "constants.h"
double precision :: A,A_dot,deRIC,tp,ts,dt
- integer :: label,j
+ integer :: j
A=PI*(dt*(J-1)-ts)/tp
A=A*A
@@ -197,7 +192,7 @@
double precision SC
- double complex CX(LX),CARG,CW,CTEMP
+ complex(selected_real_kind(15,300)) :: CX(LX),CARG,CW,CTEMP
double precision SIGNI
@@ -222,7 +217,7 @@
ISTEP=2*L
DO M=1,L
CARG=(0.0d0,1.0d0)*(PI*SIGNI*(M-1))/L
- CW=CDEXP(CARG)
+ CW=EXP(CARG)
DO I=M,LX,ISTEP
CTEMP=CW*CX(I+L)
CX(I+L)=CX(I)-CTEMP
Modified: seismo/2D/SPECFEM2D/trunk/specfem2D.F90
===================================================================
--- seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2008-05-03 17:24:53 UTC (rev 11912)
+++ seismo/2D/SPECFEM2D/trunk/specfem2D.F90 2008-05-04 00:41:04 UTC (rev 11913)
@@ -177,7 +177,7 @@
! curl in an element
real(kind=CUSTOM_REAL), dimension(NGLLX,NGLLX) :: curl_element
- integer :: i,j,k,l,it,irec,ipoin,ip,id,nbpoin,inump,n,ispec,npoin,npgeo,iglob
+ integer :: i,j,k,l,it,irec,ipoin,ip,id,n,ispec,npoin,npgeo,iglob
logical :: anyabs
double precision :: dxd,dzd,dcurld,valux,valuz,valcurl,hlagrange,rhol,cosrot,sinrot,xi,gamma,x,z
@@ -207,7 +207,7 @@
! for acoustic medium
real(kind=CUSTOM_REAL), dimension(:), allocatable :: potential_dot_dot_acoustic,potential_dot_acoustic,potential_acoustic
real(kind=CUSTOM_REAL), dimension(:), allocatable :: rmass_inverse_elastic,rmass_inverse_acoustic
- double precision, dimension(:), allocatable :: density,displread,velocread,accelread
+ double precision, dimension(:), allocatable :: density
double precision, dimension(:), allocatable :: vp_display
@@ -376,13 +376,11 @@
double precision, dimension(2) :: A_plane, B_plane, C_plane
double precision :: PP, PS, SP, SS, z0_source, x0_source, xmax, xmin, zmax, zmin, time_offset
-!over critical angle
+! over critical angle
integer , dimension(:), allocatable :: left_bound,right_bound,bot_bound
double precision , dimension(:,:), allocatable :: v0x_left,v0z_left,v0x_right,v0z_right,v0x_bot,v0z_bot
double precision , dimension(:,:), allocatable :: t0x_left,t0z_left,t0x_right,t0z_right,t0x_bot,t0z_bot
- integer count_left,count_right,count_bot,itime,ibegin,iend
-!!$ double precision :: a1,b1,a2,b2,a3,b3,sin0,cos0,Zr,Zr_dot,Zr_dot_dot,Zi,Zi_dot,Zi_dot_dot,Cte
- double precision :: X_temp,Z_temp,ignore
+ integer count_left,count_right,count_bot,ibegin,iend
logical :: over_critical_angle
!***********************************************************************
@@ -2185,11 +2183,10 @@
! first call, computation on outer elements, absorbing conditions and source
call compute_forces_acoustic(npoin,nspec,nelemabs,numat, &
- iglob_source,ispec_selected_source,is_proc_source,source_type,it,NSTEP,anyabs, &
- assign_external_model,initialfield,ibool,kmato,numabs, &
+ anyabs,assign_external_model,ibool,kmato,numabs, &
elastic,codeabs,potential_dot_dot_acoustic,potential_dot_acoustic, &
potential_acoustic,density,elastcoef,xix,xiz,gammax,gammaz,jacobian, &
- vpext,rhoext,source_time_function,hprime_xx,hprimewgll_xx, &
+ vpext,rhoext,hprime_xx,hprimewgll_xx, &
hprime_zz,hprimewgll_zz,wxgll,wzgll, &
ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
jbegin_left,jend_left,jbegin_right,jend_right, &
@@ -2279,11 +2276,10 @@
! second call, computation on inner elements
if(any_acoustic) then
call compute_forces_acoustic(npoin,nspec,nelemabs,numat, &
- iglob_source,ispec_selected_source,is_proc_source,source_type,it,NSTEP,anyabs, &
- assign_external_model,initialfield,ibool,kmato,numabs, &
+ anyabs,assign_external_model,ibool,kmato,numabs, &
elastic,codeabs,potential_dot_dot_acoustic,potential_dot_acoustic, &
potential_acoustic,density,elastcoef,xix,xiz,gammax,gammaz,jacobian, &
- vpext,rhoext,source_time_function,hprime_xx,hprimewgll_xx, &
+ vpext,rhoext,hprime_xx,hprimewgll_xx, &
hprime_zz,hprimewgll_zz,wxgll,wzgll, &
ibegin_bottom,iend_bottom,ibegin_top,iend_top, &
jbegin_left,jend_left,jbegin_right,jend_right, &
@@ -2311,6 +2307,25 @@
if(any_acoustic) then
potential_dot_dot_acoustic = potential_dot_dot_acoustic * rmass_inverse_acoustic
+
+! --- add the source
+ if(.not. initialfield) then
+! if this processor carries the source and the source element is acoustic
+ if (is_proc_source == 1 .and. .not. elastic(ispec_selected_source)) then
+! collocated force
+! beware, for acoustic medium, source is a pressure source
+! the sign is negative because pressure p = - Chi_dot_dot therefore we need
+! to add minus the source to Chi_dot_dot to get plus the source in pressure
+ if(source_type == 1) then
+ potential_dot_dot_acoustic(iglob_source) = potential_dot_dot_acoustic(iglob_source) - source_time_function(it)
+
+! moment tensor
+ else if(source_type == 2) then
+ call exit_MPI('cannot have moment tensor source in acoustic element')
+ endif
+ endif ! if this processor carries the source and the source element is acoustic
+ endif ! if not using an initial field
+
potential_dot_acoustic = potential_dot_acoustic + deltatover2*potential_dot_dot_acoustic
! free surface for an acoustic medium
@@ -2327,7 +2342,7 @@
! first call, computation on outer elements, absorbing conditions and source
if(any_elastic) &
- call compute_forces_elastic(npoin,nspec,nelemabs,numat,iglob_source, &
+ call compute_forces_elastic(npoin,nspec,nelemabs,numat, &
ispec_selected_source,is_proc_source,source_type,it,NSTEP,anyabs,assign_external_model, &
initialfield,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,angleforce,deltatcube, &
deltatfourth,twelvedeltat,fourdeltatsquare,ibool,kmato,numabs,elastic,codeabs, &
@@ -2420,9 +2435,9 @@
endif
#endif
-! second call, computation on inner elements and update of
+! second call, computation on inner elements and update
if(any_elastic) &
- call compute_forces_elastic(npoin,nspec,nelemabs,numat,iglob_source, &
+ call compute_forces_elastic(npoin,nspec,nelemabs,numat, &
ispec_selected_source,is_proc_source,source_type,it,NSTEP,anyabs,assign_external_model, &
initialfield,TURN_ATTENUATION_ON,TURN_ANISOTROPY_ON,angleforce,deltatcube, &
deltatfourth,twelvedeltat,fourdeltatsquare,ibool,kmato,numabs,elastic,codeabs, &
More information about the cig-commits
mailing list