module OW_pickbp_a

use sep
use OW_type
use OW_parms_areal
use lint2
use igrad2

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						:: z(:),coorxy(:,:)
real, allocatable,dimension(:)				:: coorz1,coorz2
integer							:: nhorz
integer							:: i,m,n,j,stat

 nhorz=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
    do n=1,m
       read(51,*,end=1) pick1(n,:)
    end do
1   close(51)
    allocate(coorz1(sog%sh%n*sog%xm%n))
    allocate(coorxy(m,2),z(m))
    coorxy=float(pick1(:,3:4))
    z=float(pick1(:,2))

    call lint2_init ( sog%sh%n,sog%xm%n,1.,1.,1.,1.,coorxy)
    call igrad2_init(sog%sh%n,sog%xm%n)  
    call solver_reg(m=coorz1, d=z,Fop=lint2_lop,Aop=igrad2_lop,stepper=cgstep,nAop=sog%sh%n*sog%xm%n*2, &
	 niter=120 ,eps=1.,verb=.true.)
    deallocate(pick1,coorxy,z)
 end if	

 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
    do n=1,m
       read(51,*,end=2) pick1(n,:)
    end do
2   close(51)
    allocate(coorz2(sog%sh%n*sog%xm%n))
    allocate(coorxy(m,2),z(m))
    coorxy=float(pick1(:,3:4))
    z=float(pick1(:,2))

    call lint2_init ( sog%sh%n,sog%xm%n,1.,1.,1.,1.,coorxy)
    call igrad2_init(sog%sh%n,sog%xm%n)  
    call solver_reg(m=coorz2, d=z,Fop=lint2_lop,Aop=igrad2_lop,stepper=cgstep,nAop=sog%sh%n*sog%xm%n*2, &
	 niter=120 ,eps=1.,verb=.true.) 
    deallocate(pick1,coorxy,z)
 end if	

 allocate(pickr(nhorz,sog%sh%n,sog%xm%n))
 if (nhorz.eq.1) then
    do i=1,sog%sh%n
       do j=1,sog%xm%n
          pickr(1,:,:)=coorz1((i-1)*sog%xm%n+j)
       end do
    end do
    deallocate(coorz1)
 else if (nhorz.eq.2) then
    do i=1,sog%sh%n
       do j=1,sog%xm%n
          pickr(1,:,:)=coorz1((i-1)*sog%xm%n+j)
          pickr(2,:,:)=coorz2((i-1)*sog%xm%n+j)
       end do
    end do
    deallocate(coorz1,coorz2)
 end if

 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
    do n=1,m
       read(51,*,end=3) pick1(n,:)
    end do
3   close(51)
    allocate(coorz1(sog%sh%n*sog%xm%n))
    allocate(coorxy(m,2),z(m))
    coorxy=float(pick1(:,3:4))
    z=float(pick1(:,2))

    call lint2_init ( sog%sh%n,sog%xm%n,1.,1.,1.,1.,coorxy)
    call igrad2_init(sog%sh%n,sog%xm%n)  
    call solver_reg(m=coorz1, d=z,Fop=lint2_lop,Aop=igrad2_lop,stepper=cgstep,nAop=sog%sh%n*sog%xm%n*2, &
	niter=120 ,eps=1.,verb=.true.)
    deallocate(pick1,coorxy,z)
 end if	

 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
    do n=1,m
       read(51,*,end=4) pick1(n,:)
    end do
4   close(51)
    allocate(coorz2(sog%sh%n*sog%xm%n))
    allocate(coorxy(m,2),z(m))
    coorxy=float(pick1(:,3:4))
    z=float(pick1(:,2))

    call lint2_init ( sog%sh%n,sog%xm%n,1.,1.,1.,1.,coorxy)
    call igrad2_init(sog%sh%n,sog%xm%n)  
    call solver_reg(m=coorz2, d=z,Fop=lint2_lop,Aop=igrad2_lop,stepper=cgstep,nAop=sog%sh%n*sog%xm%n*2, &
	 niter=120 ,eps=1.,verb=.true.)
    deallocate(pick1,coorxy,z)
 end if	

 allocate(picks(nhorz,sog%sh%n,sog%xm%n))
 if (nhorz.eq.1) then
    do i=1,sog%sh%n
       do j=1,sog%xm%n
          picks(1,i,j)=coorz1((i-1)*sog%xm%n+j)
       end do
    end do
    deallocate(coorz1)
 else if (nhorz.eq.2) then
    do i=1,sog%sh%n
       do j=1,sog%xm%n
          picks(1,i,j)=coorz1((i-1)*sog%xm%n+j)
          picks(2,i,j)=coorz2((i-1)*sog%xm%n+j)
       end do
    end do
    deallocate(coorz1,coorz2)
 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
