module_modes.f90 51.3 KB
Newer Older
capdevil's avatar
capdevil committed
1
2
3
4
5
!--------------------------------------------------------------------------------
module modes
!--------------------------------------------------------------------------------
  use def_gparam
  implicit none
6
  public :: init_modes,get_and_sum,get_and_sum_G
capdevil's avatar
capdevil committed
7
8
9
10
11
12
  private
  real(SP), dimension(:,:,:), allocatable :: source_g
  complex(DP), dimension(:,:,:), allocatable ::source_gw
  real(DP) :: fmax,fmaxt,t0_delay,wmax1,wmax2,duree,deltaf,dtf,dtf2
  integer :: logunit,coef_nbt,NBF,LmaxS,LmaxT,NBF2
  integer :: Lmax_all,Nmax_allS,Nmax_allT,Nmax_all
capdevil's avatar
   
capdevil committed
13
  real(DP), parameter :: coef_damp1=1.03_DP, coef_damp2=1.15_DP
capdevil's avatar
capdevil committed
14
15
16
17
18
19
20
21
22
23
24
25
26
27
  real(DP), dimension(:,:), allocatable  :: src_coord,Msrc
  real(DP), dimension(:,:), allocatable  :: rec_coord
  real(DP), dimension(:), allocatable  :: weight
  real(DP), dimension(1) :: rho,Vs,Vp,Qshear,radtmp
  integer, dimension(:)  ,allocatable :: LLLmax 
  integer, dimension(:,:),allocatable :: NNNmax 
  integer, dimension(:), allocatable :: NmaxS,NmaxT
!freq propre et attenuation
  real(DP), dimension(:,:,:), allocatable :: eigw,qw
!
  real(DP), dimension(  :,:), allocatable :: UU,UUp,VV,VVp,WW,WWp
  real(DP), dimension(:,:,:), allocatable :: UUs,UUps,VVs,VVps,WWs,WWps
  real(DP), dimension(:) , allocatable :: al,b0,b1,b2
!
capdevil's avatar
   
capdevil committed
28
29
  integer :: length_qnl,length_qnl_all,nbchunk
  integer, parameter :: memorymax=512
30
!  integer, parameter :: memorymax=128000000
capdevil's avatar
   
capdevil committed
31
32
33
  integer, dimension(:), allocatable :: qnldeb_chunk,qnlfin_chunk
  integer :: qnldeb,qnlfin,iqnld
!
capdevil's avatar
capdevil committed
34
35
36
37
38
39
40
41
  complex*16, dimension(:,:), allocatable :: time_serie
  integer, dimension(2), parameter :: ldeb=(/0,1/)
!
  complex(DP), parameter :: czero=(0._DP,0._DP),II=(0._DP,1._DP)
  complex(SP), parameter :: czero_sp=(0._SP,0._SP)
!
  complex(DP), dimension(:,:,:,:,:  ), allocatable :: SkN
  complex(DP), dimension(:,:,:,:,:  ), allocatable :: RkN
42
  complex(DP), dimension(:,:,:,:,:  ), allocatable :: GkN
capdevil's avatar
capdevil committed
43
!
capdevil's avatar
   
capdevil committed
44
45
46
  real(DP), parameter :: rhobar=5515._DP
  real(DP) :: scale2,scale1
  real(DP) :: scale,dscale,receiver_rad
capdevil's avatar
capdevil committed
47
48
49
50
51
52
53
54
55
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
contains
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
!-------------------------------------------------
  subroutine init_modes(logunit_)
!-------------------------------------------------
    use def_gparam
    use global_main_param, only:source_type,NBT,dt,f0,f1,f2,f3,f4,t0,get_source_dat &
                               ,get_receivers,eigenfileS,eigenfileT,RA  &
capdevil's avatar
   
capdevil committed
56
                               ,NBE,NBR,nbcomp,nmax_in, homogenization
capdevil's avatar
capdevil committed
57
58
    use time_function 
    use util_fctp
capdevil's avatar
   
capdevil committed
59
    use corr_modes
capdevil's avatar
capdevil committed
60
61
62
63
64
65
66
67
68
69
70
71
!
    implicit none
    integer, intent(in) :: logunit_
!
    real(DP), dimension(:,:), allocatable :: xdum,tmp
    real(DP), dimension(:)  , allocatable :: rad
    integer :: NBFtmp,l,n,ir,ic,unit,i
    real(DP):: t0delay
    character(len=30) :: name
    logical :: flag
    real :: bid
!
capdevil's avatar
   
capdevil committed
72
    t0delay=0.d0
capdevil's avatar
capdevil committed
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
    unit=121
    logunit=logunit_
    allocate(source_g(NBT,1,1))
    t0_delay=0.0d0
    call def_timefunc(source_type,NBT,dt,f0,f1,f2,f3,f4,t0,source_g(:,1,1),fmax)  
    fmaxt=coef_damp1*fmax
    fmax =coef_damp2*fmax
    wmax1=fmaxt*2._DP*PI
    wmax2=fmax *2._DP*PI
    duree=NBT*dt+t0delay
!pas en frequence pour la FFT
    deltaf=1./duree
!nombre de frequence  , on va a 1.1 fois Nyquist
!mais aussi nombre de pas de temps effectif
    NBFtmp   = int(2*fmax*((NBT-1)*dt+t0delay))+1
    call get_power(int(NBFtmp*1.1),NBF)
    NBF2=2*NBF
!on utilise deux fois plus long pour la convolution (zero padding sur la mointie
!du signal)
!on recupere le plus grand coef_nbt tel que NBF*coef_nbt < NBT
    coef_nbt=1
    do while (coef_nbt*2*NBF < NBT)
       coef_nbt=2*coef_nbt
    enddo
    dtf   = ((NBT-1)*dt)/real(NBF-1)
!pas de temps apres coef_nbt
    dtf2  = dtf/coef_nbt
!
!ouverture de fichiers de fonctions propres
    print*,'Opening eigenfunction files ...'
    call open_fctp(eigenfileS,eigenfileT,rad_recepteur=receiver_rad)    
capdevil's avatar
   
capdevil committed
104
105
106
107
!RA is now initialized
    scale2=1.D0/(rhobar*RA*RA*RA)
    scale1=scale2/RA/RA
!
capdevil's avatar
capdevil committed
108
109
110
111
112
113
114
115
!lecture des sources:
    print*,'Getting sources and receivers informations ...'
    allocate(src_coord(3,NBE),Msrc(6,NBE),rec_coord(2,NBR),weight(NBE))
    call get_source_dat(src_coord,Msrc)
    call get_receivers(rec_coord)
    write(logunit,*) 'There are ',NBE,' sources'
    write(logunit,*) 'There are ',NBR,' receivers'
