source: trunk/00FlowSolve_PL/PROJECTS/NEW_SRC/compact_derivs6.f90 @ 4

subroutine compact_ddx(data,deriv,flag,end_values,m1,m2,m3)
! Routine to compute the x derivative (1st dimension)
! of a 3d field stored in data(1:m1,1:m2,1:m3) assuming
! a uniformly spaced grid mesh with dx=h.
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!       inputs
!              data        m1,m2,m3 array of input values
!              m1,m2,m3    dimensions of field to be differentiated
!              h           x grid spacing
!              flag        flag which identifies periodicity
!              end_values  the end values for the derivatives
!                          (required but now obsolete)
!
!       output
!              deriv       m1,m2,m3 array containing the optimized
!                          6th order accurate d/dx(data)
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 use grid_info
 use boundary_information
use user_parameters,        only: LAPACK_FLAG
 implicit none

 integer, intent(in)                              :: flag,m1,m2,m3
 real, dimension(m1,m2,m3)                        :: data
 real, dimension(m1,m2,m3)                        :: deriv
 real, dimension(2)                               :: end_values
 real, dimension(:,:), allocatable    :: work
 integer                                          :: i,j,k,ierr,level,maxlevels
 real                                             :: h
 parameter(maxlevels=12)

 real, allocatable                                :: kx(:),in(:),out(:)
 real                                             :: pi
 integer*8                                        :: plan_f, plan_i !size must match pointer size, 8 for DEC alphas
 integer                                          :: N
 include 'fftw_kw.inc'                            !  predefined FFTW constants

 if( boundary_type2(1,1)=='periodic' ) then  ! spectral differentiation in x, use FFTW

! preliminaries, setup fftw "plans"
  pi=4.*atan(1.)
  N=m1-1                                     ! m1 includes both boundary points, N only 1
  allocate( kx(N),in(N),out(N) )

  call rfftw_f77_create_plan(plan_f,N,FFTW_REAL_TO_COMPLEX,FFTW_ESTIMATE)
  call rfftw_f77_create_plan(plan_i,N,FFTW_COMPLEX_TO_REAL,FFTW_ESTIMATE)
  do i=1,N/2+1
   kx(i)=2.*pi*(i-1.)    ! FFTW real to half-complex storage
  enddo
  do i=2,N/2
   kx(N-i+2)=kx(i)       ! FFTW real to half-complex storage
  enddo

  do j=1,m2
   do k=1,m3

    in(1:N)=data(1:N,j,k)
    call rfftw_f77_one(plan_f,in,out)   ! FORWARD TRANSFORM

!   multiply by wavenumber array & normalize by N
    out(1:N) = kx(1:N)*out(1:N)/float(N)
!   multiply by i  (fold and negate)
    do i=2,N/2
     in(i)=-out(N-i+2)  ! real <-- (-imaginary)
     in(N-i+2)=out(i)   ! imag <-- ( real )
    enddo
    in(1)=0.     ! dc
    in(N/2+1)=0. !nyquist

    call rfftw_f77_one(plan_i,in,out)   ! INVERSE TRANSFORM
    deriv(1:N,j,k)=out(1:N)             ! store results in deriv
    deriv(m1,j,k)=out(1)                ! carry 2nd periodic end point along
    
   enddo
  enddo

  call rfftw_f77_destroy_plan(plan_f)
  call rfftw_f77_destroy_plan(plan_i)
  deallocate( kx,in,out )
  return
 
 endif

!*************  nonperiodic, use compact differentiation ********************
 allocate( work(m1,m2) )
 
! determine an appropriate grid level from m1
 do level=0,maxlevels
  if( m1 == Grid(level)%n1 ) goto 999
 enddo
999 continue   ! can use ddx and piv arrays from grid level 
 h = Grid(level)%dxi

!loop through each xy horizontal plane
       
 do k=1,m3

!Create rhs vectors, store in work(:,:)
 work(1,:)=( -(274.)*data(1,:,k)+(600.)*data(2,:,k)   &
             -(600.)*data(3,:,k)+(400.)*data(4,:,k)   &
             -(150.)*data(5,:,k)+ (24.)*data(6,:,k) )/(120.*h)   
