module_modes.f90

!--------------------------------------------------------------------------------
module modes
!--------------------------------------------------------------------------------
  use def_gparam
  implicit none
  public :: init_modes,get_and_sum,get_and_sum_G
  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
  real(DP), parameter :: coef_damp1=1.03_DP, coef_damp2=1.15_DP
  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
!
  integer :: length_qnl,length_qnl_all,nbchunk
  integer, parameter :: memorymax=512
!  integer, parameter :: memorymax=128000000
  integer, dimension(:), allocatable :: qnldeb_chunk,qnlfin_chunk
  integer :: qnldeb,qnlfin,iqnld
!
  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
  complex(DP), dimension(:,:,:,:,:  ), allocatable :: GkN
!
  real(DP), parameter :: rhobar=5515._DP
  real(DP) :: scale2,scale1
  real(DP) :: scale,dscale,receiver_rad
!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
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  &
                               ,NBE,NBR,nbcomp,nmax_in, homogenization
    use time_function 
    use util_fctp
    use corr_modes
!
    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
!
    t0delay=0.d0
    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)    
!RA is now initialized
    scale2=1.D0/(rhobar*RA*RA*RA)
    scale1=scale2/RA/RA
!
!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'
!
!test
!receiver_rad=receiver_rad-500.
!receiver_rad=receiver_rad-450000.
!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)    
!
    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
!
    call split_work()
!
    print*,'init_modes done!'
!-------------------------------------------------
  end subroutine init_modes
!-------------------------------------------------

!-------------------------------------------------
  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
!-------------------------------------------------

!-------------------------------------------------
  subroutine get_and_sum(traces_)
!-------------------------------------------------
    use def_gparam
    use global_main_param, only: NBR,nbcomp,NBT,dt,comp,compr,stations_name, &
         rotation,NBE,stacking_sources,station_response,NBZ,NBP,zeros,poles &
         ,response_constant, channel,VEL
    implicit none

    complex(DP), dimension(:,:,:,:), pointer, optional, intent(out)    :: traces_
    integer :: ir,ic,i,itp,itt,is,ichunk
    complex(DP), dimension(:,:),allocatable :: RSkN
    complex(DP), dimension(NBF2)    :: tracestmp
    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
    integer :: channel_
    logical :: clean

    if (present(traces_)) then
      channel_=VEL 
      allocate(traces_(NBF2,nbcomp,NBR,NBE))
      clean=.false.
   else
      channel_=channel
      clean=.true.
    endif
    traces(:,:,:,:)=cmplx(0.d0,0.d0)
    print*,'Computing RkN and SkN terms ...'
    call init_RSkN(channel_,clean)
    do ichunk=0,nbchunk-1
       qnldeb=qnldeb_chunk(ichunk)
       qnlfin=qnlfin_chunk(ichunk)
       call init_time_serie(channel_)
       allocate(RSkN(length_qnl,nbcomp))
       do is=1,NBE
          if (NBE/=1) print*,'source :',is,NBE
          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
          enddo
       enddo
       deallocate(time_serie,RSkN)
    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
                do is=1,NBE  
                   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
          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
          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)
!write(name,'("U",a1,"_",a6)') comp(ic) ,stations_name(ir)
             endif
             open(17,file=name)
             write(17,'(f8.1,1x,e14.7)') ((i-1)*dt,trace_time(i,ic),i=1,NBT)
             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)
                if (station_response) call apply_response(NBZ(ir),zeros(:,ir) &
                     ,NBP(ir),poles(:,ir),response_constant(ir),NBT,trace_time(:,ic),dt)
             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
                   
                else
                   write(name,'("U",a1,"_",a4,"_",i3.3)') comp(ic) ,stations_name(ir),is
                endif
                open(17,file=name)
                write(17,'(f8.1,1x,e14.7)') ((i-1)*dt,trace_time(i,ic),i=1,NBT)
                close(17)
             enddo
          enddo
       enddo
    endif
    if (present(traces_)) then
       traces_=traces
       deallocate(SkN,RkN)
       return
    endif
!-------------------------------------------------
  end subroutine get_and_sum
!-------------------------------------------------
!-------------------------------------------------
  subroutine get_and_sum_G(traces_vel)
