module gradchain  
  use chain_mod
  use causint
  use igrad1
  use adj_mod
  implicit none
  integer, private, save                     :: n1, n2
  real, private, save                        :: epsx
  real, dimension(:), private, save, pointer :: wt
  contains
!---------------------------------------------------------------------------------
  subroutine gradchain_init( wt_in )
    real, dimension(:), pointer :: wt_in
    allocate( wt(size(wt_in)) )
    wt = wt_in
  end subroutine
!---------------------------------------------------------------------------------
!---------------------------------------------------------------------------------
  subroutine gradchain2d_init( n1_in, n2_in, epsx_in )
    integer                    :: n1_in, n2_in
    real                       :: epsx_in
    n1 = n1_in
    n2 = n2_in
    epsx = epsx_in
  end subroutine
!---------------------------------------------------------------------------------
!---------------------------------------------------------------------------------
  integer function diagweight_lop( adj, add, model, data)
    real, dimension(:)            :: model, data
    logical,intent(in)            :: adj,add
    call adjnull( adj, add, model, data)
    if (adj) then
      model = model + wt*data
    else
      data  = data  + wt*model
    end if 
    diagweight_lop = 0
  end function
!---------------------------------------------------------------------------------
!---------------------------------------------------------------------------------
  integer function weighted_grad_lop( adj, add, model, data)
    logical,intent(in) :: adj,add
    real, dimension(:) :: model, data
    real, pointer      :: temp(:)
    allocate( temp(size(model)) )
    call chain( diagweight_lop, igrad1_lop, adj, add, model, data, temp&
      & )
    deallocate( temp)
    weighted_grad_lop = 0
  end function
!---------------------------------------------------------------------------------
!---------------------------------------------------------------------------------
  integer function weighted_causint_lop( adj, add, model, data)
    logical,intent(in) :: adj,add
    real, dimension(:) :: model, data
    real, pointer      :: temp(:)
    allocate( temp(size(model)) )
    call chain( causint_lop, diagweight_lop, adj, add, model, data, temp&
      & )
    deallocate( temp)
    weighted_causint_lop = 0
  end function
!---------------------------------------------------------------------------------
!---------------------------------------------------------------------------------
  integer function weighted_txgrad_lop( adj, add, model, data)
    logical,intent(in) :: adj,add
    real, dimension(:) :: model, data
    real, pointer      :: temp1(:), temp2(:)
    allocate( temp1(size(model)), temp2(size(model)) )
    call chain( diagweight_lop, xgrad1_lop, tgrad1_lop, adj, add, model&
      &, data, temp1, temp2 )
    deallocate( temp1, temp2)
    weighted_txgrad_lop = 0
  end function
!---------------------------------------------------------------------------------
!---------------------------------------------------------------------------------
  integer function xgrad1_lop( adj, add, model, data)
    logical,intent(in)            :: adj,add
    real, dimension(:)            :: model, data
    real, dimension(:,:), pointer :: model2d, data2d
    integer                       :: stat, i1, i2
    call adjnull( adj, add, model, data)
    allocate( model2d(n1,n2), data2d(n1,n2) )
    model2d = reshape(model,(/n1,n2/))
    data2d = reshape(data,(/n1,n2/))
    do i1=1,n1  
      if ( adj) then
        stat = igrad1_lop( .true., add, model2d(i1,:), data2d(i1,:) )
      else
        stat = igrad1_lop( .false., add, model2d(i1,:), data2d(i1,:) )
      end if
    end do 
    model = epsx*reshape(model2d,(/n1*n2/))
    data = epsx*reshape(data2d,(/n1*n2/))
!		model = reshape(model2d,(/n1*n2/))
!		data  = reshape(data2d,(/n1*n2/))
    deallocate( model2d, data2d)
    xgrad1_lop = 0
  end function
!---------------------------------------------------------------------------------
!---------------------------------------------------------------------------------
  integer function tgrad1_lop( adj, add, model, data)
    logical,intent(in)            :: adj,add
    real, dimension(:)            :: model, data
    real, dimension(:,:), pointer :: model2d, data2d
    integer                       :: stat, i1, i2
    call adjnull( adj, add, model, data)
    allocate( model2d(n1,n2), data2d(n1,n2) )
    model2d = reshape(model,(/n1,n2/))
    data2d = reshape(data,(/n1,n2/))
    do i2=1,n2  
      if ( adj) then
        stat = igrad1_lop( .true., add, model2d(:,i2), data2d(:,i2) )
      else
        stat = igrad1_lop( .false., add, model2d(:,i2), data2d(:,i2) )
      end if
    end do 
    model = reshape(model2d,(/n1*n2/))
    data = reshape(data2d,(/n1*n2/))
    deallocate( model2d, data2d)
    tgrad1_lop = 0
  end function
!---------------------------------------------------------------------------------
end module 
