adding implcit none to jpc functions, additional testing for jpcpack/jpcunpack

src/jpcpack.F90

@@ -1,8 +1,8 @@
 !> @file
-!> @brief Pack up a data field into a JPEG2000 code stream.
+!> @brief Pack a data field into a JPEG2000 code stream.
 !> @author Stephen Gilbert @date 2002-12-17
 
-!> Pack up a data field into a JPEG2000 code stream.
+!> Pack a data field into a JPEG2000 code stream.
 !>
 !> After the data field is scaled, and the reference value is
 !> subtracted out, it is treated as a grayscale image and passed to a
@@ -28,13 +28,16 @@
 !> the packed data in bytes.
 !>
 !> @author Stephen Gilbert @date 2002-12-17
-subroutine jpcpack(fld,width,height,idrstmpl,cpack,lcpack)
-
-  integer,intent(in) :: width,height
-  real,intent(in) :: fld(width*height)
+subroutine jpcpack(fld, width, height, idrstmpl, cpack, lcpack)
+  implicit none
+
+  integer,intent(in) :: width, height
+  real,intent(in) :: fld(width * height)
   character(len=1),intent(out) :: cpack(*)
   integer,intent(inout) :: idrstmpl(*)
   integer,intent(inout) :: lcpack
+  integer :: ndpts, nbytes, nbits, maxdif, j, imin, imax
+  real :: temp, dscale, bscale
 
   interface
      function enc_jpeg2000(cin, width, height, nbits, ltype, ratio, retry, outjpc, jpclen) &