!
capdevil's avatar
   
capdevil committed
116
!test
117
!receiver_rad=receiver_rad-500.
capdevil's avatar
   
capdevil committed
118
!receiver_rad=receiver_rad-450000.
capdevil's avatar
capdevil committed
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
!calcule des LmaxS, LmaxT, NmaxS,NmaxT
    call get_Lmax(fmax,LmaxS,LmaxT)
    print*,'LmaxS, LmaxT:',LmaxS, LmaxT
    write(logunit,*) 'LmaxS, LmaxT:',LmaxS, LmaxT
    allocate(NmaxS(0:LmaxS),NmaxT(0:LmaxT))
    call get_Nmax(fmax,NmaxS,NmaxT)
    Nmax_all=max(maxval(NmaxS(:)),maxval(NmaxT(:)))
    Lmax_all=max(LmaxS,LmaxT)
    Nmax_allS=maxval(NmaxS(:))
    Nmax_allT=maxval(NmaxT(:))
    allocate(LLLmax(2))
    allocate(NNNmax(0:lmax_all,2))
    NNNmax(:,:)=0
    LLLmax(1)        =LmaxS
    LLLmax(2)        =LmaxT
    NNNmax(0:LmaxS,1)=NmaxS(:)
    NNNmax(0:LmaxT,2)=NmaxT(:)
    if (nmax_in/=-1) NNNmax(:,:)=min(nmax_in,NNNmax(:,:))
! 
!init b0 etc ...
!
    allocate(al(0:Lmax_all),b0(0:Lmax_all),b1(0:Lmax_all),b2(0:Lmax_all))
    do l=0,lmax_all
       al(l)   =dble(l*(l+1))
       b0(l)   =dsqrt(dble(2*l+1)/(4.0d0*PI))
       b1(l)   =dsqrt(al(l))*b0(l)/2.0d0
       b2(l)   =dsqrt(dble((dble(l)-1)*dble(l)*(dble(l)+1) &
            *(dble(l)+2)))*b0(l)/4.0d0
    enddo
    allocate(eigw(0:Nmax_all,0:Lmax_all,2),qw(0:Nmax_all,0:Lmax_all,2))
    call get_eigenw(eigw,qw,LmaxS,NmaxS,LmaxT,NmaxT)
    print*,'Retreving eigenfrequencies for ',1+NBE,' sources and receiver depth'
    print*,'allocating eigenfunctions (',            &
         ((Nmax_all+1)*(LmaxS+1)*4+         &
         (maxval(NmaxT(:))+1)*(LmaxT+1)*2)*8/1024./1024. &
         ,' Mb/proc)'
    allocate(UU(0:Nmax_allS,0:LmaxS),UUp(0:Nmax_allS,0:LmaxS) &
         ,VV(0:Nmax_allS,0:LmaxS),VVp(0:Nmax_allS,0:LmaxS) )
    print*,'allocating temporary eigenfunctions (',            &
         real(NBE)*((Nmax_all+1)*(LmaxS+1)*4+         &
         (maxval(NmaxT(:))+1)*(LmaxT+1)*2)*8/1024./1024. &
         ,' Mb/proc)'
    allocate(UUs(NBE,0:Nmax_allS,0:LmaxS),UUps(NBE,0:Nmax_allS,0:LmaxS) &
         ,VVs(NBE,0:Nmax_allS,0:LmaxS),VVps(NBE,0:Nmax_allS,0:LmaxS) )
    allocate(WW(0:Nmax_allT,0:LmaxT),WWp(0:Nmax_allT,0:LmaxT))
    allocate(WWs(NBE,0:Nmax_allT,0:LmaxT),WWps(NBE,0:Nmax_allT,0:LmaxT))
    UU=0.0_DP;VV=0.0_DP;UUp=0.0_DP;VVp=0.0_DP;WW=0.0_DP;WWp=0.0_DP
!
    print*,'reading and interpolating eigenfunction ...'
!
    allocate(xdum(4,NBE+1),rad(NBE+1))
!
    print*,'receiver_rad=',receiver_rad
    rad(1)=receiver_rad
    rad(2:NBE+1)=src_coord(1,:)
    !

    print*,'->for S modes ...'
    do l=0,LmaxS
       do n=0,NmaxS(l)
          call read_interpS(n,l,1+NBE,rad,xdum)
          UU( n,l) =xdum(1,1)
          UUp(n,l) =xdum(2,1)
          VV( n,l) =xdum(3,1)
          VVp(n,l) =xdum(4,1)
!
          UUs( :,n,l)=xdum(1,2:)
          UUps(:,n,l)=xdum(2,2:)
          VVs( :,n,l)=xdum(3,2:)
          VVps(:,n,l)=xdum(4,2:)
          
       enddo
    enddo
    deallocate(xdum)
    print*,'->for T modes ...'
    allocate(xdum(2,NBE+1))
    do l=0,LmaxT
       do n=0,NmaxT(l)
          call read_interpT(n,l,1+NBE,rad,xdum)
          WW( n,l) =xdum(1,1)
          WWp(n,l) =xdum(2,1)
!  
          WWs( :,n,l)=xdum(1,2:)
          WWps(:,n,l)=xdum(2,2:)
       enddo
    enddo
    deallocate(xdum,rad)    
capdevil's avatar
   
capdevil committed
206
207
208
209
210
211
!
    if (homogenization>=0) then
       print*,'First order correction of the moment tensor for homogenization requested'
       print*,'Init of module_corr_modes ....'
       call init_corr_modes(NBE,src_coord(1,:),homogenization)
    endif
capdevil's avatar
capdevil committed
212
!
capdevil's avatar
   
capdevil committed
213
    call split_work()
capdevil's avatar
capdevil committed
214
215
216
217
218
219
!
    print*,'init_modes done!'
!-------------------------------------------------
  end subroutine init_modes
!-------------------------------------------------

capdevil's avatar
   
capdevil committed
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
!-------------------------------------------------
  subroutine split_work()
!-------------------------------------------------
    use global_main_param, only: nmin
    implicit none
    real(DP) :: memall
    integer :: q,l,n,reste,pas,ip
    length_qnl_all=0  
    do q=1,2
       do l=ldeb(q),LLLmax(q)
          do n=nmin,NNNmax(l,q)
            length_qnl_all= length_qnl_all+1
          enddo
       enddo
    enddo
    memall=real(length_qnl_all)*real(NBF2)*16./1024./1024.
    nbchunk=int(memall/real(memorymax))+1
    if (nbchunk/=1) print*,'The job will be splitted in ',nbchunk,' parts of ',memorymax,' Mbytes'
    allocate(qnldeb_chunk(0:nbchunk),qnlfin_chunk(0:nbchunk))
!splitting work in case of MPI computing
    pas  =int(length_qnl_all/nbchunk)
    reste=mod(length_qnl_all,nbchunk)
!
    qnldeb_chunk(0)=1
    qnlfin_chunk(0)=pas
    if (0<reste) qnlfin_chunk(0)=qnlfin_chunk(0)+1
    do ip=1,nbchunk-1
       qnldeb_chunk(ip)=qnlfin_chunk(ip-1)+1
       qnlfin_chunk(ip)=qnldeb_chunk(ip)+pas-1
       if (ip<reste) qnlfin_chunk(ip)=qnlfin_chunk(ip)+1
    enddo
!-------------------------------------------------
  end subroutine split_work
!-------------------------------------------------

capdevil's avatar
capdevil committed
255
!-------------------------------------------------
256
  subroutine get_and_sum(traces_)
capdevil's avatar
capdevil committed
257
258
!-------------------------------------------------
    use def_gparam
capdevil's avatar
capdevil committed
259
260
    use global_main_param, only: NBR,nbcomp,NBT,dt,comp,compr,stations_name, &
         rotation,NBE,stacking_sources,station_response,NBZ,NBP,zeros,poles &
261
         ,response_constant, channel,VEL
capdevil's avatar
capdevil committed
262
    implicit none
263
264

    complex(DP), dimension(:,:,:,:), pointer, optional, intent(out)    :: traces_
capdevil's avatar
   
capdevil committed
265
266
267
    integer :: ir,ic,i,itp,itt,is,ichunk
    complex(DP), dimension(:,:),allocatable :: RSkN
    complex(DP), dimension(NBF2)    :: tracestmp
capdevil's avatar
capdevil committed
268
269
270
271
272
273
    complex(DP), dimension(NBF2,nbcomp,NBR,NBE)    :: traces
    real   (SP), dimension(NBT,nbcomp) :: trace_time
    real   (SP), dimension(NBT) :: tmptr
    complex*16 :: alpha=(1.d0,0.d0),beta=(0.d0,0.d0)
    character(len=30) :: name
    real(DP) :: t,p
274
275
    integer :: channel_
    logical :: clean
capdevil's avatar
   
capdevil committed
276

277
278
279
280
281
282
283
284
    if (present(traces_)) then
      channel_=VEL 
      allocate(traces_(NBF2,nbcomp,NBR,NBE))
      clean=.false.
   else
      channel_=channel
      clean=.true.
    endif
capdevil's avatar
   
capdevil committed
285
    traces(:,:,:,:)=cmplx(0.d0,0.d0)
capdevil's avatar
capdevil committed
286
    print*,'Computing RkN and SkN terms ...'
287
    call init_RSkN(channel_,clean)
capdevil's avatar
   
capdevil committed
288
289
290
    do ichunk=0,nbchunk-1
       qnldeb=qnldeb_chunk(ichunk)
       qnlfin=qnlfin_chunk(ichunk)
291
       call init_time_serie(channel_)
capdevil's avatar
   
capdevil committed
292
293
       allocate(RSkN(length_qnl,nbcomp))
       do is=1,NBE
capdevil's avatar
minor    
capdevil committed
294
          if (NBE/=1) print*,'source :',is,NBE
capdevil's avatar
   
capdevil committed
295
296
297
298
299
300
301
302
          do ir=1,NBR
             RSkN(:,:)=czero
             call get_RSkN(ir,is,RSkN)
             do ic=1,nbcomp
                call ZGEMV('T',length_qnl,NBF2,alpha,time_serie ,length_qnl  &
                     ,RSkN(1,ic),1,beta,tracestmp,1)  
                traces(:,ic,ir,is)=traces(:,ic,ir,is)+tracestmp(:)
             enddo
capdevil's avatar
capdevil committed
303
304
          enddo
       enddo
capdevil's avatar
   
capdevil committed
305
       deallocate(time_serie,RSkN)
capdevil's avatar
capdevil committed
306
307
308
309
310
311
312
    enddo
    print*,'converting in time and writing ascii output ...'
    if (NBE==1.or.stacking_sources) then
       do ir=1,NBR
          if (stacking_sources) then
             trace_time(:,:)=0.
             do ic=1,nbcomp
capdevil's avatar
   
capdevil committed
313
                do is=1,NBE  
capdevil's avatar
capdevil committed
314
315
316
317
318
319
320
321
322
                   call convert_f2t(traces(:,ic,ir,is),tmptr(:),is)
                   trace_time(:,ic)=trace_time(:,ic)+tmptr(:)
                enddo
             enddo
          else
             do ic=1,nbcomp
                call convert_f2t(traces(:,ic,ir,1),trace_time(:,ic),1)
             enddo
          endif
capdevil's avatar
capdevil committed
323
324
325
326
327
328
          if (station_response) then
             do ic=1,nbcomp
             call apply_response(NBZ(ir),zeros(:,ir) &
                  ,NBP(ir),poles(:,ir),response_constant(ir),NBT,trace_time(:,ic),dt)
             enddo
          endif
capdevil's avatar
capdevil committed
329
330
331
332
333
334
335
          call apply_source_filter(trace_time)
          if (rotation) call rotate(trace_time,ir,1)
          do ic=1,nbcomp
             if (rotation) then
                write(name,'("U",a1,"_",a4)') compr(ic),stations_name(ir)
             else
                write(name,'("U",a1,"_",a4)') comp(ic) ,stations_name(ir)
336
!write(name,'("U",a1,"_",a6)') comp(ic) ,stations_name(ir)
capdevil's avatar
capdevil committed
337
338
             endif
             open(17,file=name)
capdevil's avatar
   
capdevil committed
339
             write(17,'(f8.1,1x,e14.7)') ((i-1)*dt,trace_time(i,ic),i=1,NBT)
capdevil's avatar
capdevil committed
340
341
342
343
344
345
346
347
             close(17)
          enddo
       enddo
    else
       do is=1,NBE
          do ir=1,NBR
             do ic=1,nbcomp
                call convert_f2t(traces(:,ic,ir,is),trace_time(:,ic),is)
capdevil's avatar
   
capdevil committed
348
349
                if (station_response) call apply_response(NBZ(ir),zeros(:,ir) &
                     ,NBP(ir),poles(:,ir),response_constant(ir),NBT,trace_time(:,ic),dt)
capdevil's avatar
capdevil committed
350
351
352
353
354
355
             enddo
             call apply_source_filter(trace_time)
             if (rotation) call rotate(trace_time,ir,is)
             do ic=1,nbcomp
                if (rotation) then
                   write(name,'("U",a1,"_",a4,"_",i3.3)') compr(ic),stations_name(ir),is
capdevil's avatar
   
capdevil committed
356
                   