!-------------------------------------------------
    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 &
         ,response_constant,channel
    implicit none
    complex(DP), dimension(:,:,:,:), optional, intent(in)    :: traces_vel
!
    integer :: ir,iN,iNp,i,itp,itt,is,ichunk,ic,jc
    complex(DP), dimension(:,:,:),allocatable :: GSkN
    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)
    print*,'Computing GkN and SkN terms ...',nbchunk
    call init_RSkN(channel,clean=.true.)
    do ichunk=0,nbchunk-1
       qnldeb=qnldeb_chunk(ichunk)
       qnlfin=qnlfin_chunk(ichunk)
       call init_time_serie(channel)
       allocate(GSkN(length_qnl,-1:1,-1:1))
       do is=1,NBE
          if (NBE/=1) print*,'source :',is,NBE
          do ir=1,NBR
             GSkN(:,:,:)=czero
             call get_GSkN(ir,is,GSkN)
             do iN=-1,1
                do iNp=-1,1
                   call ZGEMV('T',length_qnl,NBF2,alpha,time_serie ,length_qnl  &
                        ,GSkN(1,iN,iNp),1,beta,tracestmp,1)  
                   traces(:,iN+2,iNp+2,ir,is)=traces(:,iN+2,iNp+2,ir,is)+tracestmp(:)
                enddo
             enddo
          enddo
       enddo
       deallocate(time_serie,GSkN)
    enddo
!
    call convert_gradient(traces)
    if (present(traces_vel)) call correct_spherical(traces,traces_vel)
!
    print*,'converting in time and writing ascii output ...'
    if (NBE==1.or.stacking_sources) then
       is=1
       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
                   trace_time(:,ic,jc)=0.d0
                   call convert_f2t(traces(:,ic,jc,ir,is),trace_time(:,ic,jc),1)
                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)
!gradiant:
          do jc=1,3
             do ic=1,3
                write(name,'("grad",a1,a1,"_",a4)') comp(ic),comp(jc) ,stations_name(ir)
                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
!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)
       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
             call apply_source_filter_gr(trace_time)             
             do jc=1,3
                do ic=1,3
                   write(name,'("grad",a1,a1,"_",a4,"_",i3.3)') comp(ic),comp(jc) ,stations_name(ir),is
                   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
!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)
          enddo
       enddo
    endif
!-------------------------------------------------
  end subroutine get_and_sum_G
!-------------------------------------------------
!-------------------------------------------------
  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
    use global_main_param, only: NBT,dt,nbcomp ,t_delay_sources,half_duration
    implicit none
    complex(DP), dimension(NBF2), intent(inout) :: tracef
    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
    integer :: ic,ir,i,num,if
    real(DP) :: df,w
!      
    if (half_duration(is)  /= 0.0d0) then
       df = 1._DP/(NBF2*dtf)
       do if=2,NBF+1
          w=2._DP*PI*(if-1)*df*half_duration(is)
          tracef(if)=tracef(if)*sin(w)/w
       enddo
       do if=NBF+1,NBF2
          w=-2._DP*PI*(NBF2-if+1)*df*half_duration(is)
          tracef(if)=tracef(if)*sin(w)/w
       enddo      
    endif
    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) &
         =tracef(NBF+1:NBF2)

    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
!-------------------------------------------------

!-------------------------------------------------
  subroutine init_RSkN(channel_,clean)
!-------------------------------------------------
    use def_gparam
    use global_main_param, only: nbcomp,NBE,RA,DEP,VEL,ACC,GRA
    implicit none
    integer, intent(in) :: channel_
    logical, intent(in) :: clean
    integer :: l
    real :: mega
!
    scale =1._DP/sqrt(rhobar*ra**3)
    dscale=1._DP/sqrt(rhobar*ra**3)/ra
!
    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.

    allocate(SkN(-2:2       ,0:Nmax_all,0:Lmax_all,2,NBE))
    SkN=0.0_DP
    print*,'->SkN ...',channel_
    call SkNfun(SkN)
    select case (channel_)
    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)
    case(GRA)
       print*,'->GkN ...'
       allocate(GkN(-1:1,-1:1,0:Nmax_all,0:Lmax_all,2))
       GkN=0.0_DP
       call GkNfun(GkN)
    end select
