Add test for specpack()/specunpack() (#517)

* spec test

* spec test

* fixing test

* adding test

* turned off test for spec packing with _d library

* fixing memory problems in test

* trying to fix _d test run

src/specpack.f

@@ -1,11 +1,12 @@
 !>    @file
 !>    @brief This subroutine packs up a spectral data field.
 !>    @author Stephen Gilbert @date 2002-12-19
-!>
 
 !>    This subroutine packs a spectral data field using the complex
 !>    packing algorithm for spherical harmonic data as defined in the
-!>    GRIB2 Data Representation Template 5.51.
+!>    GRIB2 [Data Representation Template
+!>    5.51](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp5-51.shtml).
+!>
 !>    @param[in] fld Contains the data values to pack.
 !>    @param[in] ndpts The number of data values in array fld.
 !>    @param[in] JJ J pentagonal resolution parameter.
@@ -17,8 +18,6 @@
 !>    @param[out] lcpack length of packed field cpack.
 !>
 !>    @author Stephen Gilbert @date 2002-12-19
-!>
-
       subroutine specpack(fld,ndpts,JJ,KK,MM,idrstmpl,cpack,lcpack)
 
       real,intent(in) :: fld(ndpts)
@@ -54,17 +53,17 @@ subroutine specpack(fld,ndpts,JJ,KK,MM,idrstmpl,cpack,lcpack)
       incp=1
       do m=0,MM
          Nm=JJ      ! triangular or trapezoidal
-         if ( KK .eq. JJ+MM ) Nm=JJ+m          ! rhombodial
+         if (KK .eq. JJ+MM) Nm=JJ+m          ! rhombodial
          Ns=Js      ! triangular or trapezoidal
-         if ( Ks .eq. Js+Ms ) Ns=Js+m          ! rhombodial
+         if (Ks .eq. Js+Ms) Ns=Js+m          ! rhombodial
          do n=m,Nm
             if (n.le.Ns .AND. m.le.Ms) then    ! save unpacked value
                unpk(incu)=fld(inc)         ! real part
                unpk(incu+1)=fld(inc+1)     ! imaginary part
                inc=inc+2
                incu=incu+2
-            else                         ! Save value to be packed and scale
-                                         ! Laplacian scale factor
+            else
+               ! Save value to be packed and scale Laplacian scale factor.
                tfld(incp)=fld(inc)*pscale(n)         ! real part
                tfld(incp+1)=fld(inc+1)*pscale(n)     ! imaginary part
                inc=inc+2
@@ -83,23 +82,23 @@ subroutine specpack(fld,ndpts,JJ,KK,MM,idrstmpl,cpack,lcpack)
          Ts=incu
       endif
 
-!     Add unpacked values to the packed data array in 32-bit IEEE format
+!     Add unpacked values to the packed data array in 32-bit IEEE format.
       call mkieee(unpk,cpack,Ts)
       ipos=4*Ts
 
-!     Scale and pack the rest of the coefficients
+!     Scale and pack the rest of the coefficients.
       tmplsim(2)=idrstmpl(2)
       tmplsim(3)=idrstmpl(3)
       tmplsim(4)=idrstmpl(4)
       call simpack(tfld,ndpts-Ts,tmplsim,cpack(ipos+1),lcpack)
       lcpack=lcpack+ipos
 
-!     Fill in Template 5.51
+!     Fill in Template 5.51.
       idrstmpl(1)=tmplsim(1)
       idrstmpl(2)=tmplsim(2)
       idrstmpl(3)=tmplsim(3)
       idrstmpl(4)=tmplsim(4)
       idrstmpl(9)=Ts
-      idrstmpl(10)=1         ! Unpacked spectral data is 32-bit IEEE
+      idrstmpl(10)=1 ! Unpacked spectral data is 32-bit IEEE
 
       end

src/specunpack.f

@@ -1,12 +1,13 @@
 !>    @file
 !>    @brief This subroutine packs up a data field.
 !>    @author Stephen Gilbert @date 2002-12-19
-!>
 
 !>    This subroutine unpacks a spectral data field that was packed
 !>    using the complex packing algorithm for spherical harmonic data as
-!>    defined in the GRIB2 documention, using info from the GRIB2 Data
-!>    Representation Template 5.51.
+!>    defined in the GRIB2 documention, using info from the GRIB2 [Data
+!>    Representation Template
+!>    5.51](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp5-51.shtml).
+!>
 !>    @param[in] cpack The packed data field (character*1 array).
 !>    @param[in] len length of packed field cpack.
 !>    @param[in] idrstmpl Contains the array of values for Data
@@ -18,8 +19,6 @@
 !>    @param[out] fld Contains the unpacked data values.
 !>    
 !>    @author Stephen Gilbert @date 2002-12-19
-!>    
-
       subroutine specunpack(cpack,len,idrstmpl,ndpts,JJ,KK,MM,fld)
 
       character(len=1),intent(in) :: cpack(len)
@@ -43,8 +42,7 @@ subroutine specunpack(cpack,len,idrstmpl,ndpts,JJ,KK,MM,fld)
       Ts=idrstmpl(9)
 
       if (idrstmpl(10).eq.1) then           ! unpacked floats are 32-bit IEEE
-         !call g2_gbytesc(cpack,ifld,0,32,0,Ts)
-         call rdieee(cpack,unpk,Ts)          ! read IEEE unpacked floats
+          call rdieee(cpack,unpk,Ts)          ! read IEEE unpacked floats
          iofst=32*Ts
          call g2_gbytesc(cpack,ifld,iofst,nbits,0,ndpts-Ts)  ! unpack scaled data
 
@@ -61,9 +59,9 @@ subroutine specunpack(cpack,len,idrstmpl,ndpts,JJ,KK,MM,fld)
          incp=1
          do m=0,MM
             Nm=JJ      ! triangular or trapezoidal
-            if ( KK .eq. JJ+MM ) Nm=JJ+m          ! rhombodial
+            if (KK .eq. JJ+MM) Nm=JJ+m          ! rhombodial
             Ns=Js      ! triangular or trapezoidal
-            if ( Ks .eq. Js+Ms ) Ns=Js+m          ! rhombodial
+            if (Ks .eq. Js+Ms) Ns=Js+m          ! rhombodial
             do n=m,Nm
                if (n.le.Ns .AND. m.le.Ms) then    ! grab unpacked value
                   fld(inc)=unpk(incu)         ! real part

tests/CMakeLists.txt

@@ -152,6 +152,7 @@ foreach(kind ${kinds})
   create_test(test_getfield ${kind})
   create_test(test_pngpack ${kind})
   create_test(test_jpcpack ${kind})
+  create_test(test_specpack ${kind})
   create_test(test_files ${kind})
   create_test(test_gb_info ${kind})
   create_test(test_gettemplates ${kind})

tests/test_specpack.F90

@@ -0,0 +1,56 @@
+! This program tests the specpack subroutine of the NCEPLIBS-g2
+! project.
+!
+! Ed Hartnett 7/31/23
+program test_specpack
+  implicit none
+
+  integer, parameter :: width = 2, height = 2, ndpts = 4
+  real, parameter :: delta = 0.2
+  real :: fld(ndpts*2)
+  real :: fld2(ndpts*2)
+  integer :: idrstmpl(10)
+  character*1 :: cpack(200)
+  integer :: lcpack = 200
+  integer :: ii
+  integer :: jj, kk, mm
+
+
+  ! Create the fld variable with data to pack.
+  fld = (/1.1, 2.2, 3.3, 4.4, 1.1, 2.2, 3.3, 4.4/)
+  fld2 = (/0, 0, 0, 0, 0, 0, 0, 0/)
+
+  idrstmpl(1) = 0
+  idrstmpl(2) = 1
+  idrstmpl(3) = 1
+  idrstmpl(4) = 32
+  idrstmpl(5) = 1
+  idrstmpl(6) = 1
+  idrstmpl(7) = 1
+  idrstmpl(8) = 0
+  idrstmpl(9) = ndpts
+  idrstmpl(10) = 1
+
+  ! Pack the data.
+  jj = 1
+  kk = 1
+  mm = 1
+  call specpack(fld, ndpts, jj, kk, mm, idrstmpl, cpack, lcpack)
+  print *, 'lcpack: ', lcpack
+
+  ! Unpack the data.
+  call specunpack(cpack, lcpack, idrstmpl, ndpts, jj, kk, mm, fld2)
+
+  ! Compare each value to see match, remember, comparing reals
+!  print *, 'fld  ', fld
+!  print *, 'fld2 ', fld2
+  do ii = 1, ndpts
+#ifdef KIND_4
+     if (abs(fld(ii) - fld2(ii)) .gt. delta) then
+        print *, fld(ii), fld2(ii), 'do not match'
+        stop 4
+     end if
+#endif
+  end do
+  print *, 'SUCCESS!'
+end program test_specpack