Skip to content
GitLab
Commit 2c849251 authored by Yann CAPDEVILLE's avatar Yann CAPDEVILLE
Browse files

On branch srcinv 

commit intermediare. Il y a beaucoup de choses qui ne vont pas.
Debut du travail sur le type nms
	modified:   nms/global_main.f90
	modified:   nms/makefile
	modified:   nms/module_modes.f90
	modified:   nms/module_srcinv.F90
	modified:   nms/srcdirect.f90
	modified:   nms/time_function.f90
parent 6b78fdd1
Hide whitespace changes
Inline Side-by-side
nms/global_main.f90
View file @ 2c849251
......@@ -119,6 +119,7 @@ contains
jj=ic+(i-1)*NBC
src%i2rad(jj)=is2rad_(i)
!rotate xyz to rtp here?
src%Mtmp(: ,jj)=0.d0
src%Mtmp(ic,jj)=1.d0
t=srccoord_(2,i)*deg2rad
p=srccoord_(3,i)*deg2rad
......@@ -130,7 +131,7 @@ contains
else
srcinv=.false.
endif
!
!
call load_nms()
call init_rcv_src()
if (station_response) call load_station_response()
......
nms/makefile
View file @ 2c849251
......@@ -10,27 +10,30 @@ LDFLAGS = $(Opt)
SRCS = nms.f90
OBJS = $(LOBJ)xerbla.o $(LOBJ)def_gparam.o $(LOBJ)lsame.o $(LOBJ)zgemv.o \
OBLAS = $(LOBJ)lsame.o $(LOBJ)zgemv.o
OBLASN =
OBJS = $(LOBJ)xerbla.o $(LOBJ)def_gparam.o $(OBLAS) \
$(LOBJ)dfour1.o $(LOBJ)nrtype.o $(LOBJ)nrutil.o $(LOBJ)module_modetypes.o $(LOBJ)global_main.o $(LOBJ)time_function.o \
$(LOBJ)module_correction_modes.o $(LOBJ)module_modele.o $(LOBJ)module_spline.o \
$(LOBJ)util_fctp.o $(LOBJ)module_modes.o $(LOBJ)nms.o
OBJSI = $(LOBJ)xerbla.o $(LOBJ)def_gparam.o $(LOBJ)lsame.o $(LOBJ)zgemv.o \
OBJSI = $(LOBJ)xerbla.o $(LOBJ)def_gparam.o $(OBLASN) \
$(LOBJ)dfour1.o $(LOBJ)nrtype.o $(LOBJ)nrutil.o $(LOBJ)module_modetypes.o $(LOBJ)global_main.o $(LOBJ)time_function.o \
$(LOBJ)module_correction_modes.o $(LOBJ)module_modele.o $(LOBJ)module_spline.o \
$(LOBJ)util_fctp.o $(LOBJ)module_modes.o $(LOBJ)speclib.o $(LOBJ)module_srcinv.o $(LOBJ)srcinv.o
OBJSD = $(LOBJ)xerbla.o $(LOBJ)def_gparam.o $(LOBJ)lsame.o $(LOBJ)zgemv.o \
OBJSD = $(LOBJ)xerbla.o $(LOBJ)def_gparam.o $(OBLAS) \
$(LOBJ)dfour1.o $(LOBJ)nrtype.o $(LOBJ)nrutil.o $(LOBJ)module_modetypes.o $(LOBJ)global_main.o $(LOBJ)time_function.o \
$(LOBJ)module_correction_modes.o $(LOBJ)module_modele.o $(LOBJ)module_spline.o \
$(LOBJ)util_fctp.o $(LOBJ)module_modes.o $(LOBJ)speclib.o $(LOBJ)module_srcinv.o $(LOBJ)srcdirect.o
OBJSC = $(LOBJ)xerbla.o $(LOBJ)def_gparam.o $(LOBJ)lsame.o $(LOBJ)zgemv.o \
OBJSC = $(LOBJ)xerbla.o $(LOBJ)def_gparam.o $(OBLAS) \
$(LOBJ)dfour1.o $(LOBJ)nrtype.o $(LOBJ)nrutil.o $(LOBJ)module_modetypes.o $(LOBJ)global_main.o $(LOBJ)time_function.o \
$(LOBJ)module_correction_modes.o $(LOBJ)module_modele.o $(LOBJ)module_spline.o \
$(LOBJ)util_fctp.o $(LOBJ)module_modes.o $(LOBJ)speclib.o $(LOBJ)module_srcinv.o $(LOBJ)cutsrc.o
OBJSH = $(LOBJ)xerbla.o $(LOBJ)def_gparam.o $(LOBJ)lsame.o $(LOBJ)zgemv.o \
OBJSH = $(LOBJ)xerbla.o $(LOBJ)def_gparam.o $(OBLAS) \
$(LOBJ)dfour1.o $(LOBJ)nrtype.o $(LOBJ)nrutil.o $(LOBJ)module_modetypes.o $(LOBJ)global_main.o $(LOBJ)time_function.o \
$(LOBJ)module_correction_modes.o $(LOBJ)module_modele.o $(LOBJ)module_spline.o \
$(LOBJ)util_fctp.o $(LOBJ)module_modes.o $(LOBJ)speclib.o $(LOBJ)module_srcinv.o $(LOBJ)module_hsrc.o $(LOBJ)hsrc.o
......
nms/module_modes.f90
View file @ 2c849251
......@@ -113,7 +113,7 @@ contains
dtf2 = dtf/coef_nbt
!
!ouverture de fichiers de fonctions propres
print*,'Opening eigenfunction files ...',eigenfileS,eigenfileT
print*,'Opening eigenfunction files ...(',TRIM(eigenfileS),' and ',TRIM(eigenfileT),')'
call open_fctp(eigenfileS,eigenfileT,rad_recepteur=receiver_rad)
!RA is now initialized
call set_src_coord()
......@@ -191,7 +191,7 @@ contains
! rad(2:NBE+1)=src_coord(1,:)
rad(1:rcv%mb%NBRR)=rcv%mb%rad(:)
rad(rcv%mb%NBRR+1:rcv%mb%NBRR+src%mb%NBRR)=src%mb%rad(:)
!
print*,'->for S modes ...'
do l=0,LmaxS
do n=0,NmaxS(l)
......@@ -268,7 +268,7 @@ contains
!!$ f3=f3_(1)
!!$ f4=f4_(1)
tmax_trace=dt*NBT
call init_time_function(sourcew,src%source_code,src%support,src%f1,src%f1,src%f2,src%f3,src%f4,src%t0,src%amp_source,tmax_trace)
call init_time_function(sourcew,src%source_code,src%support,src%f1,src%f1,src%f2,src%f3,src%f4,src%t0,src%amp_source,tmax_trace)
fmax=fh4
wmax1=fh1*2._DP*PI
wmax2=fh2 *2._DP*PI
......@@ -487,7 +487,7 @@ contains
clean=.true.
endif
traces(:,:,:,:)=cmplx(0.d0,0.d0)
print*,'Computing RkN and SkN terms ...',src%NBE,rcv%NBR
print*,'Computing RkN and SkN terms ...',rcv%NBR,src%NBE!,src%source_force
call init_RSkN(channel_,clean)
print*,'Computing green function traces ...'
do ichunk=0,nbchunk-1
......@@ -857,10 +857,10 @@ contains
scale =1._DP/sqrt(rhobar*ra**3)
dscale=1._DP/sqrt(rhobar*ra**3)/ra
!
mega=(3*rcv%mb%NBRR*nbcomp*(real(Nmax_allS+1)*(LmaxS+1)+real((Nmax_allT+1)*(LmaxT+1))) + &
5*src%mb%NBRR*(real(Nmax_allS+1)*(LmaxS+1)+real((Nmax_allT+1)*(LmaxT+1))) )/1024./1204.
mega=(3*rcv%NBR*nbcomp*(real(Nmax_allS+1)*(LmaxS+1)+real((Nmax_allT+1)*(LmaxT+1))) + &
5*src%NBE*(real(Nmax_allS+1)*(LmaxS+1)+real((Nmax_allT+1)*(LmaxT+1))) )/1024./1204.
allocate(SkN(-2:2 ,0:Nmax_all,0:Lmax_all,2,src%mb%NBRR))
allocate(SkN(-2:2 ,0:Nmax_all,0:Lmax_all,2,src%NBE))
SkN=0.0_DP
print*,'->SkN ...',channel_
call SkNfun(SkN)
......@@ -899,9 +899,8 @@ contains
real(DP) :: ff,xx,vor,dmtp,egu,egv,degu,degv,rr,x
real(DP) :: zz,wor,fw,dfw,u,v,xu,xv,vir,vit,vip
real(DP), dimension(6) :: am
integer :: is,n,l
integer :: is,n,l,isr
!
if (src%source_force) then
do l=0,LmaxS
do n=0,NmaxS(l)
......@@ -948,52 +947,54 @@ contains
!dans le cas d'une force, il n'y a pas d'intergation par partie, tout - pour tout le monde: (? ... j'ai quand meme un doute)
S=-S
else
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(l,n,is,rr,egu,egv,degu,am,ff,xx,vor,dmtp)
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(l,n,is,isr,rr,egu,egv,degu,am,ff,xx,vor,dmtp)
do l=0,LmaxS
do n=0,NmaxS(l)
do is=1,src%mb%NBRR
! rr=src_coord(1,is)/RA
rr=src%mb%rad(is)/RA
egu =src%mb%UU (is,n,l)
egv =src%mb%VV (is,n,l)
degu =src%mb%UUp(is,n,l)
degv =src%mb%VVp(is,n,l)
if (homogenization>=0) call corrmodesS(is,eigw(n,l,1),l,rr,egu,egv,degu,degv)
am(:)=src%Msrc(:,is)
!
ff =(2.0d0*egu-al(l)*egv)/rr
xx =-b1(l)*(degv+(egu-egv)/rr)
vor =b2(l)*egv/rr
dmtp=am(3)-am(2)
S( 0,n,l,1,is)=-b0(l)*(degu*am(1)+ff*0.5*(am(2)+am(3)))*dscale
S( 1,n,l,1,is)=xx*dcmplx(-am(4),am(5))*dscale
S(-1,n,l,1,is)=xx*dcmplx( am(4),am(5))*dscale
S( 2,n,l,1,is)=vor*dcmplx(dmtp, 2.0d0*am(6))*dscale
S(-2,n,l,1,is)=vor*dcmplx(dmtp,-2.0d0*am(6))*dscale
enddo
do is=1,src%NBE
! rr=src_coord(1,is)/RA
isr=src%i2rad(is)
rr=src%mb%rad(isr)/RA
egu =src%mb%UU (isr,n,l)
egv =src%mb%VV (isr,n,l)
degu =src%mb%UUp(isr,n,l)
degv =src%mb%VVp(isr,n,l)
if (homogenization>=0) call corrmodesS(isr,eigw(n,l,1),l,rr,egu,egv,degu,degv)
am(:)=src%Msrc(:,is)
!
ff =(2.0d0*egu-al(l)*egv)/rr
xx =-b1(l)*(degv+(egu-egv)/rr)
vor =b2(l)*egv/rr
dmtp=am(3)-am(2)
S( 0,n,l,1,is)=-b0(l)*(degu*am(1)+ff*0.5*(am(2)+am(3)))*dscale
S( 1,n,l,1,is)=xx*dcmplx(-am(4),am(5))*dscale
S(-1,n,l,1,is)=xx*dcmplx( am(4),am(5))*dscale
S( 2,n,l,1,is)=vor*dcmplx(dmtp, 2.0d0*am(6))*dscale
S(-2,n,l,1,is)=vor*dcmplx(dmtp,-2.0d0*am(6))*dscale
enddo
enddo
enddo
!$OMP END PARALLEL DO
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(l,n,is,rr,fw,dfw,am,zz,wor,dmtp)
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(l,n,is,isr,rr,fw,dfw,am,zz,wor,dmtp)
do l=1,LmaxT
do n=0,NmaxT(l)
do is=1,src%mb%NBRR
do is=1,src%NBE
! rr=src_coord(1,is)/RA
rr=src%mb%rad(is)/RA
fw =src%mb%WW (is,n,l)
dfw =src%mb%WWp(is,n,l)
if (homogenization>=0) call corrmodesT(is,eigw(n,l,2),l,rr,fw,dfw)
am(:) =src%Msrc(:,is)
!
zz =-b1(l)*(dfw-fw/rr)
wor = b2(l)*fw/rr
dmtp =am(2)-am(3)
S( 1,n,l,2,is)=zz *dcmplx( am(5),am(4))*dscale
S(-1,n,l,2,is)=zz *dcmplx( -am(5),am(4))*dscale
S( 2,n,l,2,is)=wor*dcmplx(2.0d0*am(6), dmtp)*dscale
S(-2,n,l,2,is)=wor*dcmplx(2.0d0*am(6),-dmtp)*dscale
enddo
enddo
isr=src%i2rad(is)
rr=src%mb%rad(isr)/RA
fw =src%mb%WW (isr,n,l)
dfw =src%mb%WWp(isr,n,l)
if (homogenization>=0) call corrmodesT(isr,eigw(n,l,2),l,rr,fw,dfw)
am(:) =src%Msrc(:,is)
!
zz =-b1(l)*(dfw-fw/rr)
wor = b2(l)*fw/rr
dmtp =am(2)-am(3)
S( 1,n,l,2,is)=zz *dcmplx( am(5),am(4))*dscale
S(-1,n,l,2,is)=zz *dcmplx( -am(5),am(4))*dscale
S( 2,n,l,2,is)=wor*dcmplx(2.0d0*am(6), dmtp)*dscale
S(-2,n,l,2,is)=wor*dcmplx(2.0d0*am(6),-dmtp)*dscale
enddo
enddo
enddo
!$OMP END PARALLEL DO
endif
......@@ -1190,6 +1191,20 @@ contains
enddo
enddo
!$OMP END PARALLEL DO
! open(11,file='source.gnu')
print*,'coucou dd',dt,NBT
do i=1,NBT
tt=tmin+dt*(i-1)
is=1
amp=time_func(sourcew,tt-src%t_delay(is))
write(344,*) tt,amp
enddo
close(344)
! open(11,file='source.gnu')
! do i=1,1000
! write(11,*) (i-1)*2.,time_func(sourcew,(i-1)*2.d0)
! enddo
! close(11)
!--------------------------------------------------
end subroutine homogenize_source_b
!--------------------------------------------------
......@@ -1438,7 +1453,7 @@ contains
do bN=-1,1
sd = czero
do bNp=-2,2
sd=sd+dsr(bN,bNp,l)*SkN(bNp,n,l,q,irs)*csd(bNp)
sd=sd+dsr(bN,bNp,l)*SkN(bNp,n,l,q,is)*csd(bNp)
enddo
rs=rs+sd*RkN(bN,ic,n,l,q,irr)*csr(bN)
enddo
......
nms/module_srcinv.F90
View file @ 2c849251
......@@ -92,7 +92,7 @@ module module_srcinv
read(11,*) geom%damp_coef
read(11,*) !time tapper length:
read(11,*) geom%timetapper
geom%save_green=.false.
geom%save_green=.true.
if (geom%timetapper>0.d0) then
geom%applymask=.true.
else
......@@ -344,8 +344,8 @@ subroutine build_green_fcts(geom)
type(meshinv), intent(inout) :: geom
!
real (DP), dimension(:,:,:,:), pointer :: green
integer :: logunit,i,j,k,ii,jj,kk,ic,ip,ipint,it,ir,ips,icp,idelta,irg,a
real(DP) :: tshift,t,ta,tb,coef,time_damp
integer :: logunit,i,j,k,ii,jj,kk,ic,ip,ipint,it,ir,ips,icp,idelta,irg,a,iploc
real(DP) :: tshift,t,ta,tb,coef,time_damp,test
real(DP), dimension(:), allocatable :: Mc
character(len=30) :: name
logunit=1327
......@@ -393,12 +393,20 @@ subroutine build_green_fcts(geom)
call init_modes(logunit)
call copy_timefunction(sourcew,geom%base)
call get_and_sum(traces_time_out=green)
do i=1,geom%ring(1)%ngllv_int
ipint=ipnumber(geom,1,1,1,1,i,1)
!print*,'coucou vv',ipint,src%i2rad(ipint),geom%mesh_rad(src%i2rad(ipint)),src%coord(1,ipint),sngl(src%Mtmp(:,ipint))
do j=1,NBT
write(2000+i,*) (j-1)*dt,green(j,1,1,ipint)
enddo
close(2000+i)
enddo
!
allocate(geom%green(geom%NBTtrace,3,geom%NBR,geom%NBP))
geom%green(:,:,:,:)=0.d0
ip=0
do irg=1,geom%nring
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(kk,jj,ii,ic,ip,ir,icp,k,j,i,ipint,Mc)
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(kk,jj,ii,ic,ip,ir,icp,k,j,i,a,ipint,Mc,iploc)
allocate(Mc(6*geom%ring(irg)%ngllh_int*geom%ring(irg)%ngllh_int*geom%ring(irg)%ngllv_int))
!$OMP DO COLLAPSE(4)
do kk=1,geom%ring(irg)%ngllv_shape
......@@ -408,19 +416,22 @@ subroutine build_green_fcts(geom)
ip=ipnumber(geom,ic,irg,ii,jj,kk,2)
call get_Mc(geom,ic,irg,ii,jj,kk,Mc)
! do ir=1,geom%NBR
!test=0.d0
iploc=0
do k=1,geom%ring(irg)%ngllv_int
do j=1,geom%ring(irg)%ngllh_int
do i=1,geom%ring(irg)%ngllh_int
do a=1,geom%NBCM
iploc=iploc+1
ipint=ipnumber(geom,a,irg,i,j,k,1)
geom%green(:,:,:,ip)=geom%green(:,:,:,ip) &
+green(:,:,:,ipint)*Mc(ipint)
+green(:,:,:,ipint)*Mc(iploc)
enddo
! geom%green(:,:,ir,ip)=geom%green(:,:,ir,ip) &
! +green(:,:,ir,ipint)*geom%ring(irg)%w_int(i,j,k)*geom%ring(irg)%shape_fcts(i,j,k,ii,jj,kk)
enddo
enddo
enddo
enddo
enddo
! enddo
enddo
enddo
......@@ -430,7 +441,14 @@ subroutine build_green_fcts(geom)
deallocate(Mc)
!$OMP END PARALLEL
enddo
do i=1,geom%ring(1)%ngllh_shape
ip=ipnumber(geom,1,1,i,1,1,2)
do j=1,NBT
write(3000+i,*) (j-1)*dt,geom%green(j,1,1,ip)
enddo
close(3000+i)
enddo
!stop
!$OMP PARALLEL DO COLLAPSE(2) DEFAULT(SHARED) PRIVATE(it,ips,ip,tshift,ic,ir)
do it=2,geom%NBTinv
do ips=1,geom%NBPs
......@@ -475,6 +493,7 @@ subroutine build_green_fcts(geom)
!$OMP END PARALLEL DO
endif
close(logunit)
print*,'coucou jj',geom%save_green
if (geom%save_green) then
print*,'writing down green functions for debug (only for the first receiver)...'
print*,geom%NBTinv,geom%NBPs
......@@ -503,7 +522,7 @@ end subroutine build_green_fcts
subroutine sum_green_fcts(geom,coefin,dir)
!----------------------------------------------------------------------------
use def_gparam
use global_main_param, only: init_global_main,NBT,dt,nbcomp,rotation,stations_name,comp,compr,src,rcv
use global_main_param, only: init_global_main,NBT,dt,nbcomp,rcv,comp,compr,src,rcv,rotation
use modes,only: sourcew,get_and_sum,init_modes,wtcoef
use time_function
implicit none
......@@ -653,9 +672,9 @@ subroutine sum_green_fcts(geom,coefin,dir)
do ir=1,geom%NBR
do ic=1,3
if (rotation) then
write(name,'("U",a1,"_",a4)') compr(ic),stations_name(ir)
write(name,'("U",a1,"_",a4)') compr(ic),rcv%name(ir)
else
write(name,'("U",a1,"_",a4)') comp(ic) ,stations_name(ir)
write(name,'("U",a1,"_",a4)') comp(ic) ,rcv%name(ir)
endif
open(17,file=add_dir(name,dir))
call fftresample(traces(:,ic,ir),0._DP,geom%dttrace,geom%NBTtrace,tmp,geom%dtout,geom%NBTout)
......@@ -684,17 +703,18 @@ subroutine get_Mc(geom,ic,irg,ii,jj,kk,Mc)
integer, intent(in) :: ic,irg,ii,jj,kk
real(DP), dimension(:), intent(out) :: Mc
!
integer :: a,i,j,k,ipint,b,irad,ijk
integer :: a,i,j,k,b,irad,ijk,iploc
real(DP), parameter:: c2=dsqrt(2.d0)/2.d0
real(DP), dimension(6) :: epsi
!
ijk=geom%ring(irg)%NBPic
do k=1,geom%ring(irg)%ngllv_int
iploc=0
do k=1,geom%ring(irg)%ngllv_int
do j=1,geom%ring(irg)%ngllh_int
do i=1,geom%ring(irg)%ngllh_int
ijk=ijk+1
ijk=ijk+1
do a=1,geom%NBCM
ipint=ipnumber(geom,a,irg,i,j,k,1)
iploc=iploc+1
epsi(:)=0.d0
select case(ic) !psi component
case (1) !x
......@@ -707,15 +727,16 @@ subroutine get_Mc(geom,ic,irg,ii,jj,kk,Mc)
epsi(6)=geom%ring(irg)%deriv(3,2,i,j,k)*geom%ring(irg)%dshape_fcts(2,i,j,k,ii,jj,kk)*c2
case (3) !z
epsi(3)=geom%ring(irg)%deriv(3,3,i,j,k)*geom%ring(irg)%dshape_fcts(3,i,j,k,ii,jj,kk)
epsi(5)=geom%ring(irg)%deriv(2,3,i,j,k)*geom%ring(irg)%dshape_fcts(1,i,j,k,ii,jj,kk)*c2
epsi(6)=geom%ring(irg)%deriv(3,3,i,j,k)*geom%ring(irg)%dshape_fcts(1,i,j,k,ii,jj,kk)*c2
epsi(5)=geom%ring(irg)%deriv(1,3,i,j,k)*geom%ring(irg)%dshape_fcts(3,i,j,k,ii,jj,kk)*c2
epsi(6)=geom%ring(irg)%deriv(2,3,i,j,k)*geom%ring(irg)%dshape_fcts(3,i,j,k,ii,jj,kk)*c2
end select
irad=geom%iglob2irad(ijk)
Mc(ipint)=0.d0
Mc(iploc)=0.d0
do b=1,6
Mc(ipint)=Mc(ipint)+geom%CIJ(b,a,irad)*epsi(b)
Mc(iploc)=Mc(iploc)+geom%CIJ(b,a,irad)*epsi(b)
enddo
Mc(ipint)=Mc(ipint)*geom%ring(irg)%w_int(i,j,k) !w_int already contains the jacobian
Mc(iploc)=Mc(iploc)*geom%ring(irg)%w_int(i,j,k) !w_int already contains the jacobian
!print*,'coucou zz',iploc,Mc(iploc)
enddo
enddo
enddo
......@@ -734,7 +755,7 @@ subroutine build_ATCA(geom)
real(DP) :: sum,maxdiag
logical :: blas=.true.
!
print*,'Allocating GtG ...',real(geom%NBP)**2*4/1024/1024,' Mbyte'
print*,'Allocating GtG ...',real(geom%NBP)**2*4./1024./1024./1024.,' Gbyte'
allocate(geom%ATCA(geom%NBP,geom%NBP),stat=ier)
if (ier/=0) STOP 'memory allocation error in build_ATCA (GtG matrix)'
geom%ATCA(:,:)=0.
......@@ -931,7 +952,7 @@ subroutine generate_invmesh(geom)
do i=1,geom%ring(irg)%ngllh_shape
eta=xgllh_shape(i)*geom%dphi/2._DP
X=tan(eta)
isdelta=1._DP/sqrt(1+X**2+Y**2)
isdelta=1._DP/sqrt(1._DP+X**2+Y**2)
geom%ring(irg)%mesh_shape(1,i,j,k)=r *isdelta
geom%ring(irg)%mesh_shape(2,i,j,k)=r*X*isdelta
geom%ring(irg)%mesh_shape(3,i,j,k)=r*Y*isdelta
......@@ -952,11 +973,11 @@ subroutine generate_invmesh(geom)
if (k==geom%ring(irg)%ngllv_int) r=r-eps !catching interfaces
do j=1,geom%ring(irg)%ngllh_int
nu=xgllh_int(j)*geom%dtheta/2._DP
Y=tan(nu)
Y=dtan(nu)
do i=1,geom%ring(irg)%ngllh_int
eta=xgllh_int(i)*geom%dphi/2._DP
X=tan(eta)
isdelta=1._DP/sqrt(1+X**2+Y**2)
X=dtan(eta)
isdelta=1._DP/dsqrt(1._DP+X**2+Y**2)
xx(1)=r *isdelta
xx(2)=r*X*isdelta
xx(3)=r*Y*isdelta
......@@ -978,7 +999,9 @@ subroutine generate_invmesh(geom)
ijk=ijk+1
call xyz2rtp(geom%ring(irg)%mesh_int(:,i,j,k),geom%mesh_int_rtp(1,ijk),geom%mesh_int_rtp(2,ijk) &
,geom%mesh_int_rtp(3,ijk),deg=.true.)
if (j==1.and.i==1) geom%mesh_rad(k+geom%ring(irg)%NBradic)=geom%mesh_int_rtp(1,ijk)
if (j==1.and.i==1) &
geom%mesh_rad(k+geom%ring(irg)%NBradic)=geom%mesh_int_rtp(1,ijk)
geom%iglob2irad(ijk)=k+geom%ring(irg)%NBradic
! print*,'hooo',ijk,geom%iglob2irad(ijk)
! geom%mesh_int_rtp(1,ijk)=(RA-geom%mesh_int_rtp(1,ijk))/1000.d0 !depth in meter
......@@ -1364,6 +1387,9 @@ subroutine read_src(name,geom,coef)
read(11,*) (coef(i),i=1,geom%NBP)
geom%tmin=max(0.d0,geom%t0-geom%time_support/2.d0)
geom%tmax=geom%t0+geom%time_support/2.d0
geom%NBCM=6
!test
!geom%ring(1)%ngllv_int=10; geom%ring(1)%ngllh_int=10
! geom%dtinv =1.d0/geom%fmax/2.2d0
! print*,'raaa:',geom%tmin,geom%tmax
! geom%NBTinv=int((geom%tmax-geom%tmin)/geom%dtinv)+1
......@@ -1422,7 +1448,7 @@ end subroutine comp_misfit
subroutine sum_green(geom,coef,dir)
!-----------------------------------------------------------
use def_gparam
use global_main_param, only:rotation,stations_name,comp,compr
use global_main_param, only:rotation,rcv,comp,compr
implicit none
type(meshinv), intent(in) :: geom
real(DP), dimension(:), intent(in) :: coef
......@@ -1441,7 +1467,7 @@ subroutine sum_green(geom,coef,dir)
traces(:,ic,ir)=traces(:,ic,ir)+coef(ip)*geom%green(:,ic,ir,ip)
enddo
enddo
enddo
enddo
#ifdef __GFORTRAN__
inquire(file=dir,exist=flag)
#else
......@@ -1454,9 +1480,9 @@ subroutine sum_green(geom,coef,dir)
do ir=1,geom%NBR
do ic=1,3
if (rotation) then
write(name,'("U",a1,"_",a4)') compr(ic),stations_name(ir)
write(name,'("U",a1,"_",a4)') compr(ic),rcv%name(ir)
else
write(name,'("U",a1,"_",a4)') comp(ic) ,stations_name(ir)
write(name,'("U",a1,"_",a4)') comp(ic) ,rcv%name(ir)
endif
open(17,file=add_dir(name,dir))
call fftresample(traces(:,ic,ir),0._DP,geom%dttrace,geom%NBTtrace,tmp,geom%dtout,geom%NBTout)
......@@ -1472,7 +1498,7 @@ subroutine sum_green(geom,coef,dir)
subroutine read_data(geom,data)
!-------------------------------------------------------------
use def_gparam
use global_main_param, only:rotation,stations_name,comp,compr
use global_main_param, only:rotation,rcv,comp,compr
implicit none
type(meshinv), intent(inout) :: geom
real(DP), dimension(:,:,:), pointer, intent(out) :: data
......@@ -1482,7 +1508,7 @@ subroutine sum_green(geom,coef,dir)
integer :: i,ir,ic,iend,NBTdata
character(len=100) :: name
!get DATA dt and NBT:
write(name,'("U",a1,"_",a4)') compr(1),stations_name(1)
write(name,'("U",a1,"_",a4)') compr(1),rcv%name(1)
open(17,file=add_dir(name,geom%datadir),status='old',action='read')
iend=0
NBTdata=0
......@@ -1499,9 +1525,9 @@ subroutine sum_green(geom,coef,dir)
do ir=1,geom%NBR
do ic=1,3
if (rotation) then
write(name,'("U",a1,"_",a4)') compr(ic),stations_name(ir)
write(name,'("U",a1,"_",a4)') compr(ic),rcv%name(ir)
else
write(name,'("U",a1,"_",a4)') comp(ic) ,stations_name(ir)
write(name,'("U",a1,"_",a4)') comp(ic) ,rcv%name(ir)
endif
allocate(tmp(NBTdata))
open(17,file=add_dir(name,geom%datadir),status='old',action='read')
......
nms/srcdirect.f90
View file @ 2c849251
......@@ -15,12 +15,15 @@ program srcdirect
read(*,'(a)') model
print*,'Reading source configuration ...'
call read_src(name,geom,coef)
print*,'fmax, fmaxb:',geom%fmax,geom%fmaxb
print*,'NBT, dt (in source) : ',geom%NBTinv, geom%dtinv
call read_modele(model,geom%modele1d)
print*,'Generating modeling mesh and green functions ...'
call generate_invmesh(geom)
geom%applymask=.true.
geom%timetapper=-1._DP !=> computed automatically
! call sum_green_fcts(geom,coef,'direct')
geom%save_green=.false.
call build_green_fcts(geom)
print*,'Summing'
call sum_green(geom,coef,'direct')
......
nms/time_function.f90
View file @ 2c849251
......@@ -146,12 +146,12 @@
allocate(fct%time(fct%nstep),fct%g(fct%nstep),fct%gp(fct%nstep))
!
select case(fct%type)
case(CRICKER)
do istep = 1,fct%nstep
case(CRICKER)
do istep = 1,fct%nstep
fct%time(istep)=fct%dt * istep
fct%g (istep) = ricker (fct%dt * istep,fct%t0,fct%f0)
fct%gp(istep) = ricker_p(fct%dt * istep,fct%t0,fct%f0)
enddo
enddo
case(CGAUSS)
do istep = 1,fct%nstep
fct%time(istep)=fct%dt * istep
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment