compress_field.f90

Go to the documentation of this file.

!-----------------------------------------------------------------
!--------------- special set of characters for RCS information
!-----------------------------------------------------------------
! $Source: /home/cvsroot/mesonh/libtools/lib/COMPRESS/src/compress.f90,v $ $Revision: 1.1.2.1 $ $Date: 2007/03/02 13:15:26 $
!-----------------------------------------------------------------
SUBROUTINE compress_field(XTAB,KX,KY,KNBTOT,KNBUSE)
USE modd_comppar
USE mode_searchgrp

#ifdef NAGf95
use,INTRINSIC :: ieee_arithmetic
#endif

IMPLICIT NONE 

REAL,PARAMETER :: PPFLOATMIN = 2.0**(-126)

INTEGER, INTENT(IN) :: KX,KY
!INTEGER, INTENT(IN) :: KNBLEV
INTEGER, INTENT(IN) :: KNBTOT
REAL(KIND=8),DIMENSION(KNBTOT),INTENT(INOUT) :: XTAB

INTEGER, INTENT(OUT) :: KNBUSE

INTEGER :: INBLEV
INTEGER,DIMENSION(:), ALLOCATABLE :: ITAB
REAL :: XMIN,XMAX
TYPE(sop_t) :: SOPRES
INTEGER :: IND1, IND2
INTEGER :: GELT,IBE
INTEGER :: ILEVNBELT
INTEGER :: NBITCOD
INTEGER :: II, JI, JJ
INTEGER :: BITOFFSET
INTEGER :: GRPIDX,GRPOFF,IDXSAVE,OFFSAVE
INTEGER :: nbgroupmod
INTEGER :: IEXTCOD
 CHARACTER(LEN=8),PARAMETER :: KEYWORD='COMPRESS'
REAL,DIMENSION(KNBTOT) :: XWORKTAB
LOGICAL :: LUPREAL,LNAN
#ifndef NAGf95
LOGICAL, EXTERNAL :: IEEE_IS_NAN
#endif

ilevnbelt = kx*ky
lupreal = .false.
lnan    = .false.

! Check for NAN and change Upper and Lower bound according to 32bits real limits.
DO ji=1,knbtot
  IF (ieee_is_nan(xtab(ji))) THEN 
    xtab(ji)=0.
    lnan = .true.
  ELSE IF (abs(xtab(ji)) > huge(1.0_4)) THEN
    xtab(ji) = sign(REAL(HUGE(1.0_4)/1.1,8),XTAB(ji))
    lupreal = .true.
  ELSEIF (abs(xtab(ji)) < tiny(1.0_4)) THEN
    xtab(ji) = 0.
  END IF
END DO

xmin=minval(xtab)
xmax=maxval(xtab)
print *,'MINVAL,MAXVAL= ',xmin,xmax
IF (lnan) print *,"==================> NAN values DETECTED : set to 0.0"
IF (lupreal) print *,"==================> OVERFLOW values DETECTED : set to ",huge(1.0_4)/1.1

! Convert 64 bits real to 32 bits real
xworktab(:) = xtab(:)
!
! BEWARE : Now XTAB is overwritten. 
!          XWORKTAB contains the 32 bits floating point data.
!
 CALL set_fillidx(0,0)
! store 8 characters header string in buffer
DO ii=1,len(keyword)
  CALL fill_bbuff(xtab,8,ichar(keyword(ii:ii)))
END DO

! is whole array XTAB64 a constant field ?

IF (xmin == xmax) THEN
  print *,"--------> CONSTANT ARRAY !"
  CALL fill_bbuff(xtab,32,jpcstencod)
  CALL fill_bbuff(xtab,32,knbtot)
  CALL fill_bbuff(xtab,32,xmin)
  CALL get_fillidx(knbuse,bitoffset)
  knbuse=knbuse+1
  RETURN
END IF


inblev = knbtot/(ilevnbelt)
IF (knbtot /= (inblev*ilevnbelt)) THEN
  print *,'Pb in COMPRESS_FIELD : KNBTOT must be a multiple of KX*KY'
  stop
END IF



ALLOCATE(itab(ilevnbelt))
 CALL ini_sopdata(sopres)

 CALL fill_bbuff(xtab,32,jpextencod)
 CALL fill_bbuff(xtab,32,knbtot)
 CALL fill_bbuff(xtab,32,kx)
 CALL fill_bbuff(xtab,32,ky)

