module allder_3d  
  use adj_mod
  implicit none
  integer, private :: nx 
  integer, private :: ny 
  integer, private :: nz 
  real , private,  dimension (:), pointer  :: pp 
  logical, private  :: dave 
  contains
  subroutine allder_3d_init ( nx_in,ny_in,nz_in,pp_in,dave_in)
    integer    :: nx_in, ny_in, nz_in
    real ,   dimension (:), pointer :: pp_in 
    logical    :: dave_in
    nx = nx_in; ny = ny_in; nz = nz_in; pp => pp_in;
    dave=dave_in
  end subroutine
  function allder_3d_lop ( adj, add, xx, yy) result(stat)
    integer            :: stat 
    logical,intent(in) :: adj,add 
    real,dimension(:)  :: xx,yy
        call       adjnull (adj,add,xx,yy )
        if(dave)then
                call allder_3d_lop3(adj,add,xx,yy )
        else
                call allder_3d_lop2(adj,add,xx,yy )
        endif
        stat=0
  end function 
  subroutine allder_3d_lop2(adj,add,xx,yy)
    logical,intent(in) :: adj,add 
    real, dimension (nx,ny,nz)  :: xx 
    real, dimension (nx,ny,nz)  :: yy 
    integer :: ix, iy,i2,iz
    real    :: a, ap, am, p
    integer, parameter :: dp = selected_real_kind(10,50)
    real(kind=dp) :: t
    do iz = 1,nz 
      do iy = 2,ny  
        do ix = 2,nx-1  
          p = pp(ix+(iy-2)*nx+nx*ny*(iz-1))            
          ap = -3.-2.*p
          am =  3.-2.*p
          a = 4.*p        
          ap=ap/12.
          am=am/12.
          a=a/12.
          if ( adj) then
            xx(ix-1,iy,iz)   =xx(ix-1,iy  ,iz) + real(yy(ix,iy-1,iz),dp) * real(am,dp)
            xx(ix-1,iy-1,iz) =xx(ix-1,iy-1,iz) - real(yy(ix,iy-1,iz),dp) * real(ap,dp)
            xx(ix,  iy-1,iz) =xx(ix,  iy-1,iz) - real(yy(ix,iy-1,iz),dp) * real(a ,dp)
            xx(ix,  iy,iz)   =xx(ix,  iy  ,iz) + yy(ix,iy,iz) * a
            xx(ix+1,iy-1,iz) =xx(ix+1,iy-1,iz) - real(yy(ix,iy-1,iz),dp) * real(am,dp)
            xx(ix+1,iy,iz)   =xx(ix+1,iy  ,iz) + real(yy(ix,iy-1,iz),dp) * real(ap,dp)
          else
            t= yy(ix,iy-1,iz)
            t=t+real((xx(ix-1,iy,iz)-xx(ix+1,iy-1,iz)),dp) *real( am,dp)
            t=t+real((xx(ix  ,iy,iz)-xx(ix  ,iy-1,iz)),dp) *real( a ,dp)
            t=t+real((xx(ix+1,iy,iz)-xx(ix-1,iy-1,iz)),dp) *real( ap,dp)
            yy(ix,iy-1,iz)=t
            !write(0,'(2I5,10(f8.5,2x))')ix,iz,t,dp,ap,real(ap,dp),real((xx(ix-1,iy,iz)-xx(ix+1,iy-1,iz)),dp)
          end if
        end do
      end do
    end do
  end subroutine 
  subroutine allder_3d_close()
  end subroutine 
  subroutine allder_3d_lop3(adj,add,xx,yy)
    logical,intent(in) :: adj,add 
    real, dimension (nx,ny,nz)  :: xx 
    real, dimension (nx,ny,nz)  :: yy 
    integer :: ix, iy,i2,iz
    real    :: a, ap, am, p, m22, p22, mp, x 
    integer, parameter :: dp = selected_real_kind(10,50)
    real(kind=dp) :: t
    do iz = 1,nz 
      do iy = 2,ny  
        do ix = 2,nx-1  
          p = pp(ix+(iy-2)*nx+nx*ny*(iz-1)); x=p
          m22 =-(1/2*x/(1+x**2)-1/2*x**3/(1+x**2)**2-1/2/(1+x**2)**2*x+1/2/(1+x**2)-x**2/(1+x**2)**2)
          p22 =-(1/2*x/(1+x**2)-1/2*x**3/(1+x**2)**2-1/2/(1+x**2)**2*x-1/2/(1+x**2)+x**2/(1+x**2)**2)
          mp  =-(1/2*x/(1+x**2)-1/2*x**3/(1+x**2)**2+1/2/(1+x**2)**2*x)
          m22 = -0.25*(((1.+p*p)*(2.*(1.+p)) - (1.+p)**2 * 2.*p) / (1.+p*p)**2)
          p22 = -0.25*(((1.+p*p)*(2.*(1.-p)) + (1.-p)**2 * 2.*p) / (1.+p*p)**2)
          mp =   p / (1.+p*p)**2
          if ( adj) then
            xx(ix-1,iy,iz)   =xx(ix-1,iy  ,iz) + real(yy(ix,iy-1,iz),dp) * real(am,dp)
            xx(ix-1,iy-1,iz) =xx(ix-1,iy-1,iz) - real(yy(ix,iy-1,iz),dp) * real(ap,dp)
            xx(ix,  iy-1,iz) =xx(ix,  iy-1,iz) - real(yy(ix,iy-1,iz),dp) * real(a ,dp)
            xx(ix,  iy,iz)   =xx(ix,  iy  ,iz) + yy(ix,iy,iz) * a
            xx(ix+1,iy-1,iz) =xx(ix+1,iy-1,iz) - real(yy(ix,iy-1,iz),dp) * real(am,dp)
            xx(ix+1,iy,iz)   =xx(ix+1,iy  ,iz) + real(yy(ix,iy-1,iz),dp) * real(ap,dp)
          else
            t= yy(ix,iy-1,iz)
            t=t+real(xx(ix-1,iy,iz),dp)*real(2.*mp,dp)-real(xx(ix+1,iy-1,iz),dp)*real(2.*p22,dp)
            t=t+real(xx(ix  ,iy,iz),dp)*real(0.   ,dp)-real(xx(ix  ,iy-1,iz),dp)*real(4.*mp ,dp)
            t=t+real(xx(ix+1,iy,iz),dp)*real(2.*mp,dp)-real(xx(ix-1,iy-1,iz),dp)*real(2.*m22,dp)
            yy(ix,iy-1,iz)=t
          end if
        end do
      end do
    end do
  end subroutine 
end module 
