module wei_slow_mod
  use wei_wavespace_type_mod
  use wei_step_mod
  implicit none

contains


  subroutine init_one_slow(space,tag,tag_ref)
    character(len=*)        :: tag,tag_ref
    type(wspace)           :: space
    type(sep3d)             :: sref

#ifdef DEBUG
    write(0,*) "IN INIT ONE SLOW"
#endif

    allocate(space%slow1)
    call init_slow(space%slow1,tag,tag_ref,space)
    space%src%slow=>space%slow1
    space%rec%slow=>space%slow1
    space%src%x=>space%slow1%xmap_s
    space%src%y=>space%slow1%ymap_s
    space%rec%x=>space%slow1%xmap_g
    space%rec%y=>space%slow1%ymap_g

#ifdef DEBUG
    write(0,*) "OUT INIT ONE SLOW"
#endif

  end subroutine

  subroutine init_two_slow(space,tag1,tag1_ref,tag2,tag2_ref)
    character(len=*)        :: tag1,tag2,tag1_ref,tag2_ref
    type(wspace)           :: space

#ifdef DEBUG
    write(0,*) "IN INIT TWO SLOW"
#endif

    allocate(space%slow1)
    allocate(space%slow2)
    call init_slow(space%slow1,tag1,tag1_ref,space)
    call init_slow(space%slow2,tag2,tag2_ref,space)
    space%src%slow=space%slow1
    space%rec%slow=space%slow2
    space%src%x=>space%slow1%xmap_s
    space%src%y=>space%slow1%ymap_s
    space%rec%x=>space%slow2%xmap_g
    space%rec%y=>space%slow2%ymap_g

#ifdef DEBUG
    write(0,*) "OUT INIT TWO SLOW"
#endif

  end subroutine


  subroutine init_slow(struct,tag,tag_ref,space)
    type(wslow)         ::   struct
    type(wspace)       ::   space
    character(len=*)    ::   tag,tag_ref
    type(sep3d)         ::   ref
    integer             ::   i1,i2,i3,ic
    real                ::   pos,pos1,pos2
    real,allocatable    ::   temp(:,:,:)
#ifdef DEBUG
    write(0,*) "IN INIT SLOW"
#endif

    struct%tag=tag
    
    call init_sep3d(tag_ref,ref,"SCRATCH")
    call init_sep3d(struct%tag,struct%struct,"SCRATCH")

    allocate(struct%buffer(struct%struct%n(1),struct%struct%n(2),space%nzb))

    if(sep3d_ndims(ref)/=3) call seperr("sref must be 3-D")
    if(sep3d_ndims(struct%struct)/=3) call seperr("slowness must be 3-D")
    if(size(space%dz)/=struct%struct%n(3)) then
      write(0,*) "SIZE OF DZ",size(space%dz)
      write(0,*) "SIZE OF SLOWNESS",struct%struct%n(3)
      call seperr("Depth in slowness and nz steps not equal")
    end if

    if(ref%n(2)/=2) call seperr("expecting slowness reference axis 2 to be 2")
    allocate(temp(ref%n(1),ref%n(2),ref%n(3)))
    call sreed(tag_ref,temp,size(temp)*4)
    allocate(struct%reference(ref%n(3)))
    do i3=1,ref%n(3)
      i2=0
      do i1=1,ref%n(1)
       if(temp(i1,1,i3) >0.) i2=i2+1
      end do
      allocate(struct%reference(i3)%sref(i2))
      allocate(struct%reference(i3)%maxs(i2))
      i2=0
      do i1=1,ref%n(1)
       if(temp(i1,1,i3) >0.)  then
         i2=i2+1
         struct%reference(i3)%sref(i2)=temp(i1,1,i3)
         struct%reference(i3)%maxs(i2)=temp(i1,2,i3)
       end if
      end do
    end do
    

    allocate(struct%xmap_s(space%npad(1),space%npad(3))) !for now 
    allocate(struct%ymap_s(space%npad(2),space%npad(4))) !for now 
    allocate(struct%xmap_g(space%npad(1),space%npad(3))) !for now 
    allocate(struct%ymap_g(space%npad(2),space%npad(4))) !for now 
    do i2=1,space%npad(3)
      do i1=1,space%npad(1)
        pos1=space%opad(1)+space%dpad(1)*(i1-1)+space%opad(3)+space%dpad(3)*(i2-1)
        struct%xmap_g(i1,i2)=min(struct%struct%n(1),max(1,nint((pos1-struct%struct%o(1))/struct%struct%d(1)+1.)))
        pos2=space%opad(1)+space%dpad(1)*(i1-1)-space%opad(3)-space%dpad(3)*(i2-1)
        struct%xmap_s(i1,i2)=min(struct%struct%n(1),max(1,nint((pos2-struct%struct%o(1))/struct%struct%d(1)+1.)))
!        write(0,*) "ASSIGNING",i1,i2,struct%xmap_s(i1,i2),struct%xmap_g(i1,i2),pos1,pos2,space%opad(3),space%dpad(3)
      end do
    end do
    do i2=1,space%npad(4)
      do i1=1,space%npad(2)
        pos=space%opad(2)+space%dpad(2)*(i1-1)+space%opad(4)+space%dpad(4)*(i2-1)
        struct%ymap_g(i1,i2)=min(struct%struct%n(2),max(1,nint((pos-struct%struct%o(2))/struct%struct%d(2)+1.)))
        pos=space%opad(2)+space%dpad(2)*(i1-1)-space%opad(4)-space%dpad(4)*(i2-1)
        struct%ymap_s(i1,i2)=min(struct%struct%n(2),max(1,nint((pos-struct%struct%o(2))/struct%struct%d(2)+1.)))
      end do
    end do
   
