Skip to content
Snippets Groups Projects
Normal view Permalink
decompress.f90 8.45 KiB
Newer Older
WAUTELET Philippe's avatar
Beginning of open source history
WAUTELET Philippe committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 
!-----------------------------------------------------------------
!--------------- special set of characters for RCS information
!-----------------------------------------------------------------
! $Source$ $Revision$ $Date$
!-----------------------------------------------------------------
SUBROUTINE GET_COMPHEADER(KTAB,SIZEKTAB,KNBELT,KTYPECOD)

INTEGER, INTENT(IN) :: SIZEKTAB
INTEGER(KIND=8), DIMENSION(SIZEKTAB), INTENT(IN) :: KTAB
INTEGER, INTENT(OUT) :: KNBELT    ! size of decompressed array
INTEGER, INTENT(OUT) :: KTYPECOD  ! code for compression type

CHARACTER(LEN=8) :: STRKEY

INTEGER :: INTCHAR
INTEGER :: JI

CALL SET_EXTRACTIDX(0,0)
! extract string header 
DO JI=1,8
  CALL EXTRACT_BBUFF(KTAB,8,INTCHAR)
  STRKEY(JI:JI) = CHAR(INTCHAR)
END DO

! Treat array if it is compressed
IF (STRKEY == 'COMPRESS') THEN
  CALL EXTRACT_BBUFF(KTAB,32,KTYPECOD)
  CALL EXTRACT_BBUFF(KTAB,32,KNBELT)
ELSE
  KNBELT    =-1
  KTYPECOD = 0
END IF

END SUBROUTINE GET_COMPHEADER

SUBROUTINE DECOMPRESS_FIELD(XTAB,NBELT,COMPTAB,NBCOMPELT,CODINGTYPE)
USE MODD_COMPPAR
USE MODE_SEARCHGRP

IMPLICIT NONE 
INTEGER,                                INTENT(IN)  :: NBELT 
INTEGER,                                INTENT(IN)  :: NBCOMPELT 
REAL   (KIND=8),DIMENSION(NBELT),TARGET,INTENT(OUT) :: XTAB
INTEGER(KIND=8),DIMENSION(NBCOMPELT),   INTENT(IN)  :: COMPTAB
INTEGER,                                INTENT(IN)  :: CODINGTYPE

INTEGER,DIMENSION(:), ALLOCATABLE  :: ITAB
LOGICAL,DIMENSION(:), ALLOCATABLE  :: GMASK

REAL :: XREF, XCOEFF
INTEGER :: INBLEV
INTEGER :: ILEVNBELT
INTEGER :: JI
INTEGER :: IND1, IND2
INTEGER :: IDIMX,IDIMY
INTEGER :: IEXTCOD
REAL(KIND=8),DIMENSION(:),POINTER  :: XPTRTAB
REAL :: XMIN,XMAX

SELECT CASE (CODINGTYPE)
CASE (JPCSTENCOD)
  CALL EXTRACT_BBUFF(COMPTAB,32,XREF)
  XTAB(:) = XREF

CASE (JPSOPENCOD)
  CALL EXTRACT_BBUFF(COMPTAB,32,IDIMX)
  CALL EXTRACT_BBUFF(COMPTAB,32,IDIMY)
  ILEVNBELT = IDIMX * IDIMY
  INBLEV = NBELT/(ILEVNBELT)
  ALLOCATE(ITAB(ILEVNBELT))
  DO JI=1,INBLEV
    IND1=(JI-1)*ILEVNBELT+1
    IND2=JI*ILEVNBELT
    XPTRTAB=>XTAB(IND1:IND2)
    IF (LPDEBUG) PRINT *,'######   Decompress(SOPENCOD) LEVEL ',JI,'######'
    CALL EXTRACT_BBUFF(COMPTAB,32,XREF)
    CALL EXTRACT_BBUFF(COMPTAB,32,XCOEFF)
    CALL EXTRACTINTARRAY(ITAB)
    CALL DECOMP_FOP(XPTRTAB,ITAB,XREF,XCOEFF)
  END DO
  