!               work(2,:)=( -(5./9.)*data(1,:,k)-(1./2.)*data(2,:,k)   &
!                           +(1.)*data(3,:,k)+(1./18.)*data(4,:,k) )/h
 work(2,:)=( -(24.)*data(1,:,k)-(130.)*data(2,:,k)   &
             +(240.)*data(3,:,k)-(120.)*data(4,:,k)   &
             +(40.)*data(5,:,k)- (6.)*data(6,:,k) )/(120.*h)   
 work(3,:)= ( (1./36.)*(data(5,:,k)-data(1,:,k))   &
             + (7./9.)*(data(4,:,k)-data(2,:,k)) )/h  

 do i=4,m1-3
  work(i,:)= ( (c1/6.)*(data(i+3,:,k)-data(i-3,:,k))   &
              + (b1/4.)*(data(i+2,:,k)-data(i-2,:,k))   & 
              + (a1/2.)*(data(i+1,:,k)-data(i-1,:,k)) )/h  
 enddo

 work(m1-2,:)= (-(1./36.)*(data(m1-4,:,k)-data(m1,:,k))   &
                - (7./9.)*(data(m1-3,:,k)-data(m1-1,:,k)) )/h  
!         work(m1-1,:)=( (5./9.)*data(m1,:,k)+(1./2.)*data(m1-1,:,k)  &
!                         -(1)*data(m1-2,:,k)-(1./18.)*data(m1-3,:,k) )/h
 work(m1-1,:)=( (24.)*data(m1,:,k)+(130.)*data(m1-1,:,k)   &
             -(240.)*data(m1-2,:,k)+(120.)*data(m1-3,:,k)   &
             -(40.)*data(m1-4,:,k)+ (6.)*data(m1-5,:,k) )/(120.*h)   
 work(m1,:)=( (274.)*data(m1,:,k)-(600.)*data(m1-1,:,k)   &
             +(600.)*data(m1-2,:,k)-(400.)*data(m1-3,:,k)   &
             +(150.)*data(m1-4,:,k)- (24.)*data(m1-5,:,k) )/(120.*h)   

!Compute the x derivatives for plane k, overwriting the array work

 if( LAPACK_FLAG == 'double' ) then
!!!  call DGBTRS('N',m1,2,2,m2,Grid(level)%ddx,7,Grid(level)%piv_dx,work(1,1),m1,ierr)
      write(0,*) 'LAPACK_FLAG = double instead of single'
      stop
 elseif( LAPACK_FLAG == 'single' ) then
  call SGBTRS('N',m1,2,2,m2,Grid(level)%ddx,7,Grid(level)%piv_dx,work(1,1),m1,ierr)
 endif

!If deriv. end conditions are 'periodic' (flag = 5) then
!overwrite the endpts with nonlocal estimate
 if( flag .eq. 5 ) then

  work(1,:) = ( 2.*data(m1-2,:,k) -16.*data(m1-1,:,k)      &
              +16.*data(2,:,k) -2.*data(3,:,k)  )/(24.*h)  
  work(m1,:)=work(1,:)

  work(2,:) = ( 2.*data(m1-1,:,k) -16.*data(m1,:,k)      &
              +16.*data(3,:,k) -2.*data(4,:,k)  )/(24.*h)  
  work(m1-1,:) = ( 2.*data(m1-3,:,k) -16.*data(m1-2,:,k)      &
              +16.*data(1,:,k) -2.*data(2,:,k)  )/(24.*h)  

  work(3,:) = ( 2.*data(m1,:,k) -16.*data(2,:,k)      &
              +16.*data(4,:,k) -2.*data(5,:,k)  )/(24.*h)  
  work(m1-2,:) = ( 2.*data(m1-4,:,k) -16.*data(m1-3,:,k)      &
              +16.*data(m1-1,:,k) -2.*data(m1,:,k)  )/(24.*h)  

 endif
      
!store the results in deriv(:,:,k)
  deriv(:,:,k)=work(:,:) 

