!!.. This module generates the cubic B-spline basis functions
module base_spline_mod
  use OW_type
  use OW_parms_areal
  implicit none

  type(axis), private :: axd1, ayd1, azd1
  type(axis), private :: axk1, ayk1, azk1

contains

!------------------------------------------------------------------------
  integer function base_spln1d_init(axd1_in,axk1_in) result(stat)
    type(axis) :: axd1_in
    type(axis) :: axk1_in
    axd1  = axd1_in
    axk1  = axk1_in
    stat = 0
  end function
!------------------------------------------------------------------------
  integer function base_spln2d_init(axd1_in,ayd1_in,axk1_in,ayk1_in) result(stat)
    type(axis) :: axd1_in,ayd1_in
    type(axis) :: axk1_in,ayk1_in
    axd1 = axd1_in; ayd1 = ayd1_in
    axk1 = axk1_in; ayk1 = ayk1_in
    stat = 0
  end function
!------------------------------------------------------------------------
  integer function base_spln3d_init(axd1_in,ayd1_in,azd1_in,axk1_in,ayk1_in,azk1_in) result(stat)
    type(axis) :: axd1_in,ayd1_in,azd1_in
    type(axis) :: axk1_in,ayk1_in,azk1_in
    axd1 = axd1_in; ayd1 = ayd1_in; azd1 = azd1_in
    axk1 = axk1_in; ayk1 = ayk1_in; azk1 = azk1_in
    stat = 0
  end function
!------------------------------------------------------------------------
  integer function base_spln1d_clos() result(stat)
    stat = 0
  end function
!------------------------------------------------------------------------
  integer function base_spln2d_clos() result(stat)
    stat = 0
  end function
!------------------------------------------------------------------------
  integer function base_spln3d_clos() result(stat)
    stat = 0
  end function
!------------------------------------------------------------------------
  integer function base_spln1d_oper(adj, add, modl, data) result(stat)
    logical                :: adj, add
    double precision, dimension(axk1%n) :: modl
    double precision, dimension(axd1%n) :: data

    integer                :: ixd, ixk
    real                   :: xd, xk
    if (.not.add) then
       if (adj) then
          modl = 0.
       else 
          data = 0.
       end if
    end if

    if (adj) then
       do ixd=1, axd1%n; xd = axd1%o+(ixd-1)*axd1%d
       do ixk=1, axk1%n; xk = axk1%o+(ixk-1)*axk1%d
          modl(ixk) = modl(ixk) + data(ixd)*spline(xd-xk,axk1%d)
       end do 
       end do
    else
       do ixd=1, axd1%n; xd = axd1%o+(ixd-1)*axd1%d
       do ixk=1, axk1%n; xk = axk1%o+(ixk-1)*axk1%d
          data(ixd) = data(ixd) + modl(ixk)*spline(xd-xk,axk1%d)
       end do
       end do
    end if
    stat = 0
  end function
!------------------------------------------------------------------------
  integer function base_spln2d_oper(adj, add, modl, data) result(stat)
    logical            :: adj, add
    double precision, dimension(:) :: modl
    double precision, dimension(:) :: data

    if (.not.add) then
       if (adj) then
          modl = 0.
       else
          data = 0.
       end if
    end if
    call base_spln2d(adj, modl, data)
    stat = 0 
  end function