! Cleaning unecessary arrays:
    if (clean) deallocate(UUs,UUps,VVs,VVps,WWs,WWps,Msrc)   
!-------------------------------------------------
  end subroutine init_RSkN
!------------------------------------------------
!-------------------------------------------------
  subroutine SkNfun(S)
!-------------------------------------------------
    use def_gparam
    use global_main_param, only:RA,NBE,source_force, homogenization
    use corr_modes
    IMPLICIT NONE 
    complex(DP), dimension(-2:,0:,0:,:,:), intent(out) :: S
    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
!

    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)
                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
                xu= b0(l)*u
                xv= b1(l)*v
                S( 0,n,l,1,is)=xu*vir*scale
                S( 1,n,l,1,is)=xv*dcmplx(-vit,vip)*scale
                S(-1,n,l,1,is)=xv*dcmplx(+vit,vip)*scale
             enddo
          enddo
       enddo
       do l=1,LmaxT
          do n=0,NmaxT(l)
             do is=1,NBE
!
                vit=Msrc(2,is)
                vip=Msrc(3,is)
                if (homogenization>=0) then
                   dfw   =WWps(is,n,l)
                   call corrmodesT(is,eigw(n,l,2),l,rr,fw,dfw)
                endif
                x= b1(l)*WWs (is,n,l)
!

                S ( 1,n,l,2,is)=x*dcmplx(+vip,vit)*scale
                S (-1,n,l,2,is)=x*dcmplx(-vip,vit)*scale
             enddo
          enddo
       enddo
!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
       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) 
                if (homogenization>=0) call corrmodesS(is,eigw(n,l,1),l,rr,egu,egv,degu,degv)
                am(:)=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
       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)
                if (homogenization>=0) call corrmodesT(is,eigw(n,l,2),l,rr,fw,dfw)
                am(:) =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
    endif
!--------------------------------------------------
  end subroutine SkNfun
!--------------------------------------------------

!--------------------------------------------------
  subroutine RkNfun(R)
!--------------------------------------------------
    use def_gparam
    use global_main_param, only: NBE,nbcomp,homogenization
    use corr_modes
    IMPLICIT NONE
    complex(DP), dimension(-1:,:,0:,0:,:), intent(out) :: R
    real(DP) :: xu,xv,u,v,vir,vip,vit,x,fw,rr,degu,degv,dfw
    integer :: ic,n,l,is
!
    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)
             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   
             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
             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
!
             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
!--------------------------------------------------
!--------------------------------------------------
  subroutine EkNfun(EE)
!strain
!--------------------------------------------------
    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) :: EE
    real(DP) :: f,x,z,rd
    integer :: n,l
!
    rd=1._DP
    do l=0,LmaxS
       do n=0,NmaxS(l)
!
          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) 
       enddo
    enddo
    do l=1,LmaxT
       do n=0,NmaxT(l)
!
          z    = WWp(n ,l )-WW(n ,l )/rd
!
          EE( :, :,n ,l,2) = 0._DP
!
          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)
       enddo
    enddo
!--------------------------------------------------
  end subroutine EkNfun
!--------------------------------------------------
!--------------------------------------------------
  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
!--------------------------------------------------
!--------------------------------------------------
  subroutine convert_gradient(tr)
!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
!--------------------------------------------------
    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(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
    integer :: is,ir,i,j,if,iN,iNp
    do is=1,NBE
       do ir=1,NBR
          ctmp(:,:,:)=0.d0
          do i=1,3
             do j=1,3
                do if=1,NBF2
                   do iN=1,3
                      do iNp=1,3
                         ctmp(if,i,j)=ctmp(if,i,j)+Cia(i,iN)*Cia(j,iNp)*tr(if,iN,iNp,ir,is)
                      enddo
                   enddo
                enddo
             enddo            
          enddo
          tr(:,:,:,ir,is)=ctmp(:,:,:)
       enddo
    enddo
!--------------------------------------------------
  end subroutine convert_gradient
!--------------------------------------------------
!--------------------------------------------------
  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
!-------------------------------------------------
!-------------------------------------------------
  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)