!====================================================================
!
      subroutine z_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,rdz2,rdz2)
!
      if( 2*promax+1 .gt. pkR-pkL+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)* rdz2
        ae  = ( 1.d0/ 1.d0)* rdz2
!
!  Internal nodes:
!
        do k = zbegin, zend
          do j = ijL, ijR
            do i = iiL, iiR
              d2f(i,j,k) = Up *  f(i,j,k+0)
     &                   + ae *( f(i,j,k+1)+f(i,j,k-1) )
            end do   
          end do
        end do
!
!  Boundary nodes:
!   Left side:
!
        if ( nproLz .eq. 1 ) then
          do j = ijL, ijR
            do i = iiL, iiR
              d2f(i,j,zbegin-1) = (bnd2A(1)*f(i,j,zbegin-1) +
     &                             bnd2A(2)*f(i,j,zbegin+0) +
     &                             bnd2A(3)*f(i,j,zbegin+1))*rdz2
            end do
          end do
        end if
!
!   Right side:
!
        if ( nproRz .eq. 1 ) then
          do j = ijL, ijR
            do i = iiL, iiR
              d2f(i,j,zend+1) = (bnd2A(1)*f(i,j,zend+1) +
     &                           bnd2A(2)*f(i,j,zend+0) +
     &                           bnd2A(3)*f(i,j,zend-1))*rdz2
            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 * rdz2
        ae  =  4.d0/ 3.d0 * rdz2
        be  = -1.d0/12.d0 * rdz2
!
!  Internal nodes:
!
        do k = zbegin, zend
          do j = ijL, ijR
            do i = iiL, iiR
              d2f(i,j,k) = Up *  f(i,j,k+0)
     &                   + ae *( f(i,j,k+1)+f(i,j,k-1) )
     &                   + be *( f(i,j,k+2)+f(i,j,k-2) )
            end do   
          end do   
        end do
!
!  Boundary nodes:
!
!   Left side:
!
!
        if ( nproLz .ge. 1 ) then
          do j = ijL, ijR
            do i = iiL, iiR
              d2f(i,j,zbegin-1) = (bnd3B(1)*f(i,j,zbegin-2) +
     &                             bnd3B(2)*f(i,j,zbegin-1) +
     &                             bnd3B(3)*f(i,j,zbegin+0) +
     &                             bnd3B(4)*f(i,j,zbegin+1) +
     &                             bnd3B(5)*f(i,j,zbegin+2))*rdz2
            end do
          end do
        end if
!
        if ( nproLz .eq. 2 ) then
          do j = ijL, ijR
            do i = iiL, iiR
              d2f(i,j,zbegin-2) = (bnd3A(1)*f(i,j,zbegin-2) +
     &                             bnd3A(2)*f(i,j,zbegin-1) +
     &                             bnd3A(3)*f(i,j,zbegin+0) +
     &                             bnd3A(4)*f(i,j,zbegin+1) +
     &                             bnd3A(5)*f(i,j,zbegin+2))*rdz2
            end do
          end do
        end if
!
!   Right side:
!
        if ( nproRz .ge. 1 ) then
          do j = ijL, ijR
            do i = iiL, iiR
              d2f(i,j,zend+1) = (bnd3B(1)*f(i,j,zend+2) +
     &                           bnd3B(2)*f(i,j,zend+1) +
     &                           bnd3B(3)*f(i,j,zend+0) +
     &                           bnd3B(4)*f(i,j,zend-1) +
     &                           bnd3B(5)*f(i,j,zend-2))*rdz2
            end do
          end do
        end if
!
        if ( nproRz .eq. 2 ) then
          do j = ijL, ijR
            do i = iiL, iiR
              d2f(i,j,zend+2) = (bnd3A(1)*f(i,j,zend+2) +
     &                           bnd3A(2)*f(i,j,zend+1) +
     &                           bnd3A(3)*f(i,j,zend+0) +
     &                           bnd3A(4)*f(i,j,zend-1) +
     &                           bnd3A(5)*f(i,j,zend-2))*rdz2
            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 * rdz2
        ae  =  3.d0/ 2.d0 * rdz2
        be  = -3.d0/20.d0 * rdz2
        ce  =  1.d0/90.d0 * rdz2
