Philippe 02/05/2016: moved all LIBTOOLS files in LIBTOOLS directory
[MNH-git_open_source-lfs.git] / LIBTOOLS / lib / COMPRESS / src / compress.f90
1 !-----------------------------------------------------------------
2 !--------------- special set of characters for RCS information
3 !-----------------------------------------------------------------
4 ! $Source$ $Revision$ $Date$
5 !-----------------------------------------------------------------
6 SUBROUTINE COMPRESS_FIELD(XTAB,KX,KY,KNBTOT,KNBUSE)
7 USE MODD_COMPPAR
8 USE MODE_SEARCHGRP
9
10 #ifdef NAGf95
11 USE,INTRINSIC :: IEEE_ARITHMETIC
12 #endif
13
14 IMPLICIT NONE 
15
16 REAL,PARAMETER :: PPFLOATMIN = 2.0**(-126)
17
18 INTEGER, INTENT(IN) :: KX,KY
19 !INTEGER, INTENT(IN) :: KNBLEV
20 INTEGER, INTENT(IN) :: KNBTOT
21 REAL(KIND=8),DIMENSION(KNBTOT),INTENT(INOUT) :: XTAB
22
23 INTEGER, INTENT(OUT) :: KNBUSE
24
25 INTEGER :: INBLEV
26 INTEGER,DIMENSION(:), ALLOCATABLE :: ITAB
27 REAL :: XMIN,XMAX
28 TYPE(SOP_t) :: SOPRES
29 INTEGER :: IND1, IND2
30 INTEGER :: GELT,IBE
31 INTEGER :: ILEVNBELT
32 INTEGER :: NBITCOD
33 INTEGER :: II, JI, JJ
34 INTEGER :: BITOFFSET
35 INTEGER :: GRPIDX,GRPOFF,IDXSAVE,OFFSAVE
36 INTEGER :: nbgroupmod
37 INTEGER :: IEXTCOD
38 CHARACTER(LEN=8),PARAMETER :: KEYWORD='COMPRESS'
39 REAL,DIMENSION(KNBTOT) :: XWORKTAB
40 LOGICAL :: LUPREAL,LNAN
41 #ifndef NAGf95
42 LOGICAL, EXTERNAL :: IEEE_IS_NAN
43 #endif
44
45 ILEVNBELT = KX*KY
46 LUPREAL = .FALSE.
47 LNAN    = .FALSE.
48
49 ! Check for NAN and change Upper and Lower bound according to 32bits real limits.
50 DO JI=1,KNBTOT
51   IF (IEEE_IS_NAN(XTAB(JI))) THEN 
52     XTAB(JI)=0.
53     LNAN = .TRUE.
54   ELSE IF (ABS(XTAB(JI)) > HUGE(1.0_4)) THEN
55     XTAB(JI) = SIGN(REAL(HUGE(1.0_4)/1.1,8),XTAB(JI))
56     LUPREAL = .TRUE.
57   ELSEIF (ABS(XTAB(JI)) < TINY(1.0_4)) THEN
58     XTAB(JI) = 0.
59   END IF
60 END DO
61
62 XMIN=MINVAL(XTAB)
63 XMAX=MAXVAL(XTAB)
64 PRINT *,'MINVAL,MAXVAL= ',XMIN,XMAX
65 IF (LNAN) PRINT *,"==================> NAN values DETECTED : set to 0.0"
66 IF (LUPREAL) PRINT *,"==================> OVERFLOW values DETECTED : set to ",HUGE(1.0_4)/1.1
67
68 ! Convert 64 bits real to 32 bits real
69 XWORKTAB(:) = XTAB(:)
70 !
71 ! BEWARE : Now XTAB is overwritten. 
72 !          XWORKTAB contains the 32 bits floating point data.
73 !
74 CALL SET_FILLIDX(0,0)
75 ! store 8 characters header string in buffer
76 DO II=1,LEN(KEYWORD)
77   CALL FILL_BBUFF(XTAB,8,ICHAR(KEYWORD(II:II)))
78 END DO
79
80 ! is whole array XTAB64 a constant field ?
81
82 IF (xmin == xmax) THEN
83   PRINT *,"--------> CONSTANT ARRAY !"
84   CALL FILL_BBUFF(XTAB,32,JPCSTENCOD)
85   CALL FILL_BBUFF(XTAB,32,KNBTOT)
86   CALL FILL_BBUFF(XTAB,32,xmin)
87   CALL GET_FILLIDX(KNBUSE,BITOFFSET)
88   KNBUSE=KNBUSE+1
89   RETURN
90 END IF
91
92
93 INBLEV = KNBTOT/(ILEVNBELT)
94 IF (KNBTOT /= (INBLEV*ILEVNBELT)) THEN
95   PRINT *,'Pb in COMPRESS_FIELD : KNBTOT must be a multiple of KX*KY'
96   STOP
97 END IF
98
99
100
101 ALLOCATE(ITAB(ILEVNBELT))
102 CALL INI_SOPDATA(SOPRES)
103
104 CALL FILL_BBUFF(XTAB,32,JPEXTENCOD)
105 CALL FILL_BBUFF(XTAB,32,KNBTOT)
106 CALL FILL_BBUFF(XTAB,32,KX)
107 CALL FILL_BBUFF(XTAB,32,KY)
108
109 DO JI=1,INBLEV
110   IND1=(JI-1)*ILEVNBELT+1
111   IND2=JI*ILEVNBELT
112   IF (LPDEBUG) PRINT *,"---- Compressing Level ",JI," ----"
113   CALL COMP_FOPEXT(XWORKTAB(IND1:IND2),ITAB,IEXTCOD)
114   IF (IEXTCOD /= JPCONST) THEN
115     CALL INVERTCOL(ITAB,KX,KY)
116     CALL RECSEARCH(ITAB,SOPRES)
117     GELT = MAXVAL(SOPRES%IEND(1:SOPRES%NBGRP)-SOPRES%IBEG(1:SOPRES%NBGRP)+1)
118     IBE = FMINBITS_IN_WORD(GELT)
119     CALL GET_FILLIDX(GRPIDX,GRPOFF) ! save the idx/offset for future NBGRP modification
120     CALL FILL_BBUFF(XTAB,32,SOPRES%NBGRP)
121     CALL FILL_BBUFF(XTAB,5,IBE)
122   
123     NBGROUPMOD = SOPRES%NBGRP
124     DO II=1,SOPRES%NBGRP
125       GELT = SOPRES%IEND(II)-SOPRES%IBEG(II)+1
126       nbitcod  = FMINBITS_IN_WORD(SOPRES%VALMAX(II)-SOPRES%VALMIN(II))
127       !    PRINT *, 'Groupe',II,'(',GELT,')',':',SOPRES%IBEG(II),SOPRES%IEND(II),&
128       !         &'MIN,MAX=',SOPRES%VALMIN(II),SOPRES%VALMAX(II),&
129       !         &'(',SOPRES%VALMAX(II)-SOPRES%VALMIN(II),'/',&
130       !         &nbitcod,')'
131       IF (nbitcod >= 16) THEN
132         PRINT *,'-----> ERREUR FATALE : Groupe',II,'codage sur ',nbitcod,'bits'
133       END IF
134       IF (GELT > 1) THEN
135         ! Plus d'un element dans le groupe
136         IF ((17*GELT) < (17+4+IBE+nbitcod*GELT)) THEN
137           ! on prefere GELT groupes de 1 elt
138           DO JJ=SOPRES%IBEG(II),SOPRES%IEND(II)
139             ! 1 seul elt par groupe
140             CALL FILL_BBUFF(XTAB,1,1)
141             CALL FILL_BBUFF(XTAB,16,ITAB(JJ))
142           END DO
143           NBGROUPMOD = NBGROUPMOD+GELT-1
144         ELSE
145           CALL FILL_BBUFF(XTAB,1,0)
146           CALL FILL_BBUFF(XTAB,16,SOPRES%VALMIN(II))
147           CALL FILL_BBUFF(XTAB,4,nbitcod)
148           CALL FILL_BBUFF(XTAB,IBE,GELT)
149           IF (nbitcod > 0) THEN
150             DO JJ=SOPRES%IBEG(II),SOPRES%IEND(II)
151               ! stockage des GELT écarts/VALMIN
152               CALL FILL_BBUFF(XTAB,nbitcod,ITAB(JJ)-SOPRES%VALMIN(II))
153             END DO
154           END IF
155         END IF
156       ELSE
157         ! 1 seul elt dans groupe
158         CALL FILL_BBUFF(XTAB,1,1)
159         CALL FILL_BBUFF(XTAB,16,SOPRES%VALMIN(II))
160       END IF
161     END DO
162     IF (NBGROUPMOD > SOPRES%NBGRP) THEN
163       ! we must change the number of elements 
164       CALL GET_FILLIDX(IDXSAVE,OFFSAVE) ! save the current idx/offset
165       CALL SET_FILLIDX(GRPIDX,GRPOFF)   
166       CALL FILL_BBUFF(XTAB,32,NBGROUPMOD)
167       CALL SET_FILLIDX(IDXSAVE,OFFSAVE) ! restore the current idx/offset
168     END IF
169   END IF
170 END DO
171
172 CALL GET_FILLIDX(IDXSAVE,OFFSAVE)
173 KNBUSE=IDXSAVE+1
174
175 DEALLOCATE(ITAB)
176
177 CONTAINS 
178
179 SUBROUTINE COMP_FOPEXT(PTAB,KTAB,KEXTCOD)
180 REAL,    DIMENSION(:), INTENT(IN) :: PTAB
181 INTEGER, DIMENSION(:), INTENT(OUT):: KTAB 
182 INTEGER,               INTENT(OUT):: KEXTCOD
183
184 LOGICAL,DIMENSION(SIZE(PTAB)) :: GMASK
185 REAL,DIMENSION(SIZE(PTAB)) :: PTABWORK
186 REAL :: XMIN1,XMAX1,XRANGE1
187 REAL :: XMIN2,XMAX2,XRANGE2
188 REAL :: XREF,XMAX,XCOEFF
189 INTEGER :: INTRANGE
190 INTEGER :: INDCOR   ! correction d'index pour la supression du min
191 LOGICAL :: GMINEXCL,GMAXEXCL,GLOG
192 INTEGER :: IEXTCOD2
193 REAL, PARAMETER :: XUNDEF = 999.
194 REAL, PARAMETER :: XUNDEFSURF =  1.E+20
195
196
197 !! G. TANGUY avril 2010 : on change la valeur indéfinie 999. a une valeur
198 !indéfinie plus grande que sera de façon certaine le max du champ s'il est
199 !present. POur ça on travaille dans le tableau de travail PTABWORK
200 PTABWORK=PTAB
201 WHERE(PTABWORK == XUNDEF)
202         PTABWORK=XUNDEFSURF
203 END WHERE
204
205 XMIN1=MINVAL(PTABWORK(:))
206 XMAX1=MAXVAL(PTABWORK(:))
207 XRANGE1=XMAX1-XMIN1
208 IF (LPDEBUG) PRINT *,"XMIN1,XMAX1,XRANGE1 = ",XMIN1,XMAX1,XRANGE1
209
210 IF (XRANGE1 > 0.) THEN
211   XMIN2=MINVAL(PTABWORK,MASK=PTABWORK>XMIN1)
212   XMAX2=MAXVAL(PTABWORK,MASK=PTABWORK<XMAX1)
213   XRANGE2 = XMAX2-XMIN2
214   IF (LPDEBUG) PRINT *,"XMIN2,XMAX2,XRANGE2 = ",XMIN2,XMAX2,XRANGE2
215   IF (XRANGE2 > 0.) THEN
216     GLOG     = .FALSE.
217     GMINEXCL = .FALSE.
218     GMAXEXCL = .FALSE.
219     GMASK(:) = .TRUE.
220     INDCOR = 0
221     KEXTCOD = JPNORM
222     INTRANGE=65535
223     XREF = XMIN1
224     XMAX = XMAX1
225
226     ! Check for range between 0 and 1 to convert to LOG values
227     IF (XMIN1 >= 0. .AND. XMAX1 < 1.) THEN
228       IF ((XMAX2/XMIN2)>10.) THEN 
229         GLOG = .TRUE.
230         KEXTCOD = JPOTHER
231         IEXTCOD2 = JPLOG
232         INTRANGE=INTRANGE-1
233         INDCOR = 1           ! On reserve la valeur 0 dans tous les cas
234         IF (XMIN1 == 0.0) THEN
235           XREF = LOG(XMIN2)
236           WHERE (PTABWORK < XMIN2)
237             KTAB  = 0
238             GMASK = .FALSE.
239           END WHERE
240         ELSE
241           XREF = LOG(XMIN1)
242         END IF
243         XMAX1 = LOG(XMAX1)
244         XMAX  = XMAX1
245         XMAX2 = LOG(XMAX2)
246         XRANGE2 = XMAX2 - XREF
247         IF (LPDEBUG) PRINT *,"EXTENCOD,  LOG conversion enabled : XMIN1, XREF, XMAX1, XMAX2 =",&
248              &XMIN1,XREF,XMAX1,XMAX2
249       END IF
250     ELSE
251       ! Check for MIN value exclusion
252       IF (XMIN1 == XUNDEFSURF .OR. (XMIN2-XMIN1) > XRANGE2) THEN
253         ! Min value excluded 
254         GMINEXCL = .TRUE.
255         XREF=XMIN2
256         INTRANGE=INTRANGE-1
257         INDCOR = 1
258         WHERE (PTABWORK < XMIN2)
259           KTAB = 0
260           GMASK = .FALSE.
261         END WHERE
262         IF (LPDEBUG) PRINT *,"EXTENCOD,     Min value isolated :",XMIN1
263         KEXTCOD = JPMINEXCL
264         IF (XMIN1 == XUNDEFSURF) THEN 
265                 XMIN1=XUNDEF
266         END IF
267       END IF
268       ! Check for MAX value exclusion
269       IF (XMAX1 == XUNDEFSURF .OR. (XMAX1-XMAX2) > XRANGE2) THEN
270         ! Max value excluded
271         GMAXEXCL = .TRUE.
272         XMAX=XMAX2
273         INTRANGE=INTRANGE-1
274         WHERE (PTABWORK > XMAX2)
275           KTAB = 65535
276           GMASK = .FALSE.
277         END WHERE
278         
279         IF (GMINEXCL) THEN
280           KEXTCOD = JPMINMAXEXCL ! Min et Max exclus
281           IF (LPDEBUG) PRINT *,"EXTENCOD, and Max value isolated :",XMAX1
282         ELSE
283           KEXTCOD = JPMAXEXCL ! Max exclus
284           IF (LPDEBUG) PRINT *,"EXTENCOD, Max value isolated :",XMAX1
285         END IF
286         ! avril 2010 : on remet la valeur indefine de mesonh 999.
287         IF (XMAX1 == XUNDEFSURF) THEN 
288                 XMAX1=XUNDEF
289         END IF
290       END IF
291     END IF
292     !
293     XCOEFF=(XMAX-XREF)/INTRANGE
294     IF (XCOEFF < PPFLOATMIN) THEN
295       XCOEFF = PPFLOATMIN
296       PRINT *, "very low range DATA : XCOEFF set to",XCOEFF
297     END IF
298     IF (LPDEBUG) PRINT *,"XCOEFF = ",XCOEFF
299     IF (GLOG) THEN
300       WHERE(GMASK)
301         KTAB = INDCOR + NINT((LOG(PTABWORK)-XREF)/XCOEFF)
302       END WHERE
303     ELSE
304       WHERE(GMASK)
305         KTAB = INDCOR + NINT((PTABWORK(:)-XREF)/XCOEFF)
306       END WHERE
307     END IF
308     IF (LPDEBUG) PRINT *,"KEXTCOD = ",KEXTCOD
309     CALL FILL_BBUFF(XTAB,3,KEXTCOD)
310     IF (GLOG)     CALL FILL_BBUFF(XTAB,3,IEXTCOD2)
311     IF (GMINEXCL) CALL FILL_BBUFF(XTAB,32,XMIN1)
312     IF (GMAXEXCL) CALL FILL_BBUFF(XTAB,32,XMAX1)
313     CALL FILL_BBUFF(XTAB,32,XREF)
314     CALL FILL_BBUFF(XTAB,32,XCOEFF)
315   ELSE
316     IF (XRANGE2 < 0.) THEN
317       ! only 2 values in PTAB array
318       !
319       ! KTAB(i)= 0 if PTAB(i)==XMIN1
320       !          1 if PTAB(i)==XMAX1
321       !
322       IF (LPDEBUG) PRINT *,"EXTENCOD, 2 values in array :",XMIN1,XMAX1
323         IF (XMAX1 == XUNDEFSURF) THEN 
324                 XMAX1=XUNDEF
325         END IF
326         IF (XMIN1 == XUNDEFSURF) THEN 
327                 XMIN1=XUNDEF
328         END IF
329       KEXTCOD = JP2VAL
330       CALL FILL_BBUFF(XTAB,3,KEXTCOD)
331       CALL FILL_BBUFF(XTAB,32,XMIN1)
332       CALL FILL_BBUFF(XTAB,32,XMAX1)
333       WHERE (PTABWORK < XMAX1)
334         KTAB = 0
335       ELSEWHERE
336         KTAB = 1
337       END WHERE
338     ELSE
339       ! XRANGE2 == 0. <==> XMIN2=XMAX2 
340       ! 3 values in PTAB array :
341       !
342       !          0 if PTAB(i)==XMIN1      ! KTAB(i)= 1 if PTAB(i)==XMIN2(=XMAX2)
343       !          2 if PTAB(i)==XMAX1
344       !
345       IF (LPDEBUG) PRINT *,"EXTENCOD, 3 values in array :",XMIN1,XMIN2,XMAX1
346         IF (XMAX1 == XUNDEFSURF) THEN 
347                 XMAX1=XUNDEF
348         END IF
349         IF (XMIN1 == XUNDEFSURF) THEN 
350                 XMIN1=XUNDEF
351         END IF
352
353       KEXTCOD = JP3VAL
354       CALL FILL_BBUFF(XTAB,3,KEXTCOD)
355       CALL FILL_BBUFF(XTAB,32,XMIN1)
356       CALL FILL_BBUFF(XTAB,32,XMIN2)
357       CALL FILL_BBUFF(XTAB,32,XMAX1)
358       WHERE (PTABWORK < XMIN2)
359         KTAB = 0
360       ELSEWHERE
361         KTAB = 1
362       END WHERE
363       WHERE (PTABWORK > XMIN2) KTAB = 2
364     END IF
365     
366   END IF
367 ELSE
368         IF (XMIN1 == XUNDEFSURF) THEN 
369                 XMIN1=XUNDEF
370         END IF
371
372   ! Constant array found : save its 32 bits real value.
373   KEXTCOD=JPCONST
374   CALL FILL_BBUFF(XTAB,3,KEXTCOD)
375   CALL FILL_BBUFF(XTAB,32,XMIN1)
376   IF (LPDEBUG) PRINT *,"EXTENCOD, constant array : ",XMIN1
377 END IF
378 END SUBROUTINE COMP_FOPEXT
379
380 END SUBROUTINE COMPRESS_FIELD