program IQM
  !! Compute image quality maps
  !! Paul Sava (paul@sep.stanford.edu)
  !! 10/14/02

  !! input:  z-ax-mx-ro (ADCIG)
  !! output: z-   mx-ro (image quality map)

  use sep
  use iqtransf_
  implicit none

  integer :: method
  logical :: verb

  real, dimension(:,:,:), allocatable :: adcig
  real, dimension(:,:),   allocatable :: iqmap

  integer :: nz,nmx,nax, nro, iz,imx,iro
  real    ::    omx,oax, oro
  real    ::    dmx,dax, dro
  integer :: xwin,zwin

  call sep_init()

  call from_param("method",method,0)
  call from_param("verb",verb,.false.)

  call from_history("n1",nz )

  call from_history("n2",nax)
  call from_history("o2",oax,0.)
  call from_history("d2",dax,1.)

  call from_history("n3",nmx)
  call from_history("o3",omx,0.)
  call from_history("d3",dmx,1.)

  call from_history("n4",nro,1)
  call from_history("o4",oro,0.)
  call from_history("d4",dro,1.)

  call   to_history("n1",nz )

  call   to_history("n2",nmx)
  call   to_history("o2",omx)
  call   to_history("d2",dmx)
  call   to_history("label2","mx")

  call   to_history("n3",nro)
  call   to_history("o3",oro)
  call   to_history("d3",dro)
  call   to_history("label3","ro")

  call   to_history("n4",1 )
  call   to_history("o4",0.)
  call   to_history("d4",1.)

  call from_param("xwin",xwin,0)
  call from_param("zwin",zwin,0)

  if(nmx==1) xwin=0

  call iqinit(nmx,nax,nz,xwin,zwin,verb)

  allocate(adcig(nz,nax,nmx)); adcig=0
  allocate(iqmap(nz,    nmx)); iqmap=0.

  write(0,*) "Method=",method
  do iro=1,nro
     if(verb) write (0,*) "ro=",iro
     call sep_read(adcig)

     select case(method)
     case(3)
        if(verb) write(0,*) "semblance"
        call to_history("title","semblance")
        call se(adcig,iqmap)
     case(2)
        if(verb) write(0,*) "differential semblance"
        call to_history("title","differential semblance")
        call ds(adcig,iqmap)
     case(1)
        if(verb) write(0,*) "stack intensity"
        call to_history("title","stack intensity")
        call si(adcig,iqmap)
     case default 
        if(verb) write(0,*) "stack power"
        call to_history("title","stack power")
        call sp(adcig,iqmap)
     end select

     do iz=1,zwin
        iqmap(iz,:) = iqmap(zwin+1,:)
     end do
     do iz=nz-zwin+1,nz
        iqmap(iz,:) = iqmap(nz-zwin,:) 
     end do

     do imx=1,xwin
        iqmap(:,imx) = iqmap(:,xwin+1)
     end do
     do imx=nmx-xwin+1,nmx
        iqmap(:,imx) = iqmap(:,nmx-xwin) 
     end do

     call sep_write(iqmap)
  end do

  call sep_close()
  deallocate(adcig,iqmap)

  call exit(0)
end program IQM













