module down_mod
  use type_mod
  use external_mod
  use util_mod
  implicit none
  type(tdown), public :: down

contains
  !----------------------------------------------------------------
  logical function down_init()
    logical::founds=.false.
    logical::absentD=.true.
    logical::absentH=.true.
    logical::absentR=.true.
    logical::absentM=.true.
    logical::absentW=.true.
    logical::absentD_d=.true.
    logical::absentH_d=.true.
    logical::absentW_d=.true.
    logical::absentN_d=.true.
    logical::badsign=.false. 
    logical::tmp
    integer::ierr
    character(len=1024)::present
#ifdef DEBUG
    character(len=128), parameter :: unit='down_init'
    call in(unit)
#endif

    down_init=F
    !OPERATION SETUP
    call from_param('adj', down%adj, .true.)
    !! set the propagation exponential sign
    down%opsign = -1
    call from_param('badsign',down%badsign,.false.)
    if (down%badsign) then
	down%opsign=down%opsign*-1
   	write(0,*) '   You changed the propagation direction since you FTd your data wrongly.'
    end if
    call from_param('add', down%add, .false.)

    !!Files expecting on the command line
    down%D    ='D'     !! receiver data
    down%D_d  ='D_d'   !! src data
    down%N_d  ='N_d'   !! DD*=source normalization
    down%M    ='M'     !! UU=SRME multiple model
    down%R    ='R'     !! image space
    down%W    ='W'     !! working space sized wavefield coming up
    down%W_d  ='W_d'   !! working space sized wavefield going down
    down%HD   ='HD'    !! trace headers
    down%HD_d ='HD_d'  !! shot headers

    !! if S= is given as parameter, assumes you have slowness
    founds=.false.
    !HANDLE CONVERTED WAVE 
    ierr=getch("S_S","s",present)+getch("S_R","s",present)
    if(2==ierr) then
       down%VS='S_S'   !!       slowness
       down%VR='S_R'   !!       slowness
       down%cwave=.true.
       down%vel=.false.
       founds=.true.
    else if( 1==ierr) then
       write(run%io,*) "Must provide S_S and S_R"
       return
    else
       ierr=getch("S","s",present)
       if(0.ne.ierr) then
          down%VS='S'
          down%VR='S'
          down%cwave=.false.
          down%vel=.false.
          founds=.true.
       end if
    end if
    !! if V= is given as parameter, assumes you have velocity
    if(.not.founds) then
       ierr=getch("V_S","s",present)+getch("V_R","s",present)
       if(2==ierr) then
          down%VS='V_S'   !!       slowness
          down%VR='V_R'   !!       slowness
          down%cwave=.true.
          down%vel=.true.
          founds=.true.
       else if( 1==ierr) then
          write(run%io,*) "Must provide V_S and V_R"
          return
       else
          ierr=getch("V","s",present)
          if(0.ne.ierr) then
             down%VS='V'
             down%VR='V'
             down%cwave=.false.
             down%vel=.true.
             founds=.true.
          end if
       end if
    end if
    if(.not.founds) call seperr("Must input slowness/velocity file.")
    !! get basic depth components from slowness/velocity file
    call hinaxis(down%VR,down%az,3);

    down%wIOd=.false.  !! flag whether I/O on wavefield.
    down%wIOu=.false.  !! flag whether I/O on wavefield.

    call from_param('scale',down%scale,1.)

    call from_param("passive",down%passive,.false.);
    call from_param("datum",down%datum,.false.);

    call from_param("offsets",tmp,.true.)
    call from_param("offset",down%offset,tmp)
    call from_param("roll",tmp,.true.)
    call from_param("rollalong",down%roll,tmp)
    ierr=getch("marine",'s',present)
    if(ierr.ne.0) then
       call from_param("marine",tmp)
       if(tmp) then
          down%offset=.true.
          down%roll=.true.
       else
          down%offset=.false.
          down%roll=.false.
       end if
       write(run%io,*)"Parameter _marine_ no longer exists.  Use _offset_ and _rollalong_"
    end if
    
    

    ierr=getch("D",'s',present)
    if(ierr.ne.0) then
       down%gathers=.true.
       absentD=.false.
    end if
    ierr=getch("R",'s',present)
    if(ierr.ne.0)then
       down%image=.true.
       absentR=.false.
    end if
    ierr=getch("M",'s',present)
    if(ierr.ne.0)then
       down%mult=.true.
       absentM=.false.
    end if
    ierr=getch("N_d",'s',present)
    if(ierr.ne.0)then
       down%normd=.true.
       absentN_d=.false.
    end if
    ierr=getch("D_d",'s',present)
    if(ierr.ne.0) then
       down%sourcefn=.true.
       absentD_d=.false.
    end if
    ierr=getch("W_d",'s',present)
    if(ierr.ne.0) then
       down%WIOd=.true.
       absentW_d=.false.
    end if
    ierr=getch("W",'s',present)
    if(ierr.ne.0) then
       down%WIOu=.true.
       absentW=.false.
    end if
    ierr=getch("HD",'s',present)
    if(0.ne.ierr)then
       down%headD=.true.
       absentD=.false.
    end if
    ierr=getch("HD_d",'s',present)
    if(0.ne.ierr)then
       down%headD_d=.true.
       absentH_d=.false.
    end if
    
    !!R=T, D=T, D_d=F, W=F, W_d=F then passsive=T
    if(down%image.and.down%gathers.and.absentD_d.and.absentW.and.absentW_d) then
       down%passive=.true.
       down%roll=.false.
       down%offset=.false.
    end if
    
    !!R=F, D=T, D_d=F, W=T, W_d=F then datum=T
    if(down%gathers.and.absentD_d.and.down%WIOu.and.absentW_d.and.absentR) then
       down%datum=.true.
       down%offset=.false.
       down%roll=.false.
    end if
    !!R=F, D=F, D_d=T, W=F, W_d=T 
    if(down%sourcefn.and.absentD.and.absentW.and.down%WIOd.and.absentR) then
       down%offset=.false.
       down%roll=.false.
    end if

    if(.not.down%passive) then
       if(down%headD.and.absentH_d) call seperr("Must have HD_d with HD.")
    end if
    if(down%headD.and..not.down%gathers) &
      call seperr("Data headers and no data is absurd.")
    if(down%headD_d.and..not.down%sourcefn) &
      call seperr("Source headers and no source function is absurd.")

    call from_param("image_real",down%image_real,.true.)
    call from_param('image_freqs',down%image_freqs,.false.)
    call from_param("firstz",down%firstz,1) !! first depth slice to write out
    call from_param("firstzW",down%firstzW,1)

    !!depth blocking info
    call from_param("autonzs",down%autonzs,.true.)
    call from_param("nzs",down%nzs,down%az%n)
    down%nzs=min(down%nzs,down%az%n)
    if(down%nzs.ne.down%az%n) down%autonzs=.false.

    if(down%image_freqs) down%autonzs=.true.
    if(down%autonzs) down%nzs=down%az%n
    !! make depth block axis
    down%zblock=down%nzs
    down%bz%o=down%az%o
    down%bz%n=down%az%n / down%nzs
    down%bz%d=down%nzs*down%az%d
    if(down%bz%n==0) call seperr("Number of depth blocks should never be 0! ")

    !! datum through water column in huge DC steps
    call from_param("watercolumn",down%awtr%o,-99999.99)
    down%awtr%j=0;down%awtr%n=0
    if (down%awtr%o.ne.-99999.99) then
       down%watercolumn=.true.
       call from_param("wtr_steps",down%awtr%n,2)
       down%awtr%d=(down%awtr%o-down%az%o)/down%awtr%n
       down%awtr%j=1.5+(down%awtr%o-down%az%o)/down%az%d

       if (down%firstz .lt. down%awtr%j) down%firstz=down%awtr%j

