external/io_grib2/g2lib/specunpack.F

   1       subroutine specunpack(cpack,len,idrstmpl,ndpts,JJ,KK,MM,fld)
   2 !$$$  SUBPROGRAM DOCUMENTATION BLOCK
   3 !                .      .    .                                       .
   4 ! SUBPROGRAM:    specunpack
   5 !   PRGMMR: Gilbert          ORG: W/NP11    DATE: 2002-12-19
   6 !
   7 ! ABSTRACT: This subroutine unpacks a spectral data field that was packed
   8 !   using the complex packing algorithm for spherical harmonic data as
   9 !   defined in the GRIB2 documention,
  10 !   using info from the GRIB2 Data Representation Template 5.51.
  11 !
  12 ! PROGRAM HISTORY LOG:
  13 ! 2002-12-19  Gilbert
  14 !
  15 ! USAGE:    CALL specunpack(cpack,len,idrstmpl,ndpts,JJ,KK,MM,fld)
  16 !   INPUT ARGUMENT LIST:
  17 !     cpack    - The packed data field (character*1 array)
  18 !     len      - length of packed field cpack().
  19 !     idrstmpl - Contains the array of values for Data Representation
  20 !                Template 5.51
  21 !     ndpts    - The number of data values to unpack
  22 !     JJ       - J - pentagonal resolution parameter
  23 !     KK       - K - pentagonal resolution parameter
  24 !     MM       - M - pentagonal resolution parameter
  25 !
  26 !   OUTPUT ARGUMENT LIST:
  27 !     fld()    - Contains the unpacked data values
  28 !
  29 ! REMARKS: None
  30 !
  31 ! ATTRIBUTES:
  32 !   LANGUAGE: XL Fortran 90
  33 !   MACHINE:  IBM SP
  34 !
  35 !$$$
  36
  37       character(len=1),intent(in) :: cpack(len)
  38       integer,intent(in) :: ndpts,len,JJ,KK,MM
  39       integer,intent(in) :: idrstmpl(*)
  40       real,intent(out) :: fld(ndpts)
  41
  42       integer :: ifld(ndpts),Ts
  43       integer(4) :: ieee
  44       real :: ref,bscale,dscale,unpk(ndpts)
  45       real,allocatable :: pscale(:)
  46
  47       ieee = idrstmpl(1)
  48       call rdieee(ieee,ref,1)
  49       bscale = 2.0**real(idrstmpl(2))
  50       dscale = 10.0**real(-idrstmpl(3))
  51       nbits = idrstmpl(4)
  52       Js=idrstmpl(6)
  53       Ks=idrstmpl(7)
  54       Ms=idrstmpl(8)
  55       Ts=idrstmpl(9)
  56
  57       if (idrstmpl(10).eq.1) then           ! unpacked floats are 32-bit IEEE
  58          !call g2lib_gbytes(cpack,ifld,0,32,0,Ts)
  59          call rdieee(cpack,unpk,Ts)          ! read IEEE unpacked floats
  60          iofst=32*Ts
  61          call g2lib_gbytes(cpack,ifld,iofst,nbits,0,ndpts-Ts)  ! unpack scaled data
  62 !
  63 !   Calculate Laplacian scaling factors for each possible wave number.
  64 !
  65          allocate(pscale(JJ+MM))
  66          tscale=real(idrstmpl(5))*1E-6
  67          do n=Js,JJ+MM
  68             pscale(n)=real(n*(n+1))**(-tscale)
  69          enddo
  70 !
  71 !   Assemble spectral coeffs back to original order.
  72 !
  73          inc=1
  74          incu=1
  75          incp=1
  76          do m=0,MM
  77             Nm=JJ      ! triangular or trapezoidal
  78             if ( KK .eq. JJ+MM ) Nm=JJ+m          ! rhombodial
  79             Ns=Js      ! triangular or trapezoidal
  80             if ( Ks .eq. Js+Ms ) Ns=Js+m          ! rhombodial
  81             do n=m,Nm
  82                if (n.le.Ns .AND. m.le.Ms) then    ! grab unpacked value
  83                   fld(inc)=unpk(incu)         ! real part
  84                   fld(inc+1)=unpk(incu+1)     ! imaginary part
  85                   inc=inc+2
  86                   incu=incu+2
  87                else                         ! Calc coeff from packed value
  88                   fld(inc)=((real(ifld(incp))*bscale)+ref)*
  89      &                      dscale*pscale(n)           ! real part
  90                   fld(inc+1)=((real(ifld(incp+1))*bscale)+ref)*
  91      &                      dscale*pscale(n)           ! imaginary part
  92                   inc=inc+2
  93                   incp=incp+2
  94                endif
  95             enddo
  96          enddo
  97
  98          deallocate(pscale)
  99
 100       else
 101          print *,'specunpack: Cannot handle 64 or 128-bit floats.'
 102          fld=0.0
 103          return
 104       endif
 105
 106       return
 107       end