DO ji=1,inblev
  ind1=(ji-1)*ilevnbelt+1
  ind2=ji*ilevnbelt
  IF (lpdebug) print *,"---- Compressing Level ",ji," ----"
  CALL comp_fopext(xworktab(ind1:ind2),itab,iextcod)
  IF (iextcod /= jpconst) THEN
    CALL invertcol(itab,kx,ky)
    CALL recsearch(itab,sopres)
    gelt = maxval(sopres%IEND(1:sopres%NBGRP)-sopres%IBEG(1:sopres%NBGRP)+1)
    ibe = fminbits_in_word(gelt)
    CALL get_fillidx(grpidx,grpoff) ! save the idx/offset for future NBGRP modification
    CALL fill_bbuff(xtab,32,sopres%NBGRP)
    CALL fill_bbuff(xtab,5,ibe)
  
    nbgroupmod = sopres%NBGRP
    DO ii=1,sopres%NBGRP
      gelt = sopres%IEND(ii)-sopres%IBEG(ii)+1
      nbitcod  = fminbits_in_word(sopres%VALMAX(ii)-sopres%VALMIN(ii))
      !    PRINT *, 'Groupe',II,'(',GELT,')',':',SOPRES%IBEG(II),SOPRES%IEND(II),&
      !         &'MIN,MAX=',SOPRES%VALMIN(II),SOPRES%VALMAX(II),&
      !         &'(',SOPRES%VALMAX(II)-SOPRES%VALMIN(II),'/',&
      !         &nbitcod,')'
      IF (nbitcod >= 16) THEN
        print *,'-----> ERREUR FATALE : Groupe',ii,'codage sur ',nbitcod,'bits'
      END IF
      IF (gelt > 1) THEN
        ! Plus d'un element dans le groupe
        IF ((17*gelt) < (17+4+ibe+nbitcod*gelt)) THEN
          ! on prefere GELT groupes de 1 elt
          DO jj=sopres%IBEG(ii),sopres%IEND(ii)
            ! 1 seul elt par groupe
            CALL fill_bbuff(xtab,1,1)
            CALL fill_bbuff(xtab,16,itab(jj))
          END DO
          nbgroupmod = nbgroupmod+gelt-1
        ELSE
          CALL fill_bbuff(xtab,1,0)
          CALL fill_bbuff(xtab,16,sopres%VALMIN(ii))
          CALL fill_bbuff(xtab,4,nbitcod)
          CALL fill_bbuff(xtab,ibe,gelt)
          IF (nbitcod > 0) THEN
            DO jj=sopres%IBEG(ii),sopres%IEND(ii)
              ! stockage des GELT écarts/VALMIN
              CALL fill_bbuff(xtab,nbitcod,itab(jj)-sopres%VALMIN(ii))
            END DO
          END IF
        END IF
      ELSE
        ! 1 seul elt dans groupe
        CALL fill_bbuff(xtab,1,1)
        CALL fill_bbuff(xtab,16,sopres%VALMIN(ii))
      END IF
    END DO
    IF (nbgroupmod > sopres%NBGRP) THEN
      ! we must change the number of elements 
      CALL get_fillidx(idxsave,offsave) ! save the current idx/offset
      CALL set_fillidx(grpidx,grpoff)   
      CALL fill_bbuff(xtab,32,nbgroupmod)
      CALL set_fillidx(idxsave,offsave) ! restore the current idx/offset
    END IF
  END IF
END DO

 CALL get_fillidx(idxsave,offsave)
knbuse=idxsave+1

DEALLOCATE(itab)

CONTAINS 

SUBROUTINE comp_fopext(PTAB,KTAB,KEXTCOD)
REAL,    DIMENSION(:), INTENT(IN) :: PTAB
INTEGER, DIMENSION(:), INTENT(OUT):: KTAB 
INTEGER,               INTENT(OUT):: KEXTCOD

LOGICAL,DIMENSION(SIZE(PTAB)) :: GMASK
REAL :: XMIN1,XMAX1,XRANGE1
REAL :: XMIN2,XMAX2,XRANGE2
REAL :: XREF,XMAX,XCOEFF
INTEGER :: INTRANGE
INTEGER :: INDCOR   ! correction d'index pour la supression du min
LOGICAL :: GMINEXCL,GMAXEXCL,GLOG
INTEGER :: IEXTCOD2

xmin1=minval(ptab(:))
xmax1=maxval(ptab(:))
xrange1=xmax1-xmin1
IF (lpdebug) print *,"XMIN1,XMAX1,XRANGE1 = ",xmin1,xmax1,xrange1