capdevil's avatar
capdevil committed
357
358
359
360
                else
                   write(name,'("U",a1,"_",a4,"_",i3.3)') comp(ic) ,stations_name(ir),is
                endif
                open(17,file=name)
capdevil's avatar
   
capdevil committed
361
                write(17,'(f8.1,1x,e14.7)') ((i-1)*dt,trace_time(i,ic),i=1,NBT)
capdevil's avatar
capdevil committed
362
363
364
365
366
                close(17)
             enddo
          enddo
       enddo
    endif
367
368
369
370
371
    if (present(traces_)) then
       traces_=traces
       deallocate(SkN,RkN)
       return
    endif
capdevil's avatar
capdevil committed
372
373
374
!-------------------------------------------------
  end subroutine get_and_sum
!-------------------------------------------------
yann's avatar
yann committed
375
!-------------------------------------------------
376
  subroutine get_and_sum_G(traces_vel)
yann's avatar
yann committed
377
378
379
380
!-------------------------------------------------
    use def_gparam
    use global_main_param, only: NBR,NBT,dt,comp,compr,stations_name, &
         rotation,NBE,stacking_sources,station_response,NBZ,NBP,zeros,poles &
381
         ,response_constant,channel
yann's avatar
yann committed
382
    implicit none
383
384
    complex(DP), dimension(:,:,:,:), optional, intent(in)    :: traces_vel
!
yann's avatar
yann committed
385
    integer :: ir,iN,iNp,i,itp,itt,is,ichunk,ic,jc
386
    complex(DP), dimension(:,:,:),allocatable :: GSkN
yann's avatar
yann committed
387
388
389
390
391
392
393
394
395
    complex(DP), dimension(NBF2)    :: tracestmp
    complex(DP), dimension(NBF2,1:3,1:3,NBR,NBE)    :: traces
    real   (SP), dimension(NBT,1:3,1:3) :: trace_time
    real   (SP), dimension(NBT) :: tmptr
    complex*16 :: alpha=(1.d0,0.d0),beta=(0.d0,0.d0)
    character(len=30) :: name
    real(DP) :: t,p

    traces(:,:,:,:,:)=cmplx(0.d0,0.d0)
396
397
    print*,'Computing GkN and SkN terms ...',nbchunk
    call init_RSkN(channel,clean=.true.)
yann's avatar
yann committed
398
399
400
    do ichunk=0,nbchunk-1
       qnldeb=qnldeb_chunk(ichunk)
       qnlfin=qnlfin_chunk(ichunk)
401
402
       call init_time_serie(channel)
       allocate(GSkN(length_qnl,-1:1,-1:1))
yann's avatar
yann committed
403
404
405
       do is=1,NBE
          if (NBE/=1) print*,'source :',is,NBE
          do ir=1,NBR
406
407
             GSkN(:,:,:)=czero
             call get_GSkN(ir,is,GSkN)
yann's avatar
yann committed
408
409
410
             do iN=-1,1
                do iNp=-1,1
                   call ZGEMV('T',length_qnl,NBF2,alpha,time_serie ,length_qnl  &
411
                        ,GSkN(1,iN,iNp),1,beta,tracestmp,1)  
yann's avatar
yann committed
412
413
414
415
416
                   traces(:,iN+2,iNp+2,ir,is)=traces(:,iN+2,iNp+2,ir,is)+tracestmp(:)
                enddo
             enddo
          enddo
       enddo
417
       deallocate(time_serie,GSkN)
yann's avatar
yann committed
418
    enddo
419
!
yann's avatar
yann committed
420
    call convert_gradient(traces)
421
422
    if (present(traces_vel)) call correct_spherical(traces,traces_vel)
!
yann's avatar
yann committed
423
424
    print*,'converting in time and writing ascii output ...'
    if (NBE==1.or.stacking_sources) then
425
       is=1
yann's avatar
yann committed
426
427
428
429
430
431
432
433
434
435
436
437
438
439
       do ir=1,NBR
          if (stacking_sources) then
             trace_time(:,:,:)=0.
             do ic=1,3
                do jc=1,3
                   do is=1,NBE  
                      call convert_f2t(traces(:,ic,jc,ir,is),tmptr(:),is)
                      trace_time(:,ic,jc)=trace_time(:,ic,jc)+tmptr(:)
                   enddo
                enddo
             enddo
          else
             do jc=1,3
                do ic=1,3
440
441
                   trace_time(:,ic,jc)=0.d0
                   call convert_f2t(traces(:,ic,jc,ir,is),trace_time(:,ic,jc),1)
yann's avatar
yann committed
442
443
444
445
446
447
448
449
450
451
452
453
                enddo
             enddo
          endif
          if (station_response) then
             do jc=1,3
                do ic=1,3
                   call apply_response(NBZ(ir),zeros(:,ir) &
                        ,NBP(ir),poles(:,ir),response_constant(ir),NBT,trace_time(:,ic,jc),dt)
                enddo
             enddo
          endif
          call apply_source_filter_gr(trace_time)
454
!gradiant:
yann's avatar
yann committed
455
456
          do jc=1,3
             do ic=1,3
457
                write(name,'("grad",a1,a1,"_",a4)') comp(ic),comp(jc) ,stations_name(ir)
yann's avatar
yann committed
458
459
460
461
462
                open(17,file=name)
                write(17,'(f8.1,1x,e14.7)') ((i-1)*dt,trace_time(i,ic,jc),i=1,NBT)
                close(17)
             enddo
          enddo
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
!rotational: 
!Z
          write(name,'("rot",a1,"_",a4)') comp(1),stations_name(ir)
          open(17,file=name)
          write(17,'(f8.1,1x,e14.7)') ((i-1)*dt,-trace_time(i,2,3)+trace_time(i,3,2),i=1,NBT)
          close(17)
!N
          write(name,'("rot",a1,"_",a4)') comp(2),stations_name(ir)
          open(17,file=name)
          write(17,'(f8.1,1x,e14.7)') ((i-1)*dt,-trace_time(i,3,1)+trace_time(i,1,3),i=1,NBT)
          close(17)
!E
          write(name,'("rot",a1,"_",a4)') comp(3),stations_name(ir)
          open(17,file=name)
          write(17,'(f8.1,1x,e14.7)') ((i-1)*dt,-trace_time(i,1,2)+trace_time(i,2,1),i=1,NBT)
          close(17)
yann's avatar
yann committed
479
480
481
482
483
484
485
486
487
488
489
       enddo
    else
       do is=1,NBE
          do ir=1,NBR
             do jc=1,3
                do ic=1,3
                   call convert_f2t(traces(:,ic,jc,ir,is),trace_time(:,ic,jc),is)
                   if (station_response) call apply_response(NBZ(ir),zeros(:,ir) &
                        ,NBP(ir),poles(:,ir),response_constant(ir),NBT,trace_time(:,ic,jc),dt)
                enddo
             enddo