!       if (down%az%o+(down%nzs-1)*down%az%d .lt. down%awtr%o ) &
!         call seperr("Water column must be less than 1 depth block (Increase nzs).")
    end if
    call from_param("tshift",down%tshift,0.)

    !! get basic job size parameters
    if(down%gathers) then
       if(down%headD_d) then
          call hinaxis(down%HD_d,down%evnt,2)
          call auxclose(down%HD_d)
       end if
       if(down%headD) then !!are there headers?
          call hinaxis(down%D,down%aw, 1)
          call hinaxis(down%D,down%trc,2)
          call hinaxis(down%D,down%asx,3)
          call hinaxis(down%D,down%asy,4)
       else
          call hinaxis(down%D,down%amx,1)
          call hinaxis(down%D,down%amy,2)
          call hinaxis(down%D,down%aw, 3)
          call hinaxis(down%D,down%asx,4)
          call hinaxis(down%D,down%asy,5)
       end if
       down%evnt%n=down%asx%n*down%asy%n
    else
       call hinaxis(down%D_d,down%amx,1)
       call hinaxis(down%D_d,down%amy,2)
       call hinaxis(down%D_d,down%aw, 3)
       call hinaxis(down%D_d,down%asx,4)
       call hinaxis(down%D_d,down%asx,5)
       down%evnt%n=down%asx%n*down%asy%n
    end if
    if(down%aw%n==1) call erexit("Only one frequency?? Rediculous")

    !!frequency blocking info
    call from_param("autonws",down%autonws,.true.)
    call from_param("nws",down%nws,down%aw%n)
    down%nws=min(down%nws,down%aw%n)
    if(down%nws.ne.down%aw%n) down%autonws=.false.
    if(down%image_freqs) down%autonws=.true. !necessary??
    if (down%autonws) then
       down%nws=down%aw%n
    end if

    !! make freq. block axis
    down%wblock=down%nws
    down%bw%n=down%aw%n/down%nws
    down%bw%d=down%nws*down%aw%d
    down%bw%o=down%aw%o
    if(down%bw%n==0) call seperr("Number of frequency blocks should never be 0! ")

    !! get width of source function
    if (down%passive.or.down%datum) then
       down%fatx%n=1; down%faty%n=1
    else
       call hinaxis('D_d',down%fatx,1)
       call hinaxis('D_d',down%faty,2)
    end if



    !! number of reference velocities
    call from_param("nref",down%nref,1)

    !! forward scattering migration
    call from_param("fwdscat",down%fwdscat,.false.)

    !! antialias limits
    call from_param("antialias",down%antion,.false.)
    call from_param("pmax_x",down%pmax_x,1000000.)
    call from_param("pmax_y",down%pmax_y,down%pmax_x)
    call from_param("pmin_x",down%pmin_x,-down%pmax_x)
    call from_param("pmin_y",down%pmin_y,-down%pmax_y)
    if(down%pmax_x.ne.1000000.) down%antion=.true.
    if(down%pmax_x.lt.0.or.down%pmin_x.gt.0) call seperr("Must have negative pmin and positive pmax(x).")
    if(down%pmax_y.lt.0.or.down%pmin_y.gt.0) call seperr("Must have negative pmin and positive pmax(y).")
    !! total shots in survey for parallel processing
    call from_param("nx_shots",down%nx_shots,1)
    call from_param("ox_shots",down%ox_shots,999999999.)
    call from_param("ny_shots",down%ny_shots,1)
    call from_param("oy_shots",down%oy_shots,999999999.)

    !! smoothing parameters for deconvolution read in Smooth.F90
    call from_param("decon",down%dcnv,.false.)
    call from_param('eps',down%eps,0.001)
    call from_param("smooth",down%smoothme,.true.)
    if(.not.down%dcnv)down%smoothme=.false.
    call from_param("rectx",down%rect1,2)
    call from_param("recty",down%rect2,down%rect1)

    !! cosine taper - parameters initialized in taper_parameters(D.F90)
    call from_param("taper",down%taperme,.true.)

    call from_param("xyreverse",down%xyr,.false.)
    call from_param("modelplane",down%modelplane,.false.)
    if(down%modelplane)down%xyr=.true.

    if(run%verb) call down_report()

    if(down%adj) then
       run%operation="migration"
    else
       run%operation="modeling"
    end if
    down_init=T
