NCEPLIBS-g2  3.4.5
pack_gp.f
Go to the documentation of this file.
1 C> @file
2 C> @brief This subroutine determines groups of variable size.
3 C> @author Harry Glahn @date 1994-02-01
4 C>
5 
6 C> This subroutine determines groups of variable size, but at least
7 C> of size minpk, the associated max(JMAX) and min(JMIN), the number
8 C> of bits necessary to hold the values in each group LBIT, the
9 C> number of values in each group NOV, the number of bits necessary
10 C> to pack the JMIN values IBIT, the number of bits necessary to pack
11 C> the LBIT values JBIT, and the number of bits necessary to pack the
12 C> NOV values KBIT. The routine is designed to determine the groups
13 C> such that a small number of bits is necessary to pack the data
14 C> without excessive computations. If all values in the group are
15 C> zero, the number of bits to use in packing is defined as zero when
16 C> there can be no missing values; when there can be missing values,
17 C> the number of bits must be at least 1 to have the capability to
18 C> recognize the missing value. However, if all values in a group are
19 C> missing, the number of bits needed is 0, and the unpacker recognizes
20 C> this. All variables are integer. even though the groups are
21 C> initially of size minpk or larger, an adjustment between two groups
22 C> (the lookback procedure) may make a group smaller than minpk. The
23 C> control on group size is that the sum of the sizes of the two
24 C> consecutive groups, each of size minpk or larger, is not decreased.
25 C> When determining the number of bits necessary for packing, the
26 C> largest value that can be accommodated in, say mbits is 2**mbits-1
27 C> this largest value (and the next smallest value) is reserved for
28 C> the missing value indicator (only) when is523 ne 0. If the
29 C> dimension NDG is not large enough to hold all the groups, the
30 C> local value of minpk is increased by 50 percent. this is repeated
31 C> until ndg will suffice. A diagnostic is printed whenever this
32 C> happens, which should be very rarely. If it happens often, NDG in
33 C> subroutine pack should be increased and a corresponding increase
34 C> in subroutine unpack made. Considerable code is provided so that
35 C> no more checking for missing values within loops is done than
36 C> necessary; the added efficiency of this is relatively minor,
37 C> but does no harm. For grib2, the reference value for the length
38 C> of groups in nov and for the number of bits necessary to pack
39 C> group values are determined, and subtracted before jbit and kbit
40 C> are determined. When 1 or more groups are large compared to the
41 C> others, the width of all groups must be as large as the largest.
42 C> A subroutine reduce breaks up large groups into 2 or more to reduce
43 C> total bits required. If reduce should abort, pack_gp will be
44 C> executed again without the call to reduce.
45 C>
46 C> PROGRAM HISTORY LOG:
47 C> - 1994-02-01 Harry Glahn tdl mos-2000.
48 C> - 1995-06-01 Harry Glahn modified for lmiss error.
49 C> - 1996-07-01 Harry Glahn added misss.
50 C> - 1997-02-01 Harry Glahn removed 4 redundant tests for missp.eq.0;
51 C> inserted a test to better handle a string of 9999's.
52 C> - 1997-02-01 Harry Glahn added loops to eliminate test for misss
53 C> when misss = 0.
54 C> - 1997-03-01 Harry Glahn corrected for secondary missing value.
55 C> - 1997-03-01 Harry Glahn corrected for use of local value of minpk.
56 C> - 1997-03-01 Harry Glahn corrected for secondary missing value.
57 C> - 1997-03-01 Harry Glahn changed calculating number of bits
58 C> through exponents to an array (improved overall packing performance
59 C> by about 35 percent). Allowed 0 bit for packing JMIN, LBIT, and NOV.
60 C> - 1997-05-01 Harry Glahn a number of changes for efficiency. mod
61 C> functions eliminated and one ifthen added. Jount removed.
62 C> Recomputation of bits not made unless necessary after moving points
63 C> from one group to another. Nendb adjusted to eliminate possibility
64 C> of very small group at the end. About 8 percent improvement in
65 C> overall packing. ISKIPA removed; There is always a group b that can
66 C> become group A. Control on size of group b (statement below 150)
67 C> added. Added adda, and use of ge and le instead of gt and lt in
68 C> loop between 150 and 160. IBITBS added to shorten trip through loop.
69 C> - 2000-03-01 Harry Glahn modified for grib2; changed name from
70 C> packgp.
71 C> - 2001-01-01 Harry Glahn Add comments; ier = 706 substituted for
72 C> stops; added return; removed statement number 110; added ier.
73 C> - 2001-11-01 Harry Glahn changed some diagnostic formats to
74 C> allow printing larger numbers
75 C> - 2001-11-01 Harry Glahn added misslx to put maximum value into JMIN
76 C> when all values missing to agree with grib standard.
77 C> - 2001-11-01 Harry Glahn changed two tests on missp and misss eq 0
78 C> to tests on is523. However, missp and misss cannot in general be 0.
79 C> - 2001-11-01 Harry Glahn added call to reduce; defined itest
80 C> before loops to reduce computation; started large group when all
81 C> same value.
82 C> - 2001-12-01 Harry Glahn modified and added a few comments.
83 C> - 2002-01-01 Harry Glahn removed loop before 150 to determine
84 C> a group of all same value.
85 C> - 2002-01-01 Harry Glahn changed mallow from 9999999 to 2**30+1,
86 C> and made it a parameter.
87 C> - 2002-03-01 Harry Glahn added non fatal ier = 716, 717; removed
88 C> nendb=nxy above 150; added iersav=0.
89 C>
90 C> @param[in] KFILDO unit number for output/print file.
91 C> @param[in] IC array to hold data for packing. The values do not
92 C> have to be positive at this point, but must be in the range
93 C> -2**30 to +2**30 (the value of mallow). These integer values
94 C> will be retained exactly through packing and unpacking.
95 C> @param[in] NXY number of values in IC. also treated as
96 C> its dimension.
97 C> @param[in] IS523 missing value management 0=data contains no
98 C> missing values: 1 data contains primary missing values; 2=data
99 C> contains primary and secondary missing values.
100 C> @param[in] MINPK the minimum size of each group, except possibly
101 C> the last one.
102 C> @param[in] INC the number of values to add to an already existing
103 C> group in determining whether or not to start a new group. Ideally,
104 C> this would be 1, but each time inc values are attempted, the max
105 C> and min of the next minpk values must be found. This is "a loop
106 C> within a loop," and a slightly larger value may give about as good
107 C> results with slightly less computational time. If inc is le 0, 1
108 C> is used, and a diagnostic is output. note: it is expected that
109 C> INC will equal 1. The code uses inc primarily in the loops
110 C> starting at statement 180. If INC were 1, there would not need
111 C> to be loops as such. However, kinc (the local value of INC) is
112 C> set ge 1 when near the end of the data to forestall a very small
113 C> group at the end.
114 C> @param[in] MISSP when missing points can be present in the data,
115 C> they will have the value missp or misss. missp is the primary
116 C> missing value and misss is the secondary missing value. These
117 C> must not be values that would occur with subtracting the minimum
118 C> (reference) value or scaling. for example, missp = 0 would not
119 C> be advisable.
120 C> @param[in] MISSS secondary missing value indicator (see missp).
121 C> @param[out] JMIN the minimum of each group (j=1,lx).
122 C> @param[out] JMAX the maximum of each group (j=1,lx). This is not
123 C> really needed, but since the max of each group must be found,
124 C> saving it here is cheap in case the user wants it.
125 C> @param[out] LBIT the number of bits necessary to pack each group
126 C> (j=1,lx). It is assumed the minimum of each group will be removed
127 C> before packing, and the values to pack will, therefore, all be
128 C> positive. However, IC does not necessarily contain all positive
129 C> values. If the overall minimum has been removed (the usual case),
130 C> then IC will contain only positive values.
131 C> @param[out] NOV the number of values in each group (j=1,lx).
132 C> @param[in] NDG the dimension of JMIN, JMAX, LBIT, and NOV.
133 C> @param[out] LX the number of groups determined.
134 C> @param[out] IBIT the number of bits necessary to pack the JMIN(j)
135 C> values, j=1,LX.
136 C> @param[out] JBIT the number of bits necessary to pack the LBIT(j)
137 C> values, j=1,LX.
138 C> @param[out] KBIT the number of bits necessary to pack the NOV(j)
139 C> values, j=1,LX.
140 C> @param[out] NOVREF reference value for NOV.
141 C> @param[out] LBITREF reference value for LBIT.
142 C> @param[out] IER error return.
143 C> - 706 value will not pack in 30 bits--fatal
144 C> - 714 error in reduce--non-fatal
145 C> - 715 ngp not large enough in reduce--non-fatal
146 C> - 716 minpk inceased--non-fatal
147 C> - 717 inc set = 1--non-fatal
148 C>
149 C> @author Harry Glahn @date 1994-02-01
150 C>
151 
152  SUBROUTINE pack_gp(KFILDO,IC,NXY,IS523,MINPK,INC,MISSP,MISSS,
153  1 JMIN,JMAX,LBIT,NOV,NDG,LX,IBIT,JBIT,KBIT,
154  2 NOVREF,LBITREF,IER)
155 
156  parameter(mallow=2**30+1)
157 C
158  CHARACTER*1 CFEED
159  LOGICAL ADDA
160 C
161  dimension ic(nxy)
162  dimension jmin(ndg),jmax(ndg),lbit(ndg),nov(ndg)
163  dimension misslx(ndg)
164 C MISSLX( ) IS AN AUTOMATIC ARRAY.
165  INTEGER, PARAMETER :: IBXX2(0:30) = (/ 1, 2, 4, 8, 16, 32, 64, &
166  & 128, 256, 512, 1024, 2048, 4096, 8192, 16384, 32768, 65536, &
167  & 131072, 262144, 524288, 1048576, 2097152, 4194304, 8388608, &
168  & 16777216, 33554432, 67108864, 134217728, 268435456, &
169  & 536870912, 1073741824 /)
170 C
171 
172  PARAMETER ifeed=12
173 C
174  ier=0
175  iersav=0
176 C CALL TIMPR(KFILDO,KFILDO,'START PACK_GP ')
177  cfeed=char(ifeed)
178 C
179  ired=0
180 C IRED IS A FLAG. WHEN ZERO, REDUCE WILL BE CALLED.
181 C IF REDUCE ABORTS, IRED = 1 AND IS NOT CALLED. IN
182 C THIS CASE PACK_GP EXECUTES AGAIN EXCEPT FOR REDUCE.
183 C
184  IF(inc.LE.0)THEN
185  iersav=717
186 C WRITE(KFILDO,101)INC
187 C101 FORMAT(/' ****INC ='I8,' NOT CORRECT IN PACK_GP. 1 IS USED.')
188  ENDIF
189 C
190 C THERE WILL BE A RESTART OF PACK_GP IF SUBROUTINE REDUCE
191 C ABORTS. THIS SHOULD NOT HAPPEN, BUT IF IT DOES, PACK_GP
192 C WILL COMPLETE WITHOUT SUBROUTINE REDUCE. A NON FATAL
193 C DIAGNOSTIC RETURN IS PROVIDED.
194 C
195  102 kinc=max(inc,1)
196  lminpk=minpk
197 C
198 C THERE WILL BE A RESTART AT 105 IS NDG IS NOT LARGE ENOUGH.
199 C A NON FATAL DIAGNOSTIC RETURN IS PROVIDED.
200 C
201  105 kstart=1
202  ktotal=0
203  lx=0
204  adda=.false.
205  lmiss=0
206  IF(is523.EQ.1)lmiss=1
207  IF(is523.EQ.2)lmiss=2
208 C
209 C *************************************
210 C
211 C THIS SECTION COMPUTES STATISTICS FOR GROUP A. GROUP A IS
212 C A GROUP OF SIZE LMINPK.
213 C
214 C *************************************
215 C
216  ibita=0
217  mina=mallow
218  maxa=-mallow
219  minak=mallow
220  maxak=-mallow
221 C
222 C FIND THE MIN AND MAX OF GROUP A. THIS WILL INITIALLY BE OF
223 C SIZE LMINPK (IF THERE ARE STILL LMINPK VALUES IN IC( )), BUT
224 C WILL INCREASE IN SIZE IN INCREMENTS OF INC UNTIL A NEW
225 C GROUP IS STARTED. THE DEFINITION OF GROUP A IS DONE HERE
226 C ONLY ONCE (UPON INITIAL ENTRY), BECAUSE A GROUP B CAN ALWAYS
227 C BECOME A NEW GROUP A AFTER A IS PACKED, EXCEPT IF LMINPK
228 C HAS TO BE INCREASED BECAUSE NDG IS TOO SMALL. THEREFORE,
229 C THE SEPARATE LOOPS FOR MISSING AND NON-MISSING HERE BUYS
230 C ALMOST NOTHING.
231 C
232  nenda=min(kstart+lminpk-1,nxy)
233  IF(nxy-nenda.LE.lminpk/2)nenda=nxy
234 C ABOVE STATEMENT GUARANTEES THE LAST GROUP IS GT LMINPK/2 BY
235 C MAKING THE ACTUAL GROUP LARGER. IF A PROVISION LIKE THIS IS
236 C NOT INCLUDED, THERE WILL MANY TIMES BE A VERY SMALL GROUP
237 C AT THE END. USE SEPARATE LOOPS FOR MISSING AND NO MISSING
238 C VALUES FOR EFFICIENCY.
239 C
240 C DETERMINE WHETHER THERE IS A LONG STRING OF THE SAME VALUE
241 C UNLESS NENDA = NXY. THIS MAY ALLOW A LARGE GROUP A TO
242 C START WITH, AS WITH MISSING VALUES. SEPARATE LOOPS FOR
243 C MISSING OPTIONS. THIS SECTION IS ONLY EXECUTED ONCE,
244 C IN DETERMINING THE FIRST GROUP. IT HELPS FOR AN ARRAY
245 C OF MOSTLY MISSING VALUES OR OF ONE VALUE, SUCH AS
246 C RADAR OR PRECIP DATA.
247 C
248  IF(nenda.NE.nxy.AND.ic(kstart).EQ.ic(kstart+1))THEN
249 C NO NEED TO EXECUTE IF FIRST TWO VALUES ARE NOT EQUAL.
250 C
251  IF(is523.EQ.0)THEN
252 C THIS LOOP IS FOR NO MISSING VALUES.
253 C
254  DO 111 k=kstart+1,nxy
255 C
256  IF(ic(k).NE.ic(kstart))THEN
257  nenda=max(nenda,k-1)
258  GO TO 114
259  ENDIF
260 C
261  111 CONTINUE
262 C
263  nenda=nxy
264 C FALL THROUGH THE LOOP MEANS ALL VALUES ARE THE SAME.
265 C
266  ELSEIF(is523.EQ.1)THEN
267 C THIS LOOP IS FOR PRIMARY MISSING VALUES ONLY.
268 C
269  DO 112 k=kstart+1,nxy
270 C
271  IF(ic(k).NE.missp)THEN
272 C
273  IF(ic(k).NE.ic(kstart))THEN
274  nenda=max(nenda,k-1)
275  GO TO 114
276  ENDIF
277 C
278  ENDIF
279 C
280  112 CONTINUE
281 C
282  nenda=nxy
283 C FALL THROUGH THE LOOP MEANS ALL VALUES ARE THE SAME.
284 C
285  ELSE
286 C THIS LOOP IS FOR PRIMARY AND SECONDARY MISSING VALUES.
287 C
288  DO 113 k=kstart+1,nxy
289 C
290  IF(ic(k).NE.missp.AND.ic(k).NE.misss)THEN
291 C
292  IF(ic(k).NE.ic(kstart))THEN
293  nenda=max(nenda,k-1)
294  GO TO 114
295  ENDIF
296 C
297  ENDIF
298 C
299  113 CONTINUE
300 C
301  nenda=nxy
302 C FALL THROUGH THE LOOP MEANS ALL VALUES ARE THE SAME.
303  ENDIF
304 C
305  ENDIF
306 C
307  114 IF(is523.EQ.0)THEN
308 C
309  DO 115 k=kstart,nenda
310  IF(ic(k).LT.mina)THEN
311  mina=ic(k)
312  minak=k
313  ENDIF
314  IF(ic(k).GT.maxa)THEN
315  maxa=ic(k)
316  maxak=k
317  ENDIF
318  115 CONTINUE
319 C
320  ELSEIF(is523.EQ.1)THEN
321 C
322  DO 117 k=kstart,nenda
323  IF(ic(k).EQ.missp)GO TO 117
324  IF(ic(k).LT.mina)THEN
325  mina=ic(k)
326  minak=k
327  ENDIF
328  IF(ic(k).GT.maxa)THEN
329  maxa=ic(k)
330  maxak=k
331  ENDIF
332  117 CONTINUE
333 C
334  ELSE
335 C
336  DO 120 k=kstart,nenda
337  IF(ic(k).EQ.missp.OR.ic(k).EQ.misss)GO TO 120
338  IF(ic(k).LT.mina)THEN
339  mina=ic(k)
340  minak=k
341  ENDIF
342  IF(ic(k).GT.maxa)THEN
343  maxa=ic(k)
344  maxak=k
345  ENDIF
346  120 CONTINUE
347 C
348  ENDIF
349 C
350  kounta=nenda-kstart+1
351 C
352 C INCREMENT KTOTAL AND FIND THE BITS NEEDED TO PACK THE A GROUP.
353 C
354  ktotal=ktotal+kounta
355  mislla=0
356  IF(mina.NE.mallow)GO TO 125
357 C ALL MISSING VALUES MUST BE ACCOMMODATED.
358  mina=0
359  maxa=0
360  mislla=1
361  ibitb=0
362  IF(is523.NE.2)GO TO 130
363 C WHEN ALL VALUES ARE MISSING AND THERE ARE NO
364 C SECONDARY MISSING VALUES, IBITA = 0.
365 C OTHERWISE, IBITA MUST BE CALCULATED.
366 C
367  125 itest=maxa-mina+lmiss
368 C
369  DO 126 ibita=0,30
370  IF(itest.LT.ibxx2(ibita))GO TO 130
371 C*** THIS TEST IS THE SAME AS:
372 C*** IF(MAXA-MINA.LT.IBXX2(IBITA)-LMISS)GO TO 130
373  126 CONTINUE
374 C
375 C WRITE(KFILDO,127)MAXA,MINA
376 C127 FORMAT(' ****ERROR IN PACK_GP. VALUE WILL NOT PACK IN 30 BITS.',
377 C 1 ' MAXA ='I13,' MINA ='I13,'. ERROR AT 127.')
378  ier=706
379  GO TO 900
380 C
381  130 CONTINUE
382 C
383 C***D WRITE(KFILDO,131)KOUNTA,KTOTAL,MINA,MAXA,IBITA,MISLLA
384 C***D131 FORMAT(' AT 130, KOUNTA ='I8,' KTOTAL ='I8,' MINA ='I8,
385 C***D 1 ' MAXA ='I8,' IBITA ='I3,' MISLLA ='I3)
386 C
387  133 IF(ktotal.GE.nxy)GO TO 200
388 C
389 C *************************************
390 C
391 C THIS SECTION COMPUTES STATISTICS FOR GROUP B. GROUP B IS A
392 C GROUP OF SIZE LMINPK IMMEDIATELY FOLLOWING GROUP A.
393 C
394 C *************************************
395 C
396  140 minb=mallow
397  maxb=-mallow
398  minbk=mallow
399  maxbk=-mallow
400  ibitbs=0
401  mstart=ktotal+1
402 C
403 C DETERMINE WHETHER THERE IS A LONG STRING OF THE SAME VALUE.
404 C THIS WORKS WHEN THERE ARE NO MISSING VALUES.
405 C
406  nendb=1
407 C
408  IF(mstart.LT.nxy)THEN
409 C
410  IF(is523.EQ.0)THEN
411 C THIS LOOP IS FOR NO MISSING VALUES.
412 C
413  DO 145 k=mstart+1,nxy
414 C
415  IF(ic(k).NE.ic(mstart))THEN
416  nendb=k-1
417  GO TO 150
418  ENDIF
419 C
420  145 CONTINUE
421 C
422  nendb=nxy
423 C FALL THROUGH THE LOOP MEANS ALL REMAINING VALUES
424 C ARE THE SAME.
425  ENDIF
426 C
427  ENDIF
428 C
429  150 nendb=max(nendb,min(ktotal+lminpk,nxy))
430 C**** 150 NENDB=MIN(KTOTAL+LMINPK,NXY)
431 C
432  IF(nxy-nendb.LE.lminpk/2)nendb=nxy
433 C ABOVE STATEMENT GUARANTEES THE LAST GROUP IS GT LMINPK/2 BY
434 C MAKING THE ACTUAL GROUP LARGER. IF A PROVISION LIKE THIS IS
435 C NOT INCLUDED, THERE WILL MANY TIMES BE A VERY SMALL GROUP
436 C AT THE END. USE SEPARATE LOOPS FOR MISSING AND NO MISSING
437 C
438 C USE SEPARATE LOOPS FOR MISSING AND NO MISSING VALUES
439 C FOR EFFICIENCY.
440 C
441  IF(is523.EQ.0)THEN
442 C
443  DO 155 k=mstart,nendb
444  IF(ic(k).LE.minb)THEN
445  minb=ic(k)
446 C NOTE LE, NOT LT. LT COULD BE USED BUT THEN A
447 C RECOMPUTE OVER THE WHOLE GROUP WOULD BE NEEDED
448 C MORE OFTEN. SAME REASONING FOR GE AND OTHER
449 C LOOPS BELOW.
450  minbk=k
451  ENDIF
452  IF(ic(k).GE.maxb)THEN
453  maxb=ic(k)
454  maxbk=k
455  ENDIF
456  155 CONTINUE
457 C
458  ELSEIF(is523.EQ.1)THEN
459 C
460  DO 157 k=mstart,nendb
461  IF(ic(k).EQ.missp)GO TO 157
462  IF(ic(k).LE.minb)THEN
463  minb=ic(k)
464  minbk=k
465  ENDIF
466  IF(ic(k).GE.maxb)THEN
467  maxb=ic(k)
468  maxbk=k
469  ENDIF
470  157 CONTINUE
471 C
472  ELSE
473 C
474  DO 160 k=mstart,nendb
475  IF(ic(k).EQ.missp.OR.ic(k).EQ.misss)GO TO 160
476  IF(ic(k).LE.minb)THEN
477  minb=ic(k)
478  minbk=k
479  ENDIF
480  IF(ic(k).GE.maxb)THEN
481  maxb=ic(k)
482  maxbk=k
483  ENDIF
484  160 CONTINUE
485 C
486  ENDIF
487 C
488  kountb=nendb-ktotal
489  misllb=0
490  IF(minb.NE.mallow)GO TO 165
491 C ALL MISSING VALUES MUST BE ACCOMMODATED.
492  minb=0
493  maxb=0
494  misllb=1
495  ibitb=0
496 C
497  IF(is523.NE.2)GO TO 170
498 C WHEN ALL VALUES ARE MISSING AND THERE ARE NO SECONDARY
499 C MISSING VALUES, IBITB = 0. OTHERWISE, IBITB MUST BE
500 C CALCULATED.
501 C
502  165 DO 166 ibitb=ibitbs,30
503  IF(maxb-minb.LT.ibxx2(ibitb)-lmiss)GO TO 170
504  166 CONTINUE
505 C
506 C WRITE(KFILDO,167)MAXB,MINB
507 C167 FORMAT(' ****ERROR IN PACK_GP. VALUE WILL NOT PACK IN 30 BITS.',
508 C 1 ' MAXB ='I13,' MINB ='I13,'. ERROR AT 167.')
509  ier=706
510  GO TO 900
511 C
512 C COMPARE THE BITS NEEDED TO PACK GROUP B WITH THOSE NEEDED
513 C TO PACK GROUP A. IF IBITB GE IBITA, TRY TO ADD TO GROUP A.
514 C IF NOT, TRY TO ADD A'S POINTS TO B, UNLESS ADDITION TO A
515 C HAS BEEN DONE. THIS LATTER IS CONTROLLED WITH ADDA.
516 C
517  170 CONTINUE
518 C
519 C***D WRITE(KFILDO,171)KOUNTA,KTOTAL,MINA,MAXA,IBITA,MISLLA,
520 C***D 1 MINB,MAXB,IBITB,MISLLB
521 C***D171 FORMAT(' AT 171, KOUNTA ='I8,' KTOTAL ='I8,' MINA ='I8,
522 C***D 1 ' MAXA ='I8,' IBITA ='I3,' MISLLA ='I3,
523 C***D 2 ' MINB ='I8,' MAXB ='I8,' IBITB ='I3,' MISLLB ='I3)
524 C
525  IF(ibitb.GE.ibita)GO TO 180
526  IF(adda)GO TO 200
527 C
528 C *************************************
529 C
530 C GROUP B REQUIRES LESS BITS THAN GROUP A. PUT AS MANY OF A'S
531 C POINTS INTO B AS POSSIBLE WITHOUT EXCEEDING THE NUMBER OF
532 C BITS NECESSARY TO PACK GROUP B.
533 C
534 C *************************************
535 C
536  kounts=kounta
537 C KOUNTA REFERS TO THE PRESENT GROUP A.
538  mintst=minb
539  maxtst=maxb
540  mintstk=minbk
541  maxtstk=maxbk
542 C
543 C USE SEPARATE LOOPS FOR MISSING AND NO MISSING VALUES
544 C FOR EFFICIENCY.
545 C
546  IF(is523.EQ.0)THEN
547 C
548  DO 1715 k=ktotal,kstart,-1
549 C START WITH THE END OF THE GROUP AND WORK BACKWARDS.
550  IF(ic(k).LT.minb)THEN
551  mintst=ic(k)
552  mintstk=k
553  ELSEIF(ic(k).GT.maxb)THEN
554  maxtst=ic(k)
555  maxtstk=k
556  ENDIF
557  IF(maxtst-mintst.GE.ibxx2(ibitb))GO TO 174
558 C NOTE THAT FOR THIS LOOP, LMISS = 0.
559  minb=mintst
560  maxb=maxtst
561  minbk=mintstk
562  maxbk=maxtstk
563  kounta=kounta-1
564 C THERE IS ONE LESS POINT NOW IN A.
565  1715 CONTINUE
566 C
567  ELSEIF(is523.EQ.1)THEN
568 C
569  DO 1719 k=ktotal,kstart,-1
570 C START WITH THE END OF THE GROUP AND WORK BACKWARDS.
571  IF(ic(k).EQ.missp)GO TO 1718
572  IF(ic(k).LT.minb)THEN
573  mintst=ic(k)
574  mintstk=k
575  ELSEIF(ic(k).GT.maxb)THEN
576  maxtst=ic(k)
577  maxtstk=k
578  ENDIF
579  IF(maxtst-mintst.GE.ibxx2(ibitb)-lmiss)GO TO 174
580 C FOR THIS LOOP, LMISS = 1.
581  minb=mintst
582  maxb=maxtst
583  minbk=mintstk
584  maxbk=maxtstk
585  misllb=0
586 C WHEN THE POINT IS NON MISSING, MISLLB SET = 0.
587  1718 kounta=kounta-1
588 C THERE IS ONE LESS POINT NOW IN A.
589  1719 CONTINUE
590 C
591  ELSE
592 C
593  DO 173 k=ktotal,kstart,-1
594 C START WITH THE END OF THE GROUP AND WORK BACKWARDS.
595  IF(ic(k).EQ.missp.OR.ic(k).EQ.misss)GO TO 1729
596  IF(ic(k).LT.minb)THEN
597  mintst=ic(k)
598  mintstk=k
599  ELSEIF(ic(k).GT.maxb)THEN
600  maxtst=ic(k)
601  maxtstk=k
602  ENDIF
603  IF(maxtst-mintst.GE.ibxx2(ibitb)-lmiss)GO TO 174
604 C FOR THIS LOOP, LMISS = 2.
605  minb=mintst
606  maxb=maxtst
607  minbk=mintstk
608  maxbk=maxtstk
609  misllb=0
610 C WHEN THE POINT IS NON MISSING, MISLLB SET = 0.
611  1729 kounta=kounta-1
612 C THERE IS ONE LESS POINT NOW IN A.
613  173 CONTINUE
614 C
615  ENDIF
616 C
617 C AT THIS POINT, KOUNTA CONTAINS THE NUMBER OF POINTS TO CLOSE
618 C OUT GROUP A WITH. GROUP B NOW STARTS WITH KSTART+KOUNTA AND
619 C ENDS WITH NENDB. MINB AND MAXB HAVE BEEN ADJUSTED AS
620 C NECESSARY TO REFLECT GROUP B (EVEN THOUGH THE NUMBER OF BITS
621 C NEEDED TO PACK GROUP B HAVE NOT INCREASED, THE END POINTS
622 C OF THE RANGE MAY HAVE).
623 C
624  174 IF(kounta.EQ.kounts)GO TO 200
625 C ON TRANSFER, GROUP A WAS NOT CHANGED. CLOSE IT OUT.
626 C
627 C ONE OR MORE POINTS WERE TAKEN OUT OF A. RANGE AND IBITA
628 C MAY HAVE TO BE RECOMPUTED; IBITA COULD BE LESS THAN
629 C ORIGINALLY COMPUTED. IN FACT, GROUP A CAN NOW CONTAIN
630 C ONLY ONE POINT AND BE PACKED WITH ZERO BITS
631 C (UNLESS MISSS NE 0).
632 C
633  nouta=kounts-kounta
634  ktotal=ktotal-nouta
635  kountb=kountb+nouta
636  IF(nenda-nouta.GT.minak.AND.nenda-nouta.GT.maxak)GO TO 200
637 C WHEN THE ABOVE TEST IS MET, THE MIN AND MAX OF THE
638 C CURRENT GROUP A WERE WITHIN THE OLD GROUP A, SO THE
639 C RANGE AND IBITA DO NOT NEED TO BE RECOMPUTED.
640 C NOTE THAT MINAK AND MAXAK ARE NO LONGER NEEDED.
641  ibita=0
642  mina=mallow
643  maxa=-mallow
644 C
645 C USE SEPARATE LOOPS FOR MISSING AND NO MISSING VALUES
646 C FOR EFFICIENCY.
647 C
648  IF(is523.EQ.0)THEN
649 C
650  DO 1742 k=kstart,nenda-nouta
651  IF(ic(k).LT.mina)THEN
652  mina=ic(k)
653  ENDIF
654  IF(ic(k).GT.maxa)THEN
655  maxa=ic(k)
656  ENDIF
657  1742 CONTINUE
658 C
659  ELSEIF(is523.EQ.1)THEN
660 C
661  DO 1744 k=kstart,nenda-nouta
662  IF(ic(k).EQ.missp)GO TO 1744
663  IF(ic(k).LT.mina)THEN
664  mina=ic(k)
665  ENDIF
666  IF(ic(k).GT.maxa)THEN
667  maxa=ic(k)
668  ENDIF
669  1744 CONTINUE
670 C
671  ELSE
672 C
673  DO 175 k=kstart,nenda-nouta
674  IF(ic(k).EQ.missp.OR.ic(k).EQ.misss)GO TO 175
675  IF(ic(k).LT.mina)THEN
676  mina=ic(k)
677  ENDIF
678  IF(ic(k).GT.maxa)THEN
679  maxa=ic(k)
680  ENDIF
681  175 CONTINUE
682 C
683  ENDIF
684 C
685  mislla=0
686  IF(mina.NE.mallow)GO TO 1750
687 C ALL MISSING VALUES MUST BE ACCOMMODATED.
688  mina=0
689  maxa=0
690  mislla=1
691  IF(is523.NE.2)GO TO 177
692 C WHEN ALL VALUES ARE MISSING AND THERE ARE NO SECONDARY
693 C MISSING VALUES IBITA = 0 AS ORIGINALLY SET. OTHERWISE,
694 C IBITA MUST BE CALCULATED.
695 C
696  1750 itest=maxa-mina+lmiss
697 C
698  DO 176 ibita=0,30
699  IF(itest.LT.ibxx2(ibita))GO TO 177
700 C*** THIS TEST IS THE SAME AS:
701 C*** IF(MAXA-MINA.LT.IBXX2(IBITA)-LMISS)GO TO 177
702  176 CONTINUE
703 C
704 C WRITE(KFILDO,1760)MAXA,MINA
705 C1760 FORMAT(' ****ERROR IN PACK_GP. VALUE WILL NOT PACK IN 30 BITS.',
706 C 1 ' MAXA ='I13,' MINA ='I13,'. ERROR AT 1760.')
707  ier=706
708  GO TO 900
709 C
710  177 CONTINUE
711  GO TO 200
712 C
713 C *************************************
714 C
715 C AT THIS POINT, GROUP B REQUIRES AS MANY BITS TO PACK AS GROUPA.
716 C THEREFORE, TRY TO ADD INC POINTS TO GROUP A WITHOUT INCREASING
717 C IBITA. THIS AUGMENTED GROUP IS CALLED GROUP C.
718 C
719 C *************************************
720 C
721  180 IF(mislla.EQ.1)THEN
722  minc=mallow
723  minck=mallow
724  maxc=-mallow
725  maxck=-mallow
726  ELSE
727  minc=mina
728  maxc=maxa
729  minck=minak
730  maxck=minak
731  ENDIF
732 C
733  nount=0
734  IF(nxy-(ktotal+kinc).LE.lminpk/2)kinc=nxy-ktotal
735 C ABOVE STATEMENT CONSTRAINS THE LAST GROUP TO BE NOT LESS THAN
736 C LMINPK/2 IN SIZE. IF A PROVISION LIKE THIS IS NOT INCLUDED,
737 C THERE WILL MANY TIMES BE A VERY SMALL GROUP AT THE END.
738 C
739 C USE SEPARATE LOOPS FOR MISSING AND NO MISSING VALUES
740 C FOR EFFICIENCY. SINCE KINC IS USUALLY 1, USING SEPARATE
741 C LOOPS HERE DOESN'T BUY MUCH. A MISSING VALUE WILL ALWAYS
742 C TRANSFER BACK TO GROUP A.
743 C
744  IF(is523.EQ.0)THEN
745 C
746  DO 185 k=ktotal+1,min(ktotal+kinc,nxy)
747  IF(ic(k).LT.minc)THEN
748  minc=ic(k)
749  minck=k
750  ENDIF
751  IF(ic(k).GT.maxc)THEN
752  maxc=ic(k)
753  maxck=k
754  ENDIF
755  nount=nount+1
756  185 CONTINUE
757 C
758  ELSEIF(is523.EQ.1)THEN
759 C
760  DO 187 k=ktotal+1,min(ktotal+kinc,nxy)
761  IF(ic(k).EQ.missp)GO TO 186
762  IF(ic(k).LT.minc)THEN
763  minc=ic(k)
764  minck=k
765  ENDIF
766  IF(ic(k).GT.maxc)THEN
767  maxc=ic(k)
768  maxck=k
769  ENDIF
770  186 nount=nount+1
771  187 CONTINUE
772 C
773  ELSE
774 C
775  DO 190 k=ktotal+1,min(ktotal+kinc,nxy)
776  IF(ic(k).EQ.missp.OR.ic(k).EQ.misss)GO TO 189
777  IF(ic(k).LT.minc)THEN
778  minc=ic(k)
779  minck=k
780  ENDIF
781  IF(ic(k).GT.maxc)THEN
782  maxc=ic(k)
783  maxck=k
784  ENDIF
785  189 nount=nount+1
786  190 CONTINUE
787 C
788  ENDIF
789 C
790 C***D WRITE(KFILDO,191)KOUNTA,KTOTAL,MINA,MAXA,IBITA,MISLLA,
791 C***D 1 MINC,MAXC,NOUNT,IC(KTOTAL),IC(KTOTAL+1)
792 C***D191 FORMAT(' AT 191, KOUNTA ='I8,' KTOTAL ='I8,' MINA ='I8,
793 C***D 1 ' MAXA ='I8,' IBITA ='I3,' MISLLA ='I3,
794 C***D 2 ' MINC ='I8,' MAXC ='I8,
795 C***D 3 ' NOUNT ='I5,' IC(KTOTAL) ='I9,' IC(KTOTAL+1) =',I9)
796 C
797 C IF THE NUMBER OF BITS NEEDED FOR GROUP C IS GT IBITA,
798 C THEN THIS GROUP A IS A GROUP TO PACK.
799 C
800  IF(minc.EQ.mallow)THEN
801  minc=mina
802  maxc=maxa
803  minck=minak
804  maxck=maxak
805  misllc=1
806  GO TO 195
807 C WHEN THE NEW VALUE(S) ARE MISSING, THEY CAN ALWAYS
808 C BE ADDED.
809 C
810  ELSE
811  misllc=0
812  ENDIF
813 C
814  IF(maxc-minc.GE.ibxx2(ibita)-lmiss) GO TO 200
815 C
816 C THE BITS NECESSARY FOR GROUP C HAS NOT INCREASED FROM THE
817 C BITS NECESSARY FOR GROUP A. ADD THIS POINT(S) TO GROUP A.
818 C COMPUTE THE NEXT GROUP B, ETC., UNLESS ALL POINTS HAVE BEEN
819 C USED.
820 C
821  195 ktotal=ktotal+nount
822  kounta=kounta+nount
823  mina=minc
824  maxa=maxc
825  minak=minck
826  maxak=maxck
827  mislla=misllc
828  adda=.true.
829  IF(ktotal.GE.nxy)GO TO 200
830 C
831  IF(minbk.GT.ktotal.AND.maxbk.GT.ktotal)THEN
832  mstart=nendb+1
833 C THE MAX AND MIN OF GROUP B WERE NOT FROM THE POINTS
834 C REMOVED, SO THE WHOLE GROUP DOES NOT HAVE TO BE LOOKED
835 C AT TO DETERMINE THE NEW MAX AND MIN. RATHER START
836 C JUST BEYOND THE OLD NENDB.
837  ibitbs=ibitb
838  nendb=1
839  GO TO 150
840  ELSE
841  GO TO 140
842  ENDIF
843 C
844 C *************************************
845 C
846 C GROUP A IS TO BE PACKED. STORE VALUES IN JMIN( ), JMAX( ),
847 C LBIT( ), AND NOV( ).
848 C
849 C *************************************
850 C
851  200 lx=lx+1
852  IF(lx.LE.ndg)GO TO 205
853  lminpk=lminpk+lminpk/2
854 C WRITE(KFILDO,201)NDG,LMINPK,LX
855 C201 FORMAT(' ****NDG ='I5,' NOT LARGE ENOUGH.',
856 C 1 ' LMINPK IS INCREASED TO 'I3,' FOR THIS FIELD.'/
857 C 2 ' LX = 'I10)
858  iersav=716
859  GO TO 105
860 C
861  205 jmin(lx)=mina
862  jmax(lx)=maxa
863  lbit(lx)=ibita
864  nov(lx)=kounta
865  kstart=ktotal+1
866 C
867  IF(mislla.EQ.0)THEN
868  misslx(lx)=mallow
869  ELSE
870  misslx(lx)=ic(ktotal)
871 C IC(KTOTAL) WAS THE LAST VALUE PROCESSED. IF MISLLA NE 0,
872 C THIS MUST BE THE MISSING VALUE FOR THIS GROUP.
873  ENDIF
874 C
875 C***D WRITE(KFILDO,206)MISLLA,IC(KTOTAL),KTOTAL,LX,JMIN(LX),JMAX(LX),
876 C***D 1 LBIT(LX),NOV(LX),MISSLX(LX)
877 C***D206 FORMAT(' AT 206, MISLLA ='I2,' IC(KTOTAL) ='I5,' KTOTAL ='I8,
878 C***D 1 ' LX ='I6,' JMIN(LX) ='I8,' JMAX(LX) ='I8,
879 C***D 2 ' LBIT(LX) ='I5,' NOV(LX) ='I8,' MISSLX(LX) =',I7)
880 C
881  IF(ktotal.GE.nxy)GO TO 209
882 C
883 C THE NEW GROUP A WILL BE THE PREVIOUS GROUP B. SET LIMITS, ETC.
884 C
885  ibita=ibitb
886  mina=minb
887  maxa=maxb
888  minak=minbk
889  maxak=maxbk
890  mislla=misllb
891  nenda=nendb
892  kounta=kountb
893  ktotal=ktotal+kounta
894  adda=.false.
895  GO TO 133
896 C
897 C *************************************
898 C
899 C CALCULATE IBIT, THE NUMBER OF BITS NEEDED TO HOLD THE GROUP
900 C MINIMUM VALUES.
901 C
902 C *************************************
903 C
904  209 ibit=0
905 C
906  DO 220 l=1,lx
907  210 IF(jmin(l).LT.ibxx2(ibit))GO TO 220
908  ibit=ibit+1
909  GO TO 210
910  220 CONTINUE
911 C
912 C INSERT THE VALUE IN JMIN( ) TO BE USED FOR ALL MISSING
913 C VALUES WHEN LBIT( ) = 0. WHEN SECONDARY MISSING
914 C VALUES CAN BE PRESENT, LBIT(L) WILL NOT = 0.
915 C
916  IF(is523.EQ.1)THEN
917 C
918  DO 226 l=1,lx
919 C
920  IF(lbit(l).EQ.0)THEN
921 C
922  IF(misslx(l).EQ.missp)THEN
923  jmin(l)=ibxx2(ibit)-1
924  ENDIF
925 C
926  ENDIF
927 C
928  226 CONTINUE
929 C
930  ENDIF
931 C
932 C *************************************
933 C
934 C CALCULATE JBIT, THE NUMBER OF BITS NEEDED TO HOLD THE BITS
935 C NEEDED TO PACK THE VALUES IN THE GROUPS. BUT FIND AND
936 C REMOVE THE REFERENCE VALUE FIRST.
937 C
938 C *************************************
939 C
940 C WRITE(KFILDO,228)CFEED,LX
941 C228 FORMAT(A1,/' *****************************************'
942 C 1 /' THE GROUP WIDTHS LBIT( ) FOR ',I8,' GROUPS'
943 C 2 /' *****************************************')
944 C WRITE(KFILDO,229) (LBIT(J),J=1,MIN(LX,100))
945 C229 FORMAT(/' '20I6)
946 C
947  lbitref=lbit(1)
948 C
949  DO 230 k=1,lx
950  IF(lbit(k).LT.lbitref)lbitref=lbit(k)
951  230 CONTINUE
952 C
953  IF(lbitref.NE.0)THEN
954 C
955  DO 240 k=1,lx
956  lbit(k)=lbit(k)-lbitref
957  240 CONTINUE
958 C
959  ENDIF
960 C
961 C WRITE(KFILDO,241)CFEED,LBITREF
962 C241 FORMAT(A1,/' *****************************************'
963 C 1 /' THE GROUP WIDTHS LBIT( ) AFTER REMOVING REFERENCE ',
964 C 2 I8,
965 C 3 /' *****************************************')
966 C WRITE(KFILDO,242) (LBIT(J),J=1,MIN(LX,100))
967 C242 FORMAT(/' '20I6)
968 C
969  jbit=0
970 C
971  DO 320 k=1,lx
972  310 IF(lbit(k).LT.ibxx2(jbit))GO TO 320
973  jbit=jbit+1
974  GO TO 310
975  320 CONTINUE
976 C
977 C *************************************
978 C
979 C CALCULATE KBIT, THE NUMBER OF BITS NEEDED TO HOLD THE NUMBER
980 C OF VALUES IN THE GROUPS. BUT FIND AND REMOVE THE
981 C REFERENCE FIRST.
982 C
983 C *************************************
984 C
985 C WRITE(KFILDO,321)CFEED,LX
986 C321 FORMAT(A1,/' *****************************************'
987 C 1 /' THE GROUP SIZES NOV( ) FOR ',I8,' GROUPS'
988 C 2 /' *****************************************')
989 C WRITE(KFILDO,322) (NOV(J),J=1,MIN(LX,100))
990 C322 FORMAT(/' '20I6)
991 C
992  novref=nov(1)
993 C
994  DO 400 k=1,lx
995  IF(nov(k).LT.novref)novref=nov(k)
996  400 CONTINUE
997 C
998  IF(novref.GT.0)THEN
999 C
1000  DO 405 k=1,lx
1001  nov(k)=nov(k)-novref
1002  405 CONTINUE
1003 C
1004  ENDIF
1005 C
1006 C WRITE(KFILDO,406)CFEED,NOVREF
1007 C406 FORMAT(A1,/' *****************************************'
1008 C 1 /' THE GROUP SIZES NOV( ) AFTER REMOVING REFERENCE ',I8,
1009 C 2 /' *****************************************')
1010 C WRITE(KFILDO,407) (NOV(J),J=1,MIN(LX,100))
1011 C407 FORMAT(/' '20I6)
1012 C WRITE(KFILDO,408)CFEED
1013 C408 FORMAT(A1,/' *****************************************'
1014 C 1 /' THE GROUP REFERENCES JMIN( )'
1015 C 2 /' *****************************************')
1016 C WRITE(KFILDO,409) (JMIN(J),J=1,MIN(LX,100))
1017 C409 FORMAT(/' '20I6)
1018 C
1019  kbit=0
1020 C
1021  DO 420 k=1,lx
1022  410 IF(nov(k).LT.ibxx2(kbit))GO TO 420
1023  kbit=kbit+1
1024  GO TO 410
1025  420 CONTINUE
1026 C
1027 C DETERMINE WHETHER THE GROUP SIZES SHOULD BE REDUCED
1028 C FOR SPACE EFFICIENCY.
1029 C
1030  IF(ired.EQ.0)THEN
1031  CALL reduce(kfildo,jmin,jmax,lbit,nov,lx,ndg,ibit,jbit,kbit,
1032  1 novref,ibxx2,ier)
1033 C
1034  IF(ier.EQ.714.OR.ier.EQ.715)THEN
1035 C REDUCE HAS ABORTED. REEXECUTE PACK_GP WITHOUT REDUCE.
1036 C PROVIDE FOR A NON FATAL RETURN FROM REDUCE.
1037  iersav=ier
1038  ired=1
1039  ier=0
1040  GO TO 102
1041  ENDIF
1042 C
1043  ENDIF
1044 C
1045 C CALL TIMPR(KFILDO,KFILDO,'END PACK_GP ')
1046  IF(iersav.NE.0)THEN
1047  ier=iersav
1048  RETURN
1049  ENDIF
1050 C
1051 C 900 IF(IER.NE.0)RETURN1
1052 C
1053  900 RETURN
1054  END
reduce
subroutine reduce(KFILDO, JMIN, JMAX, LBIT, NOV, LX, NDG, IBIT, JBIT, KBIT, NOVREF, IBXX2, IER)
This subroutine determines whether the number of groups should be increased in order to reduce the si...
Definition: reduce.f:43
pack_gp
subroutine pack_gp(KFILDO, IC, NXY, IS523, MINPK, INC, MISSP, MISSS, JMIN, JMAX, LBIT, NOV, NDG, LX, IBIT, JBIT, KBIT, NOVREF, LBITREF, IER)
This subroutine determines groups of variable size, but at least of size minpk, the associated max(JM...
Definition: pack_gp.f:155