program SRM_
  !! Stolt Residual Migration
  !! Paul Sava (paul@sep.stanford.edu)
  !! 12/1/00

  !! complex input/output
  !! input  is kz-khx-khy-kmx-kmy
  !! output is kz-khx-khy-kmx-kmy-ro

  use sep
  use utilities
  use omppar_
  use srm_storm_

  integer, parameter:: dbl=selected_int_kind(15)

  type(paxis) :: az,ahx,ahy,amx,amy,ar, aj

  integer, external :: omp_get_num_threads
  integer, external :: omp_get_thread_num

  logical     :: verb, push, scale, reserve
  real        :: ro, eps
  integer     :: stat

  !! z-hx-hy-ngathers
  complex, dimension(:,:,:,:), pointer :: data

  integer :: imx,imy,ir

  real    :: maxmem    !! max memory available  (Mb)
  real    :: gathermem !! memory for one gather (Mb)
  integer :: nth,      ith
  integer :: nblocks,  iblocks
  integer(kind=dbl) :: ngathers, igathers, allgathers

  logical :: depthmove !! remove zero offset resmig 
  logical :: depthfrez !! freeze original depth

  call sep_init()

  call hinaxis(az ,1); az%l ="z";  call houaxis(az ,1)
  call hinaxis(ahx,2); ahx%l="hx"; call houaxis(ahx,2)
  call hinaxis(ahy,3); ahy%l="hy"; call houaxis(ahy,3)
  call hinaxis(amx,4); amx%l="mx"; call houaxis(amx,4)
  call hinaxis(amy,5); amy%l="my"; call houaxis(amy,5)

  call from_param("or",ar%o,1.0)
  call from_param("dr",ar%d,0.01) 
  call from_param("nr",ar%n,0)
  ar%o = ar%o - ar%n * ar%d
  ar%n = 1 + 2*ar%n 
  ar%l="ro"
  call houaxis(ar,6)

  aj%n=1; aj%o=0.; aj%d=1.; aj%l=" "
  call houaxis(aj,7)

  call from_param("maxmem",maxmem,100.)
  call from_param("verb",verb,.false.)
  call from_param("push",push,.true.)
  call from_param("scale",scale,.false.)
  call from_param("reserve",reserve,.false.)
  call from_param("eps",eps,0.5)

  call from_param("depthmove",depthmove,.false.)
  call from_param("depthfrez",depthfrez,.false.)
  call sep_close()

  !----------------------------------------------------------------
  !! find out how many threads the environment variable actually allows
  !$OMP PARALLEL
  nth = omp_get_num_threads()
  !$OMP END PARALLEL

  !! total # of gathers
  allgathers = amy%n * amx%n

  !! memory needed for one gather
  gathermem = (1/1024./1024.) * az%n * ahx%n * ahy%n * 8 
  !! z*hx*hy*(4+8) (Mb)

  !! find the max # of gathers that fit in the allocated memory
  ngathers = maxmem/gathermem

  !! find, if possible, the max value of ngathers 
  !!       which divides allgathers
  !!      and is divided by nth
  call omppar(allgathers,ngathers,nth)

  !! set and freeze the # of threads
  call omp_set_dynamic(.false.)
  call omp_set_num_threads(nth)

  !! # of data blocks
  nblocks = allgathers/ngathers
  !----------------------------------------------------------------  

  if(verb) then
     write(0,*) "    maxmem=",maxmem,'Mb'
     write(0,*) "    effmem=",az%n*ahx%n*ahy%n*ngathers*8/1024./1024.,'Mb'
     write(0,*) "allgathers=",allgathers
     write(0,*) "  ngathers=",ngathers
     write(0,*) "       nth=",nth
     write(0,*) "   nblocks=",nblocks
     write(0,*) " depthmove=",depthmove
     write(0,*) " depthfrez=",depthfrez
  end if

  !----------------------------------------------------------------
  allocate(data(az%n,ahx%n,ahy%n,ngathers), stat=stat)
  if(stat /=0) call seperr("cant allocate data")
  data = (0.0,0.0)

  !----------------------------------------------------------------
  if(reserve) then
     !! reserve enough disk space for the entire output
     if(verb) write(0,*) "reserving disk space"
     do ir=1,ar%n
        do imy=1,amy%n
           do iblocks=1,nblocks
              call sep_write(data)
           end do
        end do
     end do
     call sseek('out',0,0)
     if(verb) write(0,*) "OK"
  end if

  !----------------------------------------------------------------
  call storm_srm_init(nth,az,ahx,ahy,amx,amy,eps,push,scale)

  !----------------------------------------------------------------
  ratio: do ir=1,ar%n
     if(verb) write(0,*) "ratio",ir,'/',ar%n
     call sseek('in',0,0)
     ro = ar%o + (ir-1) * ar%d

     if(depthmove) then
        !----------------------------------------------------------------
        do imy=1,amy%n
           do iblocks=1,nblocks
              if(verb) write(0,*) "block",iblocks,'/',nblocks
              call sep_read(data)

              !$OMP  PARALLEL
              !$OMP  DO &
              !$OMP& SCHEDULE(DYNAMIC) &
              !$OMP& PRIVATE(ith,imx,igathers)
              do igathers=1,ngathers
                 imx = igathers + (iblocks-1)*ngathers
                 ith = omp_get_thread_num()+1

                 call storm_srm_compute(ith,igathers,data,ro,az,ahx,ahy,imx,imy,depthfrez)

              end do
              !$OMP END DO
              !$OMP END PARALLEL

              call sep_write(data)
           end do
        end do
        !----------------------------------------------------------------
     else
        !----------------------------------------------------------------
        do imy=1,amy%n
           do iblocks=1,nblocks
              if(verb) write(0,*) "block",iblocks,'/',nblocks
              call sep_read(data)

              !$OMP  PARALLEL
              !$OMP  DO &
              !$OMP& SCHEDULE(DYNAMIC) &
              !$OMP& PRIVATE(ith,imx,igathers)
              do igathers=1,ngathers
                 imx = igathers + (iblocks-1)*ngathers
                 ith = omp_get_thread_num()+1

                 call storm_srm_compute(ith,igathers,data,ro,az,ahx,ahy,imx,imy)
                 call storm_srm_nomove (ith,igathers,data,ro,az,ahx,ahy,imx,imy)

              end do
              !$OMP END DO
              !$OMP END PARALLEL

              call sep_write(data)
           end do
        end do
        !----------------------------------------------------------------
     end if
  end do ratio
  !----------------------------------------------------------------
  deallocate(data, stat=stat)
  if(stat /=0) call seperr("cant deallocate data")

  call exit(0)
end program SRM_





















