# 1 "<stdin>"
# 1 "<built-in>"
# 1 "<command line>"
# 1 "<stdin>"
module pick_mod
  use sep3d_struct_mod
  implicit none
  contains
  logical function read_data(data_s,data,stack)
    type(sep3d) :: data_s
    real, pointer :: data(:,:,:),stack(:,:,:)
    real, allocatable :: datas(:,:,:)
    integer, external :: sreed
    read_data=.false.
    call init_sep3d("in",data_s,"INPUT")
    if (sep3d_ndims(data_s).ne.4) then
      write(0,*) "expecting 3-D data"
      return
    end if
    allocate(datas(data_s%n(1),data_s%n(2),data_s%n(4)))
    allocate(data(data_s%n(1),data_s%n(2),data_s%n(4)))
    allocate(stack(data_s%n(1),data_s%n(2),data_s%n(4)))
    if (4*size(datas).ne. sreed("in",datas,size(Datas)*4)) then
      write(0,*) "trouble reading data"
      return
    end if
    data=datas
!if(4*size(datas)!= sreed("stack",datas,size(Datas)*4)){
! write(0,*) "trouble reading stack";return
!}
    stack=datas
    deallocate(datas)
    read_data=.true.
  end function
  logical function calc_var(data,model,var,data_s)
    type(sep3d) :: data_s
    real :: data(:,:,:),model(:)
    real, pointer:: var(:)
    real :: gam(size(data,3)),v,mymax
    integer :: i1,i2,i,ivar,ia
    do i1=1,size(gam)
      gam(i1)=data_s%o(3)+data_s%d(3)*(i1-1)
    end do
    allocate(var(size(model)))
    i=0
    mymax=0
    do i2=1,size(data,2)
      do i1=1,size(data,1)
        i= i+1
        v=(model(i)-data_s%o(3))/data_s%d(3)+1.
        ia=v
        v=v-ia
        if (ia<1) then
          ia=1
          v=0.
        end if
        if (ia>=size(Data,3)) then
          ia=size(data,3)-1
          v=1.
        end if
        v=(1.-v)*data(i1,i2,ia)+v*data(i1,i2,ia+1)
        data(i1,i2,:)=data(i1,i2,:)/(v+.0001)
        if (sum(data(i1,i2,:)).ne.0.) then
          var(i)=sum(((model(i)-gam(:))*data(i1,i2,:))**2)/(.001+sum&
            &(data(i1,i2,:)**2))
          mymax=max(mymax,var(i))
        else
          var(i)=-1
        end if
      end do
    end do
    where(var<0.)
      var=mymax
end where
    var=sqrt(var)
    calc_var=.true.
  end function
  logical function write_files(maxp,wt,model,data_s,var)
    type(sep3d) :: data_s
    real :: maxp(:,:),model(:),wt(:,:),var(:)
    integer, external :: getch, srite
    integer ::i
    character(len=128) :: junk
    write_files=.false.
    if (1.eq.getch("wt_out","d",i)) then
      call sep_put_data_axis_par("wt_out",1,data_s%n(1),data_s%o(1)&
        &,data_s%d(1),data_s%label(1))
      call sep_put_data_axis_par("wt_out",2,data_s%n(2),data_s%o(2)&
        &,data_s%d(2),data_s%label(2))
      if (4*size(wt).ne.srite("wt_out",wt,size(wt)*4)) then
        write(0,*) "trouble writing out wt_out"
        return
      end if
    end if
    if (1.eq.getch("pmax","s",junk)) then
      call sep_put_data_axis_par("pmax",1,data_s%n(1),data_s%o(1),data_s%d&
        &(1),data_s%label(1))
      call sep_put_data_axis_par("pmax",2,data_s%n(2),data_s%o(2),data_s%d&
        &(2),data_s%label(2))
      write(0,*) "CHECK MAXP",minval(maxp),maxval(maxp),size(maxp)
      if (4*size(maxp).ne.srite("pmax",maxp,size(maxp)*4)) then
        write(0,*) "trouble writing out max p"
        return
      end if
    end if
    if (1.eq.getch("var","s",junk)) then
      call sep_put_data_axis_par("var",1,data_s%n(1),data_s%o(1),data_s%d&
        &(1),data_s%label(1))
      call sep_put_data_axis_par("var",2,data_s%n(2),data_s%o(2),data_s%d&
        &(2),data_s%label(2))
      if (4*size(var).ne.srite("var",var,size(var)*4)) then
        write(0,*) "trouble writing out var"
        return
      end if
    end if
    call sep_put_data_axis_par("out",2,data_s%n(2),data_s%o(2),data_s%d&
      &(2),data_s%label(2))
    call sep_put_data_axis_par("out",4,1,0.,1.,"none")
    if (4*size(maxp).ne.srite("out",model,size(model)*4)) then
      write(0,*) "trouble writing out "
      return
    end if
    write_files=.true.
  end function
end module