#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
  end function down_init
!----------------------------------------------------------------
  subroutine down_reset2(adj,add,mod_,dat_,dat_d_)
    logical,            intent(in) :: adj,add
    character(len=128), intent(in) :: mod_,dat_,dat_d_

    down%adj  = adj
    down%add  = add
    down%mod_ = mod_
    down%dat_ = dat_
    down%dat_d_ = dat_d_

  end subroutine down_reset2
!----------------------------------------------------------------  
  subroutine down_report()
    call line()
    write(run%io,*)              "        *Paramerter report"
    write(run%io,'(a,i4,a,i4,a,i3)')"         *  nzs=",down%nzs,"          nws=",down%nws,"           nref=",down%nref
    write(run%io,*)              "        *  RollAlong:",down%roll,"    Passive:",down%passive,"      Velocity:",down%vel
    write(run%io,*)              "        *  Offset:",down%offset,"                      Anti-alias:",down%antion
    write(run%io,*)              "        *  Forward scatter:",down%fwdscat,'         Converted wave:',down%cwave
    write(run%io,'(a,i4,a,i4)')"         *  firstZ=",down%firstz,'                      firstzW=',down%firstzw
    write(run%io,*)"        *  add:",down%add,"  Source Tshift=",down%tshift,"adj:",down%adj
    if(down%antion) then
    write(run%io,'(a,f9.7,a,f10.7)')"         *  pmax_x=",down%pmax_x,"            pmin_x=",down%pmin_x