!------------------------------------------------------------------------
  subroutine base_spln2d(adj, modl, data)
    logical                      :: adj
    double precision, dimension(axk1%n,ayk1%n) :: modl
    double precision, dimension(axd1%n,ayd1%n) :: data
    double precision, allocatable              :: modelo(:,:,:)
    double precision, allocatable              :: dados(:,:,:)

    integer                      :: ixd, ixk,inode
    integer                      :: iyd, iyk,nodes,i
    real                         :: xd, yd
    real                         :: xk, yk
    integer, external	   :: omp_get_thread_num,omp_get_num_threads

    if (node.eq.9999) then
      !$OMP PARALLEL
      nodes = omp_get_num_threads()
      !$OMP END PARALLEL
    else
      nodes=node
    end if

    if (adj) then
       allocate(modelo(axk1%n,ayk1%n,nodes));modelo=0
       !$OMP parallel do private(inode,iyd,ixd,iyk,ixk,yd,xd,yk,xk) schedule(dynamic,1) num_threads(nodes)
       do iyd=1, ayd1%n; inode=omp_get_thread_num()+1;yd = ayd1%o+(iyd-1)*ayd1%d 
       do ixd=1, axd1%n; xd = axd1%o+(ixd-1)*axd1%d
          do iyk=1, ayk1%n; yk = ayk1%o+(iyk-1)*ayk1%d
          do ixk=1, axk1%n; xk = axk1%o+(ixk-1)*axk1%d
             modelo(ixk,iyk,inode) = modelo(ixk,iyk,inode) + data(ixd,iyd)*spline(xd-xk,axk1%d)*spline(yd-yk,ayk1%d)
          end do
          end do
       end do
       end do
       !$OMP END PARALLEL DO
       do i=1,nodes
          do iyk=1, ayk1%n
             do ixk=1, axk1%n
                modl(ixk,iyk)=modl(ixk,iyk)+modelo(ixk,iyk,i)
             end do
          end do
       end do
       deallocate(modelo)
    else
       allocate(dados(axd1%n,ayd1%n,nodes));dados=0
       !$OMP parallel do private(inode,iyd,ixd,iyk,ixk,yd,xd,yk,xk) schedule(dynamic,1) num_threads(nodes)
       do iyd=1, ayd1%n; inode=omp_get_thread_num()+1;yd = ayd1%o+(iyd-1)*ayd1%d
       do ixd=1, axd1%n; xd = axd1%o+(ixd-1)*axd1%d
          do iyk=1, ayk1%n; yk = ayk1%o+(iyk-1)*ayk1%d
          do ixk=1, axk1%n; xk = axk1%o+(ixk-1)*axk1%d
             dados(ixd,iyd,inode) = dados(ixd,iyd,inode) + modl(ixk,iyk)*spline(xd-xk,axk1%d)*spline(yd-yk,ayk1%d)
          end do
          end do
       end do
       end do
       !$OMP END PARALLEL DO
       do i=1,nodes
          do iyd=1, ayd1%n
             do ixd=1, axd1%n
                data(ixd,iyd)=data(ixd,iyd)+dados(ixd,iyd,i)
             end do
          end do
       end do
       deallocate(dados)
    end if
  end subroutine
!------------------------------------------------------------------------
  integer function base_spln3d_oper(adj, add, modl, data) result(stat)
    logical            :: adj, add
    double precision, dimension(axk1%n,ayk1%n,azk1%n) :: modl
    double precision, dimension(axd1%n,ayd1%n,azd1%n) :: data

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

    call base_spln3d(adj, modl, data)

    stat = 0
  end function