CASE (JPEXTENCOD)
  CALL EXTRACT_BBUFF(COMPTAB,32,IDIMX)
  CALL EXTRACT_BBUFF(COMPTAB,32,IDIMY)
  ILEVNBELT = IDIMX * IDIMY
  INBLEV = NBELT/(ILEVNBELT)
  ALLOCATE(ITAB(ILEVNBELT))
  ALLOCATE(GMASK(ILEVNBELT))
  DO JI=1,INBLEV

    IF (LPDEBUG) PRINT *,'###### Decompress(EXTENCOD) LEVEL ',JI,'######'
    IND1=(JI-1)*ILEVNBELT+1
    IND2=JI*ILEVNBELT
    XPTRTAB=>XTAB(IND1:IND2)
    !
    CALL EXTRACT_BBUFF(COMPTAB,3,IEXTCOD)
    IF (IEXTCOD == JPOTHER) THEN
      CALL EXTRACT_BBUFF(COMPTAB,3,IEXTCOD)
      IEXTCOD = IEXTCOD + 8
    END IF
    IF (LPDEBUG) PRINT *, "IEXTCOD = ",IEXTCOD
    SELECT CASE(IEXTCOD)
    CASE(JPLOG)
      ! Conversion to log values of original data 0<=x<1
      CALL EXTRACT_BBUFF(COMPTAB,32,XREF)
      CALL EXTRACT_BBUFF(COMPTAB,32,XCOEFF)
      CALL EXTRACTINTARRAY(ITAB)
      GMASK(:) = .TRUE.
      WHERE (ITAB == 0)
        GMASK = .FALSE.
        XPTRTAB = 0.0
      END WHERE
      CALL DECOMP_FOP(XPTRTAB,ITAB,XREF,XCOEFF,GMASK,1)
      WHERE(GMASK)
        XPTRTAB = EXP(XPTRTAB)
      END WHERE
      
    CASE(JPCONST)
      ! constant value array
      CALL EXTRACT_BBUFF(COMPTAB,32,XREF)
      XPTRTAB(:) = XREF
      IF (LPDEBUG) PRINT *,"  CONST value=",XREF

    CASE(JP2VAL)
      ! 2 different values in array
      CALL EXTRACT_BBUFF(COMPTAB,32,XMIN)
      CALL EXTRACT_BBUFF(COMPTAB,32,XMAX)
      CALL EXTRACTINTARRAY(ITAB)
      WHERE (ITAB == 0)
        XPTRTAB = XMIN
      ELSEWHERE
        XPTRTAB = XMAX
      END WHERE
      IF (LPDEBUG) PRINT *,"  2 values:",XMIN,XMAX
      
    CASE(JP3VAL)
      ! 3 different values in array
      CALL EXTRACT_BBUFF(COMPTAB,32,XMIN)
      CALL EXTRACT_BBUFF(COMPTAB,32,XREF)
      CALL EXTRACT_BBUFF(COMPTAB,32,XMAX)
      CALL EXTRACTINTARRAY(ITAB)
      WHERE (ITAB == 0)
        XPTRTAB = XMIN
      ELSEWHERE
        XPTRTAB = XREF
      END WHERE
      WHERE (ITAB == 2) XPTRTAB = XMAX
      IF (LPDEBUG) PRINT *,"  3 values:",XMIN,XREF,XMAX

    CASE(JPNORM)
      ! same as JPSOPENCOD
      CALL EXTRACT_BBUFF(COMPTAB,32,XREF)
      CALL EXTRACT_BBUFF(COMPTAB,32,XCOEFF)
      CALL EXTRACTINTARRAY(ITAB)
      CALL DECOMP_FOP(XPTRTAB,ITAB,XREF,XCOEFF)
      IF (LPDEBUG) PRINT *,"  normal, XREF/XCOEFF = ",XREF,XCOEFF 

    CASE(JPMINEXCL)
      ! Min value is isolated
      CALL EXTRACT_BBUFF(COMPTAB,32,XMIN)
      CALL EXTRACT_BBUFF(COMPTAB,32,XREF)
      CALL EXTRACT_BBUFF(COMPTAB,32,XCOEFF)
      CALL EXTRACTINTARRAY(ITAB)
      GMASK(:) = .TRUE.
      WHERE (ITAB == 0)
        GMASK = .FALSE.
        XPTRTAB = XMIN
      END WHERE
      CALL DECOMP_FOP(XPTRTAB,ITAB,XREF,XCOEFF,GMASK,1)
      IF (LPDEBUG) PRINT *,"  Min exclus, MIN/XREF/XCOEFF = ",XMIN,XREF,XCOEFF

    CASE(JPMAXEXCL)
      ! Max value is isolated
      CALL EXTRACT_BBUFF(COMPTAB,32,XMAX)
      CALL EXTRACT_BBUFF(COMPTAB,32,XREF)
      CALL EXTRACT_BBUFF(COMPTAB,32,XCOEFF)
      CALL EXTRACTINTARRAY(ITAB)
      GMASK(:) = .TRUE.
      WHERE (ITAB == 65535)
        GMASK = .FALSE.
        XPTRTAB = XMAX
      END WHERE
      CALL DECOMP_FOP(XPTRTAB,ITAB,XREF,XCOEFF,GMASK,0)   
      IF (LPDEBUG) PRINT *,"  Max exclus, MAX/XREF/XCOEFF = ",XMAX,XREF,XCOEFF

    CASE(JPMINMAXEXCL)
      ! Min&Max value are isolated
      CALL EXTRACT_BBUFF(COMPTAB,32,XMIN)        
      CALL EXTRACT_BBUFF(COMPTAB,32,XMAX)
      CALL EXTRACT_BBUFF(COMPTAB,32,XREF)
      CALL EXTRACT_BBUFF(COMPTAB,32,XCOEFF)
      CALL EXTRACTINTARRAY(ITAB)
      GMASK(:) = .TRUE.
      WHERE (ITAB == 0)
        GMASK = .FALSE.
        XPTRTAB = XMIN
      END WHERE
      WHERE (ITAB == 65535)
        GMASK = .FALSE.
        XPTRTAB = XMAX
      END WHERE
      CALL DECOMP_FOP(XPTRTAB,ITAB,XREF,XCOEFF,GMASK,1)
      IF (LPDEBUG) PRINT *,"  Min et Max exclus, MIN/MAX/XREF/XCOEFF = ",&
           &XMIN,XMAX,XREF,XCOEFF
    END SELECT
  END DO
  