#ifdef DEBUG
    write(0,*) "OUT INIT SLOW"
#endif
  end subroutine

  subroutine form_iref_slowness_slice(i3in,i4in,step,i3o,i4o)
    integer ::  i3in,i4in,i3o,i4o
    type(westep) :: step
    integer      :: i1,i2,iref,i3b,i3e,i4b,i4e,i_3,i_4,i3,i4
    real         :: s
    logical      :: doit
#ifdef DEBUG
    write(0,*) "IN FORM IREF SLOWNESS SLICE"
#endif

   if(step%wave%mem_save==0) then
     if(step%iref==1 .and. i3in==1 .and. i4in==1) then
       doit=.true.
       i3b=1; i4b=1; i3e=step%wave%npad(3); i4e=step%wave%npad(4)
     else
       doit=.false.
     end if
     i3o=i3in; i4o=i4in
   else
     doit=.true. 
     i3b=i3in; i3e=i3in; i4b=i4in; i4e=i4in
     i3o=1; i4o=1
   end if
   if(doit) then
     i_4=0
     do i4=i4b,i4e
       i_4=i_4+1
       i_3=0
       do i3=i3b,i3e
         i_3=i_3+1
         do i2=1,size(step%islow,2)
           do i1=1,size(step%islow,1)
             s=step%comp%slow%buffer(step%comp%x(i1,i3),&
               step%comp%y(i2,i4),step%iz)
             iref=1
             do while(step%maxs(iref) <s)
               iref=iref+1
             end do
             step%islow(i1,i2,i_3,i_4)=iref
           end do !}
         end do!}
       end do!}
     end do!}
   end if 
#ifdef DEBUG
    write(0,*) "OUT FORM IREF SLOWNESS SLICE"
#endif

  end subroutine

  subroutine  form_slowness_slice(i3in,i4in,step,i3o,i4o)
    integer     :: i3in,i4in,i3o,i4o
    type(westep) :: step
    integer     :: i1,i2,i3,i4,i3b,i4b,i3e,i4e,i_3,i_4
    logical    :: doit

#ifdef DEBUG
    write(0,*) "IN FORM SLOWNESS SLICE"
#endif
   if(step%wave%mem_save==0) then
     if(step%iref==1 .and. i3in==1 .and. i4in==1) then
       doit=.true.
       i3b=1; i4b=1; i3e=step%wave%npad(3); i4e=step%wave%npad(4)
     else
       doit=.false.
     end if
     i3o=i3in; i4o=i4in
   else
     doit=.true. 
     i3b=i3in; i3e=i3in; i4b=i4in; i4e=i4in
     i3o=1; i4o=1
   end if
   if(doit) then
     i_4=0
     do i4=i4b,i4e
       i_4=i_4+1
       i_3=0
       do i3=i3b,i3e
         i_3=i_3+1
         do i2=1,size(step%islow,2)
           do i1=1,size(step%islow,1)
             step%slow(i1,i2,i_3,i_4)=&
               step%comp%slow%buffer(step%comp%x(i1,i3),&
               step%comp%y(i2,i4),step%iz)
           end do
         end do
       end do
     end do
   end if 
!write(0,*) "CHECK 1",i3b,i3e
!if(i3b==2) call seperr("")
#ifdef DEBUG
    write(0,*) "OUT FORM SLOWNESS SLICE"
#endif

  end subroutine
    
  subroutine read_slow(space,iz)
    type(wspace)           :: space
    integer                 :: iz     !first slowness to read
    integer                 :: n(3),f(3),j(3),i1,i2,i3
    integer,external        :: sreed_window


#ifdef DEBUG
    write(0,*) "IN READ SLOW"
#endif

    j=1; f=0; f(3)=iz-1

    n=space%src%slow%struct%n;n(3)=min(size(space%src%slow%buffer,3),space%src%slow%struct%n(3)-iz+1)
    if(0/=sreed_window(space%src%slow%tag,3,space%src%slow%struct%n,n,f,j,4,space%src%slow%buffer)) &
      call seperr("trouble reading in slowness")
    if(space%src%slow%buffer(1,1,1) >1.)  then
      do i3=1,size(space%src%slow%buffer,3)
        do i2=1,size(space%src%slow%buffer,2)
          do i1=1,size(space%src%slow%buffer,1)
            space%src%slow%buffer(i1,i2,i3)=1./space%src%slow%buffer(i1,i2,i3)
          end do
        end do
      end do
    end if
    

    if(trim(space%src%slow%tag) /= trim(space%rec%slow%tag)) then
      n=space%rec%slow%struct%n;n(3)=min(size(space%rec%slow%buffer,3),space%rec%slow%struct%n(3)-iz-1)
      if(0/=sreed_window(space%rec%slow%tag,3,space%rec%slow%struct%n,n,f,j,4,space%rec%slow%buffer)) &
        call seperr("trouble reading in slowness")
      if(space%rec%slow%buffer(1,1,1) >1.)  then
        do i3=1,size(space%rec%slow%buffer,3)
          do i2=1,size(space%rec%slow%buffer,2)
            do i1=1,size(space%rec%slow%buffer,1)
              space%rec%slow%buffer(i1,i2,i3)=1./space%rec%slow%buffer(i1,i2,i3)
            end do
          end do
        end do
      end if
    end if
    
#ifdef DEBUG
    write(0,*) "OUT READ SLOW"
#endif
  end subroutine





end module
