!====================================================================
!
      subroutine x_der2_co(f,d2f,dx,dy,dz,orderi,orderb,
     &      biL, biR, bjL, bjR, bkL, bkR,
     &      fiL, fiR, fjL, fjR, fkL, fkR,
     &      iiL, iiR, ijL, ijR, ikL, ikR,
     &     dfiL,dfiR,dfjL,dfjR,dfkL,dfkR,
     &     iperx,ipery,iperz,ierror,error)
!
      include 'der2_include.h'
!
      call der2_pre( ierror, error,
     &   xbegin,xend,ybegin,yend,zbegin,zend,
     &   bnd2A, bnd3A, bnd3B, bnd4C, bnd6D,
     &   iperx,ipery,iperz, orderi, orderb,
     &   biL, biR, bjL, bjR, bkL, bkR,
     &   fiL, fiR, fjL, fjR, fkL, fkR,
     &   iiL, iiR, ijL, ijR, ikL, ikR,
     &   piL, piR, pjL, pjR, pkL, pkR,
     &   dfiL,dfiR,dfjL,dfjR,dfkL,dfkR,
     &   nproLx,nproRx,nproLy,nproRy,nproLz,nproRz,
     &   proL,proR,promax,dx,dy,dz,rdx2,rdy2,rdz2)
!
      if( 2*promax+1 .gt. piR-piL+1) then
        ierror = -100
        error = 'DER_2: Stencil is wider than the f domain'
        return
      endif
!
! Zero-out second-derivative
!
      do k = ikL, ikR
        do j = ijL, ijR
          do i = iiL, iiR
            d2f(i,j,k) = 0.0
          end do
        end do
      end do
!
!========================================================================
! CENTERED DIFFERENCE OPERATORS
!
!-----------------------------
!  2nd-order explicit: (2-2E-2)
!
      if ( orderi .eq. 2 ) then
!
        Up  = (-2.d0/ 1.d0)* rdx2
        ae  = ( 1.d0/ 1.d0)* rdx2
!
!  Internal nodes:
!
        do k = ikL, ikR
          do j = ijL, ijR
            do i = xbegin, xend
              d2f(i,j,k) = Up *  f(i+0,j,k)
     &                   + ae *( f(i+1,j,k)+f(i-1,j,k) )
            end do   
          end do
        end do
!
!  Boundary nodes:
!   Left side:
!
        if ( nproLx .eq. 1 ) then
          do k = ikL, ikR
            do j = ijL, ijR
              d2f(xbegin-1,j,k) = (bnd2A(1)*f(xbegin-1,j,k) +
     &                             bnd2A(2)*f(xbegin+0,j,k) +
     &                             bnd2A(3)*f(xbegin+1,j,k))*rdx2
            end do
          end do
        end if
!
!   Right side:
!
        if ( nproRx .eq. 1 ) then
          do k = ikL, ikR
            do j = ijL, ijR
              d2f(xend+1,j,k) = (bnd2A(1)*f(xend+1,j,k) +
     &                           bnd2A(2)*f(xend+0,j,k) +
     &                           bnd2A(3)*f(xend-1,j,k))*rdx2
            end do
          end do
        end if
!
      end if
!
!-----------------------------
!  4th order explicit: (3,3-4E-3,3)
!
!
      if ( orderi .eq. 4 ) then

        Up  = -5.d0/ 2.d0 * rdx2
        ae  =  4.d0/ 3.d0 * rdx2
        be  = -1.d0/12.d0 * rdx2
!
!  Internal nodes:
!
        do k = ikL, ikR
          do j = ijL, ijR
            do i = xbegin, xend
              d2f(i,j,k) = Up *  f(i  ,j,k)
     &                   + ae *( f(i+1,j,k)+f(i-1,j,k) )
     &                   + be *( f(i+2,j,k)+f(i-2,j,k) )
            end do   
          end do   
        end do
