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

!------------------------------------------------------------------------------!
! Current revisions:
! -----------------
! 
!
! Former revisions:
! ---------------------
! $Id: write_compressed.f90 623 2010-12-10 08:52:17Z heinze $
!
! 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, j, k, length, nfree, npack, nsize, &
                nx, ny, nz, pos, 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

    REAL(spk), INTENT(IN) :: field(1:((nxr-nxl+3)*(nyn-nys+3)*(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
    ny     = nyn - nys + 2
    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