@@ -48,111 +51,105 @@ function enc_jpeg2000(cin, width, height, nbits, ltype, ratio, retry, outjpc, jp
      end function enc_jpeg2000
   end interface
 
-  real(4) :: ref,rmin4
-  real(8) :: rmin,rmax
+  real(4) :: ref, rmin4
+  real(8) :: rmin, rmax
   integer(4) :: iref
   integer(8) :: nsize
-  integer :: ifld(width*height),retry
-  integer,parameter :: zero=0
-  character(len=1),allocatable :: ctemp(:)
+  integer :: ifld(width * height), retry
+  integer, parameter :: zero = 0
+  character(len = 1), allocatable :: ctemp(:)
 
-  ndpts=width*height
-  bscale=2.0**real(-idrstmpl(2))
-  dscale=10.0**real(idrstmpl(3))
+  ndpts = width * height
+  bscale = 2.0**real(-idrstmpl(2))
+  dscale = 10.0**real(idrstmpl(3))
 
-  !  Find max and min values in the data
-  if(ndpts>0) then
-     rmax=fld(1)
-     rmin=fld(1)
+  ! Find max and min values in the data.
+  if (ndpts > 0) then
+     rmax = fld(1)
+     rmin = fld(1)
   else
-     rmax=1.0
-     rmin=1.0
+     rmax = 1.0
+     rmin = 1.0
   endif
-  do j=2,ndpts
-     if (fld(j).gt.rmax) rmax=fld(j)
-     if (fld(j).lt.rmin) rmin=fld(j)
+  do j = 2, ndpts
+     if (fld(j) .gt. rmax) rmax = fld(j)
+     if (fld(j) .lt. rmin) rmin = fld(j)
   enddo
-  if (idrstmpl(2).eq.0) then
-     maxdif=nint(rmax*dscale)-nint(rmin*dscale)
+  if (idrstmpl(2) .eq. 0) then
+     maxdif = nint(rmax * dscale) - nint(rmin * dscale)
   else
-     maxdif=nint((rmax-rmin)*dscale*bscale)
+     maxdif = nint((rmax - rmin) * dscale * bscale)
   endif
 
-  !  If max and min values are not equal, pack up field.
-  !  If they are equal, we have a constant field, and the reference
-  !  value (rmin) is the value for each point in the field and
-  !  set nbits to 0.
-  if (rmin.ne.rmax .AND. maxdif.ne.0) then
-
-     !  Determine which algorithm to use based on user-supplied
-     !  binary scale factor and number of bits.
-     if (idrstmpl(2).eq.0) then
-
-        !  No binary scaling and calculate minimum number of
-        !  bits in which the data will fit.
-        imin=nint(rmin*dscale)
-        imax=nint(rmax*dscale)
-        maxdif=imax-imin
-        temp=alog(real(maxdif+1))/alog(2.0)
-        nbits=ceiling(temp)
-        rmin=real(imin)
-        !   scale data
-        do j=1,ndpts
-           ifld(j)=nint(fld(j)*dscale)-imin
+  ! If max and min values are not equal, pack up field. If they are
+  ! equal, we have a constant field, and the reference value (rmin) is
+  ! the value for each point in the field and set nbits to 0.
+  if (rmin .ne. rmax .AND. maxdif .ne. 0) then
+     ! Determine which algorithm to use based on user-supplied binary
+     ! scale factor and number of bits.
+     if (idrstmpl(2) .eq. 0) then
+        ! No binary scaling and calculate minimum number of bits in
+        ! which the data will fit.
+        imin = nint(rmin * dscale)
+        imax = nint(rmax * dscale)
+        maxdif = imax - imin
+        temp = alog(real(maxdif + 1)) / alog(2.0)
+        nbits = ceiling(temp)
+        rmin = real(imin)
+        ! Scale data.
+        do j = 1, ndpts
+           ifld(j) = nint(fld(j) * dscale) - imin
         enddo
      else
-
-        !  Use binary scaling factor and calculate minimum number of
-        !  bits in which the data will fit.
-        rmin=rmin*dscale
-        rmax=rmax*dscale
-        maxdif=nint((rmax-rmin)*bscale)
-        temp=alog(real(maxdif+1))/alog(2.0)
-        nbits=ceiling(temp)
+        ! Use binary scaling factor and calculate minimum number of
+        ! bits in which the data will fit.
+        rmin = rmin * dscale
+        rmax = rmax * dscale
+        maxdif = nint((rmax - rmin) * bscale)
+        temp = alog(real(maxdif + 1)) / alog(2.0)
+        nbits = ceiling(temp)
         !   scale data
-        do j=1,ndpts
-           ifld(j)=max(0,nint(((fld(j)*dscale)-rmin)*bscale))
+        do j = 1, ndpts
+           ifld(j) = max(0, nint(((fld(j) * dscale) - rmin) * bscale))
         enddo
      endif
 
-     !  Pack data into full octets, then do JPEG2000 encode.
-     !  and calculate the length of the packed data in bytes
-     retry=0
-     nbytes=(nbits+7)/8
-     nsize=lcpack      ! needed for input to enc_jpeg2000
-     allocate(ctemp(nbytes*ndpts))
-     call g2_sbytesc(ctemp,ifld,0,nbytes*8,0,ndpts)
-     lcpack=enc_jpeg2000(ctemp,width,height,nbits,idrstmpl(6), &
-          idrstmpl(7),retry,cpack,nsize)
-     if (lcpack.le.0) then
-        print *,'jpcpack: ERROR Packing JPC=',lcpack
-        if (lcpack.eq.-3) then
-           retry=1
+     ! Pack data into full octets, then do JPEG2000 encode.  and
+     ! calculate the length of the packed data in bytes
+     retry = 0
+     nbytes = (nbits + 7) / 8
+     nsize = lcpack      ! needed for input to enc_jpeg2000
+     allocate(ctemp(nbytes * ndpts))
+     call g2_sbytesc(ctemp, ifld, 0, nbytes * 8, 0, ndpts)
+     lcpack = enc_jpeg2000(ctemp, width, height, nbits, idrstmpl(6), &
+          idrstmpl(7), retry, cpack, nsize)
+     if (lcpack .le. 0) then
+        print *,'jpcpack: ERROR Packing JPC = ',lcpack
+        if (lcpack .eq. -3) then
+           retry = 1
            print *,'jpcpack: Retrying....'
-           lcpack=enc_jpeg2000(ctemp,width,height,nbits,idrstmpl(6), &
-                idrstmpl(7),retry,cpack,nsize)
-           if (lcpack.le.0) then
+           lcpack = enc_jpeg2000(ctemp, width, height, nbits, idrstmpl(6), &
+                idrstmpl(7), retry, cpack, nsize)
+           if (lcpack .le. 0) then
               print *,'jpcpack: Retry Failed.'
            else
               print *,'jpcpack: Retry Successful.'
            endif
         endif
      endif
      deallocate(ctemp)
-
   else
-     nbits=0
-     lcpack=0
+     nbits = 0
+     lcpack = 0
   endif
 
-  !  Fill in ref value and number of bits in Template 5.0
+  ! Fill in ref value and number of bits in Template 5.0.
   rmin4 = real(rmin, 4)
-  call mkieee(rmin4,ref,1)   ! ensure reference value is IEEE format
-  !      call g2_gbytec(ref,idrstmpl(1),0,32)
-  iref=transfer(ref,iref)
-  idrstmpl(1)=iref
-  idrstmpl(4)=nbits
-  idrstmpl(5)=0         ! original data were reals
-  if (idrstmpl(6).eq.0) idrstmpl(7)=255       ! lossy not used
+  call mkieee(rmin4, ref, 1)   ! ensure reference value is IEEE format.
+  iref = transfer(ref, iref)
+  idrstmpl(1) = iref
+  idrstmpl(4) = nbits
+  idrstmpl(5) = 0         ! original data were reals.
+  if (idrstmpl(6) .eq. 0) idrstmpl(7) = 255       ! lossy not used
 
 end subroutine

src/jpcunpack.F90

@@ -16,48 +16,50 @@
 !> @param[out] fld Contains the unpacked data values.
 !>
 !> @author Stephen Gilbert @date 2002-12-17
-      subroutine jpcunpack(cpack,len,idrstmpl,ndpts,fld)
+subroutine jpcunpack(cpack,len,idrstmpl,ndpts,fld)
+  implicit none
+
+  character(len=1),intent(in) :: cpack(len)
+  integer,intent(in) :: ndpts,len
+  integer,intent(in) :: idrstmpl(*)
+  real,intent(out) :: fld(ndpts)
 
-      character(len=1),intent(in) :: cpack(len)
-      integer,intent(in) :: ndpts,len
-      integer,intent(in) :: idrstmpl(*)
-      real,intent(out) :: fld(ndpts)
+  integer :: ifld(ndpts)
+  integer(4) :: ieee
+  integer(8) :: len8
+  real :: ref,bscale,dscale
+  integer :: nbits, j, iret
 
-      integer :: ifld(ndpts)
-      integer(4) :: ieee
-      integer(8) :: len8
-      real :: ref,bscale,dscale
+  interface
+     function dec_jpeg2000(cin, len, ifld) &
+          bind(c, name="g2c_dec_jpeg2000")
+       use iso_c_binding
+       character(kind = c_char), intent(in) :: cin(*)       
+       integer(c_size_t), value, intent(in) :: len
+       integer(c_int), intent(inout) :: ifld(*)
+       integer(c_int) :: dec_jpeg2000
+     end function dec_jpeg2000
+  end interface
 
-      interface
-         function dec_jpeg2000(cin, len, ifld) &
-              bind(c, name="g2c_dec_jpeg2000")
-           use iso_c_binding
-           character(kind = c_char), intent(in) :: cin(*)       
-           integer(c_size_t), value, intent(in) :: len
-           integer(c_int), intent(inout) :: ifld(*)
-           integer(c_int) :: dec_jpeg2000
-         end function dec_jpeg2000
-      end interface
+  ieee = idrstmpl(1)
+  call rdieee(ieee,ref,1)
+  bscale = 2.0**real(idrstmpl(2))
+  dscale = 10.0**real(-idrstmpl(3))
+  nbits = idrstmpl(4)
 
-      ieee = idrstmpl(1)
-      call rdieee(ieee,ref,1)
-      bscale = 2.0**real(idrstmpl(2))
-      dscale = 10.0**real(-idrstmpl(3))
-      nbits = idrstmpl(4)
+  !  if nbits equals 0, we have a constant field where the reference value
+  !  is the data value at each gridpoint
+  if (nbits.ne.0) then
+     !         call g2_gbytesc(cpack,ifld,0,nbits,0,ndpts)
+     len8 = len
+     iret=dec_jpeg2000(cpack,len8,ifld)
+     do j=1,ndpts
+        fld(j)=((real(ifld(j))*bscale)+ref)*dscale
+     enddo
+  else
+     do j=1,ndpts
+        fld(j)=ref
+     enddo
+  endif
 
-      !  if nbits equals 0, we have a constant field where the reference value
-      !  is the data value at each gridpoint
-      if (nbits.ne.0) then
-         !         call g2_gbytesc(cpack,ifld,0,nbits,0,ndpts)
-         len8 = len
-         iret=dec_jpeg2000(cpack,len8,ifld)
-         do j=1,ndpts
-           fld(j)=((real(ifld(j))*bscale)+ref)*dscale
-         enddo
-      else
-         do j=1,ndpts
-           fld(j)=ref
-         enddo
-      endif
-
-      end
+end subroutine jpcunpack

tests/test_jpcpack.F90

@@ -2,8 +2,8 @@
 ! project.
 !
 ! Brian Curtis 2021-10-18
+! Ed Hartnett, 7/31/23
 program test_jpcpack
-
     implicit none
 
     integer, parameter :: width = 2, height = 2, ndpts = 4
@@ -27,23 +27,64 @@ program test_jpcpack
     idrstmpl(6) = 0
     idrstmpl(7) = 1
 
+    print *, 'Testing jpcpack/jpcunpack with no data...'
+
     ! Pack the data.
-    call jpcpack(fld, width, height, idrstmpl, cpack, lcpack)
+    call jpcpack(fld, 0, height, idrstmpl, cpack, lcpack)
     print *, 'lcpack: ', lcpack
+    if (lcpack .ne. 0) stop 1
+
+    ! Unpack the data.
+    call jpcunpack(cpack, lcpack, idrstmpl, 0, fld2)
+
+    ! Compare each value to see match, remember, comparing reals.
+    print *, 'ok!'
+    print *, 'Testing jpcpack/jpcunpack with data...'
+
+    ! Pack the data.
+    lcpack = 200
+    call jpcpack(fld, width, height, idrstmpl, cpack, lcpack)
+    ! The size of lcpack will depend on the version of jasper used to compress the data.
+    !print *,lcpack
+    if (lcpack .le. 0) stop 2
 
     ! Unpack the data.
     call jpcunpack(cpack, lcpack, idrstmpl, ndpts, fld2)
 
-    ! Compare each value to see match, remember, comparing reals
-    print *, fld
-    print *, fld2
+    ! Compare each value to see match, remember, comparing reals.
+    !print *, fld
+    !print *, fld2
     do ii = 1, ndpts
         if (abs(fld(ii) - fld2(ii)) .gt. delta) then
             print *, fld(ii), fld2(ii), 'do not match'
             stop 4
         end if
     end do
 
-    print *, 'SUCCESS!'
+    print *, 'ok!'
+    print *, 'Testing jpcpack/jpcunpack with data and idrstmpl(2) = 0...'
+
+    ! Pack the data.
+    idrstmpl(2) = 0
+    lcpack = 200
+    call jpcpack(fld, width, height, idrstmpl, cpack, lcpack)
+    ! The size of lcpack will depend on the version of jasper used to compress the data.
+    !print *,lcpack
+    if (lcpack .le. 0) stop 2
+
+    ! Unpack the data.
+    call jpcunpack(cpack, lcpack, idrstmpl, ndpts, fld2)
+
+    ! Compare each value to see match, remember, comparing reals.
+    !print *, fld
+    !print *, fld2
+    do ii = 1, ndpts
+        if (abs(fld(ii) - fld2(ii)) .gt. delta) then
+            print *, fld(ii), fld2(ii), 'do not match'
+            stop 4
+        end if
+    end do
 
+    print *, 'ok!'
+    print *, 'SUCCESS!'
 end program test_jpcpack