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