module OW_wei_areal use sep use OWvelo_areal use OW_type use OW_image_areal use OW_fftw_areal use OW_parms_areal implicit none contains subroutine owwei_f(taper,recv,source,sref,vel,deltaS,deltaI,GGr,GGs) complex,optional :: GGr(:,:,:),GGs(:,:,:) complex :: recv(:,:),source(:,:) integer :: iwa,iz,ikx,iref,ixw,j,inode,i,ish,nrefvel,nodes real :: sqa,taper(:),vel(:,:),sref(:,:,:),a real :: c(5),deltaS(:,:),deltaI(:,:,:) integer,allocatable :: indw(:,:) real, allocatable :: GGG(:,:,:,:),sqa0(:,:) complex, allocatable :: ws(:,:),wr(:,:),gs0(:,:),gr0(:,:) complex, allocatable :: ddws(:,:),ddwr(:,:),gs(:,:),gr(:,:) complex, allocatable :: dws(:,:),dwr(:,:),dgs0(:,:),dgr0(:,:) complex, allocatable :: dwss(:,:),dwrr(:,:),dgs(:,:),dgr(:,:) complex,allocatable :: tflds(:,:),tfldr(:,:) integer, external :: omp_get_thread_num,omp_get_num_threads if (node.eq.9999) then !$OMP PARALLEL nodes = omp_get_num_threads() !$OMP END PARALLEL else nodes=node end if write(0,*) "****************" write(0,*) "* NUM.THREADS = ",nodes write(0,*) "****************" allocate(wr(kxm%n,nodes)); wr=0 allocate(ws(kxm%n,nodes)); ws=0 allocate(ddwr(kxm%n,nodes)); ddwr=0 allocate(ddws(kxm%n,nodes)); ddws=0 allocate(gs0(kxm%n,nodes)); gs0=0 allocate(gr0(kxm%n,nodes)); gr0=0 allocate(dwr(kxm%n,nodes)); dwr=0 allocate(dws(kxm%n,nodes)); dws=0 allocate(dwrr(kxm%n,nodes)); dwrr=0 allocate(dwss(kxm%n,nodes)); dwss=0 allocate(dgs0(kxm%n,nodes)); dgs0=0 allocate(dgr0(kxm%n,nodes)); dgr0=0 allocate(sqa0(kxm%n,nodes)); sqa0=0 allocate(GGG(img%z%n,img%sh%n,kxm%n,nodes)); GGG=0 c=(/1.,0.5,0.375,0.3125,0.2734375/) allocate(indw(rec%w%n,2));indw=0 !-------------------------------------------------- ! Setting FFTW plans allocate(tflds(kxm%n,nodes));tflds=0 call owfftws(tflds,nodes) allocate(tfldr(kxm%n,nodes));tfldr=0 call owfftwr(tfldr,nodes) !$OMP parallel do private(inode,& iwa,& iz,& ikx,& iref,& ixw,& j,& nrefvel,& gr,& gs,& dgr,& dgs,& sqa) ordered schedule(dynamic,1) num_threads(nodes) do iwa = 1,rec%w%n inode=omp_get_thread_num()+1; ws(padtraces+1:sou%h%n+padtraces,inode)=source(:,iwa) wr(padtraces+1:rec%h%n+padtraces,inode)=recv(:,iwa) dws(:,inode)=0 dwr(:,inode)=0 do iz=nwbott,img%z%n !-------------------------------------------------- ! Compute vref ! nrefvel=count(mask = sref(:,1,iz).ne.-1) if (nrefvel.gt.1) then if (allocated(gs)) deallocate(gs); allocate(gs(kxm%n,nrefvel)); if (allocated(gr)) deallocate(gr); allocate(gr(kxm%n,nrefvel)); if (allocated(dgs)) deallocate(dgs); allocate(dgs(kxm%n,nrefvel)); if (allocated(dgr)) deallocate(dgr); allocate(dgr(kxm%n,nrefvel)); gs=0 gr=0 dgs=0 dgr=0 else gs0(:,inode)=0 gr0(:,inode)=0 dgs0(:,inode)=0 dgr0(:,inode)=0 end if tflds(:,inode)=ws(:,inode) call owfftexes(tflds(:,inode),inode) ws(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n) tfldr(:,inode)=wr(:,inode) call owfftexer(tfldr(:,inode),inode) wr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n) ddwr(:,inode)=0;ddws(:,inode)=0 do j=0,4 dwss(:,inode)=kxm%all**(2*j)*ws(:,inode) dwrr(:,inode)=-kxm%all**(2*j)*wr(:,inode) tflds(:,inode)=dwss(:,inode) call owifftexes(tflds(:,inode),inode) ! dwss(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)*(vel(iz-1,:)/rec%w%all(iwa))**(2*j)*c(j+1) dwss(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)*(vel(iz,:)/rec%w%all(iwa))**(2*j)*c(j+1) tfldr(:,inode)=dwrr(:,inode) call owifftexer(tfldr(:,inode),inode) ! dwrr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n)*(vel(iz-1,:)/rec%w%all(iwa))**(2*j)*c(j+1) dwrr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n)*(vel(iz,:)/rec%w%all(iwa))**(2*j)*c(j+1) ddws(:,inode)=ddws(:,inode)+dwss(:,inode) ddwr(:,inode)=ddwr(:,inode)+dwrr(:,inode) end do ! dws(:,inode)=dws(:,inode)+ddws(:,inode)*img%z%d*rec%w%all(iwa)*exp(cmplx(0,pi/2))*deltaS(iz-1,:) ! dwr(:,inode)=dwr(:,inode)+ddwr(:,inode)*img%z%d*rec%w%all(iwa)*exp(cmplx(0,pi/2))*deltaS(iz-1,:) dws(:,inode)=dws(:,inode)+ddws(:,inode)*img%z%d*rec%w%all(iwa)*exp(cmplx(0,pi/2))*deltaS(iz,:) dwr(:,inode)=dwr(:,inode)+ddwr(:,inode)*img%z%d*rec%w%all(iwa)*exp(cmplx(0,pi/2))*deltaS(iz,:) tflds(:,inode)=dws(:,inode) call owfftexes(tflds(:,inode),inode) dws(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n) tfldr(:,inode)=dwr(:,inode) call owifftexer(tfldr(:,inode),inode) dwr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n) if (nrefvel .eq. 1) then gs0(:,inode)=ws(:,inode); gr0(:,inode)=wr(:,inode); dgs0(:,inode)=dws(:,inode); dgr0(:,inode)=dwr(:,inode); do ikx = 1,kxm%n sqa = (rec%w%all(iwa) * sref(1,1,iz))**2 - kxm%all(ikx)**2 if (sqa > 0) then gs0(ikx,inode)=gs0(ikx,inode)*exp(cmplx(0,v%z%d*sqrt(sqa))) gr0(ikx,inode)=gr0(ikx,inode)*exp(cmplx(0,-v%z%d*sqrt(sqa))) dgs0(ikx,inode)=dgs0(ikx,inode)*exp(cmplx(0,v%z%d*sqrt(sqa))) dgr0(ikx,inode)=dgr0(ikx,inode)*exp(cmplx(0,-v%z%d*sqrt(sqa))) else gs0(ikx,inode)=gs0(ikx,inode)*exp(-sqrt(-sqa)*v%z%d) gr0(ikx,inode)=gr0(ikx,inode)*exp(-sqrt(-sqa)*v%z%d) dgs0(ikx,inode)=dgs0(ikx,inode)*exp(-sqrt(-sqa)*v%z%d) dgr0(ikx,inode)=dgr0(ikx,inode)*exp(-sqrt(-sqa)*v%z%d) end if end do tflds(:,inode)=gs0(:,inode) call owifftexes(tflds(:,inode),inode) ws(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n) tfldr(:,inode)=gr0(:,inode) call owifftexer(tfldr(:,inode),inode) wr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n) tflds(:,inode)=dgs0(:,inode) call owifftexes(tflds(:,inode),inode) dws(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n) tfldr(:,inode)=dgr0(:,inode) call owifftexer(tfldr(:,inode),inode) dwr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n) else do iref=1,nrefvel gs(:,iref)=ws(:,inode) gr(:,iref)=wr(:,inode) dgs(:,iref)=dws(:,inode) dgr(:,iref)=dwr(:,inode) do ikx = 1,kxm%n sqa = (rec%w%all(iwa) * sref(iref,1,iz))**2 - kxm%all(ikx)**2 if (sqa > 0) then gs(ikx,iref)=gs(ikx,iref)*exp(cmplx(0,v%z%d*sqrt(sqa))); gr(ikx,iref)=gr(ikx,iref)*exp(cmplx(0,-v%z%d*sqrt(sqa))); dgs(ikx,iref)=dgs(ikx,iref)*exp(cmplx(0,v%z%d*sqrt(sqa))); dgr(ikx,iref)=dgr(ikx,iref)*exp(cmplx(0,-v%z%d*sqrt(sqa))); else gs(ikx,iref)=gs(ikx,iref)*exp(-sqrt(-sqa)*v%z%d) gr(ikx,iref)=gr(ikx,iref)*exp(-sqrt(-sqa)*v%z%d) dgs(ikx,iref)=dgs(ikx,iref)*exp(-sqrt(-sqa)*v%z%d) dgr(ikx,iref)=dgr(ikx,iref)*exp(-sqrt(-sqa)*v%z%d) end if end do tflds(:,inode)=gs(:,iref) call owifftexes(tflds(:,inode),inode) gs(:,iref)=tflds(:,inode)/sqrt(1.*kxm%n) tfldr(:,inode)=gr(:,iref) call owifftexer(tfldr(:,inode),inode) gr(:,iref)=tfldr(:,inode)/sqrt(1.*kxm%n) tflds(:,inode)=dgs(:,iref) call owifftexes(tflds(:,inode),inode) dgs(:,iref)=tflds(:,inode)/sqrt(1.*kxm%n) tfldr(:,inode)=dgr(:,iref) call owifftexer(tfldr(:,inode),inode) dgr(:,iref)=tfldr(:,inode)/sqrt(1.*kxm%n) gs(:,iref) = gs(:,iref) * exp(cmplx(0,v%z%d*rec%w%all(iwa))*(1./vel(iz,:)-sref(iref,1,iz))); gr(:,iref) = gr(:,iref) * exp(cmplx(0,-v%z%d*rec%w%all(iwa))*(1./vel(iz,:)-sref(iref,1,iz))); dgs(:,iref) = dgs(:,iref) * exp(cmplx(0,v%z%d*rec%w%all(iwa))*(1./vel(iz,:)-sref(iref,1,iz))); dgr(:,iref) = dgr(:,iref) * exp(cmplx(0,-v%z%d*rec%w%all(iwa))*(1./vel(iz,:)-sref(iref,1,iz))); end do do iref=1,nrefvel-1 where (vel(iz,:) >= 1/sref(iref,1,iz) .and. vel(iz,:) < 1/sref(iref+1,1,iz)) ws(:,inode)=gs(:,iref) wr(:,inode)=gr(:,iref) dws(:,inode)=dgs(:,iref) dwr(:,inode)=dgr(:,iref) end where end do where (vel(iz,:) >= 1/sref(nrefvel,1,iz)) ws(:,inode)=gs(:,nrefvel) wr(:,inode)=gr(:,nrefvel) dws(:,inode)=dgs(:,nrefvel) dwr(:,inode)=dgr(:,nrefvel) end where end if ws(:,inode)=ws(:,inode)*taper wr(:,inode)=wr(:,inode)*taper dws(:,inode)=dws(:,inode)*taper dwr(:,inode)=dwr(:,inode)*taper if (exist_file("ws")) then !$OMP ORDERED indw(iwa,2)=inode GGR(iz,:,iwa)=wr(:,indw(iwa,2)) GGS(iz,:,iwa)=ws(:,indw(iwa,2)) !$OMP END ORDERED end if call owimage_f(iwa,ws(:,inode),wr(:,inode),dws(:,inode),dwr(:,inode),GGG(iz,:,:,inode)) !========================================================================|| !========================================================================|| ! call owimage0(iwa,ws,wr,GG(iz,:,:)) ! if (nwbott.ne.1) then ! call owvref(1,minval(vel(nwbott-1,:)),vel(1,:),vref,minv,maxv) ! tfld=0.; tfld=wr; call sfftw_execute(planF); wr=tfld/sqrt(1.*kxm%n) ! tfld=0.; tfld=ws; call sfftw_execute(planF); ws=tfld/sqrt(1.*kxm%n) ! ! Operations in the wavenumber domain ! ! phase shift ! do ikx = 1,kxm%n ! sqa = rec%w%all(iwa)**2 / (minv)**2 - kxm%all(ikx)**2 ! if (sqa > 0) then ! wr(ikx)=wr(ikx)*exp(cmplx(0,dwbott*sqrt(sqa))) ! ws(ikx)=ws(ikx)*exp(cmplx(0,-dwbott*sqrt(sqa))) ! else ! wr(ikx)=wr(ikx)*exp(-sqrt(-sqa)*dwbott) ! ws(ikx)=ws(ikx)*exp(-sqrt(-sqa)*dwbott) ! end if ! end do ! tfld=wr; call sfftw_execute(planI); wr=tfld/sqrt(1.*kxm%n) ! tfld=ws; call sfftw_execute(planI); ws=tfld/sqrt(1.*kxm%n) ! wr=wr*taper ! ws=ws*taper ! end if !======================================================================== !======================================================================== end do if (mod(iwa,20)==1) write(0,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",rec%x%n, & " FROM NODE: ",inode end do !$OMP END PARALLEL DO call owfftdestroy() deallocate(tfldr,tflds) a=0 do i=1,nodes !$OMP parallel do private(ixw,iz,ish) do iz=1,img%z%n do ish=1,img%sh%n do ixw=1,kxm%n a=a+abs(GGG(iz,ish,ixw,i)) end do end do end do !$OMP END PARALLEL DO end do write(0,*) "****** ",a," ******" do i=1,nodes !$OMP parallel do private(ixw,iz,ish) do iz=1,img%z%n do ish=1,img%sh%n do ixw=1,kxm%n deltaI(iz,ish,ixw)=deltaI(iz,ish,ixw)+GGG(iz,ish,ixw,i) end do end do end do !$OMP END PARALLEL DO end do 6 format(A13,i3,A8,i3,A14,i3,A12,i3) deallocate(wr,ws,ddwr,ddws,gs0,gr0) deallocate(dwr,dws,dwrr,dwss) deallocate(dgs0,dgr0,sqa0,GGG) end subroutine !========================================================================================= subroutine owwei_b(taper,recv,source,sref,vel,deltaS,deltaI) complex :: recv(:,:,:),source(:,:,:) real :: deltaI(:,:,:) integer :: iwa,iz,ikx,iref,ixw,j,inode,i,ish,ihso,nrefvel,nodes real :: sqa,a real :: taper(:),vel(:,:),sref(:,:,:),c(5),deltaS(:,:) integer,allocatable :: indw(:,:) real, allocatable :: GGG(:,:,:),sqa0(:,:) complex, allocatable :: ws(:,:),wr(:,:),gs0(:,:),gr0(:,:) complex, allocatable :: ddws(:,:),ddwr(:,:),gs(:,:),gr(:,:) complex, allocatable :: dws(:,:),dwr(:,:),dgs0(:,:),dgr0(:,:) complex, allocatable :: dwss(:,:),dwrr(:,:),dgs(:,:),dgr(:,:) complex,allocatable :: awr(:,:),aws(:,:),adwr(:,:),adws(:,:) complex,allocatable :: tflds(:,:),tfldr(:,:) integer, external :: omp_get_thread_num,omp_get_num_threads if (node.eq.9999) then !$OMP PARALLEL nodes = omp_get_num_threads() !$OMP END PARALLEL else nodes=node end if write(0,*) "****************" write(0,*) "* NUM.THREADS = ",nodes write(0,*) "****************" allocate(wr(kxm%n,nodes)); wr=0 allocate(ws(kxm%n,nodes)); ws=0 allocate(ddwr(kxm%n,nodes)); ddwr=0 allocate(ddws(kxm%n,nodes)); ddws=0 allocate(gs0(kxm%n,nodes)); gs0=0 allocate(gr0(kxm%n,nodes)); gr0=0 allocate(dwr(kxm%n,nodes)); dwr=0 allocate(dws(kxm%n,nodes)); dws=0 allocate(adwr(kxm%n,nodes)); adwr=0 allocate(adws(kxm%n,nodes)); adws=0 allocate(awr(kxm%n,nodes)); awr=0 allocate(aws(kxm%n,nodes)); aws=0 allocate(dwrr(kxm%n,nodes)); dwrr=0 allocate(dwss(kxm%n,nodes)); dwss=0 allocate(dgs0(kxm%n,nodes)); dgs0=0 allocate(dgr0(kxm%n,nodes)); dgr0=0 allocate(sqa0(kxm%n,nodes)); sqa0=0 allocate(GGG(img%z%n,kxm%n,nodes)); GGG=0 c=(/1.,0.5,0.375,0.3125,0.2734375/) !-------------------------------------------------- ! Setting FFTW plans allocate(tflds(kxm%n,nodes));tflds=0 call owfftws(tflds,nodes) allocate(tfldr(kxm%n,nodes));tfldr=0 call owfftwr(tfldr,nodes) !$OMP parallel do private(inode,& iwa,& iz,& ikx,& iref,& ixw,& j,& nrefvel,& gr,& gs,& dgr,& dgs,& sqa) ordered schedule(dynamic,1) num_threads(nodes) ! ! Freq ! do iwa = 1,rec%w%n source=conjg(source) inode=omp_get_thread_num()+1 dws(:,inode)=0. dwr(:,inode)=0. aws(:,inode)=0. awr(:,inode)=0. ddwr(:,inode)=0 ddws(:,inode)=0 dwrr(:,inode)=0 dwss(:,inode)=0 tflds(:,inode)=source(img%z%n-1,:,iwa) call owfftexes(tflds(:,inode),inode) aws(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n) tfldr(:,inode)=recv(img%z%n-1,:,iwa) call owfftexer(tfldr(:,inode),inode) awr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n) do j=0,4 dwss(:,inode)=kxm%all**(2*j)*aws(:,inode) dwrr(:,inode)=kxm%all**(2*j)*awr(:,inode) tflds(:,inode)=dwss(:,inode) call owifftexes(tflds(:,inode),inode) dwss(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)*(vel(img%z%n-1,:)/rec%w%all(iwa))**(2*j)*c(j+1) tfldr(:,inode)=dwrr(:,inode) call owifftexer(tfldr(:,inode),inode) dwrr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n)*(vel(img%z%n-1,:)/rec%w%all(iwa))**(2*j)*c(j+1) ddws(:,inode)=ddws(:,inode)+dwss(:,inode) ddwr(:,inode)=ddwr(:,inode)+dwrr(:,inode) end do ddws(:,inode)=ddws(:,inode)*img%z%d*rec%w%all(iwa)*exp(cmplx(0,pi/2)) ddwr(:,inode)=ddwr(:,inode)*img%z%d*rec%w%all(iwa)*exp(cmplx(0,pi/2)) do ixw=1,kxm%n do ish=1,img%sh%n ihso=nint(((ish-1)*img%sh%d+img%sh%o)/(img%xm%d)) ! if (ixw-ihso > 1 .and. ixw-ihso < kxm%n .and. ixw+ihso > 1 .and. ixw+ihso < kxm%n) then ! dwr(ixw,inode)=dwr(ixw,inode)+source(img%z%n,ixw-ihso,iwa)*deltaI(img%z%n,ish,ixw) ! dws(ixw,inode)=dws(ixw,inode)+recv(img%z%n,ixw+ihso,iwa)*deltaI(img%z%n,ish,ixw) if (ixw-ish > 1 .and. ixw-ish < kxm%n .and. ixw+ish > 1 .and. ixw+ish < kxm%n) then dwr(ixw,inode)=dwr(ixw,inode)+source(img%z%n,ixw-ish,iwa)*deltaI(img%z%n,ish,ixw) dws(ixw,inode)=dws(ixw,inode)+(recv(img%z%n,ixw+ish,iwa))*deltaI(img%z%n,ish,ixw) end if end do end do ! GGG(img%z%n-1,:,inode)=GGG(img%z%n-1,:,inode)+real(conjg(ddws(:,inode))*dws(:,inode)+conjg(ddwr(:,inode))*dwr(:,inode)) GGG(img%z%n-1,:,inode)=GGG(img%z%n-1,:,inode)+real((ddws(:,inode))*dws(:,inode)+(ddwr(:,inode))*dwr(:,inode)) ! call owimage_b(iwa,dws(:,inode),dwr(:,inode),ddws(:,inode),ddwr(:,inode),source(img%z%n,:,iwa),recv(img%z%n,:,iwa), & ! deltaI(img%z%n,:,:),GGG(img%z%n-1,:,inode)) ! ! Depth ! do iz=img%z%n-1,nwbott+1,-1 !-------------------------------------------------- ! Compute vref ! nrefvel=count(mask = sref(:,1,iz).ne.-1) if (nrefvel.gt.1) then if (allocated(dgs)) deallocate(dgs); allocate(dgs(kxm%n,nrefvel)); if (allocated(dgr)) deallocate(dgr); allocate(dgr(kxm%n,nrefvel)); dgs=0 dgr=0 else dgs0(:,inode)=0 dgr0(:,inode)=0 end if tflds(:,inode)=dws(:,inode) call owfftexes(tflds(:,inode),inode) dws(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n) tfldr(:,inode)=dwr(:,inode) call owfftexer(tfldr(:,inode),inode) dwr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n) if (nrefvel .eq. 1) then dgs0(:,inode)=dws(:,inode) dgr0(:,inode)=dwr(:,inode) do ikx = 1,kxm%n sqa = (rec%w%all(iwa) * sref(1,1,iz))**2 - kxm%all(ikx)**2 if (sqa > 0) then dgs0(ikx,inode)=dgs0(ikx,inode)*exp(cmplx(0,-v%z%d*sqrt(sqa))) dgr0(ikx,inode)=dgr0(ikx,inode)*exp(cmplx(0,v%z%d*sqrt(sqa))) else dgs0(ikx,inode)=dgs0(ikx,inode)*exp(-sqrt(-sqa)*v%z%d) dgr0(ikx,inode)=dgr0(ikx,inode)*exp(-sqrt(-sqa)*v%z%d) end if end do tflds(:,inode)=dgs0(:,inode) call owifftexes(tflds(:,inode),inode) dws(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n) tfldr(:,inode)=dgr0(:,inode) call owifftexer(tfldr(:,inode),inode) dwr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n) else do iref=1,nrefvel dgs(:,iref)=dws(:,inode) dgr(:,iref)=dwr(:,inode) do ikx = 1,kxm%n sqa = (rec%w%all(iwa) * sref(iref,1,iz))**2 - kxm%all(ikx)**2; if (sqa > 0) then dgs(ikx,iref)=dgs(ikx,iref)*exp(cmplx(0,-v%z%d*sqrt(sqa))); dgr(ikx,iref)=dgr(ikx,iref)*exp(cmplx(0,v%z%d*sqrt(sqa))); else dgs(ikx,iref)=dgs(ikx,iref)*exp(-sqrt(-sqa)*v%z%d) dgr(ikx,iref)=dgr(ikx,iref)*exp(-sqrt(-sqa)*v%z%d) end if end do tflds(:,inode)=dgs(:,iref) call owifftexes(tflds(:,inode),inode) dgs(:,iref)=tflds(:,inode)/sqrt(1.*kxm%n) tfldr(:,inode)=dgr(:,iref) call owifftexer(tfldr(:,inode),inode) dgr(:,iref)=tfldr(:,inode)/sqrt(1.*kxm%n) dgs(:,iref) = dgs(:,iref) * exp(cmplx(0,-v%z%d*rec%w%all(iwa))*(1./vel(iz,:)-sref(iref,1,iz))) dgr(:,iref) = dgr(:,iref) * exp(cmplx(0,v%z%d*rec%w%all(iwa))*(1./vel(iz,:)-sref(iref,1,iz))) end do do iref=1,nrefvel-1 where (vel(iz,:) >= 1/sref(iref,1,iz) .and. vel(iz,:) < 1/sref(iref+1,1,iz)) dws(:,inode)=dgs(:,iref) dwr(:,inode)=dgr(:,iref) end where end do where (vel(iz,:) >= 1/sref(nrefvel,1,iz)) dws(:,inode)=dgs(:,nrefvel) dwr(:,inode)=dgr(:,nrefvel) end where end if dws(:,inode)=dws(:,inode)*taper dwr(:,inode)=dwr(:,inode)*taper ddwr(:,inode)=0 ddws(:,inode)=0 aws(:,inode)=0. awr(:,inode)=0. tflds(:,inode)=source(iz-1,:,iwa) call owfftexes(tflds(:,inode),inode) aws(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n) tfldr(:,inode)=recv(iz-1,:,iwa) call owfftexer(tfldr(:,inode),inode) awr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n) do j=0,4 dwss(:,inode)=kxm%all**(2*j)*aws(:,inode) dwrr(:,inode)=kxm%all**(2*j)*awr(:,inode) tflds(:,inode)=dwss(:,inode) call owifftexes(tflds(:,inode),inode) dwss(:,inode)=tflds(:,inode)/sqrt(1.*kxm%n)*(vel(iz-1,:)/rec%w%all(iwa))**(2*j)*c(j+1) tfldr(:,inode)=dwrr(:,inode) call owifftexer(tfldr(:,inode),inode) dwrr(:,inode)=tfldr(:,inode)/sqrt(1.*kxm%n)*(vel(iz-1,:)/rec%w%all(iwa))**(2*j)*c(j+1) ddws(:,inode)=ddws(:,inode)+dwss(:,inode) ddwr(:,inode)=ddwr(:,inode)+dwrr(:,inode) end do ddws(:,inode)=ddws(:,inode)*img%z%d*rec%w%all(iwa)*exp(cmplx(0,pi/2)) ddwr(:,inode)=ddwr(:,inode)*img%z%d*rec%w%all(iwa)*exp(cmplx(0,pi/2)) aws(:,inode)=0. awr(:,inode)=0. do ixw=1,kxm%n do ish=1,img%sh%n ihso=nint(((ish-1)*img%sh%d+img%sh%o)/(img%xm%d)) ! if (ixw-ihso > 1 .and. ixw-ihso < kxm%n .and. ixw+ihso > 1 .and. ixw+ihso < kxm%n) then ! awr(ixw,inode)=aws(ixw,inode)+source(iz,ixw-ihso,iwa)*deltaI(iz,ish,ixw) ! aws(ixw,inode)=awr(ixw,inode)+recv(iz,ixw+ihso,iwa)*deltaI(iz,ish,ixw) if (ixw-ish > 1 .and. ixw-ish < kxm%n .and. ixw+ish > 1 .and. ixw+ish < kxm%n) then awr(ixw,inode)=aws(ixw,inode)+source(iz,ixw-ish,iwa)*deltaI(iz,ish,ixw) aws(ixw,inode)=awr(ixw,inode)+(recv(iz,ixw+ish,iwa))*deltaI(iz,ish,ixw) end if end do dws(:,inode)=dws(:,inode)+aws(:,inode) dwr(:,inode)=dwr(:,inode)+awr(:,inode) end do ! GGG(iz-1,:,inode)=GGG(iz-1,:,inode)+real(conjg(ddws(:,inode))*dws(:,inode)+conjg(ddwr(:,inode))*dwr(:,inode)) GGG(iz-1,:,inode)=GGG(iz-1,:,inode)+real((ddws(:,inode))*dws(:,inode)+(ddwr(:,inode))*dwr(:,inode)) ! call owimage_b(iwa,dws(:,inode),dwr(:,inode),ddws(:,inode),ddwr(:,inode),source(img%z%n,:,iwa),recv(img%z%n,:,iwa), & ! deltaI(img%z%n,ish,ixw),GGG(img%z%n-1,:,inode)) !========================================================================|| !========================================================================|| ! call owimage0(iwa,ws,wr,GG(iz,:,:)) ! if (nwbott.ne.1) then ! call owvref(1,minval(vel(nwbott-1,:)),vel(1,:),vref,minv,maxv) ! tfld=0.; tfld=wr; call sfftw_execute(planF); wr=tfld/sqrt(1.*kxm%n) ! tfld=0.; tfld=ws; call sfftw_execute(planF); ws=tfld/sqrt(1.*kxm%n) ! ! Operations in the wavenumber domain ! ! phase shift ! do ikx = 1,kxm%n ! sqa = rec%w%all(iwa)**2 / (minv)**2 - kxm%all(ikx)**2 ! if (sqa > 0) then ! wr(ikx)=wr(ikx)*exp(cmplx(0,dwbott*sqrt(sqa))) ! ws(ikx)=ws(ikx)*exp(cmplx(0,-dwbott*sqrt(sqa))) ! else ! wr(ikx)=wr(ikx)*exp(-sqrt(-sqa)*dwbott) ! ws(ikx)=ws(ikx)*exp(-sqrt(-sqa)*dwbott) ! end if ! end do ! tfld=wr; call sfftw_execute(planI); wr=tfld/sqrt(1.*kxm%n) ! tfld=ws; call sfftw_execute(planI); ws=tfld/sqrt(1.*kxm%n) ! wr=wr*taper ! ws=ws*taper ! end if !======================================================================== !======================================================================== end do if (mod(iwa,20)==1) write(0,6) " FREQUENCY # ",iwa," OUT OF ",rec%w%n," / SYNTH.FCT. # ",rec%x%n, & " FROM NODE: ",inode 6 format(A13,i3,A8,i3,A14,i3,A12,i3) end do !$OMP END PARALLEL DO call owfftdestroy() deallocate(tfldr,tflds) a=0 do i=1,nodes !$OMP parallel do private(ixw,iz,ish) do iz=1,img%z%n do ixw=1,kxm%n a=a+abs(GGG(iz,ixw,i)) end do end do !$OMP END PARALLEL DO end do write(0,*) "****** ",a," ******" do i=1,nodes !$OMP parallel do private(ixw,iz,ish) do iz=1,img%z%n do ixw=1,kxm%n deltaS(iz,ixw)=deltaS(iz,ixw)+GGG(iz,ixw,i) end do end do !$OMP END PARALLEL DO end do deallocate(wr,ws,ddwr,ddws,gs0,gr0) deallocate(dwr,dws,adwr,adws,dwrr,dwss) deallocate(dgs0,dgr0,sqa0,GGG) end subroutine end module