!Go on to the next k plane
 enddo

 deallocate( work  )
end subroutine compact_ddx





subroutine compact_ddy(data,deriv,flag,end_values,m1,m2,m3)
! Routine to compute the y derivative (2nd dimension)
! of a 3d field stored in data(1:m1,1:m2,1:m3) assuming
! a uniformly spaced grid mesh with dy=h.
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!       inputs
!              data        m1,m2,m3 array of input values
!              m1,m2,m3    dimensions of field to be differentiated
!              h           y grid spacing
!              flag        flag which identifies periodicity
!              end_values  the end values for the derivatives
!                          (required but obsolete)
!
!       output
!              deriv       m1,m2,m3 array containing the optimized
!                          6th order accurate d/dy(data)
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 use grid_info
 use boundary_information
use user_parameters,        only: LAPACK_FLAG
 implicit none

 integer, intent(in)                              :: flag,m1,m2,m3
 real, dimension(m1,m2,m3), intent(in)            :: data
 real, dimension(m1,m2,m3), intent(out)           :: deriv
 real, dimension(2), intent(in)                   :: end_values
!!! double precision, dimension(:,:), allocatable    :: work
 real, dimension(:,:), allocatable    :: work
 integer                                          :: i,j,k,ierr,level,maxlevels
 real                                             :: h
 parameter(maxlevels=12)

 real, allocatable                                :: ky(:),in(:),out(:)
 real                                             :: pi
 integer*8                                        :: plan_f, plan_i !size must match pointer size, 8 for DEC alphas
 integer                                          :: N
 include 'fftw_kw.inc'                            !  predefined FFTW constants

 if( boundary_type4(1,1)=='periodic' ) then  ! spectral differentiation in y, use FFTW

! preliminaries, setup fftw "plans"
  pi=4.*atan(1.)
  N=m2-1                                     ! m2 includes both boundary points, N only 1
  allocate( ky(N),in(N),out(N) )

  call rfftw_f77_create_plan(plan_f,N,FFTW_REAL_TO_COMPLEX,FFTW_ESTIMATE)
  call rfftw_f77_create_plan(plan_i,N,FFTW_COMPLEX_TO_REAL,FFTW_ESTIMATE)
  do j=1,N/2+1
   ky(j)=2.*pi*(j-1.)    ! FFTW real to half-complex storage
  enddo
  do j=2,N/2
   ky(N-j+2)=ky(j)      ! FFTW real to half-complex storage
  enddo

  do i=1,m1
   do k=1,m3

    in(1:N)=data(i,1:N,k)
    call rfftw_f77_one(plan_f,in,out)   ! FORWARD TRANSFORM

!   multiply by wavenumber array & normalize by N
    out(1:N) = ky(1:N)*out(1:N)/float(N)
!   multiply by i  (fold and negate)
    do j=2,N/2
     in(j)=-out(N-j+2)  ! real <-- (-imaginary)
     in(N-j+2)=out(j)   ! imag <-- ( real )
    enddo
    in(1)=0.     ! dc
    in(N/2+1)=0. !nyquist

    call rfftw_f77_one(plan_i,in,out)   ! INVERSE TRANSFORM
    deriv(i,1:N,k)=out(1:N)             ! store results in deriv
    deriv(i,m2,k)=out(1)                ! carry 2nd periodic end point along
    
   enddo
  enddo

  call rfftw_f77_destroy_plan(plan_f)
  call rfftw_f77_destroy_plan(plan_i)
  deallocate( ky,in,out )
  return
  
 endif

!*************  nonperiodic, use compact differentiation ********************
 allocate( work(m2,m1) )

! determine an appropriate grid level from m2
 do level=0,maxlevels-1
  if( m2 == Grid(level)%n2 ) goto 999
 enddo
999 continue   ! can use ddx and piv arrays from grid level 
 h = Grid(level)%deta

!loop through each xy horizontal plane
       
 do k=1,m3