CASE DEFAULT
  PRINT *,'Error in CODINGTYPE : program aborted'
  STOP
END SELECT

CONTAINS 

SUBROUTINE DECOMP_FOP(PTAB,KTAB,PREF,PCOEFF,OMASK,KINDCOR)
REAL(KIND=8), DIMENSION(:), INTENT(INOUT) :: PTAB 
! Attention: avec le compilateur PGF, utiliser INTENT(OUT) provoque une recopie
! complete du tableau dans PTAB (avec ecrasement possible des valeurs 
! presentes a l'appel de la procedure). Le phenomene est genant lorsque
! DECOMP_FOP ne calcule que sur une portion de PTAB (valeurs min et/ou max 
! sont presentes). En declarant PTAB en INOUT, les valeurs en entree de la routine
! sont conservees si elles n'ont pas ete modifiees.

INTEGER,      DIMENSION(:), INTENT(IN) :: KTAB 
REAL, INTENT(IN) :: PREF
REAL, INTENT(IN) :: PCOEFF
LOGICAL, DIMENSION(:),INTENT(IN),OPTIONAL :: OMASK
INTEGER,INTENT(IN),OPTIONAL  :: KINDCOR ! 1 if Min value is isolated, 0 otherwise

INTEGER :: INDCOR

IF (.NOT. PRESENT(KINDCOR)) THEN
  INDCOR = 0
ELSE
  INDCOR = KINDCOR
END IF
  
IF (PRESENT(OMASK)) THEN
  WHERE (OMASK)
    PTAB(:) = PCOEFF*(KTAB(:)-INDCOR)+PREF
  END WHERE
ELSE
  IF (PCOEFF == 0.0) THEN
    PTAB(:) = PREF
  ELSE
    PTAB(:) = PCOEFF*KTAB(:)+PREF
  END IF
END IF

END SUBROUTINE DECOMP_FOP

SUBROUTINE EXTRACTINTARRAY(KTAB)
INTEGER,DIMENSION(:),INTENT(OUT) :: KTAB
!
! COMPTAB, IDIMX and IDIMY  are defined in the calling routine
!
INTEGER :: NBGRP
INTEGER :: IBE
INTEGER :: CPT
INTEGER :: JJ
INTEGER :: ALONE
INTEGER :: NBITCOD,IMIN
INTEGER :: GELT
INTEGER :: JELT
INTEGER :: IEPS

CALL EXTRACT_BBUFF(COMPTAB,32,NBGRP)
!      PRINT *,'Nbre de groupes =',NBGRP
CALL EXTRACT_BBUFF(COMPTAB,5,IBE)
!      PRINT *,'Nbre de bits pour coder le nombre d''elements:',IBE
CPT = 1
DO JJ=1,NBGRP
  !      PRINT *,'Groupe ',JJ,' : '
  CALL EXTRACT_BBUFF(COMPTAB,1,ALONE)
  CALL EXTRACT_BBUFF(COMPTAB,16,IMIN)
  !      PRINT *,'IREF=',IMIN
  
  IF (ALONE == 1) THEN
    ! 1 seul elt dans le groupe
    !        PRINT *,'--> un seul element dans le groupe'
    KTAB(CPT)=IMIN
    CPT=CPT+1
  ELSE
    CALL EXTRACT_BBUFF(COMPTAB,4,NBITCOD)
    CALL EXTRACT_BBUFF(COMPTAB,IBE,GELT)
    !        PRINT *,'--> ',GELT,' elts, codage ecart sur ',nbitcod,'bits'
    IF (NBITCOD > 0) THEN
      DO JELT=1,GELT
        CALL EXTRACT_BBUFF(COMPTAB,NBITCOD,IEPS)
        KTAB(CPT) = IMIN+IEPS
        CPT=CPT+1
      END DO
    ELSE
      KTAB(CPT:CPT+GELT-1) = IMIN
      CPT = CPT+GELT
    END IF
  END IF
END DO
CALL INVERTCOL(KTAB,IDIMX,IDIMY)        
END SUBROUTINE EXTRACTINTARRAY

END SUBROUTINE DECOMPRESS_FIELD

Documentation | Mentions légales | CGU | Accessibilité : non conforme