source: palm/trunk/SOURCE/write_compressed.f90 @ 1315

 SUBROUTINE write_compressed( field, fid_avs, fid_fld, my_id, nxl, nxr, nyn, &
                              nys, nzb, nz_do3d, prec, nbgp )

!--------------------------------------------------------------------------------!
! This file is part of PALM.
!
! PALM is free software: you can redistribute it and/or modify it under the terms
! of the GNU General Public License as published by the Free Software Foundation,
! either version 3 of the License, or (at your option) any later version.
!
! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License along with
! PALM. If not, see <http://www.gnu.org/licenses/>.
!
! Copyright 1997-2014 Leibniz Universitaet Hannover
!--------------------------------------------------------------------------------!
!
! Current revisions:
! ------------------
! 
!
! Former revisions:
! -----------------
! $Id: write_compressed.f90 1310 2014-03-14 08:01:56Z suehring $
!
! 1092 2013-02-02 11:24:22Z raasch
! unused variables removed
!
! 1036 2012-10-22 13:43:42Z raasch
! code put under GPL (PALM 3.9)
!
! 667 2010-12-23 12:06:00Z suehring/gryschka
! Array bounds and nx, ny adapted with nbgp
!
! 622 2010-12-10 08:08:13Z raasch
! optional barriers included in order to speed up collective operations
!
! Feb. 2007
! RCS Log replace by Id keyword, revision history cleaned up
!
! Revision 1.4  2006/02/23 13:15:09  raasch
! nz_plot3d renamed nz_do3d
!
! Revision 1.1  1999/03/02 09:25:21  raasch
! Initial revision
!
!
! Description:
! ------------
! In this routine, 3D-data (to be plotted) are scaled and compressed by
! the method of bit shifting. It is designed for the use outside of PALM
! also, which is the reason why most of the data is passed by subroutine
! arguments. Nevertheless, the module pegrid is needed by MPI calls.
!
! Arguments:
! field         = data array to be compressed
! fid_avs       = file-ID of AVS-data-file
! fid_fld       = file-ID of AVS-header-file
! my_id         = ID of the calling PE
! nxl, nxr      = index bounds of the subdomain along x
! nyn, nys      = index bounds of the subdomain along y
! nzb,nz_do3d   = index bounds of the domain along z (can be smaller than
!                 the total domain)
! prec          = precision of packed data (number of digits after decimal
!                 point)
!------------------------------------------------------------------------------!

    USE pegrid         !  needed for MPI_ALLREDUCE

    IMPLICIT NONE

    INTEGER, PARAMETER   :: ip4 = SELECTED_INT_KIND ( 9 )
    INTEGER, PARAMETER   :: spk = SELECTED_REAL_KIND( 6 )

    INTEGER ::  ampl, dummy1, dummy2, factor, i, ifieldmax, ifieldmax_l, &
                ifieldmin, ifieldmin_l, ii, length, nfree, npack, nsize, &
                nx, ny, nz, startpos
    INTEGER(ip4) ::  imask (32)
    INTEGER(ip4), DIMENSION(:), ALLOCATABLE ::  ifield, packed_ifield

    INTEGER, INTENT(IN) ::  fid_avs, fid_fld, my_id, nxl, nxr, nyn, nys, nzb, &
                            nz_do3d, prec, nbgp

    REAL(spk), INTENT(IN) :: field(1:((nxr-nxl+1+2*nbgp)*(nyn-nys+1+2*nbgp)*(nz_do3d-nzb+1)))

!
!-- Initialise local variables
    ampl      = 0
    ifieldmax = 0
    ifieldmin = 0
    npack = 0
    nsize = 0
    DO  i = 1,32
       imask(i) = (2**i) - 1
    ENDDO

    nx     = nxr - nxl + 2*nbgp
    ny     = nyn - nys + 2*nbgp
    nz     = nz_do3d - nzb
    length = (nx+1) * (ny+1) * (nz+1)

!
!-- Allocate memory for integer array
    ALLOCATE ( ifield(1:length) )

!
!-- Store data on integer (in desired precision)
    factor = 10**prec
    DO  i = 1, length
       ifield(i) = NINT( field(i) * factor )
    ENDDO

!
!-- Find minimum and maximum
    ifieldmax_l = MAXVAL( ifield )
    ifieldmin_l = MINVAL( ifield )

#if defined( __parallel )
    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    CALL MPI_ALLREDUCE( ifieldmax_l, ifieldmax, 1, MPI_INTEGER, MPI_MAX, &
                        comm2d, ierr )
    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    CALL MPI_ALLREDUCE( ifieldmin_l, ifieldmin, 1, MPI_INTEGER, MPI_MIN, &
                        comm2d, ierr )
#else
    ifieldmax = ifieldmax_l
    ifieldmin = ifieldmin_l
#endif

!
!-- Minimum scaling
    ifield = ifield - ifieldmin

!
!-- Calculate the number of bits needed for each value
    ampl  = ifieldmax - ifieldmin
    nsize = 1

    DO WHILE ( imask(nsize) < ampl )
       nsize = nsize + 1
    ENDDO

! 
!-- Calculate size of the packed array
    npack = length * nsize
    IF ( MOD( npack, 32 ) /= 0 )  npack = npack + 32
    npack = npack / 32

!
!-- Start packing the data
    ALLOCATE ( packed_ifield(1:npack) )
    packed_ifield = 0

!
!-- Starting position of a word
    startpos = 0

!
!-- Starting position of the word to which data are actually written
    ii = 1

!
!-- Compress all data
    DO  i = 1, length

!
!--    Cut the significant bits from the actual grid point value (GPV)
       dummy1 = IAND( ifield(i), imask(nsize) )

!
!--    Calculate number of free bits of the actual word after packing the GPV
       nfree = 32 - startpos - nsize

       IF ( nfree > 0 )  THEN 
!
!--       GPV fits to the actual word (ii), additional bits are still free.
!--       Shift GPV to the new position
          dummy2 = ISHFT( dummy1 ,nfree )

!
!--       Append bits to the already packed data
          packed_ifield(ii) = packed_ifield(ii) + dummy2

!
!--       Calculate new starting position
          startpos = startpos + nsize

       ELSEIF ( nfree .EQ. 0 )  THEN 
!
!--       GPV fills the actual word (ii) exactly
          packed_ifield(ii) = packed_ifield(ii) + dummy1

!
!--       Activate next (new) word
          ii = ii + 1

!
!--       Reset starting position of the new word
          startpos = 0

       ELSE 
!
!--       GPV must be split up to the actual (ii) and the following (ii+1)
!--       word. Shift first part of GPV to its position.
          dummy2 = ISHFT( dummy1, nfree )

!
!--       Append bits
          packed_ifield(ii) = packed_ifield(ii) + dummy2 

!
!--       Store rest of GPV on the next word
          ii = ii + 1
          packed_ifield(ii) = ISHFT( dummy1, 32+nfree )
!
!--       Calculate starting position of the next GPV
          startpos = -nfree

       ENDIF

    ENDDO

!
!-- Write the compressed 3D array
    WRITE ( fid_avs )  packed_ifield

!
!-- Write additional informations on  FLD-file
    IF ( my_id == 0 )  WRITE ( fid_fld, 100 )  prec, ifieldmin, nsize, length

    DEALLOCATE( ifield, packed_ifield )

!
!-- Formats
100 FORMAT ('# precision = ',I4/ &
            '# feldmin   = ',I8/ &
            '# nbits     = ',I2/ &
            '# nskip     = ',I8)

END SUBROUTINE write_compressed