receivers.f90

Spectral Element Method for wave propagation and rupture dynamics.

  integer :: n,ip,ipoint,irec


  distmax = 0.d0

  if (rec%AtNode) then

    if (info) write(iout,200)
    iglob_tmp = 0
    irec = 0
    do n=1,rec%nx
      call SE_find_nearest_node(rec%coord(:,n),grid,ipoint,distance = dist)
     ! do not keep it if already in list
      if ( irec > 1 .and. any(iglob_tmp(:irec) == ipoint)) cycle
      irec = irec+1
      iglob_tmp(irec) = ipoint
      if (info) then
        write(iout,150) irec,rec%coord(:,n),grid%coord(:,ipoint),dist
        distmax = max(dist,distmax)
      endif
    enddo
    if (irec < rec%nx) then 
      rec%nx = irec
      deallocate(rec%coord)
      allocate(rec%coord(NDIME,rec%nx))
    endif
    call storearray('rec%coord',size(rec%coord),idouble)
    allocate(rec%iglob(rec%nx))
    call storearray('rec%iglob',size(rec%iglob),isngl)
    rec%iglob = iglob_tmp(:rec%nx)
    rec%coord(:,:) = grid%coord(:,rec%iglob)

  else

    if (info) write(iout,300)
    allocate(rec%interp(grid%ngll*grid%ngll,rec%nx))
    allocate(rec%einterp(rec%nx))
    do n=1,rec%nx
      call SE_find_point(rec%coord(:,n),grid,rec%einterp(n),xi,eta,new_coord)
      call SE_init_interpol(xi,eta,rec%interp(:,n),grid)
      dist = sqrt( (rec%coord(1,n)-new_coord(1))**2 + (rec%coord(2,n)-new_coord(2))**2 )
      if (info) then
        write(iout,150) n,rec%coord(:,n),new_coord,dist
        distmax = max(dist,distmax)
      endif
      rec%coord(:,n)=new_coord
    enddo

  endif

  if (info) write(iout,160) distmax

 150   format(2x,i7,5(1x,EN12.3))
 160   format(/2x,'Maximum distance between asked and real =',EN12.3)
 200  format(//1x,'R e c e i v e r s'/1x,17('=')// &
  ' Receivers have been relocated to the nearest GLL node'// &
  ' Receiver  x-requested  z-requested   x-obtained   z-obtained   distance'/)

 300  format(//1x,'R e c e i v e r s'/1x,17('=')// &
  ' Receiver            x            z   element           xi          eta'/)
 310   format(2x,i7,1x,EN12.3,1x,EN12.3,1x,i9,1x,EN12.3,EN12.3)

  end subroutine REC_posit


!=====================================================================
!
!!  Store the seismograms
  subroutine REC_store(rec,it,grid)

  use spec_grid, only : sem_grid_type

  type(rec_type), intent(inout) :: rec
  integer, intent(in) :: it
  type(sem_grid_type), intent(in) :: grid

  integer :: itsis,n,i,k,j,iglob
  double precision, allocatable :: vloc(:,:)

  if ( mod(it,rec%isamp) /= 0 ) return

  itsis = it/rec%isamp +1
  if (itsis > rec%nt) call IO_abort('receivers.REC_store: storage is full')

  if (rec%AtNode) then
    rec%sis(itsis,:,:) = rec%field(rec%iglob,:)
  else
    allocate(vloc(grid%ngll*grid%ngll,size(rec%field,2)))
    do n=1,rec%nx
      k=1
      do j=1,grid%ngll
      do i=1,grid%ngll
        iglob = grid%ibool(i,j,rec%einterp(n))
        vloc(k,:) = rec%field(iglob,:)
        k=k+1
      enddo
      enddo
      rec%sis(itsis,n,:) = matmul(rec%interp(:,n),vloc)
    enddo
    deallocate(vloc)
  endif
    
  end subroutine REC_store



!=====================================================================
!
!!  Export the seismograms
  subroutine REC_write(rec,iout)

  use stdio, only: IO_new_unit
  double precision, parameter :: factorxsu = 3.5d0

  type(rec_type), intent(inout) :: rec
  integer, intent(in) :: iout

  double precision :: xval(rec%nx),xref,zref
  integer :: iol,ounit
  character(30) :: ylabel
  character, parameter :: posvars(3) = (/ 'X','Z','D' /)
  character :: posvar

!---- binary data -------------------------------------------------------
  write(iout,*) 'Storing sismos data (SEP format) ...'
  INQUIRE( IOLENGTH=iol ) rec%sis(:,:,1)
  ounit = IO_new_unit()

  if ( size(rec%sis,3)==1) then

    open(ounit,file='Uy_sem2d.dat',status='replace',access='direct',recl=iol)
    write(ounit,rec=1) rec%sis(:,:,1)
    close(ounit)

  else

    open(ounit,file='Ux_sem2d.dat',status='replace',access='direct',recl=iol)
    write(ounit,rec=1) rec%sis(:,:,1)
    close(ounit)
  
    open(ounit,file='Uz_sem2d.dat',status='replace',access='direct',recl=iol)
    write(ounit,rec=1) rec%sis(:,:,2)
    close(ounit)

  endif

!=== Seismic Unix scripts: =================================================

  posvar = posvars(rec%irepr)

  select case(rec%irepr)
    case(1,2) ! recepteurs suivant coordonnee X ou Z
      xval = rec%coord(rec%irepr,:)
    case(3) ! recepteurs en distance
      xref = rec%coord(1,1)
      zref = rec%coord(2,1)
      xval = sqrt((rec%coord(1,:) - xref)**2 + &
              (rec%coord(2,:) - zref)**2)
  end select

  select case(rec%SeisField)
    case('D'); ylabel='Displacement (m)'
    case('V'); ylabel='Velocity (m/s)'
    case('A'); ylabel='Acceleration (m/s^2)'
  end select

