receivers.f90

Spectral Element Method for wave propagation and rupture dynamics.

! SEM2DPACK version 2.2.11 -- A Spectral Element Method for 2D wave propagation and fracture dynamics,
!                             with emphasis on computational seismology and earthquake source dynamics.
! 
! Copyright (C) 2003-2007 Jean-Paul Ampuero
! All Rights Reserved
! 
! Jean-Paul Ampuero
! 
! ETH Zurich (Swiss Federal Institute of Technology)
! Institute of Geophysics
! Seismology and Geodynamics Group
! ETH H鰊ggerberg HPP O 13.1
! CH-8093 Z黵ich
! Switzerland
! 
! [e-mail]
! [phone] (office)
! [phone] (fax)
! 
! http://www.sg.geophys.ethz.ch/geodynamics/ampuero/
! 
! 
! This software is freely available for scientific research purposes. 
! If you use this software in writing scientific papers include proper 
! attributions to its author, Jean-Paul Ampuero.
! 
! This program is free software; you can redistribute it and/or
! modify it under the terms of the GNU General Public License
! as published by the Free Software Foundation; either version 2
! of the License, or (at your option) any later version.
! 
! This program is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU General Public License for more details.
! 
! You should have received a copy of the GNU General Public License
! along with this program; if not, write to the Free Software
! Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
! 
module receivers

  use stdio, only: IO_abort
  use constants, only : NDIME

  implicit none
  private

  type rec_type
    private
    double precision, pointer :: coord(:,:),field(:,:)
    real, pointer :: sis(:,:,:)
    integer, pointer :: iglob(:)
    double precision :: tsamp
    integer :: nx,nt,isamp,irepr
    logical :: AtNode
    double precision, pointer :: interp(:,:)
    integer, pointer :: einterp(:)
    character :: SeisField
  end type rec_type

  public :: rec_type,REC_read,REC_init,REC_store,REC_write,REC_inquire

contains

!=====================================================================
!
! BEGIN INPUT BLOCK
!
! NAME   : REC_LINE
! PURPOSE: Defines a line of receivers
! SYNTAX : &REC_LINE number,isamp,field,first,last,file,AtNode,irepr /
!
! ARG: number   [int] [0] Number of stations in the line
! ARG: isamp    [int] [1] Sampling stride (in number of timesteps). Note that
!                for stability reasons the timestep can be very small.
! ARG: field    [char] ['V'] The field in the seismogram:
!                               'D'     displacement
!                               'V'     velocity
!                               'A'     acceleration
! ARG: first    [dble(2)] Receivers can be located along a line,
!                this is the position (x,z) of the first receiver
! ARG: last     [dble(2)] Position (x,z) of the last receiver,
!                other receivers will be located with regular spacing
!                between First and Last.
! ARG: file     [name] ['none'] Station positions can instead be read 
!                from an ASCII file, with 2 columns per line:
!                (1) X position (in m) and
!                (2) Z position (in m)
! ARG: AtNode	[log] [T] Relocate the stations at the nearest GLL node
! ARG: irepr    [char] ['D'] Abscissa for the seismic multitrace plot:
!                               'X' Horizontal position
!                               'Z' Depth
!                               'D' Distance to the first station
!
! NOTE   : to locate receivers at the free surface set their vertical position 
!          above the free surface and AtNode=T
!
! END INPUT BLOCK

  subroutine REC_read(rec,iin)

  use echo, only : echo_input,iout
  use stdio, only : IO_new_unit, IO_file_length

  type(rec_type), pointer :: rec
  integer, intent(in) :: iin

  double precision :: first(NDIME),last(NDIME),init_double
  integer :: i,irec,iin2,isamp,number
  logical :: AtNode
  character(50) :: file
  character :: field,irepr
  
  NAMELIST / REC_LINE / number,isamp,field,first,last,AtNode,file,irepr

  init_double = huge(init_double) ! set to an unlikely value

  number     = 0
  isamp      = 1
  field      = 'V'
  first      = init_double
  last       = init_double
  AtNode     = .true.
  file       = 'none'
  irepr      = 'D'

  rewind(iin)
  read(iin,REC_LINE,END = 200) 
  allocate(rec)

  if (number <= 0) call IO_abort('REC_read: "number" is null or missing')

  rec%isamp  = isamp
  rec%SeisField = field
  rec%AtNode = AtNode

  if ( field /='D' .and. field /='V ' .and. field /='A') &
    call IO_abort('REC_read: parameter field has wrong value [D,V,A]')

  select case(irepr)
    case('X'); rec%irepr  = 1
    case('Z'); rec%irepr  = 2
    case('D'); rec%irepr  = 3
    case default; call IO_abort('REC_read: unknown "irepr" [X,Z,D]')
  end select
  
  if (echo_input) write(iout,100) number,isamp,field,irepr

  if (file == 'none') then

    if ( any(first == init_double) ) &
      call IO_abort('REC_read: you must set "first" station coordinates')
    if ( any(last == init_double) ) &
      call IO_abort('REC_read: you must set "last" station coordinates')
    rec%nx   = number
    allocate(rec%coord(NDIME,rec%nx))
    if (number>1) then
      do i = 1,rec%nx
        rec%coord(:,i) = first + (i-1)/dble(rec%nx-1) *(last-first)
      enddo
    else
      rec%coord(:,1) = first    
    endif

  else
    rec%nx = IO_file_length(file)
    allocate(rec%coord(NDIME,rec%nx))
    iin2 = IO_new_unit()
    open(iin2,file=file,status='old')
    do i=1,rec%nx
      read(iin2 ,*) rec%coord(:,i) 
    enddo
    close(iin2)
  endif

  return

  100 format(//1x,'R e c e i v e r s', &
  /1x,17('='),//5x, &
  'Number of receivers . . . . . . . . . . . . (number) = ',I0/5x, &
  'Subsampling for seismograms recording . . . .(isamp) = ',I0/5x, &
  'Field recorded. . . . . . . . . . . . . . . .(field) = ',A/5x, &
  'Axis of the seismogram plot . . . . . . . . .(irepr) = ',A)

  200 return

  end subroutine REC_read