!
!  Internal nodes:
!
        do k = zbegin, zend
          do j = ijL, ijR
            do i = iiL, iiR
              d2f(i,j,k) = Up *  f(i,j,k+0)
     &                   + ae *( f(i,j,k+1)+f(i,j,k-1) )
     &                   + be *( f(i,j,k+2)+f(i,j,k-2) )
     &                   + ce *( f(i,j,k+3)+f(i,j,k-3) ) 
            end do    
          end do    
        end do    
!
!  Boundary nodes:
!
!   Left side:
!
        if ( nproLz .ge. 1 ) then
          do j = ijL, ijR
            do i = iiL, iiR
              d2f(i,j,zbegin-1) = (bnd4C(1)*f(i,j,zbegin-3) +
     &                             bnd4C(2)*f(i,j,zbegin-2) +
     &                             bnd4C(3)*f(i,j,zbegin-1) +
     &                             bnd4C(4)*f(i,j,zbegin+0) +
     &                             bnd4C(5)*f(i,j,zbegin+1))*rdz2
            end do
          end do
        end if
!
        if ( nproLz .ge. 2 ) then
          do j = ijL, ijR
            do i = iiL, iiR
              d2f(i,j,zbegin-2) = (bnd3B(1)*f(i,j,zbegin-3) +
     &                             bnd3B(2)*f(i,j,zbegin-2) +
     &                             bnd3B(3)*f(i,j,zbegin-1) +
     &                             bnd3B(4)*f(i,j,zbegin+0) +
     &                             bnd3B(5)*f(i,j,zbegin+1))*rdz2
            end do
          end do
        end if
!
        if ( nproLz .eq. 3 ) then
          do j = ijL, ijR
            do i = iiL, iiR
              d2f(i,j,zbegin-3) = (bnd3A(1)*f(i,j,zbegin-3) +
     &                             bnd3A(2)*f(i,j,zbegin-2) +
     &                             bnd3A(3)*f(i,j,zbegin-1) +
     &                             bnd3A(4)*f(i,j,zbegin+0) +
     &                             bnd3A(5)*f(i,j,zbegin+1))*rdz2
            end do
          end do
        end if
!
!   Right side:
!
        if ( nproRz .ge. 1 ) then
          do j = ijL, ijR
            do i = iiL, iiR
              d2f(i,j,zend+1) = (bnd4C(1)*f(i,j,zend+3) +
     &                           bnd4C(2)*f(i,j,zend+2) +
     &                           bnd4C(3)*f(i,j,zend+1) +
     &                           bnd4C(4)*f(i,j,zend+0) +
     &                           bnd4C(5)*f(i,j,zend-1))*rdz2
            end do
          end do
        end if
!
        if ( nproRz .ge. 2 ) then
          do j = ijL, ijR
            do i = iiL, iiR
              d2f(i,j,zend+2) = (bnd3B(1)*f(i,j,zend+3) +
     &                           bnd3B(2)*f(i,j,zend+2) +
     &                           bnd3B(3)*f(i,j,zend+1) +
     &                           bnd3B(4)*f(i,j,zend+0) +
     &                           bnd3B(5)*f(i,j,zend-1))*rdz2
            end do
          end do
        end if
!
        if ( nproRz .eq. 3 ) then
          do j = ijL, ijR
            do i = iiL, iiR
              d2f(i,j,zend+3) = (bnd3A(1)*f(i,j,zend+3) +
     &                           bnd3A(2)*f(i,j,zend+2) +
     &                           bnd3A(3)*f(i,j,zend+1) +
     &                           bnd3A(4)*f(i,j,zend+0) +
     &                           bnd3A(5)*f(i,j,zend-1))*rdz2
            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 * rdz2
        ae  =   8.d0/  5.d0 * rdz2
        be  =  -1.d0/  5.d0 * rdz2
        ce  =   8.d0/315.d0 * rdz2
        de  =  -1.d0/560.d0 * rdz2
!
!  Internal nodes:
!
        do k = zbegin, zend
          do j = ijL, ijR
            do i = iiL, iiR
              d2f(i,j,k) = Up *  f(i,j,k+0)
     &                   + ae *( f(i,j,k+1)+f(i,j,k-1) )
     &                   + be *( f(i,j,k+2)+f(i,j,k-2) )
     &                   + ce *( f(i,j,k+3)+f(i,j,k-3) ) 
     &                   + de *( f(i,j,k+4)+f(i,j,k-4) ) 
            end do    
          end do    
        end do    
