IMPLEMENTATION

The following Fortran 90 module implements an object-oriented design of the generic NMO program.

module nnnmo
  use interp

logical,                          private :: pul
  integer,                          private :: nt
  real,                             private :: ot, dt
  real, dimension (:), allocatable, private :: t, z, w

contains

subroutine nnnmo_init (pull, t0, z0, dz, n)
    logical, intent (in) :: pull
    real,    intent (in) :: t0, z0, dz
    integer, intent (in) :: n

integer              :: i

pul = pull ; ot = t0 ; dt = dz; nt = n 
    if (.not. allocated (z)) allocate (z (nt), t (nt), w (nt))
    if (pul) then                                     ! create a regular trace
       z = (/ (z0 + dt * (i - 1), i = 1, nt) /) ; ot = t0  
    else
       t = (/ (t0 + dt * (i - 1), i = 1, nt) /) ; ot = z0
    end if
  end subroutine nnnmo_init

subroutine nnnmo_step (type, time, ampl, x)
    integer,             intent (in) :: type
    real, dimension (:), intent (in) :: x
    interface
       function time (pull, x, z, t) result (stat)
         integer             :: stat
         logical             :: pull
         real, dimension (:) :: x, z, t
       end function time
       function ampl (pull, x, z, t, weight) result (stat)
         integer             :: stat
         logical             :: pull
         real, dimension (:) :: x, z, t, weight
       end function ampl
    end interface

integer                         :: stat1,stat2

stat1 = time (pul, x, z, t)                      ! time transform
    stat2 = ampl (pul, x, z, t, w)                   ! amplitude

if (pul) then                                    ! interpolation setup
       call interp_init (type, t, ot, dt, nt, w) 
    else
       call interp_init (type, z, ot, dt, nt, w)
    end if
 end subroutine nnnmo_step

function nnnmo_op (adj, add, mod, dat) result (stat)
    integer              :: stat
    logical, intent (in) :: adj, add
    real, dimension (:)  :: mod, dat

integer              :: stat0

if (pul) then                                   ! interpolate to a trace
       stat0 = interp_op (.not. adj, add, dat, mod)
    else
       stat0 = interp_op (      adj, add, mod, dat)
    end if

stat = 0
  end function nnnmo_op

subroutine nnnmo_close ()
    call interp_close ()
    deallocate (z, t, w)
  end subroutine nnnmo_close

end module nnnmo

The module contains the initialization routine nnnmo_init, which receives the time grid information and saves it to module variables. It also translates the regular grid into a one-dimensional array. The next subroutine, nnnmo_step, computes the NMO time stretch and amplitude scaling. Its arguments are a logical variable, which defines pull or push operation, a real array x, which may contain the trace header information, and two functions to compute the time and amplitude transformation. The user of the nnnmo module is free to specify these functions according to any particular application provided that they comply to the generic interface. Typical examples of commonly used functions are collected in module nmo.

module nmo_mod
contains

function linear_time  (pull, x, z, t) result (stat) 
    integer             :: stat
    logical             :: pull
    real, dimension (:) :: x, z, t

if (pull) then
       t = z + sum (x)          
    else
       z = t - sum (x)
    end if

stat = 0
  end function linear_time

function hyper_time (pull, x, z, t) result (stat)
    integer             :: stat
    logical             :: pull
    real, dimension (:) :: x, z, t

if (pull) then
       t = sqrt (z*z + dot_product (x,x))
    else
       z = t*t - dot_product (x,x)
       where (z > 0.) z = sqrt (z)
    end if

stat = 0
  end function hyper_time

function vofz_hyper_time (pull, x, z, t) result (stat)
    integer             :: stat
    logical             :: pull
    real, dimension (:) :: x, z, t

if (pull) then
       t = sqrt (z*z + x*x)
    else
       z = t*t - x*x 
       where (z > 0.) z = sqrt (z)
    end if

stat = 0
  end function vofz_hyper_time

function unit_ampl (pull, x, z, t, w) result (stat)
    integer             :: stat
    logical             :: pull
    real, dimension (:) :: x, z, t, w

w = 1. ; stat = 0
  end function unit_ampl

function cos2D_ampl (pull, x, z, t, w) result (stat)
    integer             :: stat
    logical             :: pull
    real, dimension (:) :: x, z, t, w

w = 0.
    if (pull) then
       where (t > 0.) w = (z / t) / sqrt (t)
    else
       where (t > 0.) w =     1.  / sqrt (t)
     end if

stat = 0
  end function cos2D_ampl

end module nmo_mod

After initialization calls to nnnmo_init and nnnmo_step, the actual operator nnnmo_op simply calls an interpolation operator to complete the job. When the pull parameter is set to .true., the forward operator is ``push'', and the adjoint is ``pull''. In this case, the time interpolation is carried out in the data space. When pull is .false., the time interpolation is done on the model trace. The interpolation is performed by module interp, which calls either bin1 or lint1 Fomel and Claerbout (1996) depending on the user-specified parameter.

module interp

use bin1
  use lint1

implicit none

integer, private :: tp

contains

subroutine interp_init (type, coord, o1, d1, n1, weight)
    real, dimension (:), intent (in) :: coord, weight
    real,                intent (in) :: o1, d1
    integer,             intent (in) :: type, n1
    optional                         :: weight

tp = type
    select case (tp)
    case (1)
       call  bin1_init  (coord, o1, d1, n1, weight)
    case default
       call  lint1_init (coord, o1, d1, n1, weight)
    end select
  end subroutine interp_init

function interp_op (adj, add, mod, ord) result (stat)
    integer              :: stat
    logical, intent (in) :: adj, add
    real, dimension (:)  :: mod, ord

integer              :: stat0

select case (tp)
    case (1)
       stat0 = bin1_op  (adj, add, mod, ord)
    case default
       stat0 = lint1_op (adj, add, mod, ord)
    end select

stat=0
  end function interp_op

subroutine interp_close ()
    select case (tp)
    case (1)
       call bin1_close  ()
    case default
       call lint1_close ()
    end select
  end subroutine interp_close

end module interp