SURFEX v7.3: /home/dasprezs/EXPORT_v7_3/src/LIB/LFI_COMPRESS/src/compress

Go to the documentation of this file.
00001 !-----------------------------------------------------------------
00002 !--------------- special set of characters for RCS information
00003 !-----------------------------------------------------------------
00004 ! $Source: /home/cvsroot/mesonh/libtools/lib/COMPRESS/src/compress.f90,v $ $Revision: 1.1.2.1 $ $Date: 2007/03/02 13:15:26 $
00005 !-----------------------------------------------------------------
00006 SUBROUTINE COMPRESS_FIELD(XTAB,KX,KY,KNBTOT,KNBUSE)
00007 USE MODD_COMPPAR
00008 USE MODE_SEARCHGRP
00009 
00010 #ifdef NAGf95
00011 USE,INTRINSIC :: IEEE_ARITHMETIC
00012 #endif
00013 
00014 IMPLICIT NONE 
00015 
00016 REAL,PARAMETER :: PPFLOATMIN = 2.0**(-126)
00017 
00018 INTEGER, INTENT(IN) :: KX,KY
00019 !INTEGER, INTENT(IN) :: KNBLEV
00020 INTEGER, INTENT(IN) :: KNBTOT
00021 REAL(KIND=8),DIMENSION(KNBTOT),INTENT(INOUT) :: XTAB
00022 
00023 INTEGER, INTENT(OUT) :: KNBUSE
00024 
00025 INTEGER :: INBLEV
00026 INTEGER,DIMENSION(:), ALLOCATABLE :: ITAB
00027 REAL :: XMIN,XMAX
00028 TYPE(SOP_t) :: SOPRES
00029 INTEGER :: IND1, IND2
00030 INTEGER :: GELT,IBE
00031 INTEGER :: ILEVNBELT
00032 INTEGER :: NBITCOD
00033 INTEGER :: II, JI, JJ
00034 INTEGER :: BITOFFSET
00035 INTEGER :: GRPIDX,GRPOFF,IDXSAVE,OFFSAVE
00036 INTEGER :: nbgroupmod
00037 INTEGER :: IEXTCOD
00038 CHARACTER(LEN=8),PARAMETER :: KEYWORD='COMPRESS'
00039 REAL,DIMENSION(KNBTOT) :: XWORKTAB
00040 LOGICAL :: LUPREAL,LNAN
00041 #ifndef NAGf95
00042 LOGICAL, EXTERNAL :: IEEE_IS_NAN
00043 #endif
00044 
00045 ILEVNBELT = KX*KY
00046 LUPREAL = .FALSE.
00047 LNAN    = .FALSE.
00048 
00049 ! Check for NAN and change Upper and Lower bound according to 32bits real limits.
00050 DO JI=1,KNBTOT
00051   IF (IEEE_IS_NAN(XTAB(JI))) THEN 
00052     XTAB(JI)=0.
00053     LNAN = .TRUE.
00054   ELSE IF (ABS(XTAB(JI)) > HUGE(1.0_4)) THEN
00055     XTAB(JI) = SIGN(REAL(HUGE(1.0_4)/1.1,8),XTAB(JI))
00056     LUPREAL = .TRUE.
00057   ELSEIF (ABS(XTAB(JI)) < TINY(1.0_4)) THEN
00058     XTAB(JI) = 0.
00059   END IF
00060 END DO
00061 
00062 XMIN=MINVAL(XTAB)
00063 XMAX=MAXVAL(XTAB)
00064 PRINT *,'MINVAL,MAXVAL= ',XMIN,XMAX
00065 IF (LNAN) PRINT *,"==================> NAN values DETECTED : set to 0.0"
00066 IF (LUPREAL) PRINT *,"==================> OVERFLOW values DETECTED : set to ",HUGE(1.0_4)/1.1
00067 
00068 ! Convert 64 bits real to 32 bits real
00069 XWORKTAB(:) = XTAB(:)
00070 !
00071 ! BEWARE : Now XTAB is overwritten. 
00072 !          XWORKTAB contains the 32 bits floating point data.
00073 !
00074 CALL SET_FILLIDX(0,0)
00075 ! store 8 characters header string in buffer
00076 DO II=1,LEN(KEYWORD)
00077   CALL FILL_BBUFF(XTAB,8,ICHAR(KEYWORD(II:II)))
00078 END DO
00079 
00080 ! is whole array XTAB64 a constant field ?
00081 
00082 IF (xmin == xmax) THEN
00083   PRINT *,"--------> CONSTANT ARRAY !"
00084   CALL FILL_BBUFF(XTAB,32,JPCSTENCOD)
00085   CALL FILL_BBUFF(XTAB,32,KNBTOT)
00086   CALL FILL_BBUFF(XTAB,32,xmin)
00087   CALL GET_FILLIDX(KNBUSE,BITOFFSET)
00088   KNBUSE=KNBUSE+1
00089   RETURN
00090 END IF
00091 
00092 
00093 INBLEV = KNBTOT/(ILEVNBELT)
00094 IF (KNBTOT /= (INBLEV*ILEVNBELT)) THEN
00095   PRINT *,'Pb in COMPRESS_FIELD : KNBTOT must be a multiple of KX*KY'
00096   STOP
00097 END IF
00098 
00099 
00100 
00101 ALLOCATE(ITAB(ILEVNBELT))
00102 CALL INI_SOPDATA(SOPRES)
00103 
00104 CALL FILL_BBUFF(XTAB,32,JPEXTENCOD)
00105 CALL FILL_BBUFF(XTAB,32,KNBTOT)
00106 CALL FILL_BBUFF(XTAB,32,KX)
00107 CALL FILL_BBUFF(XTAB,32,KY)
00108 
00109 DO JI=1,INBLEV
00110   IND1=(JI-1)*ILEVNBELT+1
00111   IND2=JI*ILEVNBELT
00112   IF (LPDEBUG) PRINT *,"---- Compressing Level ",JI," ----"
00113   CALL COMP_FOPEXT(XWORKTAB(IND1:IND2),ITAB,IEXTCOD)
00114   IF (IEXTCOD /= JPCONST) THEN
00115     CALL INVERTCOL(ITAB,KX,KY)
00116     CALL RECSEARCH(ITAB,SOPRES)
00117     GELT = MAXVAL(SOPRES%IEND(1:SOPRES%NBGRP)-SOPRES%IBEG(1:SOPRES%NBGRP)+1)
00118     IBE = FMINBITS_IN_WORD(GELT)
00119     CALL GET_FILLIDX(GRPIDX,GRPOFF) ! save the idx/offset for future NBGRP modification
00120     CALL FILL_BBUFF(XTAB,32,SOPRES%NBGRP)
00121     CALL FILL_BBUFF(XTAB,5,IBE)
00122   
00123     NBGROUPMOD = SOPRES%NBGRP
00124     DO II=1,SOPRES%NBGRP
00125       GELT = SOPRES%IEND(II)-SOPRES%IBEG(II)+1
00126       nbitcod  = FMINBITS_IN_WORD(SOPRES%VALMAX(II)-SOPRES%VALMIN(II))
00127       !    PRINT *, 'Groupe',II,'(',GELT,')',':',SOPRES%IBEG(II),SOPRES%IEND(II),&
00128       !         &'MIN,MAX=',SOPRES%VALMIN(II),SOPRES%VALMAX(II),&
00129       !         &'(',SOPRES%VALMAX(II)-SOPRES%VALMIN(II),'/',&
00130       !         &nbitcod,')'
00131       IF (nbitcod >= 16) THEN
00132         PRINT *,'-----> ERREUR FATALE : Groupe',II,'codage sur ',nbitcod,'bits'
00133       END IF
00134       IF (GELT > 1) THEN
00135         ! Plus d'un element dans le groupe
00136         IF ((17*GELT) < (17+4+IBE+nbitcod*GELT)) THEN
00137           ! on prefere GELT groupes de 1 elt
00138           DO JJ=SOPRES%IBEG(II),SOPRES%IEND(II)
00139             ! 1 seul elt par groupe
00140             CALL FILL_BBUFF(XTAB,1,1)
00141             CALL FILL_BBUFF(XTAB,16,ITAB(JJ))
00142           END DO
00143           NBGROUPMOD = NBGROUPMOD+GELT-1
00144         ELSE
00145           CALL FILL_BBUFF(XTAB,1,0)
00146           CALL FILL_BBUFF(XTAB,16,SOPRES%VALMIN(II))
00147           CALL FILL_BBUFF(XTAB,4,nbitcod)
00148           CALL FILL_BBUFF(XTAB,IBE,GELT)
00149           IF (nbitcod > 0) THEN
00150             DO JJ=SOPRES%IBEG(II),SOPRES%IEND(II)
00151               ! stockage des GELT �carts/VALMIN
00152               CALL FILL_BBUFF(XTAB,nbitcod,ITAB(JJ)-SOPRES%VALMIN(II))
00153             END DO
00154           END IF
00155         END IF
00156       ELSE
00157         ! 1 seul elt dans groupe
00158         CALL FILL_BBUFF(XTAB,1,1)
00159         CALL FILL_BBUFF(XTAB,16,SOPRES%VALMIN(II))
00160       END IF
00161     END DO
00162     IF (NBGROUPMOD > SOPRES%NBGRP) THEN
00163       ! we must change the number of elements 
00164       CALL GET_FILLIDX(IDXSAVE,OFFSAVE) ! save the current idx/offset
00165       CALL SET_FILLIDX(GRPIDX,GRPOFF)   
00166       CALL FILL_BBUFF(XTAB,32,NBGROUPMOD)
00167       CALL SET_FILLIDX(IDXSAVE,OFFSAVE) ! restore the current idx/offset
00168     END IF
00169   END IF
00170 END DO
00171 
00172 CALL GET_FILLIDX(IDXSAVE,OFFSAVE)
00173 KNBUSE=IDXSAVE+1
00174 
00175 DEALLOCATE(ITAB)
00176 
00177 CONTAINS 
00178 
00179 SUBROUTINE COMP_FOPEXT(PTAB,KTAB,KEXTCOD)
00180 REAL,    DIMENSION(:), INTENT(IN) :: PTAB
00181 INTEGER, DIMENSION(:), INTENT(OUT):: KTAB 
00182 INTEGER,               INTENT(OUT):: KEXTCOD
00183 
00184 LOGICAL,DIMENSION(SIZE(PTAB)) :: GMASK
00185 REAL :: XMIN1,XMAX1,XRANGE1
00186 REAL :: XMIN2,XMAX2,XRANGE2
00187 REAL :: XREF,XMAX,XCOEFF
00188 INTEGER :: INTRANGE
00189 INTEGER :: INDCOR   ! correction d'index pour la supression du min
00190 LOGICAL :: GMINEXCL,GMAXEXCL,GLOG
00191 INTEGER :: IEXTCOD2
00192 
00193 XMIN1=MINVAL(PTAB(:))
00194 XMAX1=MAXVAL(PTAB(:))
00195 XRANGE1=XMAX1-XMIN1
00196 IF (LPDEBUG) PRINT *,"XMIN1,XMAX1,XRANGE1 = ",XMIN1,XMAX1,XRANGE1
00197 
00198 IF (XRANGE1 > 0.) THEN
00199   XMIN2=MINVAL(PTAB,MASK=PTAB>XMIN1)
00200   XMAX2=MAXVAL(PTAB,MASK=PTAB<XMAX1)
00201   XRANGE2 = XMAX2-XMIN2
00202   IF (LPDEBUG) PRINT *,"XMIN2,XMAX2,XRANGE2 = ",XMIN2,XMAX2,XRANGE2
00203   IF (XRANGE2 > 0.) THEN
00204     GLOG     = .FALSE.
00205     GMINEXCL = .FALSE.
00206     GMAXEXCL = .FALSE.
00207     GMASK(:) = .TRUE.
00208     INDCOR = 0
00209     KEXTCOD = JPNORM
00210     INTRANGE=65535
00211     XREF = XMIN1
00212     XMAX = XMAX1
00213 
00214     ! Check for range between 0 and 1 to convert to LOG values
00215     IF (XMIN1 >= 0. .AND. XMAX1 < 1.) THEN
00216       IF ((XMAX2/XMIN2)>10.) THEN 
00217         GLOG = .TRUE.
00218         KEXTCOD = JPOTHER
00219         IEXTCOD2 = JPLOG
00220         INTRANGE=INTRANGE-1
00221         INDCOR = 1           ! On reserve la valeur 0 dans tous les cas
00222         IF (XMIN1 == 0.0) THEN
00223           XREF = LOG(XMIN2)
00224           WHERE (PTAB < XMIN2)
00225             KTAB  = 0
00226             GMASK = .FALSE.
00227           END WHERE
00228         ELSE
00229           XREF = LOG(XMIN1)
00230         END IF
00231         XMAX1 = LOG(XMAX1)
00232         XMAX  = XMAX1
00233         XMAX2 = LOG(XMAX2)
00234         XRANGE2 = XMAX2 - XREF
00235         IF (LPDEBUG) PRINT *,"EXTENCOD,  LOG conversion enabled : XMIN1, XREF, XMAX1, XMAX2 =",&
00236              &XMIN1,XREF,XMAX1,XMAX2
00237       END IF
00238     ELSE
00239       ! Check for MIN value exclusion
00240       IF ((XMIN2-XMIN1) > XRANGE2) THEN
00241         ! Min value excluded 
00242         GMINEXCL = .TRUE.
00243         XREF=XMIN2
00244         INTRANGE=INTRANGE-1
00245         INDCOR = 1
00246         WHERE (PTAB < XMIN2)
00247           KTAB = 0
00248           GMASK = .FALSE.
00249         END WHERE
00250         IF (LPDEBUG) PRINT *,"EXTENCOD,     Min value isolated :",XMIN1
00251         KEXTCOD = JPMINEXCL
00252       END IF
00253       ! Check for MAX value exclusion
00254       IF ((XMAX1-XMAX2) > XRANGE2) THEN
00255         ! Max value excluded
00256         GMAXEXCL = .TRUE.
00257         XMAX=XMAX2
00258         INTRANGE=INTRANGE-1
00259         WHERE (PTAB > XMAX2)
00260           KTAB = 65535
00261           GMASK = .FALSE.
00262         END WHERE
00263         
00264         IF (GMINEXCL) THEN
00265           KEXTCOD = JPMINMAXEXCL ! Min et Max exclus
00266           IF (LPDEBUG) PRINT *,"EXTENCOD, and Max value isolated :",XMAX1
00267         ELSE
00268           KEXTCOD = JPMAXEXCL ! Max exclus
00269           IF (LPDEBUG) PRINT *,"EXTENCOD, Max value isolated :",XMAX1
00270         END IF
00271       END IF
00272     END IF
00273     !
00274     XCOEFF=(XMAX-XREF)/INTRANGE
00275     IF (XCOEFF < PPFLOATMIN) THEN
00276       XCOEFF = PPFLOATMIN
00277       PRINT *, "very low range DATA : XCOEFF set to",XCOEFF
00278     END IF
00279     IF (LPDEBUG) PRINT *,"XCOEFF = ",XCOEFF
00280     IF (GLOG) THEN
00281       WHERE(GMASK)
00282         KTAB = INDCOR + NINT((LOG(PTAB)-XREF)/XCOEFF)
00283       END WHERE
00284     ELSE
00285       WHERE(GMASK)
00286         KTAB = INDCOR + NINT((PTAB(:)-XREF)/XCOEFF)
00287       END WHERE
00288     END IF
00289     IF (LPDEBUG) PRINT *,"KEXTCOD = ",KEXTCOD
00290     CALL FILL_BBUFF(XTAB,3,KEXTCOD)
00291     IF (GLOG)     CALL FILL_BBUFF(XTAB,3,IEXTCOD2)
00292     IF (GMINEXCL) CALL FILL_BBUFF(XTAB,32,XMIN1)
00293     IF (GMAXEXCL) CALL FILL_BBUFF(XTAB,32,XMAX1)
00294     CALL FILL_BBUFF(XTAB,32,XREF)
00295     CALL FILL_BBUFF(XTAB,32,XCOEFF)
00296   ELSE
00297     IF (XRANGE2 < 0.) THEN
00298       ! only 2 values in PTAB array
00299       !
00300       ! KTAB(i)= 0 if PTAB(i)==XMIN1
00301       !          1 if PTAB(i)==XMAX1
00302       !
00303       IF (LPDEBUG) PRINT *,"EXTENCOD, 2 values in array :",XMIN1,XMAX1
00304       KEXTCOD = JP2VAL
00305       CALL FILL_BBUFF(XTAB,3,KEXTCOD)
00306       CALL FILL_BBUFF(XTAB,32,XMIN1)
00307       CALL FILL_BBUFF(XTAB,32,XMAX1)
00308       WHERE (PTAB < XMAX1) 
00309         KTAB = 0
00310       ELSEWHERE
00311         KTAB = 1
00312       END WHERE
00313     ELSE
00314       ! XRANGE2 == 0. <==> XMIN2=XMAX2 
00315       ! 3 values in PTAB array :
00316       !
00317       !          0 if PTAB(i)==XMIN1      ! KTAB(i)= 1 if PTAB(i)==XMIN2(=XMAX2)
00318       !          2 if PTAB(i)==XMAX1
00319       !
00320       IF (LPDEBUG) PRINT *,"EXTENCOD, 3 values in array :",XMIN1,XMIN2,XMAX1
00321       KEXTCOD = JP3VAL
00322       CALL FILL_BBUFF(XTAB,3,KEXTCOD)
00323       CALL FILL_BBUFF(XTAB,32,XMIN1)
00324       CALL FILL_BBUFF(XTAB,32,XMIN2)
00325       CALL FILL_BBUFF(XTAB,32,XMAX1)
00326       WHERE (PTAB < XMIN2)
00327         KTAB = 0
00328       ELSEWHERE
00329         KTAB = 1
00330       END WHERE
00331       WHERE (PTAB > XMIN2) KTAB = 2
00332     END IF
00333     
00334   END IF
00335 ELSE
00336   ! Constant array found : save its 32 bits real value.
00337   KEXTCOD=JPCONST
00338   CALL FILL_BBUFF(XTAB,3,KEXTCOD)
00339   CALL FILL_BBUFF(XTAB,32,XMIN1)
00340   IF (LPDEBUG) PRINT *,"EXTENCOD, constant array : ",XMIN1
00341 END IF
00342 END SUBROUTINE COMP_FOPEXT
00343 
00344 END SUBROUTINE COMPRESS_FIELD