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

 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 1321 2014-03-20 09:40:40Z raasch $
!
! 1320 2014-03-20 08:40:49Z raasch
! kind-parameters added to all INTEGER and REAL declaration statements, 
! kinds are defined in new module kinds, 
! old module precision_kind is removed,
! revision history before 2012 removed,
! comment fields (!:) to be used for variable explanations added to
! all variable declaration statements 
!
! 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)
!
! 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 kinds
    
    USE pegrid         !  needed for MPI_ALLREDUCE

    IMPLICIT NONE

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

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

    REAL(sp), 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