! station "coordinates" for plots
  open(ounit,file='x2_sem2d.tab',status='unknown')
  write(ounit,'(1000(f0.2:","))') xval
               !this number > max nb stations ever
  close(ounit)

!-- Xwindow -------------------------------------------------------------------
  open(ounit,file='Xline_sem2d.csh',status='unknown')
  write(ounit,110)
  write(ounit,*) 'set x2 = `cat x2_sem2d.tab`'
  write(ounit,100) 'xwigb',factorxsu,rec%nt,rec%tsamp,posvar,rec%nx,&
    ' title="Ux component" < Ux_sem2d.dat'
  write(ounit,100) 'xwigb',factorxsu,rec%nt,rec%tsamp,posvar,rec%nx,&
    ' title="Uz component" < Uz_sem2d.dat'
  close(ounit)

!-- PostScript ---------------------------------------------------------------
  open(ounit,file='PSline_sem2d.csh',status='unknown')
  write(ounit,110)
  write(ounit,*) 'set x2 = `cat x2_sem2d.tab`'
  write(ounit,100) 'pswigp',factorxsu,rec%nt,rec%tsamp,posvar,rec%nx,&
    ' title="Ux component" < Ux_sem2d.dat >! UxPoly_sem2d.ps'
  write(ounit,100) 'pswigp',factorxsu,rec%nt,rec%tsamp,posvar,rec%nx,&
    ' title="Uz component" < Uz_sem2d.dat >! UzPoly_sem2d.ps'
  close(ounit)

!-- one trace for Xwindow --------------------------------------------

  open(ounit,file='Xtrace_sem2d.csh',status='unknown')
  write(ounit,110)
  write(ounit,*) 'set nt = ',rec%nt
  write(ounit,*) 'set dt = ',real(rec%tsamp)
  write(ounit,*) '@ trace=10000'
  write(ounit,*) 'while ($trace > -1)'
  write(ounit,*) 'echo Trace number ?'
  write(ounit,*) 'set rep=$<'
  write(ounit,*) '@ trace = $rep'
  write(ounit,*) 'echo Trace asked : $trace'
  write(ounit,*) '# traces commencent a zero dans format SEP'
  write(ounit,*) '@ septrace = $trace - 1'
  write(ounit,*) 'subset < Ux_sem2d.dat n1=$nt', &
   ' if2s=$septrace n2s=1 | xgraph -geometry 1085x272', &
   ' label1="Time (s)" label2="',trim(ylabel),'" n=$nt style=normal d1=$dt', &
   ' title="Ux component trace "$trace &'
  write(ounit,*) 'subset < Uz_sem2d.dat n1=$nt', &
   ' if2s=$septrace n2s=1 | xgraph -geometry 1085x272', &
   ' label1="Time (s)" label2="',trim(ylabel),'" n=$nt style=normal d1=$dt', &
   ' title="Uz component trace "$trace &'
  write(ounit,*) 'end'
  close(ounit)

!-- one trace for PostScript --------------------------------------------

  open(ounit,file='PStrace_sem2d.csh',status='unknown')
  write(ounit,110)
  write(ounit,*) 'set nt = ',rec%nt
  write(ounit,*) 'set dt = ',real(rec%tsamp)
  write(ounit,*) '@ trace=10000'
  write(ounit,*) 'while ($trace > -1)'
  write(ounit,*) 'echo Trace number ?'
  write(ounit,*) 'set rep=$<'
  write(ounit,*) '@ trace = $rep'
  write(ounit,*) 'echo Trace asked : $trace'
  write(ounit,*) '# traces commencent a zero dans format SEP'
  write(ounit,*) '@ septrace = $trace - 1'
  write(ounit,*) 'rm -f UxTrace{$trace}_sem2d.ps UzTrace{$trace}_sem2d.ps'
  write(ounit,*) 'subset < Ux_sem2d.dat n1=$nt', &
   '" if2s=$septrace n2s=1 | psgraph label1="Time (s)" label2="',trim(ylabel),&
   ' n=$nt style=normal d1=$dt hbox=4.0 wbox=6.0', &
   ' title="Ux component trace "$trace > UxTrace{$trace}_sem2d.ps'
  write(ounit,*) 'subset < Uz_sem2d.dat n1=$nt', &
   '" if2s=$septrace n2s=1 | psgraph label1="Time (s)" label2="',trim(ylabel),&
   ' n=$nt style=normal d1=$dt hbox=4.0 wbox=6.0', &
   ' title="Uz component trace "$trace > UzTrace{$trace}_sem2d.ps'
  write(ounit,*) 'end'
  close(ounit)

 100  format(A,' xcur=',F0.2,' n1=',I0,' d1=',F0.8,' label1="Time (s)" label2="',A, &
             ' (m)" n2=',I0,' x2=$x2',A )
 110  format('#!/bin/csh -f')
  
  end subroutine REC_write


!=====================================================================
!
  subroutine REC_inquire(rec,coord,isamp)

  type(rec_type), intent(in) :: rec
  double precision, pointer, optional :: coord(:,:)
  integer, optional :: isamp

  if (present(coord))  coord => rec%coord
  if (present(isamp))  isamp = rec%isamp

  end subroutine REC_inquire

end module receivers