IF (xrange1 > 0.) THEN
  xmin2=minval(ptab,mask=ptab>xmin1)
  xmax2=maxval(ptab,mask=ptab<xmax1)
  xrange2 = xmax2-xmin2
  IF (lpdebug) print *,"XMIN2,XMAX2,XRANGE2 = ",xmin2,xmax2,xrange2
  IF (xrange2 > 0.) THEN
    glog     = .false.
    gminexcl = .false.
    gmaxexcl = .false.
    gmask(:) = .true.
    indcor = 0
    kextcod = jpnorm
    intrange=65535
    xref = xmin1
    xmax = xmax1

    ! Check for range between 0 and 1 to convert to LOG values
    IF (xmin1 >= 0. .AND. xmax1 < 1.) THEN
      IF ((xmax2/xmin2)>10.) THEN 
        glog = .true.
        kextcod = jpother
        iextcod2 = jplog
        intrange=intrange-1
        indcor = 1           ! On reserve la valeur 0 dans tous les cas
        IF (xmin1 == 0.0) THEN
          xref = log(xmin2)
          WHERE (ptab < xmin2)
            ktab  = 0
            gmask = .false.
          END WHERE
        ELSE
          xref = log(xmin1)
        END IF
        xmax1 = log(xmax1)
        xmax  = xmax1
        xmax2 = log(xmax2)
        xrange2 = xmax2 - xref
        IF (lpdebug) print *,"EXTENCOD,  LOG conversion enabled : XMIN1, XREF, XMAX1, XMAX2 =",&
             &xmin1,xref,xmax1,xmax2
      END IF
    ELSE
      ! Check for MIN value exclusion
      IF ((xmin2-xmin1) > xrange2) THEN
        ! Min value excluded 
        gminexcl = .true.
        xref=xmin2
        intrange=intrange-1
        indcor = 1
        WHERE (ptab < xmin2)
          ktab = 0
          gmask = .false.
        END WHERE
        IF (lpdebug) print *,"EXTENCOD,     Min value isolated :",xmin1
        kextcod = jpminexcl
      END IF
      ! Check for MAX value exclusion
      IF ((xmax1-xmax2) > xrange2) THEN
        ! Max value excluded
        gmaxexcl = .true.
        xmax=xmax2
        intrange=intrange-1
        WHERE (ptab > xmax2)
          ktab = 65535
          gmask = .false.
        END WHERE
        
        IF (gminexcl) THEN
          kextcod = jpminmaxexcl ! Min et Max exclus
          IF (lpdebug) print *,"EXTENCOD, and Max value isolated :",xmax1
        ELSE
          kextcod = jpmaxexcl ! Max exclus
          IF (lpdebug) print *,"EXTENCOD, Max value isolated :",xmax1
        END IF
      END IF
    END IF
    !
    xcoeff=(xmax-xref)/intrange
    IF (xcoeff < ppfloatmin) THEN
      xcoeff = ppfloatmin
      print *, "very low range DATA : XCOEFF set to",xcoeff
    END IF
    IF (lpdebug) print *,"XCOEFF = ",xcoeff
    IF (glog) THEN
      WHERE(gmask)
        ktab = indcor + nint((log(ptab)-xref)/xcoeff)
      END WHERE
    ELSE
      WHERE(gmask)
        ktab = indcor + nint((ptab(:)-xref)/xcoeff)
      END WHERE
    END IF
    IF (lpdebug) print *,"KEXTCOD = ",kextcod
    CALL fill_bbuff(xtab,3,kextcod)
    IF (glog)     CALL fill_bbuff(xtab,3,iextcod2)
    IF (gminexcl) CALL fill_bbuff(xtab,32,xmin1)
    IF (gmaxexcl) CALL fill_bbuff(xtab,32,xmax1)
    CALL fill_bbuff(xtab,32,xref)
    CALL fill_bbuff(xtab,32,xcoeff)
  ELSE
    IF (xrange2 < 0.) THEN
      ! only 2 values in PTAB array
      !
      ! KTAB(i)= 0 if PTAB(i)==XMIN1
      !          1 if PTAB(i)==XMAX1
      !
      IF (lpdebug) print *,"EXTENCOD, 2 values in array :",xmin1,xmax1
      kextcod = jp2val
      CALL fill_bbuff(xtab,3,kextcod)
      CALL fill_bbuff(xtab,32,xmin1)
      CALL fill_bbuff(xtab,32,xmax1)
      WHERE (ptab < xmax1) 
        ktab = 0
      ELSEWHERE
        ktab = 1
      END WHERE
    ELSE
      ! XRANGE2 == 0. <==> XMIN2=XMAX2 
      ! 3 values in PTAB array :
      !
      !          0 if PTAB(i)==XMIN1      ! KTAB(i)= 1 if PTAB(i)==XMIN2(=XMAX2)
      !          2 if PTAB(i)==XMAX1
      !
      IF (lpdebug) print *,"EXTENCOD, 3 values in array :",xmin1,xmin2,xmax1
      kextcod = jp3val
      CALL fill_bbuff(xtab,3,kextcod)
      CALL fill_bbuff(xtab,32,xmin1)
      CALL fill_bbuff(xtab,32,xmin2)
      CALL fill_bbuff(xtab,32,xmax1)
      WHERE (ptab < xmin2)
        ktab = 0
      ELSEWHERE
        ktab = 1
      END WHERE
      WHERE (ptab > xmin2) ktab = 2
    END IF
    
  END IF
ELSE
  ! Constant array found : save its 32 bits real value.
  kextcod=jpconst
  CALL fill_bbuff(xtab,3,kextcod)
  CALL fill_bbuff(xtab,32,xmin1)
  IF (lpdebug) print *,"EXTENCOD, constant array : ",xmin1
END IF
END SUBROUTINE comp_fopext

END SUBROUTINE compress_field