!
!  Boundary nodes:
!
!   Left side:
!
        if ( nproLz .ge. 1 ) then
          do j = ijL, ijR
            do i = iiL, iiR
              d2f(i,j,zbegin-1) = (bnd6D(1)*f(i,j,zbegin-4) +
     &                             bnd6D(2)*f(i,j,zbegin-3) +
     &                             bnd6D(3)*f(i,j,zbegin-2) +
     &                             bnd6D(4)*f(i,j,zbegin-1) +
     &                             bnd6D(5)*f(i,j,zbegin+0) +
     &                             bnd6D(6)*f(i,j,zbegin+1) +
     &                             bnd6D(7)*f(i,j,zbegin+2))*rdz2
            end do
          end do
        end if

        if ( nproLz .ge. 2 ) then
          do j = ijL, ijR
            do i = iiL, iiR
              d2f(i,j,zbegin-2) = (bnd4C(1)*f(i,j,zbegin-4) +
     &                             bnd4C(2)*f(i,j,zbegin-3) +
     &                             bnd4C(3)*f(i,j,zbegin-2) +
     &                             bnd4C(4)*f(i,j,zbegin-1) +
     &                             bnd4C(5)*f(i,j,zbegin+0))*rdz2
            end do
          end do
        end if
!
        if ( nproLz .ge. 3 ) then
          do j = ijL, ijR
            do i = iiL, iiR
              d2f(i,j,zbegin-3) = (bnd3B(1)*f(i,j,zbegin-4) +
     &                             bnd3B(2)*f(i,j,zbegin-3) +
     &                             bnd3B(3)*f(i,j,zbegin-2) +
     &                             bnd3B(4)*f(i,j,zbegin-1) +
     &                             bnd3B(5)*f(i,j,zbegin+0))*rdz2
            end do
          end do
        end if
!
        if ( nproLz .eq. 4 ) then
          do j = ijL, ijR
            do i = iiL, iiR
              d2f(i,j,zbegin-4) = (bnd3A(1)*f(i,j,zbegin-4) +
     &                             bnd3A(2)*f(i,j,zbegin-3) +
     &                             bnd3A(3)*f(i,j,zbegin-2) +
     &                             bnd3A(4)*f(i,j,zbegin-1) +
     &                             bnd3A(5)*f(i,j,zbegin+0))*rdz2
            end do
          end do
        end if
!
!   Right side:
!
        if ( nproRz .ge. 1 ) then
          do j = ijL, ijR
            do i = iiL, iiR
              d2f(i,j,zend+1) = (bnd6D(1)*f(i,j,zend+4) +
     &                           bnd6D(2)*f(i,j,zend+3) +
     &                           bnd6D(3)*f(i,j,zend+2) +
     &                           bnd6D(4)*f(i,j,zend+1) +
     &                           bnd6D(5)*f(i,j,zend+0) +
     &                           bnd6D(6)*f(i,j,zend-1) +
     &                           bnd6D(7)*f(i,j,zend-2))*rdz2
            end do
          end do
        end if
!
        if ( nproRz .ge. 2 ) then
          do j = ijL, ijR
            do i = iiL, iiR
              d2f(i,j,zend+2) = (bnd4C(1)*f(i,j,zend+4) +
     &                           bnd4C(2)*f(i,j,zend+3) +
     &                           bnd4C(3)*f(i,j,zend+2) +
     &                           bnd4C(4)*f(i,j,zend+1) +
     &                           bnd4C(5)*f(i,j,zend+0))*rdz2
            end do
          end do
        end if
!
        if ( nproRz .ge. 3 ) then
          do j = ijL, ijR
            do i = iiL, iiR
              d2f(i,j,zend+3) = (bnd3B(1)*f(i,j,zend+4) +
     &                           bnd3B(2)*f(i,j,zend+3) +
     &                           bnd3B(3)*f(i,j,zend+2) +
     &                           bnd3B(4)*f(i,j,zend+1) +
     &                           bnd3B(5)*f(i,j,zend+0))*rdz2
            end do
          end do
        end if
!
        if ( nproRz .eq. 4 ) then
          do j = ijL, ijR
            do i = iiL, iiR
              d2f(i,j,zend+4) = (bnd3A(1)*f(i,j,zend+4) +
     &                           bnd3A(2)*f(i,j,zend+3) +
     &                           bnd3A(3)*f(i,j,zend+2) +
     &                           bnd3A(4)*f(i,j,zend+1) +
     &                           bnd3A(5)*f(i,j,zend+0))*rdz2
            end do
          end do
        end if
!
      end if
!
! End Centered Difference Operators
!
      return
      end
!====================================================================