module OW_pickbp1

use sep
use OW_type
use OW_parms_areal
use lint1
use lint
use igrad1

implicit none

integer, public, save, allocatable,dimension(:,:,:)		:: pickr,picks
integer, public, save					:: wavz,nhorz

contains

subroutine owpick()

character(len=128) 						:: hor1
integer,pointer,dimension(:,:)				:: pick1
real,pointer,dimension(:)					:: z,x,weight,sh,sh1
real, allocatable,dimension(:,:)				:: coorz1,coorz1r,coorz2r,coorz1s,coorz2s
integer							:: nhorz
integer							:: i,m,n,j,k,stat,ns,ms,nh,nhh(999)

 nhorz=0
 allocate(sh1(sog%sh%n));sh1=0.
 if (exist_file("hor1rec")) then
    nhorz=1
    call from_param("wavz",wavz,7)
    call from_param('hor1rec',hor1)
    open(51,file=hor1,status='old')
    read(51,*) m
    allocate(pick1(m,4));pick1=0
    ms=9999
    nh=0
    j=0
    do n=1,m
       read(51,*,end=1) pick1(n,:)
       if(pick1(n,3).ne.ms) then
         ms=pick1(n,3)
         j=j+1
         sh1(j)=float(ms)
         nh=nh+1
         i=1
       else
         i=i+1
         nhh(nh)=i
       end if
    end do
1   close(51)
    allocate(sh(j));sh=sh1(1:j)
    deallocate(sh1)
    allocate(coorz1(nh,sog%xm%n))

    n=0; m=1
    do i=1,nh
       n=n+nhh(i)
       allocate(z(n-m+1));z=float(pick1(m:n,2))
       allocate(x(n-m+1));x=float(pick1(m:n,4))
       allocate(weight(n-m+1)); weight=1.
       call lint_init (.true.,float(pick1(m,4)),1.,sog%xm%n,x,weight)
       stat = lint_lop (.true.,.false., coorz1(i,:), z)
!       call lint1_init (float(pick1(m,4)),1.,x)
!       call solver_reg( m=coorz1(i,:), d=z, Fop=lint1_lop, Aop=igrad1_lop, &
!			stepper=cgstep, nAop=3, niter=120,eps=1.)
       coorz1(i,int(x(n-m+1)):sog%xm%n)=9999.
       m=m+nhh(i)
       deallocate(z,x,weight)
    end do

    allocate(coorz1r(sog%sh%n,sog%xm%n))
    allocate(weight(j)); weight=1.
    call lint1_init (1.,1.,sh)
    do i=1,sog%xm%n
       call solver_reg( m=coorz1r(:,i), d=coorz1(:,i), Fop=lint1_lop, Aop=igrad1_lop, &
			stepper=cgstep, nAop=sog%sh%n, niter=20,eps=0.1)
       coorz1r(1:int(sh(1)-1),i)=9999.
       coorz1r(int(sh(j)+1):sog%sh%n,i)=9999.
    end do
    deallocate(weight,coorz1,pick1)

 end if	

 allocate(sh1(sog%sh%n));sh1=0.
 if (exist_file("hor2rec")) then
    nhorz=nhorz+1
    call from_param('hor2rec',hor1)
    open(51,file=hor1,status='old')
    read(51,*) m
    allocate(pick1(m,4));pick1=0
    ms=9999
    nh=0
    j=0
    do n=1,m
       read(51,*,end=2) pick1(n,:)
       if(pick1(n,3).ne.ms) then
         ms=pick1(n,3)
         j=j+1
         sh1(j)=float(ms)
         nh=nh+1
         i=1
       else
         i=i+1
         nhh(nh)=i
       end if
    end do
2   close(51)
    allocate(sh(j));sh=sh1(1:j)
    deallocate(sh1)
    allocate(coorz1(nh,sog%xm%n))

    n=0; m=1
    do i=1,nh
       n=n+nhh(i)
       allocate(z(n-m+1));z=float(pick1(m:n,2))
       allocate(x(n-m+1));x=float(pick1(m:n,4))
       allocate(weight(n-m+1)); weight=1.
       call lint_init (.true.,float(pick1(m,4)),1.,sog%xm%n,x,weight)
       stat = lint_lop (.true.,.false., coorz1(i,:), z)
!       call lint1_init (float(pick1(m,4)),1.,x)
!       call solver_reg( m=coorz1(i,:), d=z, Fop=lint1_lop, Aop=igrad1_lop, &
!			stepper=cgstep, nAop=sog%xm%n, niter=20,eps=0.1)
       coorz1(i,int(x(n-m+1)):sog%xm%n)=9999.
       m=m+nhh(i)
       deallocate(z,x,weight)
    end do

    allocate(coorz2r(sog%sh%n,sog%xm%n))
    allocate(weight(j)); weight=1.
    call lint1_init (1.,1.,sh)
    do i=1,sog%xm%n
       call solver_reg( m=coorz2r(:,i), d=coorz1(:,i), Fop=lint1_lop, Aop=igrad1_lop, &
			stepper=cgstep, nAop=sog%sh%n, niter=20,eps=0.1)
       coorz2r(1:int(sh(1)-1),i)=9999.
       coorz2r(int(sh(j)+1):sog%sh%n,i)=9999.
    end do
    deallocate(weight,coorz1,pick1)

 end if	

 allocate(pickr(nhorz,sog%sh%n,sog%xm%n))
 if (nhorz.eq.1) then
    pickr(1,:,:)=coorz1r
    deallocate(coorz1r)
 else if (nhorz.eq.2) then
    pickr(1,:,:)=coorz1r
    pickr(2,:,:)=coorz2r
    deallocate(coorz1r,coorz2r)
 end if

 allocate(sh1(sog%sh%n));sh1=0.
 if (exist_file("hor1sou")) then
    call from_param('hor1sou',hor1)
    open(51,file=hor1,status='old')
    read(51,*) m
    allocate(pick1(m,4));pick1=0
    ms=9999
    nh=0
    j=0
    do n=1,m
       read(51,*,end=3) pick1(n,:)
       if(pick1(n,3).ne.ms) then
         ms=pick1(n,3)
         j=j+1
         sh1(j)=float(ms)
         nh=nh+1
         i=1
       else
         i=i+1
         nhh(nh)=i
       end if
    end do