!
!  Boundary nodes:
!
!   Left side:
!
!
        if ( nproLx .ge. 1 ) then
          do k = ikL, ikR
            do j = ijL, ijR
              d2f(xbegin-1,j,k) = (bnd3B(1)*f(xbegin-2,j,k) +
     &                             bnd3B(2)*f(xbegin-1,j,k) +
     &                             bnd3B(3)*f(xbegin+0,j,k) +
     &                             bnd3B(4)*f(xbegin+1,j,k) +
     &                             bnd3B(5)*f(xbegin+2,j,k))*rdx2
            end do
          end do
        end if
!
        if ( nproLx .eq. 2 ) then
          do k = ikL, ikR
            do j = ijL, ijR
              d2f(xbegin-2,j,k) = (bnd3A(1)*f(xbegin-2,j,k) +
     &                             bnd3A(2)*f(xbegin-1,j,k) +
     &                             bnd3A(3)*f(xbegin+0,j,k) +
     &                             bnd3A(4)*f(xbegin+1,j,k) +
     &                             bnd3A(5)*f(xbegin+2,j,k))*rdx2
            end do
          end do
        end if
!
!   Right side:
!
        if ( nproRx .ge. 1 ) then
          do k = ikL, ikR
            do j = ijL, ijR
              d2f(xend+1,j,k) = (bnd3B(1)*f(xend+2,j,k) +
     &                           bnd3B(2)*f(xend+1,j,k) +
     &                           bnd3B(3)*f(xend+0,j,k) +
     &                           bnd3B(4)*f(xend-1,j,k) +
     &                           bnd3B(5)*f(xend-2,j,k))*rdx2
            end do
          end do
        end if
!
        if ( nproRx .eq. 2 ) then
          do k = ikL, ikR
            do j = ijL, ijR
              d2f(xend+2,j,k) = (bnd3A(1)*f(xend+2,j,k) +
     &                           bnd3A(2)*f(xend+1,j,k) +
     &                           bnd3A(3)*f(xend+0,j,k) +
     &                           bnd3A(4)*f(xend-1,j,k) +
     &                           bnd3A(5)*f(xend-2,j,k))*rdx2
            end do
          end do
        end if
!
      end if
!
!-----------------------------
!  6th order explicit: (3,3,4-6E-4,3,3)
!
!
      if ( orderi .eq. 6 ) then
!
        Up  =-49.d0/18.d0 * rdx2
        ae  =  3.d0/ 2.d0 * rdx2
        be  = -3.d0/20.d0 * rdx2
        ce  =  1.d0/90.d0 * rdx2
!
!  Internal nodes:
!
        do k = ikL, ikR
          do j = ijL, ijR
            do i = xbegin, xend
              d2f(i,j,k) = Up *  f(i  ,j,k)
     &                   + ae *( f(i+1,j,k)+f(i-1,j,k) )
     &                   + be *( f(i+2,j,k)+f(i-2,j,k) )
     &                   + ce *( f(i+3,j,k)+f(i-3,j,k) ) 
            end do    
          end do    
        end do    
!
!  Boundary nodes:
!
!   Left side:
!
        if ( nproLx .ge. 1 ) then
          do k = ikL, ikR
            do j = ijL, ijR
              d2f(xbegin-1,j,k) = (bnd4C(1)*f(xbegin-3,j,k) +
     &                             bnd4C(2)*f(xbegin-2,j,k) +
     &                             bnd4C(3)*f(xbegin-1,j,k) +
     &                             bnd4C(4)*f(xbegin+0,j,k) +
     &                             bnd4C(5)*f(xbegin+1,j,k))*rdx2
            end do
          end do
        end if
!
        if ( nproLx .ge. 2 ) then
          do k = ikL, ikR
            do j = ijL, ijR
              d2f(xbegin-2,j,k) = (bnd3B(1)*f(xbegin-3,j,k) +
     &                             bnd3B(2)*f(xbegin-2,j,k) +
     &                             bnd3B(3)*f(xbegin-1,j,k) +
     &                             bnd3B(4)*f(xbegin+0,j,k) +
     &                             bnd3B(5)*f(xbegin+1,j,k))*rdx2
            end do
          end do
        end if