!      create rhs vectors, store in work(:,:,1)
 work(1,:)=( -(274.)*data(:,1,k)+(600.)*data(:,2,k)   &
             -(600.)*data(:,3,k)+(400.)*data(:,4,k)   &
             -(150.)*data(:,5,k)+ (24.)*data(:,6,k) )/(120.*h)   
!         work(2,:)=( -(5./9.)*data(:,1,k)-(1./2.)*data(:,2,k)   &
!                       +(1.)*data(:,3,k)+(1./18.)*data(:,4,k) )/h
 work(2,:)=( -(24.)*data(:,1,k)-(130.)*data(:,2,k)   &
             +(240.)*data(:,3,k)-(120.)*data(:,4,k)   &
             +(40.)*data(:,5,k)- (6.)*data(:,6,k) )/(120.*h)   
 work(3,:) = ( (1./36.)*( data(:,5,k) - data(:,1,k) )   &
              + (7./9.)*( data(:,4,k) - data(:,2,k) ) )/h  
 do j=4,m2-3
  work(j,:) = (  (c1/6.)*( data(:,j+3,k) - data(:,j-3,k) )   &
               + (b1/4.)*( data(:,j+2,k) - data(:,j-2,k) )   &
               + (a1/2.)*( data(:,j+1,k) - data(:,j-1,k) ) )/h  
 enddo

 work(m2-2,:) = (-(1./36.)*( data(:,m2-4,k) - data(:,m2,k) )   &
              - (7./9.)*( data(:,m2-3,k) - data(:,m2-1,k) ) )/h  
!         work(m2-1,:)= ( (5./9.)*data(:,m2,k)+(1./2.)*data(:,m2-1,k)   &
!                          -(1.)*data(:,m2-2,k)-(1./18.)*data(:,m2-3,k) )/h
 work(m2-1,:)=( (24.)*data(:,m2,k)+(130.)*data(:,m2-1,k)   &
             -(240.)*data(:,m2-2,k)+(120.)*data(:,m2-3,k)   &
             -(40.)*data(:,m2-4,k)+ (6.)*data(:,m2-5,k) )/(120.*h)   
 work(m2,:)=(  (274.)*data(:,m2,k)-(600.)*data(:,m2-1,k)   &
             +(600.)*data(:,m2-2,k)-(400.)*data(:,m2-3,k)   &
             +(150.)*data(:,m2-4,k)- (24.)*data(:,m2-5,k) )/(120.*h)   

!Compute the y derivatives for plane k, overwriting the array work
 if( LAPACK_FLAG == 'double' ) then
!!!  call DGBTRS('N',m2,2,2,m1,Grid(level)%ddy,7,Grid(level)%piv_dy,work(1,1),m2,ierr)
      write(0,*) 'LAPACK_FLAG = double instead of single'
      stop
 elseif( LAPACK_FLAG == 'single' ) then
  call SGBTRS('N',m2,2,2,m1,Grid(level)%ddy,7,Grid(level)%piv_dy,work(1,1),m2,ierr)
 endif

!If deriv. end conditions are 'periodic' (flag = 5) then
!overwrite the endpts with nonlocal estimate
 if( flag .eq. 5 ) then
  work(1,:) = ( 2.*data(:,m2-2,k) -16.*data(:,m2-1,k)      &
              +16.*data(:,2,k) -2.*data(:,3,k)  )/(24.*h)  
  work(m2,:)=work(1,:)

  work(2,:) = ( 2.*data(:,m2-1,k) -16.*data(:,m2,k)      &
              +16.*data(:,3,k) -2.*data(:,4,k)  )/(24.*h)  
  work(m2-1,:) = ( 2.*data(:,m2-3,k) -16.*data(:,m2-2,k)      &
              +16.*data(:,1,k) -2.*data(:,2,k)  )/(24.*h)  

  work(3,:) = ( 2.*data(:,m2,k) -16.*data(:,2,k)      &
              +16.*data(:,4,k) -2.*data(:,5,k)  )/(24.*h)  
  work(m2-2,:) = ( 2.*data(:,m2-4,k) -16.*data(:,m2-3,k)      &
              +16.*data(:,m2-1,k) -2.*data(:,m2,k)  )/(24.*h)  

 endif
      