490
             call apply_source_filter_gr(trace_time)             
yann's avatar
yann committed
491
492
             do jc=1,3
                do ic=1,3
493
                   write(name,'("grad",a1,a1,"_",a4,"_",i3.3)') comp(ic),comp(jc) ,stations_name(ir),is
yann's avatar
yann committed
494
495
496
497
498
                   open(17,file=name)
                   write(17,'(f8.1,1x,e14.7)') ((i-1)*dt,trace_time(i,ic,jc),i=1,NBT)
                   close(17)
                enddo
             enddo
499
500
501
502
503
504
505
506
507
508
509
510
511
!rotational: 
             write(name,'("rot",a1,"_",a4,"_",i3.3)') comp(1),stations_name(ir),is
             open(17,file=name)
             write(17,'(f8.1,1x,e14.7)') ((i-1)*dt,-trace_time(i,2,3)+trace_time(i,3,2),i=1,NBT)
             close(17)
             write(name,'("rot",a1,"_",a4,"_",i3.3)') comp(2),stations_name(ir),is
             open(17,file=name)
             write(17,'(f8.1,1x,e14.7)') ((i-1)*dt,-trace_time(i,3,1)+trace_time(i,1,3),i=1,NBT)
             close(17)
             write(name,'("rot",a1,"_",a4,"_",i3.3)') comp(3),stations_name(ir),is
             open(17,file=name)
             write(17,'(f8.1,1x,e14.7)') ((i-1)*dt,-trace_time(i,1,2)+trace_time(i,2,1),i=1,NBT)
             close(17)
yann's avatar
yann committed
512
513
514
515
          enddo
       enddo
    endif
!-------------------------------------------------
516
  end subroutine get_and_sum_G
yann's avatar
yann committed
517
!-------------------------------------------------
capdevil's avatar
capdevil committed
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
!-------------------------------------------------
  subroutine rotate(trace,ir,is)
!-------------------------------------------------
    use def_gparam
    use util_fctp
    use global_main_param, only: NBT,stations_name
    implicit none
    real(SP), dimension(:,:), intent(out) :: trace
    integer, intent(in) :: ir,is
!
    real(DP) :: tr,pr,ts,ps,asr,bsr,gsr,bazm
    real(SP), dimension(NBT,2) :: tmp
!    

    tr=rec_coord(1,ir)
    pr=rec_coord(2,ir)
    ts=src_coord(2,is)
    ps=src_coord(3,is)
    call euler(ts,ps,tr,pr,asr,bsr,gsr)
    bazm=-asr
    print*,stations_name(ir),', bazm= ',bazm*180./3.1415926,', delta=',bsr*180./3.1415926
    tmp(:,:)=trace(:,2:3)
    trace(:,2)=-tmp(:,1)*cos(bazm)+tmp(:,2)*sin(bazm)
    trace(:,3)=-tmp(:,1)*sin(bazm)-tmp(:,2)*cos(bazm)
!-------------------------------------------------
  end subroutine rotate
!-------------------------------------------------
!-------------------------------------------------
  subroutine convert_f2t(tracef,tracet,is)
!-------------------------------------------------
    use def_gparam
    use util_fctp
capdevil's avatar
   
capdevil committed
550
    use global_main_param, only: NBT,dt,nbcomp ,t_delay_sources,half_duration
capdevil's avatar
capdevil committed
551
    implicit none
capdevil's avatar
   
capdevil committed
552
    complex(DP), dimension(NBF2), intent(inout) :: tracef
capdevil's avatar
capdevil committed
553
554
555
556
557
558
559
560
    real, dimension(NBT), intent(out) :: tracet
    integer, intent(in) :: is
!
!
    complex(DP), dimension(NBT) :: trace_test
    real(DP), dimension(NBF*coef_nbt) :: trace_testr
    real, dimension(NBT,1,1) :: trace_r
    complex(DP), dimension(NBF2*coef_nbt) :: traces_testf
capdevil's avatar
   
capdevil committed
561
562
    integer :: ic,ir,i,num,if
    real(DP) :: df,w
capdevil's avatar
capdevil committed
563
!      
capdevil's avatar
   
capdevil committed
564
565
566
    if (half_duration(is)  /= 0.0d0) then
       df = 1._DP/(NBF2*dtf)
       do if=2,NBF+1
capdevil's avatar
   
capdevil committed
567
          w=2._DP*PI*(if-1)*df*half_duration(is)
capdevil's avatar
   
capdevil committed
568
569
570
          tracef(if)=tracef(if)*sin(w)/w
       enddo
       do if=NBF+1,NBF2
capdevil's avatar
   
capdevil committed
571
          w=-2._DP*PI*(NBF2-if+1)*df*half_duration(is)
capdevil's avatar
   
capdevil committed
572
573
574
          tracef(if)=tracef(if)*sin(w)/w
       enddo      
    endif
capdevil's avatar
capdevil committed
575
576
577
    traces_testf(:)=(0.d0,0.d0)
    traces_testf(1:NBF+1)=tracef(1:NBF+1)
    traces_testf((2*coef_nbt*NBF)-NBF+1:coef_nbt*NBF2) &
capdevil's avatar
   
capdevil committed
578
579
         =tracef(NBF+1:NBF2)

capdevil's avatar
capdevil committed
580
581
582
583
584
585
586
587
588
589
    call dfour1(traces_testf,NBF2*coef_nbt,-1)
!SP + * frequency step
    trace_testr(1:NBF*coef_nbt)=real(traces_testf(1:NBF*coef_nbt)/NBF2/dtf,DP)
    call resample(trace_testr,t_delay_sources(is),dtf2,NBF*coef_nbt &
         ,tracet,dt,NBT)
!-------------------------------------------------
  end subroutine convert_f2t
!-------------------------------------------------

!-------------------------------------------------
590
  subroutine init_RSkN(channel_,clean)
capdevil's avatar
capdevil committed
591
592
!-------------------------------------------------
    use def_gparam
593
    use global_main_param, only: nbcomp,NBE,RA,DEP,VEL,ACC,GRA
capdevil's avatar
capdevil committed
594
    implicit none
595
596
    integer, intent(in) :: channel_
    logical, intent(in) :: clean
capdevil's avatar
capdevil committed
597
598
599
    integer :: l
    real :: mega
!
capdevil's avatar
   
capdevil committed
600
601
    scale =1._DP/sqrt(rhobar*ra**3)
    dscale=1._DP/sqrt(rhobar*ra**3)/ra