!
        if ( nproLx .eq. 3 ) then
          do k = ikL, ikR
            do j = ijL, ijR
              d2f(xbegin-3,j,k) = (bnd3A(1)*f(xbegin-3,j,k) +
     &                             bnd3A(2)*f(xbegin-2,j,k) +
     &                             bnd3A(3)*f(xbegin-1,j,k) +
     &                             bnd3A(4)*f(xbegin+0,j,k) +
     &                             bnd3A(5)*f(xbegin+1,j,k))*rdx2
            end do
          end do
        end if
!
!   Right side:
!
        if ( nproRx .ge. 1 ) then
          do k = ikL, ikR
            do j = ijL, ijR
              d2f(xend+1,j,k) = (bnd4C(1)*f(xend+3,j,k) +
     &                           bnd4C(2)*f(xend+2,j,k) +
     &                           bnd4C(3)*f(xend+1,j,k) +
     &                           bnd4C(4)*f(xend+0,j,k) +
     &                           bnd4C(5)*f(xend-1,j,k))*rdx2
            end do
          end do
        end if
!
        if ( nproRx .ge. 2 ) then
          do k = ikL, ikR
            do j = ijL, ijR
              d2f(xend+2,j,k) = (bnd3B(1)*f(xend+3,j,k) +
     &                           bnd3B(2)*f(xend+2,j,k) +
     &                           bnd3B(3)*f(xend+1,j,k) +
     &                           bnd3B(4)*f(xend+0,j,k) +
     &                           bnd3B(5)*f(xend-1,j,k))*rdx2
            end do
          end do
        end if
!
        if ( nproRx .eq. 3 ) then
          do k = ikL, ikR
            do j = ijL, ijR
              d2f(xend+3,j,k) = (bnd3A(1)*f(xend+3,j,k) +
     &                           bnd3A(2)*f(xend+2,j,k) +
     &                           bnd3A(3)*f(xend+1,j,k) +
     &                           bnd3A(4)*f(xend+0,j,k) +
     &                           bnd3A(5)*f(xend-1,j,k))*rdx2
            end do
          end do
        end if
      end if
!-----------------------------
!  8th order explicit: (3,3,4,6-8E-6,4,3,3)
!
!
      if ( orderi .eq. 8 ) then
!
        Up  =-205.d0/ 72.d0 * rdx2
        ae  =   8.d0/  5.d0 * rdx2
        be  =  -1.d0/  5.d0 * rdx2
        ce  =   8.d0/315.d0 * rdx2
        de  =  -1.d0/560.d0 * rdx2
!
!  Internal nodes:
!
        do k = ikL, ikR
          do j = ijL, ijR
            do i = xbegin, xend
              d2f(i,j,k) = Up *  f(i  ,j,k)
     &                   + ae *( f(i+1,j,k)+f(i-1,j,k) )
     &                   + be *( f(i+2,j,k)+f(i-2,j,k) )
     &                   + ce *( f(i+3,j,k)+f(i-3,j,k) ) 
     &                   + de *( f(i+4,j,k)+f(i-4,j,k) ) 
            end do    
          end do    
        end do    
!
!  Boundary nodes:
!
!   Left side:
!
        if ( nproLx .ge. 1 ) then
          do k = ikL, ikR
            do j = ijL, ijR
              d2f(xbegin-1,j,k) = (bnd6D(1)*f(xbegin-4,j,k) +
     &                             bnd6D(2)*f(xbegin-3,j,k) +
     &                             bnd6D(3)*f(xbegin-2,j,k) +
     &                             bnd6D(4)*f(xbegin-1,j,k) +
     &                             bnd6D(5)*f(xbegin+0,j,k) +
     &                             bnd6D(6)*f(xbegin+1,j,k) +
     &                             bnd6D(7)*f(xbegin+2,j,k))*rdx2
            end do
          end do
        end if

        if ( nproLx .ge. 2 ) then
          do k = ikL, ikR
            do j = ijL, ijR
              d2f(xbegin-2,j,k) = (bnd4C(1)*f(xbegin-4,j,k) +
     &                             bnd4C(2)*f(xbegin-3,j,k) +
     &                             bnd4C(3)*f(xbegin-2,j,k) +
     &                             bnd4C(4)*f(xbegin-1,j,k) +
     &                             bnd4C(5)*f(xbegin+0,j,k))*rdx2
            end do
          end do
        end if