!Store the results in deriv(:,:,k)
 deriv(1:m1,1:m2,k)=transpose( work(1:m2,1:m1)  )

!Go on to the next k plane
 enddo

 deallocate( work )
end subroutine compact_ddy




subroutine compact_ddz(data,deriv,flag,end_values,m1,m2,m3)
! Routine to compute the x derivative (3rd dimension)
! of a 3d field stored in data(1:m1,1:m2,1:m3) assuming
! a uniformly spaced grid mesh with dz=h.
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
!       inputs
!              data        m1,m2,m3 array of input values
!              m1,m2,m3    dimensions of field to be differentiated
!              h           z grid spacing
!              flag        flag which identifies end point scheme
!              end_values  the end values for the derivatives
!                          (required but obsolete)
!
!       output
!              deriv       m1,m2,m3 array containing the optimized
!                          6th order accurate d/dz(data)
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 use grid_info
 use user_parameters,  only: LAPACK_FLAG
 implicit none

 integer, intent(in)                              :: flag,m1,m2,m3
 real, dimension(m1,m2,m3), intent(in)            :: data
 real, dimension(m1,m2,m3), intent(out)           :: deriv
 real, dimension(2), intent(in)                   :: end_values
!!! double precision, dimension(:,:), allocatable    :: work
 real, dimension(:,:), allocatable    :: work
 integer                                          :: i,j,k,ierr,level,maxlevels
 real                                             :: h
 parameter(maxlevels=12)

 allocate( work(m3,m1) )

! determine an appropriate grid level from m3
 do level=0,maxlevels-1
  if( m3 == Grid(level)%n3 ) goto 999
 enddo
999 continue   ! can use ddx and piv arrays from grid level 
 h = Grid(level)%dzeta

!loop through each xz plane
       
 do j=1,m2

!create rhs vectors, store in work(:,:,1)
 work(1,:)=( -(274.)*data(:,j,1)+(600.)*data(:,j,2)   &
             -(600.)*data(:,j,3)+(400.)*data(:,j,4)   &
             -(150.)*data(:,j,5)+(24.)*data(:,j,6) )/(120.*h)   
!    work(2,:)=(-(5./9.)*data(:,j,1)-(1./2.)*data(:,j,2)   &
!                 +(1.)*data(:,j,3)+(1./18.)*data(:,j,4) )/h
 work(2,:)=( -(24.)*data(:,j,1)-(130.)*data(:,j,2)   &
             +(240.)*data(:,j,3)-(120.)*data(:,j,4)   &
             +(40.)*data(:,j,5)-(6.)*data(:,j,6) )/(120.*h)   
 work(3,:)=( (1./36.)* ( data(:,j,5) - data(:,j,1) )  &
            + (7./9.)* ( data(:,j,4) - data(:,j,2) ) )/h

 do k=4,m3-3
  work(k,:)=( (c1/6.)* ( data(:,j,k+3) - data(:,j,k-3) )  &
             + (b1/4.)* ( data(:,j,k+2) - data(:,j,k-2) )  &
             + (a1/2.)* ( data(:,j,k+1) - data(:,j,k-1) ) )/h
 enddo

 work(m3-2,:)=(-(1./36.)* ( data(:,j,m3-4) - data(:,j,m3) )  &
            - (7./9.)* ( data(:,j,m3-3) - data(:,j,m3-1) ) )/h
!    work(m3-1,:)= ( (5./9.)*data(:,j,m3)+(1./2.)*data(:,j,m3-1)   &
!                      -(1.)*data(:,j,m3-2)-(1./18.)*data(:,j,m3-3) )/h
 work(m3-1,:)=( (24.)*data(:,j,m3)+(130.)*data(:,j,m3-1)   &
             -(240.)*data(:,j,m3-2)+(120.)*data(:,j,m3-3)   &
             -(40.)*data(:,j,m3-4)+(6.)*data(:,j,m3-5) )/(120.*h)   
 work(m3,:)=( (274.)*data(:,j,m3)-(600.)*data(:,j,m3-1)   &
             +(600.)*data(:,j,m3-2)-(400.)*data(:,j,m3-3)   &
             +(150.)*data(:,j,m3-4)-(24.)*data(:,j,m3-5) )/(120.*h)   