capdevil's avatar
capdevil committed
602
603
604
!
    mega=(3*nbcomp*(real(Nmax_allS+1)*(LmaxS+1)+real((Nmax_allT+1)*(LmaxT+1))) + &
          5*NBE*(real(Nmax_allS+1)*(LmaxS+1)+real((Nmax_allT+1)*(LmaxT+1))) )/1024./1204.
605

yann's avatar
yann committed
606
    allocate(SkN(-2:2       ,0:Nmax_all,0:Lmax_all,2,NBE))
capdevil's avatar
capdevil committed
607
    SkN=0.0_DP
608
    print*,'->SkN ...',channel_
capdevil's avatar
capdevil committed
609
    call SkNfun(SkN)
610
    select case (channel_)
yann's avatar
yann committed
611
612
613
614
615
    case (DEP,VEL,ACC)
       print*,'->RkN ...'
       allocate(RkN(-1:1,nbcomp,0:Nmax_all,0:Lmax_all,2))
       RkN=0.0_DP
       call RkNfun(RkN)
616
617
618
619
620
    case(GRA)
       print*,'->GkN ...'
       allocate(GkN(-1:1,-1:1,0:Nmax_all,0:Lmax_all,2))
       GkN=0.0_DP
       call GkNfun(GkN)
yann's avatar
yann committed
621
    end select
capdevil's avatar
capdevil committed
622
! Cleaning unecessary arrays:
623
    if (clean) deallocate(UUs,UUps,VVs,VVps,WWs,WWps,Msrc)   
capdevil's avatar
capdevil committed
624
625
!-------------------------------------------------
  end subroutine init_RSkN
yann's avatar
yann committed
626
!------------------------------------------------
capdevil's avatar
capdevil committed
627
628
629
630
!-------------------------------------------------
  subroutine SkNfun(S)
!-------------------------------------------------
    use def_gparam
capdevil's avatar
   
capdevil committed
631
632
    use global_main_param, only:RA,NBE,source_force, homogenization
    use corr_modes
capdevil's avatar
capdevil committed
633
634
    IMPLICIT NONE 
    complex(DP), dimension(-2:,0:,0:,:,:), intent(out) :: S
capdevil's avatar
   
capdevil committed
635
636
    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
capdevil's avatar
capdevil committed
637
638
639
    real(DP), dimension(6) :: am
    integer :: is,n,l
!
capdevil's avatar
   
capdevil committed
640

capdevil's avatar
   
capdevil committed
641
642
643
644
645
646
647
648
649
650
    if (source_force) then
       do l=0,LmaxS
          do n=0,NmaxS(l)
             do is=1,NBE
                vir=Msrc(1,is)
                vit=Msrc(2,is)
                vip=Msrc(3,is)
                
                u =UUs (is,n,l)
                v =VVs (is,n,l)
capdevil's avatar
   
capdevil committed
651
652
653
654
655
656
                if (homogenization>=0) then
                   rr=src_coord(1,is)/RA
                   degu =UUps(is,n,l)
                   degv =VVps(is,n,l) 
                   call corrmodesS(is,eigw(n,l,1),l,rr,u,v,degu,degv)
                endif
capdevil's avatar
   
capdevil committed
657
                xu= b0(l)*u
capdevil's avatar
   
capdevil committed
658
                xv= b1(l)*v
capdevil's avatar
   
capdevil committed
659
                S( 0,n,l,1,is)=xu*vir*scale
capdevil's avatar
   
capdevil committed
660
661
                S( 1,n,l,1,is)=xv*dcmplx(-vit,vip)*scale
                S(-1,n,l,1,is)=xv*dcmplx(+vit,vip)*scale
capdevil's avatar
   
capdevil committed
662
663
664
665
666
667
668
669
670
             enddo
          enddo
       enddo
       do l=1,LmaxT
          do n=0,NmaxT(l)
             do is=1,NBE
!
                vit=Msrc(2,is)
                vip=Msrc(3,is)
capdevil's avatar
   
capdevil committed
671
672
673
674
                if (homogenization>=0) then
                   dfw   =WWps(is,n,l)
                   call corrmodesT(is,eigw(n,l,2),l,rr,fw,dfw)
                endif
capdevil's avatar
   
capdevil committed
675
                x= b1(l)*WWs (is,n,l)
capdevil's avatar
   
capdevil committed
676
677
!

capdevil's avatar
   
capdevil committed
678
679
                S ( 1,n,l,2,is)=x*dcmplx(+vip,vit)*scale
                S (-1,n,l,2,is)=x*dcmplx(-vip,vit)*scale
capdevil's avatar
   
capdevil committed
680
681
682
             enddo
          enddo
       enddo
capdevil's avatar
   
capdevil committed
683
684
!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
capdevil's avatar
   
capdevil committed
685
686
687
688
689
690
691
692
693
    else
       do l=0,LmaxS
          do n=0,NmaxS(l)
             do is=1,NBE
                rr=src_coord(1,is)/RA
                egu  =UUs (is,n,l)
                egv  =VVs (is,n,l)
                degu =UUps(is,n,l)
                degv =VVps(is,n,l) 
capdevil's avatar
   
capdevil committed
694
                if (homogenization>=0) call corrmodesS(is,eigw(n,l,1),l,rr,egu,egv,degu,degv)
capdevil's avatar
   
capdevil committed
695
                am(:)=Msrc(:,is)
capdevil's avatar
capdevil committed
696
!    
capdevil's avatar
   
capdevil committed
697
698
699
700
701
702
703
704
705
706
                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
capdevil's avatar
capdevil committed
707
708
          enddo
       enddo
capdevil's avatar
   
capdevil committed
709
710
711
712
713
714
       do l=1,LmaxT
          do n=0,NmaxT(l)
             do is=1,NBE
                rr=src_coord(1,is)/RA
                fw    =WWs (is,n,l)
                dfw   =WWps(is,n,l)
capdevil's avatar
   
capdevil committed
715
                if (homogenization>=0) call corrmodesT(is,eigw(n,l,2),l,rr,fw,dfw)
capdevil's avatar
   
capdevil committed
716
                am(:) =Msrc(:,is)
capdevil's avatar
capdevil committed
717
!
capdevil's avatar
   
capdevil committed
718
719
720
721
722
723
724
725
                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
capdevil's avatar
capdevil committed
726
727
          enddo
       enddo
capdevil's avatar
   
capdevil committed
728
    endif
capdevil's avatar
capdevil committed
729
730
731
732
733
734
735
736
!--------------------------------------------------
  end subroutine SkNfun
!--------------------------------------------------

!--------------------------------------------------
  subroutine RkNfun(R)
!--------------------------------------------------
    use def_gparam
capdevil's avatar
   
capdevil committed
737
738
    use global_main_param, only: NBE,nbcomp,homogenization
    use corr_modes