!
        if ( nproLx .ge. 3 ) then
          do k = ikL, ikR
            do j = ijL, ijR
              d2f(xbegin-3,j,k) = (bnd3B(1)*f(xbegin-4,j,k) +
     &                             bnd3B(2)*f(xbegin-3,j,k) +
     &                             bnd3B(3)*f(xbegin-2,j,k) +
     &                             bnd3B(4)*f(xbegin-1,j,k) +
     &                             bnd3B(5)*f(xbegin+0,j,k))*rdx2
            end do
          end do
        end if
!
        if ( nproLx .eq. 4 ) then
          do k = ikL, ikR
            do j = ijL, ijR
              d2f(xbegin-4,j,k) = (bnd3A(1)*f(xbegin-4,j,k) +
     &                             bnd3A(2)*f(xbegin-3,j,k) +
     &                             bnd3A(3)*f(xbegin-2,j,k) +
     &                             bnd3A(4)*f(xbegin-1,j,k) +
     &                             bnd3A(5)*f(xbegin+0,j,k))*rdx2
            end do
          end do
        end if
!
!   Right side:
!
        if ( nproRx .ge. 1 ) then
          do k = ikL, ikR
            do j = ijL, ijR
              d2f(xend+1,j,k) = (bnd6D(1)*f(xend+4,j,k) +
     &                           bnd6D(2)*f(xend+3,j,k) +
     &                           bnd6D(3)*f(xend+2,j,k) +
     &                           bnd6D(4)*f(xend+1,j,k) +
     &                           bnd6D(5)*f(xend+0,j,k) +
     &                           bnd6D(6)*f(xend-1,j,k) +
     &                           bnd6D(7)*f(xend-2,j,k))*rdx2
            end do
          end do
        end if
!
        if ( nproRx .ge. 2 ) then
          do k = ikL, ikR
            do j = ijL, ijR
              d2f(xend+2,j,k) = (bnd4C(1)*f(xend+4,j,k) +
     &                           bnd4C(2)*f(xend+3,j,k) +
     &                           bnd4C(3)*f(xend+2,j,k) +
     &                           bnd4C(4)*f(xend+1,j,k) +
     &                           bnd4C(5)*f(xend+0,j,k))*rdx2
            end do
          end do
        end if
!
        if ( nproRx .ge. 3 ) then
          do k = ikL, ikR
            do j = ijL, ijR
              d2f(xend+3,j,k) = (bnd3B(1)*f(xend+4,j,k) +
     &                           bnd3B(2)*f(xend+3,j,k) +
     &                           bnd3B(3)*f(xend+2,j,k) +
     &                           bnd3B(4)*f(xend+1,j,k) +
     &                           bnd3B(5)*f(xend+0,j,k))*rdx2
            end do
          end do
        end if
!
        if ( nproRx .eq. 4 ) then
          do k = ikL, ikR
            do j = ijL, ijR
              d2f(xend+4,j,k) = (bnd3A(1)*f(xend+4,j,k) +
     &                           bnd3A(2)*f(xend+3,j,k) +
     &                           bnd3A(3)*f(xend+2,j,k) +
     &                           bnd3A(4)*f(xend+1,j,k) +
     &                           bnd3A(5)*f(xend+0,j,k))*rdx2
            end do
          end do
        end if
!
      end if
!
! End Centered Difference Operators
!
      return
      end
!====================================================================