!=====================================================================
!
!!  Initialize
  subroutine REC_init(rec,grid,time,fields)

    use fields_class, only : fields_type
    use time_evol, only : timescheme_type
    use spec_grid, only : sem_grid_type
    use memory_info
    use echo, only: iout,info=>echo_init
    use stdio, only: IO_new_unit

    type(rec_type), intent(inout) :: rec
    type(sem_grid_type), intent(in) :: grid
    type(fields_type), intent(in) :: fields
    type(timescheme_type), intent(in) :: time

    integer :: ounit,irec

    call REC_posit(rec,grid)

    rec%tsamp = time%dt*rec%isamp

    rec%nt = time%nt/rec%isamp +1
    if (rec%nt == 0) call IO_abort('REC_init: zero samples')
      
    allocate(rec%sis(rec%nt,rec%nx,fields%ndof))
    rec%sis = 0.
    call storearray('rec%sis',size(rec%sis),idouble)

    select case(rec%SeisField)
      case('D'); rec%field => fields%displ
      case('V'); rec%field => fields%veloc
      case('A'); rec%field => fields%accel
      case default
        call IO_abort('REC_init: Wrong field code for seismograms display')
    end select

  if (info) then
    write(iout,*)
    write(iout,100) 'Sampling rate (Hz)        = ',1d0/rec%tsamp
    write(iout,100) 'Sampling timestep (secs)  = ',rec%tsamp
    write(iout,200) 'Total number of samples   = ',rec%nt
    write(iout,200) 'Number of receivers       = ',rec%nx
    write(iout,*)
  endif

! Save headers
  ounit = IO_new_unit()
  open(ounit,file='SeisHeader_sem2d.hdr',status='replace')
  write(ounit,*) 'DT NSAMP NSTA'
  write(ounit,*) real(rec%tsamp),rec%nt,rec%nx
  write(ounit,*) 'XSTA ZSTA'
  do irec=1,rec%nx
    write(ounit,*) rec%coord(:,irec)
  enddo
  close(ounit)

100 format(2X,A,EN12.3)
200 format(2X,A,I0)

  end subroutine REC_init

!=====================================================================
!
!!  Real receivers position
  subroutine REC_posit(rec,grid)

  use echo, only: iout,info=>echo_init
  use spec_grid, only : sem_grid_type,SE_find_nearest_node,SE_init_interpol,SE_find_point
  use memory_info

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

  integer :: iglob_tmp(rec%nx)
  double precision :: distmax,xs,zs,xp,zp,dist, xi,eta, new_coord(NDIME)