capdevil's avatar
capdevil committed
739
740
    IMPLICIT NONE
    complex(DP), dimension(-1:,:,0:,0:,:), intent(out) :: R
capdevil's avatar
   
capdevil committed
741
742
    real(DP) :: xu,xv,u,v,vir,vip,vit,x,fw,rr,degu,degv,dfw
    integer :: ic,n,l,is
capdevil's avatar
capdevil committed
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
!
    do l=0,LmaxS
       do n=0,NmaxS(l)
          do ic=1,nbcomp !boucle sur les composantes
             select case(ic)
             case(1) ! verticale
                vir=1.0_DP
                vit=0.0_DP
                vip=0.0_DP
             case(2) !Nord=-theta
                vir=0.0_DP
                vit=-1.0_DP
                vip=0.0_DP
             case(3) !Est=phi
                vir=0.0_DP
                vit=0.0_DP
                vip=1.0_DP
             end select
             u =UU (n,l)
capdevil's avatar
   
capdevil committed
762
763
764
765
766
767
768
769
             v =VV (n,l) 
             if (homogenization>=0) then
                rr=1.d0
                is=NBE+1
                degu =UUp(n,l)
                degv =VVp(n,l) 
                call corrmodesS(is,eigw(n,l,1),l,rr,u,v,degu,degv)
             endif   
capdevil's avatar
capdevil committed
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
             xu= b0(l)*u
             xv=-b1(l)*v
             R( 0,ic,n,l,1)=xu*vir*scale
             R( 1,ic,n,l,1)=xv*dcmplx( vit,vip)*scale
             R(-1,ic,n,l,1)=xv*dcmplx(-vit,vip)*scale
          enddo
       enddo
    enddo
    do l=1,LmaxT
       do n=0,NmaxT(l)
          do ic=1,nbcomp !boucle sur les composantes
             select case(ic)
             case(1) ! verticale
                vir=1.0_DP
                vit=0.0_DP
                vip=0.0_DP
             case(2) !theta
                vir=0.0_DP
                vit=-1.0_DP
                vip=0.0_DP
             case(3) !phi
                vir=0.0_DP
                vit=0.0_DP
                vip=1.0_DP
             end select
capdevil's avatar
   
capdevil committed
795
796
797
798
799
800
801
802
             fw    =WW (n,l)
             if (homogenization>=0) then
                rr=1.d0
                is=NBE+1
                dfw   =WWp(n,l)
                call corrmodesT(is,eigw(n,l,2),l,rr,fw,dfw)
             endif
             x=-b1(l)*fw
capdevil's avatar
capdevil committed
803
804
805
806
807
808
809
810
811
!
             R ( 1,ic,n,l,2)=x*dcmplx(-vip,vit)*scale
             R (-1,ic,n,l,2)=x*dcmplx( vip,vit)*scale
          enddo
       enddo
    enddo
!--------------------------------------------------
  end subroutine RkNfun
!--------------------------------------------------
Yann CAPDEVILLE's avatar
Yann CAPDEVILLE committed
812
!--------------------------------------------------
yann's avatar
yann committed
813
  subroutine EkNfun(EE)
814
!strain
Yann CAPDEVILLE's avatar
Yann CAPDEVILLE committed
815
816
!--------------------------------------------------
    use def_gparam
yann's avatar
yann committed
817
    use global_main_param, only: NBE,nbcomp,homogenization,s2
Yann CAPDEVILLE's avatar
Yann CAPDEVILLE committed
818
819
    use corr_modes
    IMPLICIT NONE
yann's avatar
yann committed
820
821
822
823
824
    complex(DP), dimension(-1:,-1:,0:,0:,:), intent(out) :: EE
    real(DP) :: f,x,z,rd
    integer :: n,l
!
    rd=1._DP
Yann CAPDEVILLE's avatar
Yann CAPDEVILLE committed
825
826
    do l=0,LmaxS
       do n=0,NmaxS(l)
yann's avatar
yann committed
827
828
829
830
831
832
833
834
835
836
837
838
839
840
!
          f    = (2.0d0*UU(n ,l )-al(l )*VV(n ,l ))/rd
          x    = VVp(n ,l )+(UU(n ,l )-VV(n ,l ))/rd
!
          EE( 0, 0,n ,l,1) =       b0(l)*UUp (n ,l )
          EE( 1, 0,n ,l,1) =    s2*b1(l)*x/2._DP
          EE( 1, 1,n ,l,1) = 2._DP*b2(l)*VV (n ,l )/rd
          EE(-1, 1,n ,l,1) =      -b0(l)*f/2._DP
!
          EE( 1,-1,n ,l,1) = EE(-1, 1,n ,l,1) 
          EE(-1,-1,n ,l,1) = EE( 1, 1,n ,l,1)
          EE( 0, 1,n ,l,1) = EE( 1, 0,n ,l,1) 
          EE(-1, 0,n ,l,1) = EE( 1, 0,n ,l,1) 
          EE( 0,-1,n ,l,1) = EE( 1, 0,n ,l,1) 
Yann CAPDEVILLE's avatar
Yann CAPDEVILLE committed
841
842
843
844
845
       enddo
    enddo
    do l=1,LmaxT
       do n=0,NmaxT(l)
!
yann's avatar
yann committed
846
          z    = WWp(n ,l )-WW(n ,l )/rd
Yann CAPDEVILLE's avatar
Yann CAPDEVILLE committed
847
!
yann's avatar
yann committed
848
          EE( :, :,n ,l,2) = 0._DP
Yann CAPDEVILLE's avatar
Yann CAPDEVILLE committed
849
!
yann's avatar
yann committed
850
851
852
853
854
855
856
          EE( 1, 1,n ,l,2) = 2._DP*b2(l)*WW (n ,l )/rd
          EE( 0, 1,n ,l,2) = s2*b1(l)*z/2._DP
!
          EE(-1,-1,n ,l,2) =-EE( 1, 1,n ,l,2)
          EE( 1, 0,n ,l,2) = EE( 0, 1,n ,l,2) 
          EE( 0,-1,n ,l,2) =-EE( 0, 1,n ,l,2)  
          EE(-1, 0,n ,l,2) =-EE( 0, 1,n ,l,2)
Yann CAPDEVILLE's avatar
Yann CAPDEVILLE committed
857
858
859
860
861
       enddo
    enddo
!--------------------------------------------------
  end subroutine EkNfun
!--------------------------------------------------
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
!--------------------------------------------------
  subroutine GkNfun(GG)
!Gradient
!--------------------------------------------------
    use def_gparam
    use global_main_param, only: NBE,nbcomp,homogenization,s2
    use corr_modes
    IMPLICIT NONE
    complex(DP), dimension(-1:,-1:,0:,0:,:), intent(out) :: GG
    real(DP) :: f,x,z,rd
    integer :: n,l