!------------------------------------------------------------------------
  subroutine base_spln3d(adj, modl, data)
    logical                            :: adj
    double precision, dimension(axk1%n,ayk1%n,azk1%n) :: modl
    double precision, dimension(axd1%n,ayd1%n,azd1%n) :: data
    double precision, allocatable                     :: modelo(:,:,:,:)
    double precision, allocatable                     :: dados(:,:,:,:)
    real                         :: xd, yd, zd
    real                         :: xk, yk, zk
    integer                      :: ixd, ixk,inode
    integer                      :: izd, izk
    integer                      :: iyd, iyk,nodes,i
    integer, external	   :: omp_get_thread_num,omp_get_num_threads

    if (node.eq.9999) then
      !$OMP PARALLEL
      nodes = omp_get_num_threads()
      !$OMP END PARALLEL
    else
      nodes=node
    end if

    if (adj) then
       allocate(modelo(axk1%n,ayk1%n,azk1%n,nodes));modelo=0
       !$OMP parallel do private(inode,izd,iyd,ixd,izk,iyk,ixk,zd,yd,xd,zk,yk,xk) schedule(dynamic,1) num_threads(nodes)
       do izd=1, azd1%n; inode=omp_get_thread_num()+1; zd = azd1%o+(izd-1)*azd1%d!;write(0,*) "passei 1 ",izd,inode
       do iyd=1, ayd1%n; yd = ayd1%o+(iyd-1)*ayd1%d
       do ixd=1, axd1%n; xd = axd1%o+(ixd-1)*axd1%d
          do izk=1, azk1%n; zk = azk1%o+(izk-1)*azk1%d
          do iyk=1, ayk1%n; yk = ayk1%o+(iyk-1)*ayk1%d
          do ixk=1, axk1%n; xk = axk1%o+(ixk-1)*axk1%d
             modelo(ixk,iyk,izk,inode) = modelo(ixk,iyk,izk,inode) + data(ixd,iyd,izd)*spline(xd-xk,axk1%d)*spline(yd-yk,ayk1%d)*spline(zd-zk,azk1%d)
          end do
          end do
          end do
       end do
       end do
       end do
       !$OMP END PARALLEL DO
       do i=1,nodes
          do izk=1, azk1%n
             do iyk=1, ayk1%n
                do ixk=1, axk1%n
                   modl(ixk,iyk,izk)=modl(ixk,iyk,izk)+modelo(ixk,iyk,izk,i)
                end do
             end do
          end do
       end do
       deallocate(modelo)
    else
       allocate(dados(axd1%n,ayd1%n,azd1%n,nodes));dados=0
       !$OMP parallel do private(inode,izd,iyd,ixd,izk,iyk,ixk,zd,yd,xd,zk,yk,xk) schedule(dynamic,1) num_threads(nodes)
       do izd=1, azd1%n; inode=omp_get_thread_num()+1; zd = azd1%o+(izd-1)*azd1%d!;write(0,*) "passei 2 ",izd,inode
       do iyd=1, ayd1%n; yd = ayd1%o+(iyd-1)*ayd1%d
       do ixd=1, axd1%n; xd = axd1%o+(ixd-1)*axd1%d
          do izk=1, azk1%n; zk = azk1%o+(izk-1)*azk1%d
          do iyk=1, ayk1%n; yk = ayk1%o+(iyk-1)*ayk1%d
          do ixk=1, axk1%n; xk = axk1%o+(ixk-1)*axk1%d
             dados(ixd,iyd,izd,inode) = dados(ixd,iyd,izd,inode) + modl(ixk,iyk,izk)*spline(xd-xk,axk1%d)*spline(yd-yk,ayk1%d)*spline(zd-zk,azk1%d)
          end do
          end do
          end do
       end do
       end do
       end do
       !$OMP END PARALLEL DO
       do i=1,nodes
          do izd=1, azd1%n
             do iyd=1, ayd1%n
                do ixd=1, axd1%n
                   data(ixd,iyd,izd)=data(ixd,iyd,izd)+dados(ixd,iyd,izd,i)
                end do
             end do
          end do
       end do
       deallocate(dados)
    end if
  end subroutine
!------------------------------------------------------------------------
  real function spline(x, dx)
    real :: x, dx
    if (x<-2*dx) then
       spline = 0.
    else if (x>=-2*dx .and. x<-dx) then
       spline = (2*dx+x)**3/6./dx**3
    else if (x>=-dx .and. x<0.) then
       spline = (-x)**3/2./dx**3-(-x)**2/dx**2+2./3.
    else if (x>=0 .and. x<dx) then
       spline = (+x)**3/2./dx**3-(+x)**2/dx**2+2./3.
    else if (x>=dx .and. x<=2*dx) then
       spline = (2*dx-x)**3/6./dx**3
    else
       spline = 0.
    end if
  end function
!------------------------------------------------------------------------
end module
