!<
!semb_drive
!
!Usage:
!semb_drive.x  <in.H >out.H
!
!Input Parameters:
!
!
!Output Parameters:
!
!
!
!Description:
!
!
!
!>
!-------------------------------------------------
!
!Author: Robert Clapp, ESMB 463, 7230253
!
!Date Created:Thu Jan 30 12:41:38 2003
!
!Purpose:
!
program semb_drive 
  use sep3d_struct_mod
  use sep
  use scan_mod
  use verb_mod
  implicit none
  integer n1,n2,i1,i2,n3,n4,ngam,igam,i,ntot,ij,i3,iax,i4
  real, allocatable,dimension(:,:) :: semb
  real, allocatable,dimension(:,:) :: gather
  type(sep3d) ::  mig_s,semb_s
  real :: ogam,dgam,pos(2)
  logical :: verb
  character(len=128) :: intag,outtag
  call initpar()
  call doc(SOURCE)
  call sep_begin_prog()
  outtag="outtag"
  intag="intag"
  call init_sep3d(intag,mig_s,"INPUT")
  call set_alternate_putch("outtag")

  call from_param('ogam',ogam,.9)
  call from_param('dgam',dgam,.01)
  call from_param('ngam',ngam,21)
  call from_param("verb",verb,.false.)
  call init_sep3d(mig_s,semb_s,"OUTPUT")
  if(sep3d_ndims(mig_s)==4) call sep3d_reshape(mig_s,(/1,2,3,4,4/))
  semb_s%o(2)=ogam
  semb_s%d(2)=dgam
  semb_s%n(2)=ngam
  call sep3d_set_sep3d(semb_s)
  call sep3d_write_description(outtag,semb_s)
  if (sep3d_ndims(mig_s).ne.5 .and.sep3d_ndims(mig_s).ne.4) then
    call seperr("Only understand input as 5 or 4 dimenson")
  end if

  if (.not. init_verbose()) then
    call seperr("trouble initializing verbosity")
  end if

  if (.not. init_print_progress(0,mig_s%n(3)*mig_s%n(4))) then
    call seperr("trouble initializing progress")
  end if
  if (.not.init_scan(mig_s)) then
    write(0,*) "trouble initializing scan"
    call seperr("")
  end if
  allocate(semb(semb_s%n(1),ngam),gather(mig_s%n(1),mig_s%n(2)))
  i=0

  do i4=1,mig_s%n(5)
    pos(2)=mig_s%o(5)+mig_s%d(5)*(i4-1)
    do i3=1,mig_s%n(4)
      call sep3d_grab_headers(intag,mig_s,i,nwind=(/mig_s%n(2),1,1,1/)&
        &,fwind=(/0,0,i3-1,i4-1/))
      if (i>0) then
        call print_progress("begining pct ",1)
        pos(1)=mig_s%o(4)+mig_s%d(4)*(i3-1)
        if (.not. sep3d_read_data(intag,mig_s,gather)) then
          call seperr("trouble reading data")
        end if
        do igam=1,ngam 
          if (.not. compute_semb(pos,ogam+dgam*(igam-1),gather,semb(:&
            &,igam))) then
            write(0,*) i3,i4
            call seperr("trouble computing semb")
          end if
        end do 
        if (.not. sep3d_write_data(outtag,semb_s,semb,nwind=(/semb_s%n&
          &(1),semb_s%n(2),1,1,1/),fwind=(/0,0,0,i3-1,i4-1/))) then
        end if
      end if
    end do
  end do 
  call sep_end_prog()
end program  
