module base_reglar_mod
  implicit none

  integer, private :: ndim
  integer, private :: nx, ny, nz
  integer, private :: regtyp
contains

  integer function base_reglar_init(nx_in, ny_in, nz_in, regtyp_in) result(stat)
    integer :: nx_in, ny_in, nz_in
    integer :: regtyp_in
    nx     = nx_in
    ny     = ny_in
    nz     = nz_in
    regtyp = regtyp_in
    ndim   = nx*ny*nz

    stat   = 0 
  end function

  integer function base_reglar_clos() result(stat)
    stat = 0
  end function

  integer function base_reglar_oper(adj, add, modl, data) result(stat)
    logical            :: adj, add
    real, dimension(:) :: modl
    real, dimension(:) :: data

    if (.not. add) then
       if (adj) then
          modl = 0.
       else
          data = 0.
       end if 
    end if

    select case (regtyp)
      case (1)
        call reglar_damping(adj, modl, data)
      case (2)
        call reglar_horlapl(adj, modl, data)
    end select  
    stat = 0
  end function

  !!.. assuming Gaussian a prior infomation
  subroutine reglar_damping(adj, modl, data)
    logical               :: adj
    real, dimension(ndim) :: modl
    real, dimension(ndim) :: data

    if (adj) then
       modl = modl + data
    else
       data = data + modl
    end if
    return
  end subroutine

  subroutine reglar_gradt1d(adj, modl, data)
  end subroutine

  !!.. horizontal laplacian
  subroutine reglar_horlapl(adj, modl, data)
    logical                   :: adj
    real, dimension(nx,ny,nz) :: modl
    real, dimension(nx,ny,nz) :: data
    integer                   :: ix, iy, iz
    if (ny > 1) then
       do iy=2, ny-1
       do ix=2, nx-1
          if (adj) then
             modl(ix  ,iy  ,:) = modl(ix  ,iy  ,:) - 4*data(ix,iy,:)
             modl(ix+1,iy  ,:) = modl(ix+1,iy  ,:) +   data(ix,iy,:)
             modl(ix-1,iy  ,:) = modl(ix-1,iy  ,:) +   data(ix,iy,:)
             modl(ix  ,iy+1,:) = modl(ix  ,iy+1,:) +   data(ix,iy,:)
             modl(ix  ,iy-1,:) = modl(ix  ,iy-1,:) +   data(ix,iy,:)
          else
             data(ix,iy,:) = data(ix,iy,:) + modl(ix+1,iy,:)+modl(ix-1,iy,:)+modl(ix,iy+1,:)+modl(ix,iy-1,:)-4*modl(ix,iy,:)
          end if
       end do
       end do
    else
       do ix=2, nx-1
          if (adj) then
             modl(ix  ,:,:) = modl(ix  ,:,:) - 2*data(ix,:,:)
             modl(ix+1,:,:) = modl(ix+1,:,:) +   data(ix,:,:)
             modl(ix-1,:,:) = modl(ix-1,:,:) +   data(ix,:,:)
          else
             data(ix,:,:) = data(ix,:,:) + modl(ix+1,:,:)+modl(ix-1,:,:)-2*modl(ix,:,:)
          end if
       end do
       if (adj) then
          modl(1,:,:) = 0.; modl(nx,:,:) = 0.
       else
          data(1,:,:) = 0.; data(nx,:,:) = 0.
       end if  
    end if
  end subroutine
  !!.. Laplacian regularization
!  subroutine reglar_laplace(adj, modl, data)
    
!  end subroutine

  !!.. vector projection regularization
  subroutine reglar_vecproj(adj, modl, data)
    logical


    

  end subroutine
 

end module