!
    rd=1._DP
    do l=0,LmaxS
       do n=0,NmaxS(l)
!
          GG( 0, 0,n ,l,1) =       b0(l)*UUp (n ,l )
          GG( 1, 1,n ,l,1) = 2._DP*b2(l)*VV  (n ,l )/rd
          GG(-1,-1,n ,l,1) = GG( 1, 1,n ,l,1)
!         
          GG( 1, 0,n ,l,1) =    s2*b1(l)*VVp(n,l)
          GG(-1, 0,n ,l,1) = GG( 1, 0,n ,l,1) 
!
          GG( 0, 1,n ,l,1) = s2*b1(l)/rd*(UU(n,l)-VV(n,l))
          GG( 0,-1,n ,l,1) = GG( 0, 1,n ,l,1) 
!
          GG(-1, 1,n ,l,1) =  b0(l)*(-2._DP*UU(n,l)+al(l)*VV(n,l))/2._DP/rd
          GG( 1,-1,n ,l,1) = GG(-1, 1,n ,l,1) 
!
       enddo
    enddo
    do l=1,LmaxT
       do n=0,NmaxT(l)
!
          z    = WWp(n ,l )-WW(n ,l )/rd
!
!
          GG( 0, 0,n ,l,2) = 0._DP
          GG( 1, 1,n ,l,2) = 2._DP*b2(l)*WW (n ,l )/rd
          GG(-1,-1,n ,l,2) =-GG( 1, 1,n ,l,2)

          GG( 1, 0,n ,l,2) =    s2*b1(l)*WWp(n,l)
          GG(-1, 0,n ,l,2) = -GG( 1, 0,n ,l,2) 
!
          GG( 0, 1,n ,l,2) = s2*b1(l)/rd*(-WW(n,l))
          GG( 0,-1,n ,l,2) = -GG( 0, 1,n ,l,2) 
!
          GG(-1, 1,n ,l,2) =  -b0(l)*al(l)*WW(n,l)/2._DP/rd
          GG( 1,-1,n ,l,2) = -GG(-1, 1,n ,l,2) 
       enddo
    enddo
    GG(:,:,:,:,:)=GG(:,:,:,:,:)*dscale
!--------------------------------------------------
  end subroutine GkNfun
!--------------------------------------------------
yann's avatar
yann committed
917
918
!--------------------------------------------------
  subroutine convert_gradient(tr)
919
920
921
!goes from {-1,0,1}^2 to {Z,N,E}^2
!attention, Nord=-theta. Du coup, il y a un changement de signe dans la premiere ligne
!du Cia par rapport  Phinney & Burridge
yann's avatar
yann committed
922
923
924
925
926
927
!--------------------------------------------------
    use def_gparam
    use global_main_param, only: NBR,NBE,nbcomp,NBT,s2,is2,II,czero,cone
    implicit none
    complex(DP), dimension (:,:,:,:,:), intent(inout) :: tr
!
928
929
930
    complex(DP), dimension(1:3,1:3), parameter :: Cia=reshape( (/czero,-is2*cone,-II*is2,cone,czero,czero,czero,is2*cone, &
         -II*is2/),(/3,3/))
    complex(DP), dimension(NBF2,3,3) :: ctmp
yann's avatar
yann committed
931
932
933
    integer :: is,ir,i,j,if,iN,iNp
    do is=1,NBE
       do ir=1,NBR
934
          ctmp(:,:,:)=0.d0
yann's avatar
yann committed
935
936
937
938
939
          do i=1,3
             do j=1,3
                do if=1,NBF2
                   do iN=1,3
                      do iNp=1,3
940
                         ctmp(if,i,j)=ctmp(if,i,j)+Cia(i,iN)*Cia(j,iNp)*tr(if,iN,iNp,ir,is)
yann's avatar
yann committed
941
942
943
944
945
                      enddo
                   enddo
                enddo
             enddo            
          enddo
946
          tr(:,:,:,ir,is)=ctmp(:,:,:)
yann's avatar
yann committed
947
948
949
950
951
       enddo
    enddo
!--------------------------------------------------
  end subroutine convert_gradient
!--------------------------------------------------
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
!--------------------------------------------------
  subroutine correct_spherical(tr,trv)
!warning: component 2 of velocity (try) is the North component (and therefore -theta!)
!--------------------------------------------------
    use def_gparam
    use global_main_param, only: NBR,NBE,nbcomp,NBT,s2,is2,II,czero,cone
    implicit none
    complex(DP), dimension (:,:,:,:,:), intent(inout) :: tr
    complex(DP), dimension (:,:,:,:), intent(in) :: trv
!
    real(DP) :: thr,str,cotr
    integer :: is,ir
    do is=1,NBE
       do ir=1,NBR   
          thr=rec_coord(1,ir)
          str=sin(thr)
          cotr=cos(thr)/sin(thr)
!
          tr(:,2,2 ,ir,is)=tr(:,2,2 ,ir,is)-trv(:,1,ir,is)/receiver_rad !checked
          tr(:,3,3 ,ir,is)=tr(:,3,3 ,ir,is)-(trv(:,1,ir,is)-cotr*trv(:,2,ir,is))/receiver_rad !checked
!
          tr(:,1,2 ,ir,is)=tr(:,1,2 ,ir,is)+trv(:,2,ir,is)/receiver_rad !checked
          tr(:,1,3 ,ir,is)=tr(:,1,3 ,ir,is)+str*trv(:,3,ir,is)/receiver_rad !checked
          tr(:,2,3 ,ir,is)=tr(:,2,3 ,ir,is)-cotr*trv(:,3,ir,is)/receiver_rad !checked
          
       enddo
    enddo
!-------------------------------------------------
  end subroutine correct_spherical
!-------------------------------------------------
capdevil's avatar
capdevil committed
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
!-------------------------------------------------
  subroutine get_RSkN(ir,is,RSkN)
!-------------------------------------------------
    use def_gparam 
    use util_fctp
    use global_main_param, only: nbcomp,source_delay,t_delay_sources,nmin
    implicit none
    integer, intent(in) :: ir,is
    complex(DP), dimension(:,:), intent(out) :: RSkN
!
    integer     :: iqnl,bNp,qdeb,qfin,l,n,q,bN,ird,ic
    complex(DP) :: sd,wq,eiwtd,rs
    real(DP) :: ts,ps,tr,pr,asr,bsr,gsr,x
    complex(DP), dimension(-2:2) :: csd
    complex(DP), dimension(-1:1) :: csr
    real(DP), dimension(-2:2,-2:2,0:Lmax_all) :: dsr
!
    tr=rec_coord(1,ir)
    pr=rec_coord(2,ir)