3   close(51)
    allocate(sh(j));sh=sh1(1:j)
    deallocate(sh1)
    allocate(coorz1(nh,sog%xm%n))

    n=0; m=1
    do i=1,nh
       n=n+nhh(i)
       allocate(z(n-m+1));z=float(pick1(m:n,2))
       allocate(x(n-m+1));x=float(pick1(m:n,4))
       allocate(weight(n-m+1)); weight=1.
       call lint_init (.true.,float(pick1(m,4)),1.,sog%xm%n,x,weight)
       stat = lint_lop (.true.,.false., coorz1(i,:), z)
!       call lint1_init (float(pick1(m,4)),1.,x)
!       call solver_reg( m=coorz1(i,:), d=z, Fop=lint1_lop, Aop=igrad1_lop, &
!			stepper=cgstep, nAop=sog%xm%n, niter=20,eps=0.1)
       coorz1(i,int(x(n-m+1)):sog%xm%n)=9999.
       m=m+nhh(i)
       deallocate(z,x,weight)
    end do

    allocate(coorz1s(sog%sh%n,sog%xm%n))
    allocate(weight(j)); weight=1.
    call lint1_init (1.,1.,sh)
    do i=1,sog%xm%n
       call solver_reg( m=coorz1s(:,i), d=coorz1(:,i), Fop=lint1_lop, Aop=igrad1_lop, &
			stepper=cgstep, nAop=sog%sh%n, niter=20,eps=0.1)
       coorz1s(1:int(sh(1)-1),i)=9999.
       coorz1s(int(sh(j)+1):sog%sh%n,i)=9999.
    end do
    deallocate(weight,coorz1,pick1)

 end if	

 allocate(sh1(sog%sh%n));sh1=0.
 if (exist_file("hor2sou")) then
    call from_param('hor2sou',hor1)
    open(51,file=hor1,status='old')
    read(51,*) m
    allocate(pick1(m,4));pick1=0
    ms=9999
    nh=0
    j=0
    do n=1,m
       read(51,*,end=4) pick1(n,:)
       if(pick1(n,3).ne.ms) then
         ms=pick1(n,3)
         j=j+1
         sh1(j)=float(ms)
         nh=nh+1
         i=1
       else
         i=i+1
         nhh(nh)=i
       end if
    end do
4   close(51)
    allocate(sh(j));sh=sh1(1:j)
    deallocate(sh1)
    allocate(coorz1(nh,sog%xm%n))

    n=0; m=1
    do i=1,nh
       n=n+nhh(i)
       allocate(z(n-m+1));z=float(pick1(m:n,2))
       allocate(x(n-m+1));x=float(pick1(m:n,4))
       allocate(weight(n-m+1)); weight=1.
       call lint_init (.true.,float(pick1(m,4)),1.,sog%xm%n,x,weight)
       stat = lint_lop (.true.,.false., coorz1(i,:), z)
!       call lint1_init (float(pick1(m,4)),1.,x)
!       call solver_reg( m=coorz1(i,:), d=z, Fop=lint1_lop, Aop=igrad1_lop, &
!			stepper=cgstep, nAop=sog%xm%n, niter=20,eps=0.1)
       coorz1(i,int(x(n-m+1)):sog%xm%n)=9999.
       m=m+nhh(i)
       deallocate(z,x,weight)
    end do

    allocate(coorz2s(sog%sh%n,sog%xm%n))
    allocate(weight(j)); weight=1.
    call lint1_init (1.,1.,sh)
    do i=1,sog%xm%n
       call solver_reg( m=coorz2s(:,i), d=coorz1(:,i), Fop=lint1_lop, Aop=igrad1_lop, &
			stepper=cgstep, nAop=3, niter=120,eps=10.)
!       call solver_reg( m=coorz2s(:,i), d=coorz1(:,i), Fop=lint1_lop, Aop=igrad1_lop, &
!			stepper=cgstep, nAop=sog%sh%n, niter=20,eps=0.1)
       coorz2s(1:int(sh(1)-1),i)=9999.
       coorz2s(int(sh(j)+1):sog%sh%n,i)=9999.
    end do
    deallocate(weight,coorz1,pick1)

 end if	

 allocate(picks(nhorz,sog%sh%n,sog%xm%n))
 if (nhorz.eq.1) then
    picks(1,:,:)=coorz1s
    deallocate(coorz1s)
 else if (nhorz.eq.2) then
    picks(1,:,:)=coorz1s
    picks(2,:,:)=coorz2s
    deallocate(coorz1s,coorz2s)
 end if

do i=1,nhorz
   do n=1,sog%xm%n
      do m=1,sog%sh%n
         write(5,*) i,m,n,picks(i,m,n)
      end do
   end do
end do

end subroutine
end module