!    if(down%amy%n.gt.1) &
    write(run%io,'(a,f9.7,a,f10.7)')"         *  pmax_y=",down%pmax_y,"            pmin_y=",down%pmin_y
    end if
    if(down%smoothme.and.down%dcnv) then
    write(run%io,'(a,i2,a,i2)')"         *  Deconvolve: T      rectx=",down%rect1,"          recty=",down%rect2
    else
    write(run%io,*)"        *  Deconvolve:",down%dcnv,"                      Smooth:",down%smoothme
    end if
  end subroutine down_report
  !----------------------------------------------------------------
  subroutine process_radians()
    down%aw%o = 2*pi * down%aw%o
    down%aw%d = 2*pi * down%aw%d
    down%bw%o = 2*pi * down%bw%o
    down%bw%d = 2*pi * down%bw%d
  end subroutine process_radians
!----------------------------------------------------------------
!----------------------------------------------------------------
! Function for initializing the run enviro
  logical function run_init()
    character(len=128) :: operation

#ifdef DEBUG
    character(len=128), parameter :: unit='run_mpi_parallel'
    if(run%debug) call in(unit)
#endif DEBUG
  
    call from_param('verb', run%verb, T)
    call from_param('debug',run%debug,F)
    call from_param('depar',run%depar,F)

    if(      run%debug) run%verb =T
    if(.not. run%debug) run%depar=T

    run_init=F
    call from_param('operation',operation,'migration')
    run%operation = operation
    call from_param("io_file",run%file_out,"SCREEN")
    if(run%file_out(1:6)=="SCREEN") then
       run%io=0
    else!if(run%file_out(1:4)=="NONE") then
       run%io=6
       open(run%io,file=run%file_out)
    end if

    if(.not. run_omp_parallel()) then
      write(run%io,*) "trouble initializing omp enviro"
      return
    end if

   run_init=T

   if(run%verb) call run_report()
#ifdef DEBUG
    if(run%debug) call out(unit)
#endif
  end function run_init
!----------------------------------------------------------------
! Subroutine to initialize OMP run environment
! Sets number of omp threads
  logical function run_omp_parallel()
#ifdef DEBUG
    character(len=128), parameter :: unit='omp_run_parallel'
    if(run%debug) call in(unit)
#endif

#ifdef SEP_OMP
!$OMP PARALLEL
    run%o_nth = omp_get_num_threads()
    run%o_ith=1+omp_get_thread_num()
!$OMP END PARALLEL
    run%omp=T
#else
    run%o_nth=1
    run%o_ith=1
    run%omp=F
#endif

   run_omp_parallel=.true.
#ifdef DEBUG
   if(run%debug) call out(unit)
#endif
  end function  run_omp_parallel

!----------------------------------------------------------------
  subroutine run_report()
    write(run%io,*) '        SP *********************************************'
    if(run%verb)then
       write(run%io,'(a,i3)') "         *     Verbose: T              OMP threads=",run%o_nth
    else
       write(run%io,'(a,i3)') "         *     Verbose: F              OMP threads=",run%o_nth
    end if
    write(run%io,*) "        *     Debug:",run%debug,"             Parallel debug:", run%depar
  end subroutine run_report
!----------------------------------------------------------------
end module down_mod