!Compute the z derivatives for plane j, overwriting the array work
  if( LAPACK_FLAG == 'double' ) then
!!!   call DGBTRS('N',m3,2,2,m1,Grid(level)%ddz,7,Grid(level)%piv_dz,work(1,1),m3,ierr)
      write(0,*) 'LAPACK_FLAG = double instead of single'
      stop
  elseif( LAPACK_FLAG == 'single' ) then
!!!      write(0,*) '********** ddz = ', Grid(level)%ddz
!!!      write(0,*) '********** piv_dz = ', Grid(level)%piv_dz
   call SGBTRS('N',m3,2,2,m1,Grid(level)%ddz,7,Grid(level)%piv_dz,work(1,1),m3,ierr)
  endif

!!!  do i=1,m1
!!!  do k=1,m3
!!!     if(abs(work(k,i)) .lt. 1.e-20) work(k,i) = 0.0
!!!  enddo
!!!  enddo

!Store the results in deriv(:,:,k)
 do i=1,m1
  deriv(i,j,:)=work(:,i)
 enddo

!If deriv. end conditions are 'periodic' (flag = 5) then
!overwrite the endpts w/ the avg of the endpts
 if( flag .eq. 5 ) then
  work(1,:)=(deriv(:,j,1)+deriv(:,j,m3))/2.
  deriv(:,j,1)=work(1,:)
  deriv(:,j,m3)=work(1,:)
 endif
      
!Go on to the next j plane
 enddo

 deallocate ( work )
end subroutine compact_ddz





subroutine compact_ddz_mpi(data,deriv,flag,end_values,m1,m2,m3)
!     This version computes the z derivative of a block of 3d data.
!     There are three ways in which this routine can be used, depending
!     on the values of nz, m3 and nzderiv. The cases are denoted CASE1
!     CASE2 and CASE3 as described below.
!
!     CASE1: nprocs=1 and nzderiv=m3=nz
!            Serial execution. No communication calls necessary.
!            Derivatives implemented with single call to
!            compact_ddz with a vector length of nzderiv.

!     CASE2: nprocs>1 and nzderiv=m3 < nz
!            Parallel execution but no data exchanges between processors
!            required. Each processor computes z derivatives over its own
!            data space, i.e. over vector lengths of m3.
!            Derivatives implemented by having each processor execute
!            compact_ddz with a vector length of nzderiv.
!
!     CASE3: nprocs>1 and nzderiv=nz > m3. 
!            No processor possesses the required data locally and so
!            data exchanges between processors are required both before
!            and after differentiation. The algorithm is described below.            
!     On entry, each process possesses a "slice" of the total data space,
!     i.e. a block with dimensions (m1,m2,m3). Since the compact schemes
!     are nonlocal, i.e. all nz elements are required to compute a z derivative.
!     In this implementation, each process will compute z derivatives within
!     a "block of columns", i.e. a block of data of size (m1/nprocs,m2,nz)
!     or ( (m1-1)/nprocs,m2,n2 ) if ffts are being used in x&y
!     Since the basic data distribution consists of "blocks of rows", exchange
!     of data between processes is required. 
!
!     Each process lacks (nprocs-1) blocks of size (locnx,m2,m3)
!     needed to make up its "block of columns". On entry, each process
!     has (nprocs-1)  data "blocks" of size (locnx,m2,m3) that are
!     not needed locally but required by other processes to complete their
!     "block of columns".
!
!     This algorithm thus proceeds as follows:
!     (1) ship data not needed locally and receive data required to
!           complete a "block of columns".
!     (2) compute z derivatives over the processes "block of columns".
!     (3) send results to other processes, receive their results
!     (4) store the results in appropriate block locations of the original
!           "slice" of the total data space.
!*******************************************************
 use mpi_parameters
 use boundary_information
 use columnslab

 implicit none

 integer, intent(in)                              :: flag,m1,m2,m3
 real, dimension(m1,m2,m3), intent(in)            :: data
 real, dimension(m1,m2,m3), intent(out)           :: deriv
 real, dimension(2), intent(in)                   :: end_values
 integer                                          :: i,j,k,nzderiv,locnx

 include 'mpif.h'
 integer          :: tag,ierr,status_array(MPI_STATUS_SIZE,2*nprocs),reqs(2*nprocs)
 integer          :: dest,source,istart,iend,kstart,kend,nz
 integer          :: sizeofreal
 integer          :: iproc,numwords
 real, allocatable                      :: block_of_cols(:,:,:)

 if( boundary_type2(1,1)=='periodic' .and. boundary_type4(1,1)=='periodic' ) then
  locnx = (m1-1)/nprocs
 else
  locnx = m1/nprocs
 endif


 nz=m3*nprocs   ! global # of points
 nzderiv=nz     ! always do global differentiation

