external/io_grib2/g2lib/specpack.F

   1       subroutine specpack(fld,ndpts,JJ,KK,MM,idrstmpl,cpack,lcpack)
   2 !$$$  SUBPROGRAM DOCUMENTATION BLOCK
   3 !                .      .    .                                       .
   4 ! SUBPROGRAM:    specpack
   5 !   PRGMMR: Gilbert          ORG: W/NP11    DATE: 2002-12-19
   6 !
   7 ! ABSTRACT: This subroutine packs a spectral data field using the complex
   8 !   packing algorithm for spherical harmonic data as
   9 !   defined in the GRIB2 Data Representation Template 5.51.
  10 !
  11 ! PROGRAM HISTORY LOG:
  12 ! 2002-12-19  Gilbert
  13 !
  14 ! USAGE:    CALL specpack(fld,ndpts,JJ,KK,MM,idrstmpl,cpack,lcpack)
  15 !   INPUT ARGUMENT LIST:
  16 !     fld()    - Contains the packed data values
  17 !     ndpts    - The number of data values to pack
  18 !     JJ       - J - pentagonal resolution parameter
  19 !     KK       - K - pentagonal resolution parameter
  20 !     MM       - M - pentagonal resolution parameter
  21 !     idrstmpl - Contains the array of values for Data Representation
  22 !                Template 5.51
  23 !
  24 !   OUTPUT ARGUMENT LIST:
  25 !     cpack    - The packed data field (character*1 array)
  26 !     lcpack   - length of packed field cpack().
  27 !
  28 ! REMARKS: None
  29 !
  30 ! ATTRIBUTES:
  31 !   LANGUAGE: XL Fortran 90
  32 !   MACHINE:  IBM SP
  33 !
  34 !$$$
  35
  36       real,intent(in) :: fld(ndpts)
  37       integer,intent(in) :: ndpts,JJ,KK,MM
  38       integer,intent(inout) :: idrstmpl(*)
  39       character(len=1),intent(out) :: cpack(*)
  40       integer,intent(out) :: lcpack
  41
  42       integer :: ifld(ndpts),Ts,tmplsim(5)
  43       real :: bscale,dscale,unpk(ndpts),tfld(ndpts)
  44       real,allocatable :: pscale(:)
  45
  46       bscale = 2.0**real(-idrstmpl(2))
  47       dscale = 10.0**real(idrstmpl(3))
  48       nbits = idrstmpl(4)
  49       Js=idrstmpl(6)
  50       Ks=idrstmpl(7)
  51       Ms=idrstmpl(8)
  52       Ts=idrstmpl(9)
  53
  54 !
  55 !   Calculate Laplacian scaling factors for each possible wave number.
  56 !
  57       allocate(pscale(JJ+MM))
  58       tscale=real(idrstmpl(5))*1E-6
  59       do n=Js,JJ+MM
  60          pscale(n)=real(n*(n+1))**(tscale)
  61       enddo
  62 !
  63 !   Separate spectral coeffs into two lists; one to contain unpacked
  64 !   values within the sub-spectrum Js, Ks, Ms, and the other with values
  65 !   outside of the sub-spectrum to be packed.
  66 !
  67       inc=1
  68       incu=1
  69       incp=1
  70       do m=0,MM
  71          Nm=JJ      ! triangular or trapezoidal
  72          if ( KK .eq. JJ+MM ) Nm=JJ+m          ! rhombodial
  73          Ns=Js      ! triangular or trapezoidal
  74          if ( Ks .eq. Js+Ms ) Ns=Js+m          ! rhombodial
  75          do n=m,Nm
  76             if (n.le.Ns .AND. m.le.Ms) then    ! save unpacked value
  77                unpk(incu)=fld(inc)         ! real part
  78                unpk(incu+1)=fld(inc+1)     ! imaginary part
  79                inc=inc+2
  80                incu=incu+2
  81             else                         ! Save value to be packed and scale
  82                                          ! Laplacian scale factor
  83                tfld(incp)=fld(inc)*pscale(n)         ! real part
  84                tfld(incp+1)=fld(inc+1)*pscale(n)     ! imaginary part
  85                inc=inc+2
  86                incp=incp+2
  87             endif
  88          enddo
  89       enddo
  90
  91       deallocate(pscale)
  92
  93       incu=incu-1
  94       if (incu .ne. Ts) then
  95          print *,'specpack: Incorrect number of unpacked values ',
  96      &           'given:',Ts
  97          print *,'specpack: Resetting idrstmpl(9) to ',incu
  98          Ts=incu
  99       endif
 100 !
 101 !  Add unpacked values to the packed data array in 32-bit IEEE format
 102 !
 103       call mkieee(unpk,cpack,Ts)
 104       ipos=4*Ts
 105 !
 106 !  Scale and pack the rest of the coefficients
 107 !
 108       tmplsim(2)=idrstmpl(2)
 109       tmplsim(3)=idrstmpl(3)
 110       tmplsim(4)=idrstmpl(4)
 111       call simpack(tfld,ndpts-Ts,tmplsim,cpack(ipos+1),lcpack)
 112       lcpack=lcpack+ipos
 113 !
 114 !  Fill in Template 5.51
 115 !
 116       idrstmpl(1)=tmplsim(1)
 117       idrstmpl(2)=tmplsim(2)
 118       idrstmpl(3)=tmplsim(3)
 119       idrstmpl(4)=tmplsim(4)
 120       idrstmpl(9)=Ts
 121       idrstmpl(10)=1         ! Unpacked spectral data is 32-bit IEEE
 122
 123       return
 124       end