-
Notifications
You must be signed in to change notification settings - Fork 0
/
wavefield_decomp2.f90
126 lines (108 loc) · 4.06 KB
/
wavefield_decomp2.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
!! wavefield decomposition code
!! apply globally
subroutine wavefield_decomp2(p,p_hb,sou_p,sou_p_hb,&
image_du,taper_x,taper_z,&
nx,nz,lv,nxfft,nzfft,flag)
use fftw_module_f
implicit none
include "fftw3.f"
integer::nx,nz,nxfft,nzfft,lv,flag
!! up/down && left/right
real,dimension(-lv:nx+lv,-lv:nz+lv)::p,p_hb,sou_p,sou_p_hb
!!
real,dimension(-lv:nx+lv,-lv:nz+lv)::p_lu,p_ld,p_ru,p_rd
real,dimension(1:nx,1:nz)::image_du
complex,dimension(1:nx,1:nz)::sou_p_hb_t_spa,p_hb_t_spa
! complex,dimension(1:nzfft)::tmp_in_ud,tmp_out_ud,&
! tmp_in_dw,tmp_out_dw,&
! tmp_in_up,tmp_out_up
!
! complex,dimension(1:nxfft)::tmp_in_lr,tmp_out_lr,&
! tmp_in_lf,tmp_out_lf,&
! tmp_in_rt,tmp_out_rt
!
! complex,dimension(1:nxfft,1:nzfft)::tmp_in_bin,tmp_out_bin,&
! tmp_in_lu,tmp_out_lu,tmp_in_ru,tmp_out_ru,&
! tmp_in_ld,tmp_out_ld,tmp_in_rd,tmp_out_rd
real,dimension(1:nxfft)::taper_x
real,dimension(1:nzfft)::taper_z
integer::i,j
real::pi,tp
if(flag.eq.1) then
!! begin do i=1,nx
do i=1,nx
tmp_in_ud(:)=cmplx(0.0,0.0)
tmp_out_ud(:)=cmplx(0.0,0.0)
do j=1,nz
tmp_in_ud(j)=cmplx(p(i,j),p_hb(i,j))
enddo
call sfftw_execute(plan_ud,tmp_in_ud,tmp_out_ud)
do j=2,nzfft/2
tmp_out_ud(j)=2.0*cmplx(real(tmp_out_ud(j)),aimag(tmp_out_ud(j)))
enddo
do j=nzfft/2+2,nzfft
tmp_out_ud(j)=cmplx(0.0,0.0)
enddo
if(.false.) then
do j=nzfft/2+2,nzfft
tmp_out_ud(j)=2.0*cmplx(real(tmp_out_ud(j)),aimag(tmp_out_ud(j)))
enddo
do j=2,nzfft/2
tmp_out_ud(j)=cmplx(0.0,0.0)
enddo
! do j=1,nzfft/2
endif
! tmp_out_ud(j)=cmplx(aimag(tmp_out_ud(j)),-real(tmp_out_ud(j)))
! enddo
! do j=nzfft/2+1,nzfft
! tmp_out_ud(j)=cmplx(-aimag(tmp_out_ud(j)),real(tmp_out_ud(j)))
! enddo
tmp_in_ud=cmplx(0.0,0.0)
call sfftw_execute(plan_ud_inv,tmp_out_ud,tmp_in_ud)
do j=1,nz
p_hb_t_spa(i,j)=tmp_in_ud(j)/real(nzfft)
enddo
!! enddo i=1,nx
enddo
!! begin do i=1,nx
do i=1,nx
tmp_in_ud(:)=cmplx(0.0,0.0)
tmp_out_ud(:)=cmplx(0.0,0.0)
do j=1,nz
tmp_in_ud(j)=cmplx(sou_p(i,j),sou_p_hb(i,j))
enddo
call sfftw_execute(plan_ud,tmp_in_ud,tmp_out_ud)
! do j=1,nzfft/2
! tmp_out_ud(j)=cmplx(aimag(tmp_out_ud(j)),-real(tmp_out_ud(j)))
! enddo
! do j=nzfft/2+1,nzfft
! tmp_out_ud(j)=cmplx(-aimag(tmp_out_ud(j)),real(tmp_out_ud(j)))
! enddo
if(.false.) then
do j=2,nzfft/2
tmp_out_ud(j)=2.0*cmplx(real(tmp_out_ud(j)),aimag(tmp_out_ud(j)))
enddo
do j=nzfft/2+2,nzfft
tmp_out_ud(j)=cmplx(0.0,0.0)
enddo
endif
do j=2,nzfft/2
tmp_out_ud(j)=cmplx(0.0,0.0)
enddo
do j=nzfft/2+2,nzfft
tmp_out_ud(j)=2.0*cmplx(real(tmp_out_ud(j)),aimag(tmp_out_ud(j)))
enddo
tmp_in_ud=cmplx(0.0,0.0)
call sfftw_execute(plan_ud_inv,tmp_out_ud,tmp_in_ud)
do j=1,nz
sou_p_hb_t_spa(i,j)=tmp_in_ud(j)/real(nzfft)
enddo
enddo
!! enddo i=1,nx
endif
do j=1,nz
do i=1,nx
image_du(i,j)=image_du(i,j)+real(sou_p_hb_t_spa(i,j)*p_hb_t_spa(i,j))
enddo
enddo
end subroutine wavefield_decomp2