! Both CASE1 and CASE2 require simple calls to compact_ddz
  if( nprocs .eq. 1 .and. m3 .eq. nz .and. nzderiv .eq. nz ) then
   goto 999   ! CASE1
  elseif(nprocs .gt. 1 .and. nzderiv .eq. m3 .and. m3 .lt. nz) then
   goto 999   ! CASE2
  elseif(nprocs .gt. 1 .and. nzderiv .eq. nz .and. nz .gt. m3) then
!  CASE3 logic

 allocate ( block_of_cols(locnx,m2,nprocs*m3) )
 block_of_cols=0.

! Create a derived data type for ease of swaps
    call MPI_TYPE_EXTENT(MPI_REAL,sizeofreal,ierr)
! data type for 1d section of basic data slice: m1/nprocs
    call MPI_TYPE_VECTOR(locnx,1,1,MPI_REAL,oneslice,ierr)
! data type for 2d section of basic data slice: m2=ny 1d sections
    call MPI_TYPE_HVECTOR(m2,1,m1*sizeofreal,oneslice,twoslice,ierr)
! data type for 3d section of basic data slice:locnz=nz/numprocs 2d sections
    call MPI_TYPE_HVECTOR(m3,1,m1*m2*sizeofreal,twoslice,subslice,ierr)
! declare data type
    call MPI_TYPE_COMMIT(subslice,ierr)

!!*******************************************************************************
!!*************   DO THE DATA EXCHANGES  ****************************************
!!*******************************************************************************
!! (1)  Ship portions of my slice needed by other processes and
!!      receive data needed to fill in "block_of_cols"

  call slabs_to_columns(data,block_of_cols)

!! (2)  Compute z derivatives over my "block of columns".
!!       Store results in the array "col_of_derivs".
 
  call compact_ddz(block_of_cols,block_of_cols,flag,end_values,locnx,m2,nz)

!! (3) Send derivative results to other processes
!!     and receive their results

  call columns_to_slabs(block_of_cols,deriv)

!!*******************************************************************************
!!*************     END DATA EXCHANGES   ****************************************
!!*******************************************************************************
   
 if( boundary_type2(1,1)=='periodic' .and. boundary_type4(1,1)=='periodic' ) then
  deriv(m1,:,:)=deriv(1,:,:)  ! we didn't compute the extra bdry point
 endif


 call MPI_BARRIER(comm,ierr)

!!! free data type
 call MPI_TYPE_FREE(oneslice,ierr)
 call MPI_TYPE_FREE(twoslice,ierr)
 call MPI_TYPE_FREE(subslice,ierr)

 deallocate(block_of_cols)
return    ! end CASE3 

 else
  write(0,*) ' size logic error using compact_ddz_mpi '
  write(0,*) ' nz : ',nz
  write(0,*) ' locnz : ',m3
  write(0,*) ' nzderiv : ',nzderiv
  write(0,*) ' nprocs : ',nprocs
 endif

999   continue   ! CASE1 & CASE2
 call compact_ddz(data,deriv,flag,end_values,m1,m2,m3)

end subroutine compact_ddz_mpi