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

Last change on this file since 1036 was 1036, checked in by raasch, 9 years ago

code has been put under the GNU General Public License (v3)

  • Property svn:keywords set to Id
File size: 6.5 KB
Line 
1 SUBROUTINE write_compressed( field, fid_avs, fid_fld, my_id, nxl, nxr, nyn, &
2                              nys, nzb, nz_do3d, prec, nbgp )
3
4!--------------------------------------------------------------------------------!
5! This file is part of PALM.
6!
7! PALM is free software: you can redistribute it and/or modify it under the terms
8! of the GNU General Public License as published by the Free Software Foundation,
9! either version 3 of the License, or (at your option) any later version.
10!
11! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
12! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
13! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
14!
15! You should have received a copy of the GNU General Public License along with
16! PALM. If not, see <http://www.gnu.org/licenses/>.
17!
18! Copyright 1997-2012  Leibniz University Hannover
19!--------------------------------------------------------------------------------!
20!
21! Current revisions:
22! -----------------
23!
24! Former revisions:
25! ---------------------
26! $Id: write_compressed.f90 1036 2012-10-22 13:43:42Z raasch $
27!
28! 667 2010-12-23 12:06:00Z suehring/gryschka
29! Array bounds and nx, ny adapted with nbgp
30!
31! 622 2010-12-10 08:08:13Z raasch
32! optional barriers included in order to speed up collective operations
33!
34! Feb. 2007
35! RCS Log replace by Id keyword, revision history cleaned up
36!
37! Revision 1.4  2006/02/23 13:15:09  raasch
38! nz_plot3d renamed nz_do3d
39!
40! Revision 1.1  1999/03/02 09:25:21  raasch
41! Initial revision
42!
43!
44! Description:
45! ------------
46! In this routine, 3D-data (to be plotted) are scaled and compressed by
47! the method of bit shifting. It is designed for the use outside of PALM
48! also, which is the reason why most of the data is passed by subroutine
49! arguments. Nevertheless, the module pegrid is needed by MPI calls.
50!
51! Arguments:
52! field         = data array to be compressed
53! fid_avs       = file-ID of AVS-data-file
54! fid_fld       = file-ID of AVS-header-file
55! my_id         = ID of the calling PE
56! nxl, nxr      = index bounds of the subdomain along x
57! nyn, nys      = index bounds of the subdomain along y
58! nzb,nz_do3d   = index bounds of the domain along z (can be smaller than
59!                 the total domain)
60! prec          = precision of packed data (number of digits after decimal
61!                 point)
62!------------------------------------------------------------------------------!
63
64    USE pegrid         !  needed for MPI_ALLREDUCE
65
66    IMPLICIT NONE
67
68    INTEGER, PARAMETER   :: ip4 = SELECTED_INT_KIND ( 9 )
69    INTEGER, PARAMETER   :: spk = SELECTED_REAL_KIND( 6 )
70
71    INTEGER ::  ampl, dummy1, dummy2, factor, i, ifieldmax, ifieldmax_l, &
72                ifieldmin, ifieldmin_l, ii, j, k, length, nfree, npack, nsize, &
73                nx, ny, nz, pos, startpos
74    INTEGER(ip4) ::  imask (32)
75    INTEGER(ip4), DIMENSION(:), ALLOCATABLE ::  ifield, packed_ifield
76
77    INTEGER, INTENT(IN) ::  fid_avs, fid_fld, my_id, nxl, nxr, nyn, nys, nzb, &
78                            nz_do3d, prec, nbgp
79
80    REAL(spk), INTENT(IN) :: field(1:((nxr-nxl+1+2*nbgp)*(nyn-nys+1+2*nbgp)*(nz_do3d-nzb+1)))
81
82!
83!-- Initialise local variables
84    ampl      = 0
85    ifieldmax = 0
86    ifieldmin = 0
87    npack = 0
88    nsize = 0
89    DO  i = 1,32
90       imask(i) = (2**i) - 1
91    ENDDO
92
93    nx     = nxr - nxl + 2*nbgp
94    ny     = nyn - nys + 2*nbgp
95    nz     = nz_do3d - nzb
96    length = (nx+1) * (ny+1) * (nz+1)
97
98!
99!-- Allocate memory for integer array
100    ALLOCATE ( ifield(1:length) )
101
102!
103!-- Store data on integer (in desired precision)
104    factor = 10**prec
105    DO  i = 1, length
106       ifield(i) = NINT( field(i) * factor )
107    ENDDO
108
109!
110!-- Find minimum and maximum
111    ifieldmax_l = MAXVAL( ifield )
112    ifieldmin_l = MINVAL( ifield )
113
114#if defined( __parallel )
115    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
116    CALL MPI_ALLREDUCE( ifieldmax_l, ifieldmax, 1, MPI_INTEGER, MPI_MAX, &
117                        comm2d, ierr )
118    IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
119    CALL MPI_ALLREDUCE( ifieldmin_l, ifieldmin, 1, MPI_INTEGER, MPI_MIN, &
120                        comm2d, ierr )
121#else
122    ifieldmax = ifieldmax_l
123    ifieldmin = ifieldmin_l
124#endif
125
126!
127!-- Minimum scaling
128    ifield = ifield - ifieldmin
129
130!
131!-- Calculate the number of bits needed for each value
132    ampl  = ifieldmax - ifieldmin
133    nsize = 1
134
135    DO WHILE ( imask(nsize) < ampl )
136       nsize = nsize + 1
137    ENDDO
138
139!
140!-- Calculate size of the packed array
141    npack = length * nsize
142    IF ( MOD( npack, 32 ) /= 0 )  npack = npack + 32
143    npack = npack / 32
144
145!
146!-- Start packing the data
147    ALLOCATE ( packed_ifield(1:npack) )
148    packed_ifield = 0
149
150!
151!-- Starting position of a word
152    startpos = 0
153
154!
155!-- Starting position of the word to which data are actually written
156    ii = 1
157
158!
159!-- Compress all data
160    DO  i = 1, length
161
162!
163!--    Cut the significant bits from the actual grid point value (GPV)
164       dummy1 = IAND( ifield(i), imask(nsize) )
165
166!
167!--    Calculate number of free bits of the actual word after packing the GPV
168       nfree = 32 - startpos - nsize
169
170       IF ( nfree > 0 )  THEN 
171!
172!--       GPV fits to the actual word (ii), additional bits are still free.
173!--       Shift GPV to the new position
174          dummy2 = ISHFT( dummy1 ,nfree )
175
176!
177!--       Append bits to the already packed data
178          packed_ifield(ii) = packed_ifield(ii) + dummy2
179
180!
181!--       Calculate new starting position
182          startpos = startpos + nsize
183
184       ELSEIF ( nfree .EQ. 0 )  THEN 
185!
186!--       GPV fills the actual word (ii) exactly
187          packed_ifield(ii) = packed_ifield(ii) + dummy1
188
189!
190!--       Activate next (new) word
191          ii = ii + 1
192
193!
194!--       Reset starting position of the new word
195          startpos = 0
196
197       ELSE 
198!
199!--       GPV must be split up to the actual (ii) and the following (ii+1)
200!--       word. Shift first part of GPV to its position.
201          dummy2 = ISHFT( dummy1, nfree )
202
203!
204!--       Append bits
205          packed_ifield(ii) = packed_ifield(ii) + dummy2
206
207!
208!--       Store rest of GPV on the next word
209          ii = ii + 1
210          packed_ifield(ii) = ISHFT( dummy1, 32+nfree )
211!
212!--       Calculate starting position of the next GPV
213          startpos = -nfree
214
215       ENDIF
216
217    ENDDO
218
219!
220!-- Write the compressed 3D array
221    WRITE ( fid_avs )  packed_ifield
222
223!
224!-- Write additional informations on  FLD-file
225    IF ( my_id == 0 )  WRITE ( fid_fld, 100 )  prec, ifieldmin, nsize, length
226
227    DEALLOCATE( ifield, packed_ifield )
228
229!
230!-- Formats
231100 FORMAT ('# precision = ',I4/ &
232            '# feldmin   = ',I8/ &
233            '# nbits     = ',I2/ &
234            '# nskip     = ',I8)
235
236END SUBROUTINE write_compressed
Note: See TracBrowser for help on using the repository browser.