From e347b4cef5f9e7cd651ee8aaf1f6920b9c088c0b Mon Sep 17 00:00:00 2001 From: Edward Hartnett <38856240+edwardhartnett@users.noreply.github.com> Date: Wed, 31 Jan 2024 10:31:03 -0700 Subject: [PATCH] converted some .f to .F90 and re-indented (#592) * converted realloc.f to realloc.F90 and re-indented * converting other files to F90 * converting other files to F90 * converting other files to F90 * converting other files to F90 * fixing cmake * fixing cmake * fixing cmake * fixing cmake * more F90 conversion * more F90 conversion * more F90 conversion --- src/CMakeLists.txt | 21 +- src/addfield.F90 | 455 ++++++++++++++++++++++++++ src/addfield.f | 459 -------------------------- src/{cmplxpack.f => cmplxpack.F90} | 34 +- src/compack.F90 | 447 +++++++++++++++++++++++++ src/compack.f | 446 ------------------------- src/comunpack.F90 | 319 ++++++++++++++++++ src/comunpack.f | 319 ------------------ src/getdim.F90 | 93 ++++++ src/getdim.f | 93 ------ src/getpoly.F90 | 68 ++++ src/getpoly.f | 68 ---- src/intmath.F90 | 268 +++++++++++++++ src/intmath.f | 301 ----------------- src/misspack.F90 | 492 ++++++++++++++++++++++++++++ src/misspack.f | 502 ----------------------------- src/realloc.F90 | 158 +++++++++ src/realloc.f | 158 --------- src/simpack.F90 | 154 +++++++++ src/simpack.f | 154 --------- tests/test_intmath.F90 | 34 ++ 21 files changed, 2515 insertions(+), 2528 deletions(-) create mode 100644 src/addfield.F90 delete mode 100644 src/addfield.f rename src/{cmplxpack.f => cmplxpack.F90} (69%) create mode 100644 src/compack.F90 delete mode 100644 src/compack.f create mode 100644 src/comunpack.F90 delete mode 100644 src/comunpack.f create mode 100644 src/getdim.F90 delete mode 100644 src/getdim.f create mode 100644 src/getpoly.F90 delete mode 100644 src/getpoly.f create mode 100644 src/intmath.F90 delete mode 100644 src/intmath.f create mode 100644 src/misspack.F90 delete mode 100644 src/misspack.f create mode 100644 src/realloc.F90 delete mode 100644 src/realloc.f create mode 100644 src/simpack.F90 delete mode 100644 src/simpack.f diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 3e5538e1..5a0d91cb 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -4,18 +4,19 @@ # Mark Potts, Kyle Gerheiser, Ed Hartnett # These are the fortran source files. -set(fortran_src addfield.f addgrid.F90 addlocal.F90 cmplxpack.f compack.f -comunpack.f drstemplates.F90 g2_gbytesc.F90 g2grids.F90 gb_info.F90 -getdim.f getfield.F90 getg2i.F90 getg2ir.F90 getgb2.F90 getgb2l.F90 getgb2p.F90 -getgb2r.F90 getgb2rp.F90 getgb2s.F90 getidx.F90 getlocal.F90 getpoly.f -gettemplates.F90 gf_free.F90 gf_getfld.F90 gf_unpack1.F90 gf_unpack2.F90 -gf_unpack3.F90 gf_unpack4.F90 gf_unpack5.F90 gf_unpack6.F90 gf_unpack7.F90 +set(fortran_src addfield.F90 addgrid.F90 addlocal.F90 cmplxpack.F90 +compack.F90 comunpack.F90 drstemplates.F90 g2_gbytesc.F90 g2grids.F90 +gb_info.F90 getdim.F90 getfield.F90 getg2i.F90 getg2ir.F90 getgb2.F90 +getgb2l.F90 getgb2p.F90 getgb2r.F90 getgb2rp.F90 getgb2s.F90 +getidx.F90 getlocal.F90 getpoly.F90 gettemplates.F90 gf_free.F90 +gf_getfld.F90 gf_unpack1.F90 gf_unpack2.F90 gf_unpack3.F90 +gf_unpack4.F90 gf_unpack5.F90 gf_unpack6.F90 gf_unpack7.F90 gribcreate.F90 gribend.F90 gribinfo.F90 -${CMAKE_CURRENT_BINARY_DIR}/gribmod.F90 gridtemplates.F90 intmath.f -ixgb2.F90 jpcpack.F90 jpcunpack.F90 misspack.f mkieee.F90 pack_gp.f +${CMAKE_CURRENT_BINARY_DIR}/gribmod.F90 gridtemplates.F90 intmath.F90 +ixgb2.F90 jpcpack.F90 jpcunpack.F90 misspack.F90 mkieee.F90 pack_gp.f params_ecmwf.F90 params.F90 pdstemplates.F90 pngpack.F90 pngunpack.F90 -putgb2.F90 rdieee.F90 realloc.f reduce.f simpack.f simunpack.F90 skgb.F90 -specpack.F90 specunpack.F90) +putgb2.F90 rdieee.F90 realloc.F90 reduce.f simpack.F90 simunpack.F90 +skgb.F90 specpack.F90 specunpack.F90) # This function calls NCEPLIBS-w3emc. if (BUILD_WITH_W3EMC) diff --git a/src/addfield.F90 b/src/addfield.F90 new file mode 100644 index 00000000..6bebf1f0 --- /dev/null +++ b/src/addfield.F90 @@ -0,0 +1,455 @@ +!> @file +!> @brief Pack up Sections 4 through 7 for a field and add them +!> to a GRIB2 message. +!> @author Stephen Gilbert @date 2000-05-02 + +!> Pack up Sections 4 through 7 for a field and add them to a +!> GRIB2 message. +!> +!> They are the Product Definition Section, Data Representation +!> Section, Bit-Map Section and Data Sections. +!> +!> This routine is used with routines gribcreate(), addlocal(), +!> addgrid(), and gribend() to create a complete GRIB2 +!> message. Subroutine gribcreate() must be called first to initialize +!> a new GRIB2 message. Subroutine addgrid() must be called after +!> gribcreate() and before this routine to add the appropriate grid +!> description to the GRIB2 message. A call to gribend() is required +!> to complete GRIB2 message after all fields have been added. +!> +!> @param[inout] cgrib Character array to contain the GRIB2 message. +!> @param[in] lcgrib Maximum length (bytes) of array cgrib. +!> @param[in] ipdsnum Product Definition Template Number (see [Code +!> Table 4.0] +!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table4-0.shtml)). +!> @param[in] ipdstmpl Contains the data values for the +!> Product Definition Template specified by ipdsnum. +!> @param[in] ipdstmplen Max dimension of ipdstmpl. +!> @param[out] coordlist Array containg floating point values intended to +!> document the vertical discretisation associated to model data on hybrid +!> coordinate vertical levels (part of Section 4). +!> @param[in] numcoord - number of values in array coordlist. +!> @param[in] idrsnum - Data Representation Template Number (see +!> [Code Table 5.0] +!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table5-0.shtml)). +!> @param[in] idrstmpl Contains the data values for the Data +!> Representation Template specified by idrsnum. Note that some +!> values in this template (eg. reference values, number of bits, +!> etc...) may be changed by the data packing algorithms. Use this +!> to specify scaling factors and order of spatial differencing, if +!> desired. +!> @param[in] idrstmplen Max dimension of idrstmpl. This must be at +!> least as large as the length of the selected PDS template. +!> @param[in] fld Array of data points to pack. +!> @param[out] ngrdpts Number of data points in grid. i.e. size of +!> fld and bmap. +!> @param[out] ibmap Bitmap indicator (see [Code Table +!> 6.0](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table6-0.shtml)). +!> - 0 bitmap applies and is included in Section 6. +!> - 1-253 Predefined bitmap applies +!> - 254 Previously defined bitmap applies to this field +!> - 255 Bit map does not apply to this product. +!> @param[out] bmap Logical*1 array containing bitmap to be added +!> (if ibmap=0 or ibmap=254). +!> @param[out] ierr Error return code. +!> - 0 no error. +!> - 1 GRIB message was not initialized. Need to call +!> routine gribcreate first. +!> - 2 GRIB message already complete. Cannot add new section. +!> - 3 Sum of Section byte counts does not add to total +!> byte count. +!> - 4 Previous Section was not 3 or 7. +!> - 5 Could not find requested Product Definition Template. +!> - 6 Section 3 (GDS) not previously defined in message. +!> - 7 Tried to use unsupported Data Representationi Template. +!> - 8 Specified use of a previously defined bitmap, but one +!> does not exist in the GRIB message. +!> - 9 GDT of one of 5.50 through 5.53 required to pack +!> using DRT 5.51. +!> - 10 Error packing data field. +!> +!> @author Stephen Gilbert @date 2000-05-02 +subroutine addfield(cgrib,lcgrib,ipdsnum,ipdstmpl,ipdstmplen, & + coordlist,numcoord,idrsnum,idrstmpl, & + idrstmplen,fld,ngrdpts,ibmap,bmap,ierr) + use pdstemplates + use drstemplates + implicit none + + logical :: match + character(len=1),intent(inout) :: cgrib(lcgrib) + integer,intent(in) :: ipdsnum,ipdstmpl(*) + integer,intent(in) :: idrsnum,numcoord,ipdstmplen,idrstmplen + integer,intent(in) :: lcgrib,ngrdpts,ibmap + real,intent(in) :: coordlist(numcoord) + real(kind = 4) :: coordlist_4(numcoord) + real,target,intent(in) :: fld(ngrdpts) + integer,intent(out) :: ierr + integer,intent(inout) :: idrstmpl(*) + logical*1,intent(in) :: bmap(ngrdpts) + + character(len=4),parameter :: grib='GRIB',c7777='7777' + character(len=4):: ctemp + character(len=1),allocatable :: cpack(:) + real,pointer,dimension(:) :: pfld + real(4) :: coordieee(numcoord),re00 + integer(4) :: ire00,allones + integer :: mappds(ipdstmplen),intbmap(ngrdpts),mapdrs(idrstmplen) + integer,parameter :: zero=0,one=1,four=4,five=5,six=6,seven=7 + integer,parameter :: minsize=50000 + integer iofst,ibeg,lencurr,len,mappdslen,mapdrslen,lpos3 + integer width,height,ndpts + integer lensec3,lensec4,lensec5,lensec6,lensec7 + logical issec3,needext,isprevbmap + integer :: nbits, newlen, nsize, lcpack, left + integer :: ibmprev, ilen, ioctet, iscan, isecnum, itemp + integer :: i, jj, kk, mm + integer :: iret, istat + + allones = int(Z'FFFFFFFF') + ierr=0 + + ! Check to see if beginning of GRIB message exists + match=.true. + do i=1,4 + if(cgrib(i) /= grib(i:i)) then + match=.false. + endif + enddo + if (.not. match) then + print *,'addfield: GRIB not found in given message.' + print *,'addfield: Call to routine gribcreate required to initialize GRIB messge.' + ierr=1 + return + endif + + ! Get current length of GRIB message + call g2_gbytec(cgrib,lencurr,96,32) + + ! Check to see if GRIB message is already complete + ctemp=cgrib(lencurr-3)//cgrib(lencurr-2)//cgrib(lencurr-1) //cgrib(lencurr) + if (ctemp.eq.c7777) then + print *,'addfield: GRIB message already complete. Cannot add new section.' + ierr=2 + return + endif + + ! Loop through all current sections of the GRIB message to + ! find the last section number. + issec3=.false. + isprevbmap=.false. + len=16 ! length of Section 0 + do + ! Get number and length of next section + iofst=len*8 + call g2_gbytec(cgrib,ilen,iofst,32) + iofst=iofst+32 + call g2_gbytec(cgrib,isecnum,iofst,8) + iofst=iofst+8 + ! Check if previous Section 3 exists and save location of + ! the section 3 in case needed later. + if (isecnum.eq.3) then + issec3=.true. + lpos3=len+1 + lensec3=ilen + endif + ! Check if a previous defined bitmap exists + if (isecnum.eq.6) then + call g2_gbytec(cgrib,ibmprev,iofst,8) + iofst=iofst+8 + if ((ibmprev.ge.0).and.(ibmprev.le.253)) isprevbmap=.true. + endif + len=len+ilen + ! Exit loop if last section reached + if (len.eq.lencurr) exit + ! If byte count for each section does not match current + ! total length, then there is a problem. + if (len.gt.lencurr) then + print *,'addfield: Section byte counts don''t add to total.' + print *,'addfield: Sum of section byte counts = ',len + print *,'addfield: Total byte count in Section 0 = ',lencurr + ierr=3 + return + endif + enddo + + ! Sections 4 through 7 can only be added after section 3 or 7. + if ((isecnum.ne.3) .and. (isecnum.ne.7)) then + print *,'addfield: Sections 4-7 can only be added after', & + ' Section 3 or 7.' + print *,'addfield: Section ',isecnum,' was the last found in', & + ' given GRIB message.' + ierr=4 + return + + ! Sections 4 through 7 can only be added if section 3 was previously defined. + elseif (.not.issec3) then + print *,'addfield: Sections 4-7 can only be added if Section', & + ' 3 was previously included.' + print *,'addfield: Section 3 was not found in given GRIB message.' + print *,'addfield: Call to routine addgrid required', & + ' to specify Grid definition.' + ierr=6 + return + endif + + ! Add Section 4 - Product Definition Section + ibeg=lencurr*8 ! Calculate offset for beginning of section 4 + iofst=ibeg+32 ! leave space for length of section + call g2_sbytec(cgrib,four,iofst,8) ! Store section number (4) + iofst=iofst+8 + call g2_sbytec(cgrib,numcoord,iofst,16) ! Store num of coordinate values + iofst=iofst+16 + call g2_sbytec(cgrib,ipdsnum,iofst,16) ! Store Prod Def Template num. + iofst=iofst+16 + + ! Get Product Definition Template + call getpdstemplate(ipdsnum,mappdslen,mappds,needext,iret) + if (iret.ne.0) then + ierr=5 + return + endif + + ! Extend the Product Definition Template, if necessary. + ! The number of values in a specific template may vary + ! depending on data specified in the "static" part of the + ! template. + if (needext) then + call extpdstemplate(ipdsnum,ipdstmpl,mappdslen,mappds) + endif + + ! Pack up each input value in array ipdstmpl into the + ! the appropriate number of octets, which are specified in + ! corresponding entries in array mappds. + do i=1,mappdslen + nbits=iabs(mappds(i))*8 + if ((mappds(i).ge.0).or.(ipdstmpl(i).ge.0)) then + call g2_sbytec(cgrib,ipdstmpl(i),iofst,nbits) + else + call g2_sbytec(cgrib,one,iofst,1) + call g2_sbytec(cgrib,iabs(ipdstmpl(i)),iofst+1,nbits-1) + endif + iofst=iofst+nbits + enddo + + ! Add Optional list of vertical coordinate values + ! after the Product Definition Template, if necessary. + if (numcoord .ne. 0) then + do i = 1, numcoord + coordlist_4(i) = real(coordlist(i), 4) + end do + call mkieee(coordlist_4, coordieee, numcoord) + call g2_sbytesc(cgrib,coordieee,iofst,32,0,numcoord) + iofst=iofst+(32*numcoord) + endif + + ! Calculate length of section 4 and store it in octets + ! 1-4 of section 4. + lensec4=(iofst-ibeg)/8 + call g2_sbytec(cgrib,lensec4,ibeg,32) + + ! Pack Data using appropriate algorithm + + ! Get Data Representation Template + call getdrstemplate(idrsnum,mapdrslen,mapdrs,needext,iret) + if (iret.ne.0) then + ierr=5 + return + endif + + ! contract data field, removing data at invalid grid points, + ! if bit-map is provided with field. + if (ibmap.eq.0 .OR. ibmap.eq.254) then + allocate(pfld(max(2,ngrdpts))) + ndpts=0; + do jj=1,ngrdpts + intbmap(jj)=0 + if (bmap(jj)) then + intbmap(jj)=1 + ndpts=ndpts+1 + pfld(ndpts)=fld(jj); + endif + enddo + if(ndpts==0 .and. ngrdpts>0) then + pfld(1)=0 + endif + else + ndpts=ngrdpts; + pfld=>fld; + endif + lcpack=0 + nsize=ndpts*4 + if (nsize .lt. minsize) nsize=minsize + allocate(cpack(nsize),stat=istat) + if (idrsnum.eq.0) then ! Simple Packing + call simpack(pfld,ndpts,idrstmpl,cpack,lcpack) + elseif (idrsnum.eq.2.or.idrsnum.eq.3) then ! Complex Packing + call cmplxpack(pfld,ndpts,idrsnum,idrstmpl,cpack,lcpack) + elseif (idrsnum.eq.50) then ! Sperical Harmonic Simple Packing + call simpack(pfld(2),ndpts-1,idrstmpl,cpack,lcpack) + call mkieee(real(pfld(1)),re00,1) ! ensure RE(0,0) value is IEEE format + !call g2_gbytec(re00,idrstmpl(5),0,32) + ire00=transfer(re00,ire00) + idrstmpl(5)=ire00 + elseif (idrsnum.eq.51) then ! Sperical Harmonic Complex Packing + call getpoly(cgrib(lpos3),lensec3,jj,kk,mm) + if (jj.ne.0 .AND. kk.ne.0 .AND. mm.ne.0) then + call specpack(pfld,ndpts,jj,kk,mm,idrstmpl,cpack,lcpack) + else + print *,'addfield: Cannot pack DRT 5.51.' + ierr=9 + return + endif + + elseif (idrsnum.eq.40 .OR. idrsnum.eq.40000) then ! JPEG2000 encoding + if (ibmap.eq.255) then + call getdim(cgrib(lpos3),lensec3,width,height,iscan) + if (width.eq.0 .OR. height.eq.0) then + width=ndpts + height=1 + elseif (width.eq.allones .OR. height.eq.allones) then + width=ndpts + height=1 + elseif (ibits(iscan,5,1) .eq. 1) then ! Scanning mode: bit 3 + itemp=width + width=height + height=itemp + endif + else + width=ndpts + height=1 + endif + if(width<1 .or. height<1) then + ! Special case: bitmask off everywhere. + write(0,*) 'Warning: bitmask off everywhere.' + write(0,*) ' Pretend one point in jpcpack to avoid crash.' + width=1 + height=1 + endif + lcpack=nsize + !print *,'w,h=',width,height + call jpcpack(pfld,width,height,idrstmpl,cpack,lcpack) + + elseif (idrsnum.eq.41 .OR. idrsnum.eq.40010) then ! PNG encoding + if (ibmap.eq.255) then + call getdim(cgrib(lpos3),lensec3,width,height,iscan) + if (width.eq.0 .OR. height.eq.0) then + width=ndpts + height=1 + elseif (width.eq.allones .OR. height.eq.allones) then + width=ndpts + height=1 + elseif (ibits(iscan,5,1) .eq. 1) then ! Scanning mode: bit 3 + itemp=width + width=height + height=itemp + endif + else + width=ndpts + height=1 + endif + !print *,'png size ',width,height + call pngpack(pfld,width,height,idrstmpl,cpack,lcpack) + !print *,'png packed' + else + print *,'addfield: Data Representation Template 5.',idrsnum, & + ' not yet implemented.' + ierr=7 + return + endif + if (ibmap.eq.0 .OR. ibmap.eq.254) then + deallocate(pfld) + endif + if (lcpack .lt. 0) then + if(allocated(cpack))deallocate(cpack) + ierr=10 + return + endif + + ! Add Section 5 - Data Representation Section + ibeg=iofst ! Calculate offset for beginning of section 5 + iofst=ibeg+32 ! leave space for length of section + call g2_sbytec(cgrib,five,iofst,8) ! Store section number (5) + iofst=iofst+8 + call g2_sbytec(cgrib,ndpts,iofst,32) ! Store num of actual data points + iofst=iofst+32 + call g2_sbytec(cgrib,idrsnum,iofst,16) ! Store Data Repr. Template num. + iofst=iofst+16 + + ! Pack up each input value in array idrstmpl into the + ! the appropriate number of octets, which are specified in + ! corresponding entries in array mapdrs. + do i=1,mapdrslen + nbits=iabs(mapdrs(i))*8 + if ((mapdrs(i).ge.0).or.(idrstmpl(i).ge.0)) then + call g2_sbytec(cgrib,idrstmpl(i),iofst,nbits) + else + call g2_sbytec(cgrib,one,iofst,1) + call g2_sbytec(cgrib,iabs(idrstmpl(i)),iofst+1,nbits-1) + endif + iofst=iofst+nbits + enddo + + ! Calculate length of section 5 and store it in octets + ! 1-4 of section 5. + lensec5=(iofst-ibeg)/8 + call g2_sbytec(cgrib,lensec5,ibeg,32) + + ! Add Section 6 - Bit-Map Section + ibeg=iofst ! Calculate offset for beginning of section 6 + iofst=ibeg+32 ! leave space for length of section + call g2_sbytec(cgrib,six,iofst,8) ! Store section number (6) + iofst=iofst+8 + call g2_sbytec(cgrib,ibmap,iofst,8) ! Store Bit Map indicator + iofst=iofst+8 + + ! Store bitmap, if supplied + if (ibmap.eq.0) then + call g2_sbytesc(cgrib,intbmap,iofst,1,0,ngrdpts) ! Store BitMap + iofst=iofst+ngrdpts + endif + + ! If specifying a previously defined bit-map, make sure + ! one already exists in the current GRIB message. + if ((ibmap.eq.254).and.(.not.isprevbmap)) then + print *,'addfield: Requested previously defined bitmap, ', & + ' but one does not exist in the current GRIB message.' + ierr=8 + return + endif + + ! Calculate length of section 6 and store it in octets + ! 1-4 of section 6. Pad to end of octect, if necessary. + left=8-mod(iofst,8) + if (left.ne.8) then + call g2_sbytec(cgrib,zero,iofst,left) ! Pad with zeros to fill Octet + iofst=iofst+left + endif + lensec6=(iofst-ibeg)/8 + call g2_sbytec(cgrib,lensec6,ibeg,32) + + ! Add Section 7 - Data Section + ibeg=iofst ! Calculate offset for beginning of section 7 + iofst=ibeg+32 ! leave space for length of section + call g2_sbytec(cgrib,seven,iofst,8) ! Store section number (7) + iofst=iofst+8 + ! Store Packed Binary Data values, if non-constant field + if (lcpack.ne.0) then + ioctet=iofst/8 + cgrib(ioctet+1:ioctet+lcpack)=cpack(1:lcpack) + iofst=iofst+(8*lcpack) + endif + + ! Calculate length of section 7 and store it in octets + ! 1-4 of section 7. + lensec7=(iofst-ibeg)/8 + call g2_sbytec(cgrib,lensec7,ibeg,32) + + if(allocated(cpack) )deallocate(cpack) + + ! Update current byte total of message in Section 0 + newlen=lencurr+lensec4+lensec5+lensec6+lensec7 + call g2_sbytec(cgrib,newlen,96,32) + + return +end subroutine addfield diff --git a/src/addfield.f b/src/addfield.f deleted file mode 100644 index bb5d6632..00000000 --- a/src/addfield.f +++ /dev/null @@ -1,459 +0,0 @@ -!> @file -!> @brief Pack up Sections 4 through 7 for a field and add them -!> to a GRIB2 message. -!> @author Stephen Gilbert @date 2000-05-02 - -!> Pack up Sections 4 through 7 for a field and add them to a -!> GRIB2 message. -!> -!> They are the Product Definition Section, Data Representation -!> Section, Bit-Map Section and Data Sections. -!> -!> This routine is used with routines gribcreate(), addlocal(), -!> addgrid(), and gribend() to create a complete GRIB2 -!> message. Subroutine gribcreate() must be called first to initialize -!> a new GRIB2 message. Subroutine addgrid() must be called after -!> gribcreate() and before this routine to add the appropriate grid -!> description to the GRIB2 message. A call to gribend() is required -!> to complete GRIB2 message after all fields have been added. -!> -!> @param[inout] cgrib Character array to contain the GRIB2 message. -!> @param[in] lcgrib Maximum length (bytes) of array cgrib. -!> @param[in] ipdsnum Product Definition Template Number (see [Code -!> Table 4.0] -!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table4-0.shtml)). -!> @param[in] ipdstmpl Contains the data values for the -!> Product Definition Template specified by ipdsnum. -!> @param[in] ipdstmplen Max dimension of ipdstmpl. -!> @param[out] coordlist Array containg floating point values intended to -!> document the vertical discretisation associated to model data on hybrid -!> coordinate vertical levels (part of Section 4). -!> @param[in] numcoord - number of values in array coordlist. -!> @param[in] idrsnum - Data Representation Template Number (see -!> [Code Table 5.0] -!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table5-0.shtml)). -!> @param[in] idrstmpl Contains the data values for the Data -!> Representation Template specified by idrsnum. Note that some -!> values in this template (eg. reference values, number of bits, -!> etc...) may be changed by the data packing algorithms. Use this -!> to specify scaling factors and order of spatial differencing, if -!> desired. -!> @param[in] idrstmplen Max dimension of idrstmpl. This must be at -!> least as large as the length of the selected PDS template. -!> @param[in] fld Array of data points to pack. -!> @param[out] ngrdpts Number of data points in grid. i.e. size of -!> fld and bmap. -!> @param[out] ibmap Bitmap indicator (see [Code Table -!> 6.0](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table6-0.shtml)). -!> - 0 bitmap applies and is included in Section 6. -!> - 1-253 Predefined bitmap applies -!> - 254 Previously defined bitmap applies to this field -!> - 255 Bit map does not apply to this product. -!> @param[out] bmap Logical*1 array containing bitmap to be added -!> (if ibmap=0 or ibmap=254). -!> @param[out] ierr Error return code. -!> - 0 no error. -!> - 1 GRIB message was not initialized. Need to call -!> routine gribcreate first. -!> - 2 GRIB message already complete. Cannot add new section. -!> - 3 Sum of Section byte counts does not add to total -!> byte count. -!> - 4 Previous Section was not 3 or 7. -!> - 5 Could not find requested Product Definition Template. -!> - 6 Section 3 (GDS) not previously defined in message. -!> - 7 Tried to use unsupported Data Representationi Template. -!> - 8 Specified use of a previously defined bitmap, but one -!> does not exist in the GRIB message. -!> - 9 GDT of one of 5.50 through 5.53 required to pack -!> using DRT 5.51. -!> - 10 Error packing data field. -!> -!> @author Stephen Gilbert @date 2000-05-02 - subroutine addfield(cgrib,lcgrib,ipdsnum,ipdstmpl,ipdstmplen, - & coordlist,numcoord,idrsnum,idrstmpl, - & idrstmplen,fld,ngrdpts,ibmap,bmap,ierr) - use pdstemplates - use drstemplates - implicit none - - logical :: match - character(len=1),intent(inout) :: cgrib(lcgrib) - integer,intent(in) :: ipdsnum,ipdstmpl(*) - integer,intent(in) :: idrsnum,numcoord,ipdstmplen,idrstmplen - integer,intent(in) :: lcgrib,ngrdpts,ibmap - real,intent(in) :: coordlist(numcoord) - real(kind = 4) :: coordlist_4(numcoord) - real,target,intent(in) :: fld(ngrdpts) - integer,intent(out) :: ierr - integer,intent(inout) :: idrstmpl(*) - logical*1,intent(in) :: bmap(ngrdpts) - - character(len=4),parameter :: grib='GRIB',c7777='7777' - character(len=4):: ctemp - character(len=1),allocatable :: cpack(:) - real,pointer,dimension(:) :: pfld - real(4) :: coordieee(numcoord),re00 - integer(4) :: ire00,allones - integer :: mappds(ipdstmplen),intbmap(ngrdpts),mapdrs(idrstmplen) - integer,parameter :: zero=0,one=1,four=4,five=5,six=6,seven=7 - integer,parameter :: minsize=50000 - integer iofst,ibeg,lencurr,len,mappdslen,mapdrslen,lpos3 - integer width,height,ndpts - integer lensec3,lensec4,lensec5,lensec6,lensec7 - logical issec3,needext,isprevbmap - integer :: nbits, newlen, nsize, lcpack, left - integer :: ibmprev, ilen, ioctet, iscan, isecnum, itemp - integer :: i, jj, kk, mm - integer :: iret, istat - - allones = int(Z'FFFFFFFF') - ierr=0 - -! Check to see if beginning of GRIB message exists - match=.true. - do i=1,4 - if(cgrib(i) /= grib(i:i)) then - match=.false. - endif - enddo - if (.not. match) then - print *,'addfield: GRIB not found in given message.' - print *,'addfield: Call to routine gribcreate required', - & ' to initialize GRIB messge.' - ierr=1 - return - endif - -! Get current length of GRIB message - call g2_gbytec(cgrib,lencurr,96,32) - -! Check to see if GRIB message is already complete - ctemp=cgrib(lencurr-3)//cgrib(lencurr-2)//cgrib(lencurr-1) - & //cgrib(lencurr) - if (ctemp.eq.c7777) then - print *,'addfield: GRIB message already complete. Cannot', - & ' add new section.' - ierr=2 - return - endif - -! Loop through all current sections of the GRIB message to -! find the last section number. - issec3=.false. - isprevbmap=.false. - len=16 ! length of Section 0 - do - ! Get number and length of next section - iofst=len*8 - call g2_gbytec(cgrib,ilen,iofst,32) - iofst=iofst+32 - call g2_gbytec(cgrib,isecnum,iofst,8) - iofst=iofst+8 - ! Check if previous Section 3 exists and save location of - ! the section 3 in case needed later. - if (isecnum.eq.3) then - issec3=.true. - lpos3=len+1 - lensec3=ilen - endif - ! Check if a previous defined bitmap exists - if (isecnum.eq.6) then - call g2_gbytec(cgrib,ibmprev,iofst,8) - iofst=iofst+8 - if ((ibmprev.ge.0).and.(ibmprev.le.253)) isprevbmap=.true. - endif - len=len+ilen - ! Exit loop if last section reached - if (len.eq.lencurr) exit - ! If byte count for each section does not match current - ! total length, then there is a problem. - if (len.gt.lencurr) then - print *,'addfield: Section byte counts don''t add to total.' - print *,'addfield: Sum of section byte counts = ',len - print *,'addfield: Total byte count in Section 0 = ',lencurr - ierr=3 - return - endif - enddo - -! Sections 4 through 7 can only be added after section 3 or 7. - if ((isecnum.ne.3) .and. (isecnum.ne.7)) then - print *,'addfield: Sections 4-7 can only be added after', - & ' Section 3 or 7.' - print *,'addfield: Section ',isecnum,' was the last found in', - & ' given GRIB message.' - ierr=4 - return - -! Sections 4 through 7 can only be added if section 3 was previously defined. - elseif (.not.issec3) then - print *,'addfield: Sections 4-7 can only be added if Section', - & ' 3 was previously included.' - print *,'addfield: Section 3 was not found in', - & ' given GRIB message.' - print *,'addfield: Call to routine addgrid required', - & ' to specify Grid definition.' - ierr=6 - return - endif - -! Add Section 4 - Product Definition Section - ibeg=lencurr*8 ! Calculate offset for beginning of section 4 - iofst=ibeg+32 ! leave space for length of section - call g2_sbytec(cgrib,four,iofst,8) ! Store section number (4) - iofst=iofst+8 - call g2_sbytec(cgrib,numcoord,iofst,16) ! Store num of coordinate values - iofst=iofst+16 - call g2_sbytec(cgrib,ipdsnum,iofst,16) ! Store Prod Def Template num. - iofst=iofst+16 - - ! Get Product Definition Template - call getpdstemplate(ipdsnum,mappdslen,mappds,needext,iret) - if (iret.ne.0) then - ierr=5 - return - endif - - ! Extend the Product Definition Template, if necessary. - ! The number of values in a specific template may vary - ! depending on data specified in the "static" part of the - ! template. - if (needext) then - call extpdstemplate(ipdsnum,ipdstmpl,mappdslen,mappds) - endif - - ! Pack up each input value in array ipdstmpl into the - ! the appropriate number of octets, which are specified in - ! corresponding entries in array mappds. - do i=1,mappdslen - nbits=iabs(mappds(i))*8 - if ((mappds(i).ge.0).or.(ipdstmpl(i).ge.0)) then - call g2_sbytec(cgrib,ipdstmpl(i),iofst,nbits) - else - call g2_sbytec(cgrib,one,iofst,1) - call g2_sbytec(cgrib,iabs(ipdstmpl(i)),iofst+1,nbits-1) - endif - iofst=iofst+nbits - enddo - - ! Add Optional list of vertical coordinate values - ! after the Product Definition Template, if necessary. - if (numcoord .ne. 0) then - do i = 1, numcoord - coordlist_4(i) = real(coordlist(i), 4) - end do - call mkieee(coordlist_4, coordieee, numcoord) - call g2_sbytesc(cgrib,coordieee,iofst,32,0,numcoord) - iofst=iofst+(32*numcoord) - endif - - ! Calculate length of section 4 and store it in octets - ! 1-4 of section 4. - lensec4=(iofst-ibeg)/8 - call g2_sbytec(cgrib,lensec4,ibeg,32) - -! Pack Data using appropriate algorithm - - ! Get Data Representation Template - call getdrstemplate(idrsnum,mapdrslen,mapdrs,needext,iret) - if (iret.ne.0) then - ierr=5 - return - endif - - ! contract data field, removing data at invalid grid points, - ! if bit-map is provided with field. - if (ibmap.eq.0 .OR. ibmap.eq.254) then - allocate(pfld(max(2,ngrdpts))) - ndpts=0; - do jj=1,ngrdpts - intbmap(jj)=0 - if (bmap(jj)) then - intbmap(jj)=1 - ndpts=ndpts+1 - pfld(ndpts)=fld(jj); - endif - enddo - if(ndpts==0 .and. ngrdpts>0) then - pfld(1)=0 - endif - else - ndpts=ngrdpts; - pfld=>fld; - endif - lcpack=0 - nsize=ndpts*4 - if (nsize .lt. minsize) nsize=minsize - allocate(cpack(nsize),stat=istat) - if (idrsnum.eq.0) then ! Simple Packing - call simpack(pfld,ndpts,idrstmpl,cpack,lcpack) - elseif (idrsnum.eq.2.or.idrsnum.eq.3) then ! Complex Packing - call cmplxpack(pfld,ndpts,idrsnum,idrstmpl,cpack,lcpack) - elseif (idrsnum.eq.50) then ! Sperical Harmonic Simple Packing - call simpack(pfld(2),ndpts-1,idrstmpl,cpack,lcpack) - call mkieee(real(pfld(1)),re00,1) ! ensure RE(0,0) value is IEEE format - !call g2_gbytec(re00,idrstmpl(5),0,32) - ire00=transfer(re00,ire00) - idrstmpl(5)=ire00 - elseif (idrsnum.eq.51) then ! Sperical Harmonic Complex Packing - call getpoly(cgrib(lpos3),lensec3,jj,kk,mm) - if (jj.ne.0 .AND. kk.ne.0 .AND. mm.ne.0) then - call specpack(pfld,ndpts,jj,kk,mm,idrstmpl,cpack,lcpack) - else - print *,'addfield: Cannot pack DRT 5.51.' - ierr=9 - return - endif - - elseif (idrsnum.eq.40 .OR. idrsnum.eq.40000) then ! JPEG2000 encoding - if (ibmap.eq.255) then - call getdim(cgrib(lpos3),lensec3,width,height,iscan) - if (width.eq.0 .OR. height.eq.0) then - width=ndpts - height=1 - elseif (width.eq.allones .OR. height.eq.allones) then - width=ndpts - height=1 - elseif (ibits(iscan,5,1) .eq. 1) then ! Scanning mode: bit 3 - itemp=width - width=height - height=itemp - endif - else - width=ndpts - height=1 - endif - if(width<1 .or. height<1) then - ! Special case: bitmask off everywhere. - write(0,*) 'Warning: bitmask off everywhere.' - write(0,*) ' Pretend one point in jpcpack to avoid crash.' - width=1 - height=1 - endif - lcpack=nsize - !print *,'w,h=',width,height - call jpcpack(pfld,width,height,idrstmpl,cpack,lcpack) - - elseif (idrsnum.eq.41 .OR. idrsnum.eq.40010) then ! PNG encoding - if (ibmap.eq.255) then - call getdim(cgrib(lpos3),lensec3,width,height,iscan) - if (width.eq.0 .OR. height.eq.0) then - width=ndpts - height=1 - elseif (width.eq.allones .OR. height.eq.allones) then - width=ndpts - height=1 - elseif (ibits(iscan,5,1) .eq. 1) then ! Scanning mode: bit 3 - itemp=width - width=height - height=itemp - endif - else - width=ndpts - height=1 - endif - !print *,'png size ',width,height - call pngpack(pfld,width,height,idrstmpl,cpack,lcpack) - !print *,'png packed' - else - print *,'addfield: Data Representation Template 5.',idrsnum, - * ' not yet implemented.' - ierr=7 - return - endif - if (ibmap.eq.0 .OR. ibmap.eq.254) then - deallocate(pfld) - endif - if (lcpack .lt. 0) then - if(allocated(cpack))deallocate(cpack) - ierr=10 - return - endif - -! Add Section 5 - Data Representation Section - ibeg=iofst ! Calculate offset for beginning of section 5 - iofst=ibeg+32 ! leave space for length of section - call g2_sbytec(cgrib,five,iofst,8) ! Store section number (5) - iofst=iofst+8 - call g2_sbytec(cgrib,ndpts,iofst,32) ! Store num of actual data points - iofst=iofst+32 - call g2_sbytec(cgrib,idrsnum,iofst,16) ! Store Data Repr. Template num. - iofst=iofst+16 - - ! Pack up each input value in array idrstmpl into the - ! the appropriate number of octets, which are specified in - ! corresponding entries in array mapdrs. - do i=1,mapdrslen - nbits=iabs(mapdrs(i))*8 - if ((mapdrs(i).ge.0).or.(idrstmpl(i).ge.0)) then - call g2_sbytec(cgrib,idrstmpl(i),iofst,nbits) - else - call g2_sbytec(cgrib,one,iofst,1) - call g2_sbytec(cgrib,iabs(idrstmpl(i)),iofst+1,nbits-1) - endif - iofst=iofst+nbits - enddo - - ! Calculate length of section 5 and store it in octets - ! 1-4 of section 5. - lensec5=(iofst-ibeg)/8 - call g2_sbytec(cgrib,lensec5,ibeg,32) - -! Add Section 6 - Bit-Map Section - ibeg=iofst ! Calculate offset for beginning of section 6 - iofst=ibeg+32 ! leave space for length of section - call g2_sbytec(cgrib,six,iofst,8) ! Store section number (6) - iofst=iofst+8 - call g2_sbytec(cgrib,ibmap,iofst,8) ! Store Bit Map indicator - iofst=iofst+8 - - ! Store bitmap, if supplied - if (ibmap.eq.0) then - call g2_sbytesc(cgrib,intbmap,iofst,1,0,ngrdpts) ! Store BitMap - iofst=iofst+ngrdpts - endif - - ! If specifying a previously defined bit-map, make sure - ! one already exists in the current GRIB message. - if ((ibmap.eq.254).and.(.not.isprevbmap)) then - print *,'addfield: Requested previously defined bitmap, ', - & ' but one does not exist in the current GRIB message.' - ierr=8 - return - endif - - ! Calculate length of section 6 and store it in octets - ! 1-4 of section 6. Pad to end of octect, if necessary. - left=8-mod(iofst,8) - if (left.ne.8) then - call g2_sbytec(cgrib,zero,iofst,left) ! Pad with zeros to fill Octet - iofst=iofst+left - endif - lensec6=(iofst-ibeg)/8 - call g2_sbytec(cgrib,lensec6,ibeg,32) - -! Add Section 7 - Data Section - ibeg=iofst ! Calculate offset for beginning of section 7 - iofst=ibeg+32 ! leave space for length of section - call g2_sbytec(cgrib,seven,iofst,8) ! Store section number (7) - iofst=iofst+8 - ! Store Packed Binary Data values, if non-constant field - if (lcpack.ne.0) then - ioctet=iofst/8 - cgrib(ioctet+1:ioctet+lcpack)=cpack(1:lcpack) - iofst=iofst+(8*lcpack) - endif - - ! Calculate length of section 7 and store it in octets - ! 1-4 of section 7. - lensec7=(iofst-ibeg)/8 - call g2_sbytec(cgrib,lensec7,ibeg,32) - - if(allocated(cpack) )deallocate(cpack) - -! Update current byte total of message in Section 0 - newlen=lencurr+lensec4+lensec5+lensec6+lensec7 - call g2_sbytec(cgrib,newlen,96,32) - - return - end diff --git a/src/cmplxpack.f b/src/cmplxpack.F90 similarity index 69% rename from src/cmplxpack.f rename to src/cmplxpack.F90 index abacb991..9d25f9a4 100644 --- a/src/cmplxpack.f +++ b/src/cmplxpack.F90 @@ -34,24 +34,22 @@ !> idrstmpl(7) is not set correctly. !> !> @author Stephen Gilbert @date 2004-08-27 - subroutine cmplxpack(fld,ndpts,idrsnum,idrstmpl,cpack,lcpack) +subroutine cmplxpack(fld,ndpts,idrsnum,idrstmpl,cpack,lcpack) - integer,intent(in) :: ndpts,idrsnum - real,intent(in) :: fld(ndpts) - character(len=1),intent(out) :: cpack(*) - integer,intent(inout) :: idrstmpl(*) - integer,intent(out) :: lcpack + integer,intent(in) :: ndpts,idrsnum + real,intent(in) :: fld(ndpts) + character(len=1),intent(out) :: cpack(*) + integer,intent(inout) :: idrstmpl(*) + integer,intent(out) :: lcpack - + if ( idrstmpl(7) .eq. 0 ) then ! No internal missing values + call compack(fld,ndpts,idrsnum,idrstmpl,cpack,lcpack) + elseif ( idrstmpl(7).eq.1 .OR. idrstmpl(7).eq.2) then + call misspack(fld,ndpts,idrsnum,idrstmpl,cpack,lcpack) + else + print *,'cmplxpack: Do not recognize Missing value option.' + lcpack=-1 + endif - if ( idrstmpl(7) .eq. 0 ) then ! No internal missing values - call compack(fld,ndpts,idrsnum,idrstmpl,cpack,lcpack) - elseif ( idrstmpl(7).eq.1 .OR. idrstmpl(7).eq.2) then - call misspack(fld,ndpts,idrsnum,idrstmpl,cpack,lcpack) - else - print *,'cmplxpack: Do not recognize Missing value option.' - lcpack=-1 - endif - - return - end + return +end subroutine cmplxpack diff --git a/src/compack.F90 b/src/compack.F90 new file mode 100644 index 00000000..7934befa --- /dev/null +++ b/src/compack.F90 @@ -0,0 +1,447 @@ +!> @file +!> @brief Pack a data field with complex packing with or without +!> spatial differences defined in [Data Representation +!> Template5.2](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp5-2.shtml) +!> and [Data Representation Template +!> 5.3](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp5-3.shtml). +!> @author Stephen Gilbert @date 2000-06-21 + +!> Pack a data field with complex packing with or without spatial +!> differences, defined in [Data Representation Template +!> 5.2](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp5-2.shtml) +!> and [Data Representation Template +!> 5.3](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp5-3.shtml). +!> +!> It also fills in Data Representation Template 5.2 or 5.3 with +!> the appropriate values. +!> +!> @param[in] fld The data values to pack. +!> @param[in] ndpts The number of data values in array fld. +!> @param[in] idrsnum Data Representation Template number - must equal +!> 2 or 3. +!> @param[inout] idrstmpl The array of values for Data +!> Representation Template 5.2 or 5.3. +!> - 1 Reference value - ignored on input +!> - 2 Binary Scale Factor +!> - 3 Decimal Scale Factor +!> - 7 Missing value management +!> - 8 Primary missing value +!> - 9 Secondary missing value +!> - 17 Order of Spatial Differencing (1 or 2) +!> @param[out] cpack The packed data field (character*1 array). +!> @param[out] lcpack length of packed field cpack. +!> +!> @author Stephen Gilbert @date 2000-06-21 +subroutine compack(fld,ndpts,idrsnum,idrstmpl,cpack,lcpack) + + use intmath + implicit none + integer,intent(in) :: ndpts,idrsnum + real,intent(in) :: fld(ndpts) + character(len=1),intent(out) :: cpack(*) + integer,intent(inout) :: idrstmpl(*) + integer,intent(out) :: lcpack + + real :: bscale,dscale + integer :: j,iofst,imin,ival1,ival2,minsd,nbitsd,n + integer :: igmax,nbitsgref,left,iwmax,i,ilmax,kk,ij + integer :: ngwidthref,nbitsgwidth,nglenref,nglenlast + integer :: maxorig,nbitorig,isd,ngroups,itemp,minpk + integer :: kfildo,inc,maxgrps,missopt,miss1,miss2,lg + integer :: ibit,jbit,kbit,novref,lbitref,ier,ng,imax + integer :: nbitsglen + real(4) :: ref,rmin4 + real(8) :: rmin,rmax + + integer(4) :: iref + integer,allocatable :: ifld(:) + integer,allocatable :: jmin(:),jmax(:),lbit(:) + integer,parameter :: zero=0 + integer,allocatable :: gref(:),gwidth(:),glen(:) + integer :: glength,grpwidth + logical :: simple_alg + + simple_alg = .false. + + 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) + else + 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) + enddo + ! + ! 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. + ! + multival: if (rmin.ne.rmax) then + iofst=0 + allocate(ifld(ndpts)) + allocate(gref(ndpts)) + allocate(gwidth(ndpts)) + allocate(glen(ndpts)) + ! + ! Scale original data + ! + if (idrstmpl(2).eq.0) then ! No binary scaling + imin=nint(rmin*dscale) + !imax=nint(rmax*dscale) + rmin=real(imin) + do j=1,ndpts + ifld(j)=max(0,nint(fld(j)*dscale)-imin) + enddo + else ! Use binary scaling factor + rmin=rmin*dscale + !rmax=rmax*dscale + do j=1,ndpts + ifld(j)=max(0,nint(((fld(j)*dscale)-rmin)*bscale)) + enddo + endif + ! + ! Calculate Spatial differences, if using DRS Template 5.3 + ! + alg3: if (idrsnum.eq.3) then ! spatial differences + if (idrstmpl(17).ne.1.and.idrstmpl(17).ne.2) idrstmpl(17)=2 + if (idrstmpl(17).eq.1) then ! first order + ival1=ifld(1) + if(ival1<0) then + print *,'G2: negative ival1',ival1 + stop 101 + endif + do j=ndpts,2,-1 + ifld(j)=ifld(j)-ifld(j-1) + enddo + ifld(1)=0 + elseif (idrstmpl(17).eq.2) then ! second order + ival1=ifld(1) + ival2=ifld(2) + if(ival1<0) then + print *,'G2: negative ival1',ival1 + stop 11 + endif + if(ival2<0) then + print *,'G2: negative ival2',ival2 + stop 12 + endif + do j=ndpts,3,-1 + ifld(j)=ifld(j)-(2*ifld(j-1))+ifld(j-2) + enddo + ifld(1)=0 + ifld(2)=0 + endif + ! + ! subtract min value from spatial diff field + ! + isd=idrstmpl(17)+1 + minsd=minval(ifld(isd:ndpts)) + do j=isd,ndpts + ifld(j)=ifld(j)-minsd + enddo + ! + ! find num of bits need to store minsd and add 1 extra bit + ! to indicate sign + ! + nbitsd=i1log2(abs(minsd))+1 + ! + ! find num of bits need to store ifld(1) ( and ifld(2) + ! if using 2nd order differencing ) + ! + maxorig=ival1 + if (idrstmpl(17).eq.2.and.ival2.gt.ival1) maxorig=ival2 + nbitorig=i1log2(maxorig)+1 + if (nbitorig.gt.nbitsd) nbitsd=nbitorig + ! increase number of bits to even multiple of 8 ( octet ) + if (mod(nbitsd,8).ne.0) nbitsd=nbitsd+(8-mod(nbitsd,8)) + ! + ! Store extra spatial differencing info into the packed + ! data section. + ! + if (nbitsd.ne.0) then + ! pack first original value + if (ival1.ge.0) then + call g2_sbytec(cpack,ival1,iofst,nbitsd) + iofst=iofst+nbitsd + else + call g2_sbytec(cpack,1,iofst,1) + iofst=iofst+1 + call g2_sbytec(cpack,iabs(ival1),iofst,nbitsd-1) + iofst=iofst+nbitsd-1 + endif + if (idrstmpl(17).eq.2) then + ! pack second original value + if (ival2.ge.0) then + call g2_sbytec(cpack,ival2,iofst,nbitsd) + iofst=iofst+nbitsd + else + call g2_sbytec(cpack,1,iofst,1) + iofst=iofst+1 + call g2_sbytec(cpack,iabs(ival2),iofst,nbitsd-1) + iofst=iofst+nbitsd-1 + endif + endif + ! pack overall min of spatial differences + if (minsd.ge.0) then + call g2_sbytec(cpack,minsd,iofst,nbitsd) + iofst=iofst+nbitsd + else + call g2_sbytec(cpack,1,iofst,1) + iofst=iofst+1 + call g2_sbytec(cpack,iabs(minsd),iofst,nbitsd-1) + iofst=iofst+nbitsd-1 + endif + endif + !print *,'SDp ',ival1,ival2,minsd,nbitsd + endif alg3 ! end of spatial diff section + ! + ! Determine Groups to be used. + ! + simplealg: if ( simple_alg ) then + ! set group length to 10 : calculate number of groups + ! and length of last group + print *,'G2: use simple algorithm' + ngroups=ndpts/10 + glen(1:ngroups)=10 + itemp=mod(ndpts,10) + if (itemp.ne.0) then + ngroups=ngroups+1 + glen(ngroups)=itemp + endif + else + ! Use Dr. Glahn's algorithm for determining grouping. + ! + kfildo=6 + minpk=10 + inc=1 + maxgrps=((ndpts+minpk-1)/minpk) + allocate(jmin(maxgrps)) + allocate(jmax(maxgrps)) + allocate(lbit(maxgrps)) + missopt=0 + call pack_gp(kfildo,ifld,ndpts,missopt,minpk,inc,miss1,miss2, & + jmin,jmax,lbit,glen,maxgrps,ngroups,ibit,jbit, & + kbit,novref,lbitref,ier) + if(ier/=0) then + ! Dr. Glahn's algorithm failed; use simple packing method instead. +1099 format('G2: fall back to simple algorithm (glahn ier=',I0,')') + print 1099,ier + ngroups=ndpts/10 + glen(1:ngroups)=10 + itemp=mod(ndpts,10) + if (itemp.ne.0) then + ngroups=ngroups+1 + glen(ngroups)=itemp + endif + elseif(ngroups<1) then + ! Dr. Glahn's algorithm failed; use simple packing method instead. + print *,'Glahn algorithm failed; use simple packing' + ngroups=ndpts/10 + glen(1:ngroups)=10 + itemp=mod(ndpts,10) + if (itemp.ne.0) then + ngroups=ngroups+1 + glen(ngroups)=itemp + endif + else + !print *,'SAGier = ',ier,ibit,jbit,kbit,novref,lbitref + do ng=1,ngroups + glen(ng)=glen(ng)+novref + enddo + deallocate(jmin) + deallocate(jmax) + deallocate(lbit) + endif + endif simplealg + ! + ! For each group, find the group's reference value + ! and the number of bits needed to hold the remaining values + ! + n=1 + do ng=1,ngroups + ! find max and min values of group + gref(ng)=ifld(n) + imax=ifld(n) + j=n+1 + do lg=2,glen(ng) + if (ifld(j).lt.gref(ng)) gref(ng)=ifld(j) + if (ifld(j).gt.imax) imax=ifld(j) + j=j+1 + enddo + ! calc num of bits needed to hold data + if ( gref(ng).ne.imax ) then + gwidth(ng)=i1log2(imax-gref(ng)) + else + gwidth(ng)=0 + endif + ! Subtract min from data + j=n + do lg=1,glen(ng) + ifld(j)=ifld(j)-gref(ng) + j=j+1 + enddo + ! increment fld array counter + n=n+glen(ng) + enddo + ! + ! Find max of the group references and calc num of bits needed + ! to pack each groups reference value, then + ! pack up group reference values + ! + !write(77,*)'GREFS: ',(gref(j),j=1,ngroups) + igmax=maxval(gref(1:ngroups)) + if (igmax.ne.0) then + nbitsgref=i1log2(igmax) + call g2_sbytesc(cpack,gref,iofst,nbitsgref,0,ngroups) + itemp=nbitsgref*ngroups + iofst=iofst+itemp + if (mod(itemp,8).ne.0) then + left=8-mod(itemp,8) + call g2_sbytec(cpack,zero,iofst,left) + iofst=iofst+left + endif + else + nbitsgref=0 + endif + ! + ! Find max/min of the group widths and calc num of bits needed + ! to pack each groups width value, then + ! pack up group width values + ! + !write(77,*)'GWIDTHS: ',(gwidth(j),j=1,ngroups) + iwmax=maxval(gwidth(1:ngroups)) + ngwidthref=minval(gwidth(1:ngroups)) + if (iwmax.ne.ngwidthref) then + nbitsgwidth=i1log2(iwmax-ngwidthref) + do i=1,ngroups + gwidth(i)=gwidth(i)-ngwidthref + if(gwidth(i)<0) then + write(0,*) 'i,gw,ngw=',i,gwidth(i),ngwidthref + stop 9 + endif + enddo + call g2_sbytesc(cpack,gwidth,iofst,nbitsgwidth,0,ngroups) + itemp=nbitsgwidth*ngroups + iofst=iofst+itemp + ! Pad last octet with Zeros, if necessary, + if (mod(itemp,8).ne.0) then + left=8-mod(itemp,8) + call g2_sbytec(cpack,zero,iofst,left) + iofst=iofst+left + endif + else + nbitsgwidth=0 + gwidth(1:ngroups)=0 + endif + ! + ! Find max/min of the group lengths and calc num of bits needed + ! to pack each groups length value, then + ! pack up group length values + ! + !write(77,*)'GLENS: ',(glen(j),j=1,ngroups) + ilmax=maxval(glen(1:ngroups-1)) + nglenref=minval(glen(1:ngroups-1)) + nglenlast=glen(ngroups) + if (ilmax.ne.nglenref) then + nbitsglen=i1log2(ilmax-nglenref) + do i=1,ngroups-1 + glen(i)=glen(i)-nglenref + if(glen(i)<0) then + write(0,*) 'i,glen(i) = ',i,glen(i) + stop 23 + endif + enddo + call g2_sbytesc(cpack,glen,iofst,nbitsglen,0,ngroups) + itemp=nbitsglen*ngroups + iofst=iofst+itemp + ! Pad last octet with Zeros, if necessary, + if (mod(itemp,8).ne.0) then + left=8-mod(itemp,8) + call g2_sbytec(cpack,zero,iofst,left) + iofst=iofst+left + endif + else + nbitsglen=0 + glen(1:ngroups)=0 + endif + ! + ! For each group, pack data values + ! + !write(77,*)'IFLDS: ',(ifld(j),j=1,ndpts) + n=1 + ij=0 + do ng=1,ngroups + glength=glen(ng)+nglenref + if (ng.eq.ngroups ) glength=nglenlast + grpwidth=gwidth(ng)+ngwidthref + !write(77,*)'NGP ',ng,grpwidth,glength,gref(ng) + if ( grpwidth.ne.0 ) then + call g2_sbytesc(cpack,ifld(n),iofst,grpwidth,0,glength) + iofst=iofst+(grpwidth*glength) + endif + do kk=1,glength + ij=ij+1 + !write(77,*)'SAG ',ij,fld(ij),ifld(ij),gref(ng),bscale,rmin,dscale + enddo + n=n+glength + enddo + ! Pad last octet with Zeros, if necessary, + if (mod(iofst,8).ne.0) then + left=8-mod(iofst,8) + call g2_sbytec(cpack,zero,iofst,left) + iofst=iofst+left + endif + lcpack=iofst/8 + ! + if ( allocated(ifld) ) deallocate(ifld) + if ( allocated(gref) ) deallocate(gref) + if ( allocated(gwidth) ) deallocate(gwidth) + if ( allocated(glen) ) deallocate(glen) + else ! Constant field ( max = min ) + lcpack=0 + nbitsgref=0 + ngroups=0 + ngwidthref=0 + nbitsgwidth=0 + nglenref=0 + nglenlast=0 + nbitsglen=0 + nbitsd=0 + endif multival + + ! + ! Fill in ref value and number of bits in Template 5.2 + ! + 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)=nbitsgref + idrstmpl(5)=0 ! original data were reals + idrstmpl(6)=1 ! general group splitting + idrstmpl(7)=0 ! No internal missing values + idrstmpl(8)=0 ! Primary missing value + idrstmpl(9)=0 ! secondary missing value + idrstmpl(10)=ngroups ! Number of groups + idrstmpl(11)=ngwidthref ! reference for group widths + idrstmpl(12)=nbitsgwidth ! num bits used for group widths + idrstmpl(13)=nglenref ! Reference for group lengths + idrstmpl(14)=1 ! length increment for group lengths + idrstmpl(15)=nglenlast ! True length of last group + idrstmpl(16)=nbitsglen ! num bits used for group lengths + if (idrsnum.eq.3) then + idrstmpl(18)=nbitsd/8 ! num bits used for extra spatial + ! differencing values + endif + +end subroutine compack diff --git a/src/compack.f b/src/compack.f deleted file mode 100644 index 4a10e103..00000000 --- a/src/compack.f +++ /dev/null @@ -1,446 +0,0 @@ -!> @file -!> @brief Pack a data field with complex packing with -!> or without spatial differences defined in [Data Representation -!> Template5.2](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp5-2.shtml) -!> and [Data Representation Template 5.3](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp5-3.shtml). -!> @author Stephen Gilbert @date 2000-06-21 - -!> Pack a data field with complex packing with or -!> without spatial differences, defined in [Data Representation -!> Template 5.2](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp5-2.shtml) -!> and [Data Representation Template 5.3](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_temp5-3.shtml). -!> -!> It also fills in Data Representation Template 5.2 or 5.3 with -!> the appropriate values. -!> -!> @param[in] fld The data values to pack. -!> @param[in] ndpts The number of data values in array fld. -!> @param[in] idrsnum Data Representation Template number - must equal -!> 2 or 3. -!> @param[inout] idrstmpl The array of values for Data -!> Representation Template 5.2 or 5.3. -!> - 1 Reference value - ignored on input -!> - 2 Binary Scale Factor -!> - 3 Decimal Scale Factor -!> - 7 Missing value management -!> - 8 Primary missing value -!> - 9 Secondary missing value -!> - 17 Order of Spatial Differencing (1 or 2) -!> @param[out] cpack The packed data field (character*1 array). -!> @param[out] lcpack length of packed field cpack. -!> -!> @author Stephen Gilbert @date 2000-06-21 - subroutine compack(fld,ndpts,idrsnum,idrstmpl,cpack,lcpack) - - use intmath - implicit none - integer,intent(in) :: ndpts,idrsnum - real,intent(in) :: fld(ndpts) - character(len=1),intent(out) :: cpack(*) - integer,intent(inout) :: idrstmpl(*) - integer,intent(out) :: lcpack - - real :: bscale,dscale - integer :: j,iofst,imin,ival1,ival2,minsd,nbitsd,n - integer :: igmax,nbitsgref,left,iwmax,i,ilmax,kk,ij - integer :: ngwidthref,nbitsgwidth,nglenref,nglenlast - integer :: maxorig,nbitorig,isd,ngroups,itemp,minpk - integer :: kfildo,inc,maxgrps,missopt,miss1,miss2,lg - integer :: ibit,jbit,kbit,novref,lbitref,ier,ng,imax - integer :: nbitsglen - real(4) :: ref,rmin4 - real(8) :: rmin,rmax - - integer(4) :: iref - integer,allocatable :: ifld(:) - integer,allocatable :: jmin(:),jmax(:),lbit(:) - integer,parameter :: zero=0 - integer,allocatable :: gref(:),gwidth(:),glen(:) - integer :: glength,grpwidth - logical :: simple_alg - - simple_alg = .false. - - 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) - else - 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) - enddo -! -! 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. -! - multival: if (rmin.ne.rmax) then - iofst=0 - allocate(ifld(ndpts)) - allocate(gref(ndpts)) - allocate(gwidth(ndpts)) - allocate(glen(ndpts)) - ! - ! Scale original data - ! - if (idrstmpl(2).eq.0) then ! No binary scaling - imin=nint(rmin*dscale) - !imax=nint(rmax*dscale) - rmin=real(imin) - do j=1,ndpts - ifld(j)=max(0,nint(fld(j)*dscale)-imin) - enddo - else ! Use binary scaling factor - rmin=rmin*dscale - !rmax=rmax*dscale - do j=1,ndpts - ifld(j)=max(0,nint(((fld(j)*dscale)-rmin)*bscale)) - enddo - endif - ! - ! Calculate Spatial differences, if using DRS Template 5.3 - ! - alg3: if (idrsnum.eq.3) then ! spatial differences - if (idrstmpl(17).ne.1.and.idrstmpl(17).ne.2) idrstmpl(17)=2 - if (idrstmpl(17).eq.1) then ! first order - ival1=ifld(1) - if(ival1<0) then - print *,'G2: negative ival1',ival1 - stop 101 - endif - do j=ndpts,2,-1 - ifld(j)=ifld(j)-ifld(j-1) - enddo - ifld(1)=0 - elseif (idrstmpl(17).eq.2) then ! second order - ival1=ifld(1) - ival2=ifld(2) - if(ival1<0) then - print *,'G2: negative ival1',ival1 - stop 11 - endif - if(ival2<0) then - print *,'G2: negative ival2',ival2 - stop 12 - endif - do j=ndpts,3,-1 - ifld(j)=ifld(j)-(2*ifld(j-1))+ifld(j-2) - enddo - ifld(1)=0 - ifld(2)=0 - endif - ! - ! subtract min value from spatial diff field - ! - isd=idrstmpl(17)+1 - minsd=minval(ifld(isd:ndpts)) - do j=isd,ndpts - ifld(j)=ifld(j)-minsd - enddo - ! - ! find num of bits need to store minsd and add 1 extra bit - ! to indicate sign - ! - nbitsd=i1log2(abs(minsd))+1 - ! - ! find num of bits need to store ifld(1) ( and ifld(2) - ! if using 2nd order differencing ) - ! - maxorig=ival1 - if (idrstmpl(17).eq.2.and.ival2.gt.ival1) maxorig=ival2 - nbitorig=i1log2(maxorig)+1 - if (nbitorig.gt.nbitsd) nbitsd=nbitorig - ! increase number of bits to even multiple of 8 ( octet ) - if (mod(nbitsd,8).ne.0) nbitsd=nbitsd+(8-mod(nbitsd,8)) - ! - ! Store extra spatial differencing info into the packed - ! data section. - ! - if (nbitsd.ne.0) then - ! pack first original value - if (ival1.ge.0) then - call g2_sbytec(cpack,ival1,iofst,nbitsd) - iofst=iofst+nbitsd - else - call g2_sbytec(cpack,1,iofst,1) - iofst=iofst+1 - call g2_sbytec(cpack,iabs(ival1),iofst,nbitsd-1) - iofst=iofst+nbitsd-1 - endif - if (idrstmpl(17).eq.2) then - ! pack second original value - if (ival2.ge.0) then - call g2_sbytec(cpack,ival2,iofst,nbitsd) - iofst=iofst+nbitsd - else - call g2_sbytec(cpack,1,iofst,1) - iofst=iofst+1 - call g2_sbytec(cpack,iabs(ival2),iofst,nbitsd-1) - iofst=iofst+nbitsd-1 - endif - endif - ! pack overall min of spatial differences - if (minsd.ge.0) then - call g2_sbytec(cpack,minsd,iofst,nbitsd) - iofst=iofst+nbitsd - else - call g2_sbytec(cpack,1,iofst,1) - iofst=iofst+1 - call g2_sbytec(cpack,iabs(minsd),iofst,nbitsd-1) - iofst=iofst+nbitsd-1 - endif - endif - !print *,'SDp ',ival1,ival2,minsd,nbitsd - endif alg3 ! end of spatial diff section - ! - ! Determine Groups to be used. - ! - simplealg: if ( simple_alg ) then - ! set group length to 10 : calculate number of groups - ! and length of last group - print *,'G2: use simple algorithm' - ngroups=ndpts/10 - glen(1:ngroups)=10 - itemp=mod(ndpts,10) - if (itemp.ne.0) then - ngroups=ngroups+1 - glen(ngroups)=itemp - endif - else - ! Use Dr. Glahn's algorithm for determining grouping. - ! - kfildo=6 - minpk=10 - inc=1 - maxgrps=((ndpts+minpk-1)/minpk) - allocate(jmin(maxgrps)) - allocate(jmax(maxgrps)) - allocate(lbit(maxgrps)) - missopt=0 - call pack_gp(kfildo,ifld,ndpts,missopt,minpk,inc,miss1,miss2, - & jmin,jmax,lbit,glen,maxgrps,ngroups,ibit,jbit, - & kbit,novref,lbitref,ier) - if(ier/=0) then - ! Dr. Glahn's algorithm failed; use simple packing method instead. - 1099 format('G2: fall back to simple algorithm (glahn ier=',I0,& - & ')') - print 1099,ier - ngroups=ndpts/10 - glen(1:ngroups)=10 - itemp=mod(ndpts,10) - if (itemp.ne.0) then - ngroups=ngroups+1 - glen(ngroups)=itemp - endif - elseif(ngroups<1) then - ! Dr. Glahn's algorithm failed; use simple packing method instead. - print *,'Glahn algorithm failed; use simple packing' - ngroups=ndpts/10 - glen(1:ngroups)=10 - itemp=mod(ndpts,10) - if (itemp.ne.0) then - ngroups=ngroups+1 - glen(ngroups)=itemp - endif - else -!print *,'SAGier = ',ier,ibit,jbit,kbit,novref,lbitref - do ng=1,ngroups - glen(ng)=glen(ng)+novref - enddo - deallocate(jmin) - deallocate(jmax) - deallocate(lbit) - endif - endif simplealg - ! - ! For each group, find the group's reference value - ! and the number of bits needed to hold the remaining values - ! - n=1 - do ng=1,ngroups - ! find max and min values of group - gref(ng)=ifld(n) - imax=ifld(n) - j=n+1 - do lg=2,glen(ng) - if (ifld(j).lt.gref(ng)) gref(ng)=ifld(j) - if (ifld(j).gt.imax) imax=ifld(j) - j=j+1 - enddo - ! calc num of bits needed to hold data - if ( gref(ng).ne.imax ) then - gwidth(ng)=i1log2(imax-gref(ng)) - else - gwidth(ng)=0 - endif - ! Subtract min from data - j=n - do lg=1,glen(ng) - ifld(j)=ifld(j)-gref(ng) - j=j+1 - enddo - ! increment fld array counter - n=n+glen(ng) - enddo - ! - ! Find max of the group references and calc num of bits needed - ! to pack each groups reference value, then - ! pack up group reference values - ! - !write(77,*)'GREFS: ',(gref(j),j=1,ngroups) - igmax=maxval(gref(1:ngroups)) - if (igmax.ne.0) then - nbitsgref=i1log2(igmax) - call g2_sbytesc(cpack,gref,iofst,nbitsgref,0,ngroups) - itemp=nbitsgref*ngroups - iofst=iofst+itemp - if (mod(itemp,8).ne.0) then - left=8-mod(itemp,8) - call g2_sbytec(cpack,zero,iofst,left) - iofst=iofst+left - endif - else - nbitsgref=0 - endif - ! - ! Find max/min of the group widths and calc num of bits needed - ! to pack each groups width value, then - ! pack up group width values - ! - !write(77,*)'GWIDTHS: ',(gwidth(j),j=1,ngroups) - iwmax=maxval(gwidth(1:ngroups)) - ngwidthref=minval(gwidth(1:ngroups)) - if (iwmax.ne.ngwidthref) then - nbitsgwidth=i1log2(iwmax-ngwidthref) - do i=1,ngroups - gwidth(i)=gwidth(i)-ngwidthref - if(gwidth(i)<0) then - write(0,*) 'i,gw,ngw=',i,gwidth(i),ngwidthref - stop 9 - endif - enddo - call g2_sbytesc(cpack,gwidth,iofst,nbitsgwidth,0,ngroups) - itemp=nbitsgwidth*ngroups - iofst=iofst+itemp - ! Pad last octet with Zeros, if necessary, - if (mod(itemp,8).ne.0) then - left=8-mod(itemp,8) - call g2_sbytec(cpack,zero,iofst,left) - iofst=iofst+left - endif - else - nbitsgwidth=0 - gwidth(1:ngroups)=0 - endif - ! - ! Find max/min of the group lengths and calc num of bits needed - ! to pack each groups length value, then - ! pack up group length values - ! - !write(77,*)'GLENS: ',(glen(j),j=1,ngroups) - ilmax=maxval(glen(1:ngroups-1)) - nglenref=minval(glen(1:ngroups-1)) - nglenlast=glen(ngroups) - if (ilmax.ne.nglenref) then - nbitsglen=i1log2(ilmax-nglenref) - do i=1,ngroups-1 - glen(i)=glen(i)-nglenref - if(glen(i)<0) then - write(0,*) 'i,glen(i) = ',i,glen(i) - stop 23 - endif - enddo - call g2_sbytesc(cpack,glen,iofst,nbitsglen,0,ngroups) - itemp=nbitsglen*ngroups - iofst=iofst+itemp - ! Pad last octet with Zeros, if necessary, - if (mod(itemp,8).ne.0) then - left=8-mod(itemp,8) - call g2_sbytec(cpack,zero,iofst,left) - iofst=iofst+left - endif - else - nbitsglen=0 - glen(1:ngroups)=0 - endif - ! - ! For each group, pack data values - ! - !write(77,*)'IFLDS: ',(ifld(j),j=1,ndpts) - n=1 - ij=0 - do ng=1,ngroups - glength=glen(ng)+nglenref - if (ng.eq.ngroups ) glength=nglenlast - grpwidth=gwidth(ng)+ngwidthref - !write(77,*)'NGP ',ng,grpwidth,glength,gref(ng) - if ( grpwidth.ne.0 ) then - call g2_sbytesc(cpack,ifld(n),iofst,grpwidth,0,glength) - iofst=iofst+(grpwidth*glength) - endif - do kk=1,glength - ij=ij+1 - !write(77,*)'SAG ',ij,fld(ij),ifld(ij),gref(ng),bscale,rmin,dscale - enddo - n=n+glength - enddo - ! Pad last octet with Zeros, if necessary, - if (mod(iofst,8).ne.0) then - left=8-mod(iofst,8) - call g2_sbytec(cpack,zero,iofst,left) - iofst=iofst+left - endif - lcpack=iofst/8 - ! - if ( allocated(ifld) ) deallocate(ifld) - if ( allocated(gref) ) deallocate(gref) - if ( allocated(gwidth) ) deallocate(gwidth) - if ( allocated(glen) ) deallocate(glen) - else ! Constant field ( max = min ) - lcpack=0 - nbitsgref=0 - ngroups=0 - ngwidthref=0 - nbitsgwidth=0 - nglenref=0 - nglenlast=0 - nbitsglen=0 - nbitsd=0 - endif multival - -! -! Fill in ref value and number of bits in Template 5.2 -! - 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)=nbitsgref - idrstmpl(5)=0 ! original data were reals - idrstmpl(6)=1 ! general group splitting - idrstmpl(7)=0 ! No internal missing values - idrstmpl(8)=0 ! Primary missing value - idrstmpl(9)=0 ! secondary missing value - idrstmpl(10)=ngroups ! Number of groups - idrstmpl(11)=ngwidthref ! reference for group widths - idrstmpl(12)=nbitsgwidth ! num bits used for group widths - idrstmpl(13)=nglenref ! Reference for group lengths - idrstmpl(14)=1 ! length increment for group lengths - idrstmpl(15)=nglenlast ! True length of last group - idrstmpl(16)=nbitsglen ! num bits used for group lengths - if (idrsnum.eq.3) then - idrstmpl(18)=nbitsd/8 ! num bits used for extra spatial - ! differencing values - endif - - end diff --git a/src/comunpack.F90 b/src/comunpack.F90 new file mode 100644 index 00000000..2a3e457d --- /dev/null +++ b/src/comunpack.F90 @@ -0,0 +1,319 @@ +!> @file +!> @brief Unpack a data field that was packed using a complex packing +!> algorithm as defined in the GRIB2 documention. +!> @author Stephen Gilbert @date 2000-06-21 + +!> Unpack a data field that was packed using a complex packing +!> algorithm as defined in the GRIB2 documention. +!> +!> This subroutine Supports GRIB2 complex packing templates with or +!> without spatial differences: It supports GRIB2 complex packing +!> templates with or without spatial differences: Data Representation +!> Tables [5.2] +!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table5-2.shtml) +!> and +!> [5.3](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table5-3.shtml)). +!> +!> @param[in] cpack The packed data field (character*1 array). +!> @param[in] len Length of packed field cpack. +!> @param[in] lensec Length of section 7 (used for error checking). +!> @param[in] idrsnum Data Representation Template number 5.N. Must +!> equal 2 or 3. +!> @param[in] idrstmpl The array of values for Data Representation +!> Template 5.2 or 5.3. +!> @param[in] ndpts The number of data values to unpack. +!> @param[out] fld Contains the unpacked data values. +!> @param[out] ier Error return: +!> - 0 = No error. +!> - 1 = Problem - inconsistent group lengths of widths. +!> +!> @author Stephen Gilbert @date 2000-06-21 +subroutine comunpack(cpack,len,lensec,idrsnum,idrstmpl,ndpts, & + fld,ier) + + character(len=1),intent(in) :: cpack(len) + integer,intent(in) :: ndpts,len + integer,intent(in) :: idrstmpl(*) + real,intent(out) :: fld(ndpts) + + integer,allocatable :: ifld(:),ifldmiss(:) + integer(4) :: ieee + integer,allocatable :: gref(:),gwidth(:),glen(:) + real :: ref,bscale,dscale,rmiss1,rmiss2 + ! real :: fldo(6045) + integer :: totBit, totLen + + ier=0 + !print *,'IDRSTMPL: ',(idrstmpl(j),j=1,16) + ieee = idrstmpl(1) + call rdieee(ieee,ref,1) + bscale = 2.0**real(idrstmpl(2)) + dscale = 10.0**real(-idrstmpl(3)) + nbitsgref = idrstmpl(4) + itype = idrstmpl(5) + ngroups = idrstmpl(10) + nbitsgwidth = idrstmpl(12) + nbitsglen = idrstmpl(16) + if (idrsnum.eq.3) then + nbitsd=idrstmpl(18)*8 + endif + + ! Constant field + + if (ngroups.eq.0) then + do j=1,ndpts + fld(j)=ref + enddo + return + endif + + iofst=0 + allocate(ifld(ndpts),stat=is) + !print *,'ALLOC ifld: ',is,ndpts + allocate(gref(ngroups),stat=is) + !print *,'ALLOC gref: ',is + allocate(gwidth(ngroups),stat=is) + !print *,'ALLOC gwidth: ',is + ! + ! Get missing values, if supplied + ! + if ( idrstmpl(7).eq.1 ) then + if (itype.eq.0) then + call rdieee(idrstmpl(8),rmiss1,1) + else + rmiss1=real(idrstmpl(8)) + endif + elseif ( idrstmpl(7).eq.2 ) then + if (itype.eq.0) then + call rdieee(idrstmpl(8),rmiss1,1) + call rdieee(idrstmpl(9),rmiss2,1) + else + rmiss1=real(idrstmpl(8)) + rmiss2=real(idrstmpl(9)) + endif + endif + !print *,'RMISSs: ',rmiss1,rmiss2,ref + ! + ! Extract Spatial differencing values, if using DRS Template 5.3 + ! + if (idrsnum.eq.3) then + if (nbitsd.ne.0) then + call g2_gbytec(cpack,ival1,iofst,nbitsd) + iofst=iofst+nbitsd + if (idrstmpl(17).eq.2) then + call g2_gbytec(cpack,ival2,iofst,nbitsd) + iofst=iofst+nbitsd + endif + call g2_gbytec(cpack,isign,iofst,1) + iofst=iofst+1 + call g2_gbytec(cpack,minsd,iofst,nbitsd-1) + iofst=iofst+nbitsd-1 + if (isign.eq.1) minsd=-minsd + else + ival1=0 + ival2=0 + minsd=0 + endif + !print *,'SDu ',ival1,ival2,minsd,nbitsd + endif + ! + ! Extract Each Group's reference value + ! + !print *,'SAG1: ',nbitsgref,ngroups,iofst + if (nbitsgref.ne.0) then + call g2_gbytesc(cpack,gref,iofst,nbitsgref,0,ngroups) + itemp=nbitsgref*ngroups + iofst=iofst+(itemp) + if (mod(itemp,8).ne.0) iofst=iofst+(8-mod(itemp,8)) + else + gref(1:ngroups)=0 + endif + !write(78,*)'GREFs: ',(gref(j),j=1,ngroups) + ! + ! Extract Each Group's bit width + ! + !print *,'SAG2: ',nbitsgwidth,ngroups,iofst,idrstmpl(11) + if (nbitsgwidth.ne.0) then + call g2_gbytesc(cpack,gwidth,iofst,nbitsgwidth,0,ngroups) + itemp=nbitsgwidth*ngroups + iofst=iofst+(itemp) + if (mod(itemp,8).ne.0) iofst=iofst+(8-mod(itemp,8)) + else + gwidth(1:ngroups)=0 + endif + do j=1,ngroups + gwidth(j)=gwidth(j)+idrstmpl(11) + enddo + !write(78,*)'GWIDTHs: ',(gwidth(j),j=1,ngroups) + ! + ! Extract Each Group's length (number of values in each group) + ! + allocate(glen(ngroups),stat=is) + !print *,'ALLOC glen: ',is + !print *,'SAG3: ',nbitsglen,ngroups,iofst,idrstmpl(14),idrstmpl(13) + if (nbitsglen.ne.0) then + call g2_gbytesc(cpack,glen,iofst,nbitsglen,0,ngroups) + itemp=nbitsglen*ngroups + iofst=iofst+(itemp) + if (mod(itemp,8).ne.0) iofst=iofst+(8-mod(itemp,8)) + else + glen(1:ngroups)=0 + endif + do j=1,ngroups + glen(j)=(glen(j)*idrstmpl(14))+idrstmpl(13) + enddo + glen(ngroups)=idrstmpl(15) + !write(78,*)'GLENs: ',(glen(j),j=1,ngroups) + !print *,'GLENsum: ',sum(glen) + ! + ! Test to see if the group widths and lengths are consistent with number of + ! values, and length of section 7. + ! + totBit = 0 + totLen = 0 + do j=1,ngroups + totBit = totBit + (gwidth(j)*glen(j)); + totLen = totLen + glen(j); + enddo + if (totLen .NE. ndpts) then + ier=1 + return + endif + if ( (totBit/8) .GT. lensec) then + ier=1 + return + endif + ! + ! For each group, unpack data values + ! + if ( idrstmpl(7).eq.0 ) then ! no missing values + n=1 + do j=1,ngroups + !write(78,*)'NGP ',j,gwidth(j),glen(j),gref(j) + if (gwidth(j).ne.0) then + call g2_gbytesc(cpack,ifld(n),iofst,gwidth(j),0,glen(j)) + do k=1,glen(j) + ifld(n)=ifld(n)+gref(j) + n=n+1 + enddo + else + ifld(n:n+glen(j)-1)=gref(j) + n=n+glen(j) + endif + iofst=iofst+(gwidth(j)*glen(j)) + enddo + elseif ( idrstmpl(7).eq.1.OR.idrstmpl(7).eq.2 ) then + ! missing values included + allocate(ifldmiss(ndpts)) + !ifldmiss=0 + n=1 + non=1 + do j=1,ngroups + !print *,'SAGNGP ',j,gwidth(j),glen(j),gref(j) + if (gwidth(j).ne.0) then + msng1=(2**gwidth(j))-1 + msng2=msng1-1 + call g2_gbytesc(cpack,ifld(n),iofst,gwidth(j),0,glen(j)) + iofst=iofst+(gwidth(j)*glen(j)) + do k=1,glen(j) + if (ifld(n).eq.msng1) then + ifldmiss(n)=1 + elseif (idrstmpl(7).eq.2.AND.ifld(n).eq.msng2) then + ifldmiss(n)=2 + else + ifldmiss(n)=0 + ifld(non)=ifld(n)+gref(j) + non=non+1 + endif + n=n+1 + enddo + else + msng1=(2**nbitsgref)-1 + msng2=msng1-1 + if (gref(j).eq.msng1) then + ifldmiss(n:n+glen(j)-1)=1 + !ifld(n:n+glen(j)-1)=0 + elseif (idrstmpl(7).eq.2.AND.gref(j).eq.msng2) then + ifldmiss(n:n+glen(j)-1)=2 + !ifld(n:n+glen(j)-1)=0 + else + ifldmiss(n:n+glen(j)-1)=0 + ifld(non:non+glen(j)-1)=gref(j) + non=non+glen(j) + endif + n=n+glen(j) + endif + enddo + endif + !write(78,*)'IFLDs: ',(ifld(j),j=1,ndpts) + + if ( allocated(gref) ) deallocate(gref) + if ( allocated(gwidth) ) deallocate(gwidth) + if ( allocated(glen) ) deallocate(glen) + ! + ! If using spatial differences, add overall min value, and + ! sum up recursively + ! + if (idrsnum.eq.3) then ! spatial differencing + if (idrstmpl(17).eq.1) then ! first order + ifld(1)=ival1 + if ( idrstmpl(7).eq.0 ) then ! no missing values + itemp=ndpts + else + itemp=non-1 + endif + do n=2,itemp + ifld(n)=ifld(n)+minsd + ifld(n)=ifld(n)+ifld(n-1) + enddo + elseif (idrstmpl(17).eq.2) then ! second order + ifld(1)=ival1 + ifld(2)=ival2 + if ( idrstmpl(7).eq.0 ) then ! no missing values + itemp=ndpts + else + itemp=non-1 + endif + do n=3,itemp + ifld(n)=ifld(n)+minsd + ifld(n)=ifld(n)+(2*ifld(n-1))-ifld(n-2) + enddo + endif + !write(78,*)'IFLDs: ',(ifld(j),j=1,ndpts) + endif + ! + ! Scale data back to original form + ! + !print *,'SAGT: ',ref,bscale,dscale + if ( idrstmpl(7).eq.0 ) then ! no missing values + do n=1,ndpts + fld(n)=((real(ifld(n))*bscale)+ref)*dscale + !write(78,*)'SAG ',n,fld(n),ifld(n),bscale,ref,1./dscale + enddo + elseif ( idrstmpl(7).eq.1.OR.idrstmpl(7).eq.2 ) then + ! missing values included + non=1 + do n=1,ndpts + if ( ifldmiss(n).eq.0 ) then + fld(n)=((real(ifld(non))*bscale)+ref)*dscale + !print *,'SAG ',n,fld(n),ifld(non),bscale,ref,dscale + non=non+1 + elseif ( ifldmiss(n).eq.1 ) then + fld(n)=rmiss1 + elseif ( ifldmiss(n).eq.2 ) then + fld(n)=rmiss2 + endif + enddo + if ( allocated(ifldmiss) ) deallocate(ifldmiss) + endif + + if ( allocated(ifld) ) deallocate(ifld) + + !open(10,form='unformatted',recl=24180,access='direct') + !read(10,rec=1) (fldo(k),k=1,6045) + !do i =1,6045 + ! print *,i,fldo(i),fld(i),fldo(i)-fld(i) + !enddo + + return +end subroutine comunpack diff --git a/src/comunpack.f b/src/comunpack.f deleted file mode 100644 index 8e7a4d55..00000000 --- a/src/comunpack.f +++ /dev/null @@ -1,319 +0,0 @@ -!> @file -!> @brief Unpack a data field that was packed using a complex packing -!> algorithm as defined in the GRIB2 documention. -!> @author Stephen Gilbert @date 2000-06-21 - -!> Unpack a data field that was packed using a complex packing -!> algorithm as defined in the GRIB2 documention. -!> -!> This subroutine Supports GRIB2 complex packing templates with or -!> without spatial differences: It supports GRIB2 complex packing -!> templates with or without spatial differences: Data Representation -!> Tables [5.2] -!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table5-2.shtml) -!> and -!> [5.3](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table5-3.shtml)). -!> -!> @param[in] cpack The packed data field (character*1 array). -!> @param[in] len Length of packed field cpack. -!> @param[in] lensec Length of section 7 (used for error checking). -!> @param[in] idrsnum Data Representation Template number 5.N. Must -!> equal 2 or 3. -!> @param[in] idrstmpl The array of values for Data Representation -!> Template 5.2 or 5.3. -!> @param[in] ndpts The number of data values to unpack. -!> @param[out] fld Contains the unpacked data values. -!> @param[out] ier Error return: -!> - 0 = No error. -!> - 1 = Problem - inconsistent group lengths of widths. -!> -!> @author Stephen Gilbert @date 2000-06-21 - subroutine comunpack(cpack,len,lensec,idrsnum,idrstmpl,ndpts, - & fld,ier) - - character(len=1),intent(in) :: cpack(len) - integer,intent(in) :: ndpts,len - integer,intent(in) :: idrstmpl(*) - real,intent(out) :: fld(ndpts) - - integer,allocatable :: ifld(:),ifldmiss(:) - integer(4) :: ieee - integer,allocatable :: gref(:),gwidth(:),glen(:) - real :: ref,bscale,dscale,rmiss1,rmiss2 -! real :: fldo(6045) - integer :: totBit, totLen - - ier=0 - !print *,'IDRSTMPL: ',(idrstmpl(j),j=1,16) - ieee = idrstmpl(1) - call rdieee(ieee,ref,1) - bscale = 2.0**real(idrstmpl(2)) - dscale = 10.0**real(-idrstmpl(3)) - nbitsgref = idrstmpl(4) - itype = idrstmpl(5) - ngroups = idrstmpl(10) - nbitsgwidth = idrstmpl(12) - nbitsglen = idrstmpl(16) - if (idrsnum.eq.3) then - nbitsd=idrstmpl(18)*8 - endif - - ! Constant field - - if (ngroups.eq.0) then - do j=1,ndpts - fld(j)=ref - enddo - return - endif - - iofst=0 - allocate(ifld(ndpts),stat=is) - !print *,'ALLOC ifld: ',is,ndpts - allocate(gref(ngroups),stat=is) - !print *,'ALLOC gref: ',is - allocate(gwidth(ngroups),stat=is) - !print *,'ALLOC gwidth: ',is -! -! Get missing values, if supplied -! - if ( idrstmpl(7).eq.1 ) then - if (itype.eq.0) then - call rdieee(idrstmpl(8),rmiss1,1) - else - rmiss1=real(idrstmpl(8)) - endif - elseif ( idrstmpl(7).eq.2 ) then - if (itype.eq.0) then - call rdieee(idrstmpl(8),rmiss1,1) - call rdieee(idrstmpl(9),rmiss2,1) - else - rmiss1=real(idrstmpl(8)) - rmiss2=real(idrstmpl(9)) - endif - endif - !print *,'RMISSs: ',rmiss1,rmiss2,ref -! -! Extract Spatial differencing values, if using DRS Template 5.3 -! - if (idrsnum.eq.3) then - if (nbitsd.ne.0) then - call g2_gbytec(cpack,ival1,iofst,nbitsd) - iofst=iofst+nbitsd - if (idrstmpl(17).eq.2) then - call g2_gbytec(cpack,ival2,iofst,nbitsd) - iofst=iofst+nbitsd - endif - call g2_gbytec(cpack,isign,iofst,1) - iofst=iofst+1 - call g2_gbytec(cpack,minsd,iofst,nbitsd-1) - iofst=iofst+nbitsd-1 - if (isign.eq.1) minsd=-minsd - else - ival1=0 - ival2=0 - minsd=0 - endif - !print *,'SDu ',ival1,ival2,minsd,nbitsd - endif -! -! Extract Each Group's reference value -! - !print *,'SAG1: ',nbitsgref,ngroups,iofst - if (nbitsgref.ne.0) then - call g2_gbytesc(cpack,gref,iofst,nbitsgref,0,ngroups) - itemp=nbitsgref*ngroups - iofst=iofst+(itemp) - if (mod(itemp,8).ne.0) iofst=iofst+(8-mod(itemp,8)) - else - gref(1:ngroups)=0 - endif - !write(78,*)'GREFs: ',(gref(j),j=1,ngroups) -! -! Extract Each Group's bit width -! - !print *,'SAG2: ',nbitsgwidth,ngroups,iofst,idrstmpl(11) - if (nbitsgwidth.ne.0) then - call g2_gbytesc(cpack,gwidth,iofst,nbitsgwidth,0,ngroups) - itemp=nbitsgwidth*ngroups - iofst=iofst+(itemp) - if (mod(itemp,8).ne.0) iofst=iofst+(8-mod(itemp,8)) - else - gwidth(1:ngroups)=0 - endif - do j=1,ngroups - gwidth(j)=gwidth(j)+idrstmpl(11) - enddo - !write(78,*)'GWIDTHs: ',(gwidth(j),j=1,ngroups) -! -! Extract Each Group's length (number of values in each group) -! - allocate(glen(ngroups),stat=is) - !print *,'ALLOC glen: ',is - !print *,'SAG3: ',nbitsglen,ngroups,iofst,idrstmpl(14),idrstmpl(13) - if (nbitsglen.ne.0) then - call g2_gbytesc(cpack,glen,iofst,nbitsglen,0,ngroups) - itemp=nbitsglen*ngroups - iofst=iofst+(itemp) - if (mod(itemp,8).ne.0) iofst=iofst+(8-mod(itemp,8)) - else - glen(1:ngroups)=0 - endif - do j=1,ngroups - glen(j)=(glen(j)*idrstmpl(14))+idrstmpl(13) - enddo - glen(ngroups)=idrstmpl(15) - !write(78,*)'GLENs: ',(glen(j),j=1,ngroups) - !print *,'GLENsum: ',sum(glen) -! -! Test to see if the group widths and lengths are consistent with number of -! values, and length of section 7. -! - totBit = 0 - totLen = 0 - do j=1,ngroups - totBit = totBit + (gwidth(j)*glen(j)); - totLen = totLen + glen(j); - enddo - if (totLen .NE. ndpts) then - ier=1 - return - endif - if ( (totBit/8) .GT. lensec) then - ier=1 - return - endif -! -! For each group, unpack data values -! - if ( idrstmpl(7).eq.0 ) then ! no missing values - n=1 - do j=1,ngroups - !write(78,*)'NGP ',j,gwidth(j),glen(j),gref(j) - if (gwidth(j).ne.0) then - call g2_gbytesc(cpack,ifld(n),iofst,gwidth(j),0,glen(j)) - do k=1,glen(j) - ifld(n)=ifld(n)+gref(j) - n=n+1 - enddo - else - ifld(n:n+glen(j)-1)=gref(j) - n=n+glen(j) - endif - iofst=iofst+(gwidth(j)*glen(j)) - enddo - elseif ( idrstmpl(7).eq.1.OR.idrstmpl(7).eq.2 ) then - ! missing values included - allocate(ifldmiss(ndpts)) - !ifldmiss=0 - n=1 - non=1 - do j=1,ngroups - !print *,'SAGNGP ',j,gwidth(j),glen(j),gref(j) - if (gwidth(j).ne.0) then - msng1=(2**gwidth(j))-1 - msng2=msng1-1 - call g2_gbytesc(cpack,ifld(n),iofst,gwidth(j),0,glen(j)) - iofst=iofst+(gwidth(j)*glen(j)) - do k=1,glen(j) - if (ifld(n).eq.msng1) then - ifldmiss(n)=1 - elseif (idrstmpl(7).eq.2.AND.ifld(n).eq.msng2) then - ifldmiss(n)=2 - else - ifldmiss(n)=0 - ifld(non)=ifld(n)+gref(j) - non=non+1 - endif - n=n+1 - enddo - else - msng1=(2**nbitsgref)-1 - msng2=msng1-1 - if (gref(j).eq.msng1) then - ifldmiss(n:n+glen(j)-1)=1 - !ifld(n:n+glen(j)-1)=0 - elseif (idrstmpl(7).eq.2.AND.gref(j).eq.msng2) then - ifldmiss(n:n+glen(j)-1)=2 - !ifld(n:n+glen(j)-1)=0 - else - ifldmiss(n:n+glen(j)-1)=0 - ifld(non:non+glen(j)-1)=gref(j) - non=non+glen(j) - endif - n=n+glen(j) - endif - enddo - endif - !write(78,*)'IFLDs: ',(ifld(j),j=1,ndpts) - - if ( allocated(gref) ) deallocate(gref) - if ( allocated(gwidth) ) deallocate(gwidth) - if ( allocated(glen) ) deallocate(glen) -! -! If using spatial differences, add overall min value, and -! sum up recursively -! - if (idrsnum.eq.3) then ! spatial differencing - if (idrstmpl(17).eq.1) then ! first order - ifld(1)=ival1 - if ( idrstmpl(7).eq.0 ) then ! no missing values - itemp=ndpts - else - itemp=non-1 - endif - do n=2,itemp - ifld(n)=ifld(n)+minsd - ifld(n)=ifld(n)+ifld(n-1) - enddo - elseif (idrstmpl(17).eq.2) then ! second order - ifld(1)=ival1 - ifld(2)=ival2 - if ( idrstmpl(7).eq.0 ) then ! no missing values - itemp=ndpts - else - itemp=non-1 - endif - do n=3,itemp - ifld(n)=ifld(n)+minsd - ifld(n)=ifld(n)+(2*ifld(n-1))-ifld(n-2) - enddo - endif - !write(78,*)'IFLDs: ',(ifld(j),j=1,ndpts) - endif -! -! Scale data back to original form -! - !print *,'SAGT: ',ref,bscale,dscale - if ( idrstmpl(7).eq.0 ) then ! no missing values - do n=1,ndpts - fld(n)=((real(ifld(n))*bscale)+ref)*dscale - !write(78,*)'SAG ',n,fld(n),ifld(n),bscale,ref,1./dscale - enddo - elseif ( idrstmpl(7).eq.1.OR.idrstmpl(7).eq.2 ) then - ! missing values included - non=1 - do n=1,ndpts - if ( ifldmiss(n).eq.0 ) then - fld(n)=((real(ifld(non))*bscale)+ref)*dscale - !print *,'SAG ',n,fld(n),ifld(non),bscale,ref,dscale - non=non+1 - elseif ( ifldmiss(n).eq.1 ) then - fld(n)=rmiss1 - elseif ( ifldmiss(n).eq.2 ) then - fld(n)=rmiss2 - endif - enddo - if ( allocated(ifldmiss) ) deallocate(ifldmiss) - endif - - if ( allocated(ifld) ) deallocate(ifld) - - !open(10,form='unformatted',recl=24180,access='direct') - !read(10,rec=1) (fldo(k),k=1,6045) - !do i =1,6045 - ! print *,i,fldo(i),fld(i),fldo(i)-fld(i) - !enddo - - return - end diff --git a/src/getdim.F90 b/src/getdim.F90 new file mode 100644 index 00000000..38b3db88 --- /dev/null +++ b/src/getdim.F90 @@ -0,0 +1,93 @@ +!> @file +!> @brief Return info from section 3, the GRIB2 Grid Definition +!> Section. +!> @author Stephen Gilbert @date 2002-12-11 + +!> Return the dimensions and scanning mode of a grid +!> definition packed in the Grid Definition Section. +!> +!> @param[in] csec3 Character array that contains the packed GRIB2 +!> GDS. +!> @param[in] lcsec3 Length (in octets) of section 3. +!> @param[out] width x (or i) dimension of the grid. +!> @param[out] height y (or j) dimension of the grid. +!> @param[out] iscan Scanning mode (see [Code Table 3.4] +!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table3-4.shtml)). +!> +!> @note Returns width and height set to zero, if grid template not +!> recognized. +!> +!> @author Stephen Gilbert @date 2002-12-11 +subroutine getdim(csec3,lcsec3,width,height,iscan) + implicit none + + character(len=1),intent(in) :: csec3(*) + integer,intent(in) :: lcsec3 + integer,intent(out) :: width,height,iscan + + integer,pointer,dimension(:) :: igdstmpl,list_opt + integer :: igds(5) + integer iofst,igdtlen,num_opt,jerr + + interface + subroutine gf_unpack3(cgrib,lcgrib,iofst,igds,igdstmpl, & + mapgridlen,ideflist,idefnum,ierr) + character(len=1),intent(in) :: cgrib(lcgrib) + integer,intent(in) :: lcgrib + integer,intent(inout) :: iofst + integer,pointer,dimension(:) :: igdstmpl,ideflist + integer,intent(out) :: igds(5) + integer,intent(out) :: ierr,idefnum + end subroutine gf_unpack3 + end interface + + nullify(igdstmpl,list_opt) + + iofst=0 ! set offset to beginning of section + call gf_unpack3(csec3,lcsec3,iofst,igds,igdstmpl, & + igdtlen,list_opt,num_opt,jerr) + if (jerr.eq.0) then + selectcase( igds(5) ) ! Template number + case (0:3) ! Lat/Lon + width=igdstmpl(8) + height=igdstmpl(9) + iscan=igdstmpl(19) + case (10) ! Mercator + width=igdstmpl(8) + height=igdstmpl(9) + iscan=igdstmpl(16) + case (20) ! Polar Stereographic + width=igdstmpl(8) + height=igdstmpl(9) + iscan=igdstmpl(18) + case (30) ! Lambert Conformal + width=igdstmpl(8) + height=igdstmpl(9) + iscan=igdstmpl(18) + case (40:43) ! Gaussian + width=igdstmpl(8) + height=igdstmpl(9) + iscan=igdstmpl(19) + case (90) ! Space View/Orthographic + width=igdstmpl(8) + height=igdstmpl(9) + iscan=igdstmpl(17) + case (110) ! Equatorial Azimuthal + width=igdstmpl(8) + height=igdstmpl(9) + iscan=igdstmpl(16) + case default + width=0 + height=0 + iscan=0 + end select + else + width=0 + height=0 + endif + + if (associated(igdstmpl)) deallocate(igdstmpl) + if (associated(list_opt)) deallocate(list_opt) + + return +end subroutine getdim diff --git a/src/getdim.f b/src/getdim.f deleted file mode 100644 index 5d0d23a6..00000000 --- a/src/getdim.f +++ /dev/null @@ -1,93 +0,0 @@ -!> @file -!> @brief Return info from section 3, the GRIB2 Grid Definition -!> Section. -!> @author Stephen Gilbert @date 2002-12-11 - -!> Return the dimensions and scanning mode of a grid -!> definition packed in the Grid Definition Section. -!> -!> @param[in] csec3 Character array that contains the packed GRIB2 -!> GDS. -!> @param[in] lcsec3 Length (in octets) of section 3. -!> @param[out] width x (or i) dimension of the grid. -!> @param[out] height y (or j) dimension of the grid. -!> @param[out] iscan Scanning mode (see [Code Table 3.4] -!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table3-4.shtml)). -!> -!> @note Returns width and height set to zero, if grid template not -!> recognized. -!> -!> @author Stephen Gilbert @date 2002-12-11 - subroutine getdim(csec3,lcsec3,width,height,iscan) - implicit none - - character(len=1),intent(in) :: csec3(*) - integer,intent(in) :: lcsec3 - integer,intent(out) :: width,height,iscan - - integer,pointer,dimension(:) :: igdstmpl,list_opt - integer :: igds(5) - integer iofst,igdtlen,num_opt,jerr - - interface - subroutine gf_unpack3(cgrib,lcgrib,iofst,igds,igdstmpl, - & mapgridlen,ideflist,idefnum,ierr) - character(len=1),intent(in) :: cgrib(lcgrib) - integer,intent(in) :: lcgrib - integer,intent(inout) :: iofst - integer,pointer,dimension(:) :: igdstmpl,ideflist - integer,intent(out) :: igds(5) - integer,intent(out) :: ierr,idefnum - end subroutine gf_unpack3 - end interface - - nullify(igdstmpl,list_opt) - ! - iofst=0 ! set offset to beginning of section - call gf_unpack3(csec3,lcsec3,iofst,igds,igdstmpl, - & igdtlen,list_opt,num_opt,jerr) - if (jerr.eq.0) then - selectcase( igds(5) ) ! Template number - case (0:3) ! Lat/Lon - width=igdstmpl(8) - height=igdstmpl(9) - iscan=igdstmpl(19) - case (10) ! Mercator - width=igdstmpl(8) - height=igdstmpl(9) - iscan=igdstmpl(16) - case (20) ! Polar Stereographic - width=igdstmpl(8) - height=igdstmpl(9) - iscan=igdstmpl(18) - case (30) ! Lambert Conformal - width=igdstmpl(8) - height=igdstmpl(9) - iscan=igdstmpl(18) - case (40:43) ! Gaussian - width=igdstmpl(8) - height=igdstmpl(9) - iscan=igdstmpl(19) - case (90) ! Space View/Orthographic - width=igdstmpl(8) - height=igdstmpl(9) - iscan=igdstmpl(17) - case (110) ! Equatorial Azimuthal - width=igdstmpl(8) - height=igdstmpl(9) - iscan=igdstmpl(16) - case default - width=0 - height=0 - iscan=0 - end select - else - width=0 - height=0 - endif - ! - if (associated(igdstmpl)) deallocate(igdstmpl) - if (associated(list_opt)) deallocate(list_opt) - - return - end diff --git a/src/getpoly.F90 b/src/getpoly.F90 new file mode 100644 index 00000000..0acff415 --- /dev/null +++ b/src/getpoly.F90 @@ -0,0 +1,68 @@ +!> @file +!> @brief Return the J, K, and M pentagonal resolution parameters +!> specified in a GRIB2 Grid Definition Section. +!> @author Stephen Gilbert @date 2002-12-11 + +!> Return the J, K, and M pentagonal resolution +!> parameters specified in a GRIB2 Grid Definition Section used +!> spherical harmonic coefficients using GDT 5.50 through 5.53. +!> +!> @param[in] csec3 Character array containing the packed GRIB2 GDS +!> @param[in] lcsec3 Length (in octets) of section 3 +!> @param[out] JJ =J pentagonal resolution parameter +!> @param[out] KK =K pentagonal resolution parameter +!> @param[out] MM =M pentagonal resolution parameter +!> +!> @note Returns JJ, KK, and MM set to zero, if grid template not +!> recognized. +!> +!> @author Stephen Gilbert @date 2002-12-11 +subroutine getpoly(csec3,lcsec3,jj,kk,mm) + implicit none + + character(len=1),intent(in) :: csec3(*) + integer,intent(in) :: lcsec3 + integer,intent(out) :: jj,kk,mm + + integer,pointer,dimension(:) :: igdstmpl,list_opt + integer :: igds(5) + integer iofst,igdtlen,num_opt,jerr + + interface + subroutine gf_unpack3(cgrib,lcgrib,iofst,igds,igdstmpl, & + mapgridlen,ideflist,idefnum,ierr) + character(len=1),intent(in) :: cgrib(lcgrib) + integer,intent(in) :: lcgrib + integer,intent(inout) :: iofst + integer,pointer,dimension(:) :: igdstmpl,ideflist + integer,intent(out) :: igds(5) + integer,intent(out) :: ierr,idefnum + end subroutine gf_unpack3 + end interface + + nullify(igdstmpl,list_opt) + + iofst=0 ! set offset to beginning of section + call gf_unpack3(csec3,lcsec3,iofst,igds,igdstmpl, & + igdtlen,list_opt,num_opt,jerr) + if (jerr.eq.0) then + selectcase( igds(5) ) ! Template number + case (50:53) ! Spherical harmonic coefficients + jj=igdstmpl(1) + kk=igdstmpl(2) + mm=igdstmpl(3) + case default + jj=0 + kk=0 + mm=0 + end select + else + jj=0 + kk=0 + mm=0 + endif + + if (associated(igdstmpl)) deallocate(igdstmpl) + if (associated(list_opt)) deallocate(list_opt) + +end subroutine getpoly diff --git a/src/getpoly.f b/src/getpoly.f deleted file mode 100644 index a74d444c..00000000 --- a/src/getpoly.f +++ /dev/null @@ -1,68 +0,0 @@ -!> @file -!> @brief Return the J, K, and M pentagonal resolution parameters -!> specified in a GRIB2 Grid Definition Section. -!> @author Stephen Gilbert @date 2002-12-11 - -!> Return the J, K, and M pentagonal resolution -!> parameters specified in a GRIB2 Grid Definition Section used -!> spherical harmonic coefficients using GDT 5.50 through 5.53. -!> -!> @param[in] csec3 Character array containing the packed GRIB2 GDS -!> @param[in] lcsec3 Length (in octets) of section 3 -!> @param[out] JJ =J pentagonal resolution parameter -!> @param[out] KK =K pentagonal resolution parameter -!> @param[out] MM =M pentagonal resolution parameter -!> -!> @note Returns JJ, KK, and MM set to zero, if grid template not -!> recognized. -!> -!> @author Stephen Gilbert @date 2002-12-11 - subroutine getpoly(csec3,lcsec3,jj,kk,mm) - implicit none - - character(len=1),intent(in) :: csec3(*) - integer,intent(in) :: lcsec3 - integer,intent(out) :: jj,kk,mm - - integer,pointer,dimension(:) :: igdstmpl,list_opt - integer :: igds(5) - integer iofst,igdtlen,num_opt,jerr - - interface - subroutine gf_unpack3(cgrib,lcgrib,iofst,igds,igdstmpl, - & mapgridlen,ideflist,idefnum,ierr) - character(len=1),intent(in) :: cgrib(lcgrib) - integer,intent(in) :: lcgrib - integer,intent(inout) :: iofst - integer,pointer,dimension(:) :: igdstmpl,ideflist - integer,intent(out) :: igds(5) - integer,intent(out) :: ierr,idefnum - end subroutine gf_unpack3 - end interface - - nullify(igdstmpl,list_opt) - - iofst=0 ! set offset to beginning of section - call gf_unpack3(csec3,lcsec3,iofst,igds,igdstmpl, - & igdtlen,list_opt,num_opt,jerr) - if (jerr.eq.0) then - selectcase( igds(5) ) ! Template number - case (50:53) ! Spherical harmonic coefficients - jj=igdstmpl(1) - kk=igdstmpl(2) - mm=igdstmpl(3) - case default - jj=0 - kk=0 - mm=0 - end select - else - jj=0 - kk=0 - mm=0 - endif - - if (associated(igdstmpl)) deallocate(igdstmpl) - if (associated(list_opt)) deallocate(list_opt) - - end diff --git a/src/intmath.F90 b/src/intmath.F90 new file mode 100644 index 00000000..e3aeeb4d --- /dev/null +++ b/src/intmath.F90 @@ -0,0 +1,268 @@ +!> @file +!> @brief Define math functions used by compack(), +!> simpack(), and misspack(). +!> @author Stephen Gilbert @date 2000-06-21 + +!> @brief Define math functions used by compack(), +!> simpack(), and misspack(). +!> +!> This module includes the following functions: +!> - ilog2 Calculate log(x)/log(2). +!> - ilog2_8 for 8 bit integer numbers. +!> - ilog2_4 for 4 bit integer numbers. +!> - ilog2_2 for 2 bit integer numbers. +!> - ilog2_1 for 1 bit integer numbers. +!> - i1log2 Calculate log(x+1)/log(2) unless x=maxint, in which case log(x)/log(2). +!> - i1log2_8 for 8 bit integer numbers. +!> - i1log2_4 for 4 bit integer numbers. +!> - i1log2_2 for 2 bit integer numbers. +!> - i1log2_1 for 1 bit integer numbers. +!> +!> @author Stephen Gilbert @date 2000-06-21 +module intmath + implicit none + + interface ilog2 + ! log(x)/log(2) + module procedure ilog2_8 + module procedure ilog2_4 + module procedure ilog2_2 + module procedure ilog2_1 + end interface ilog2 + + interface i1log2 + ! log(x+1)/log(2) unless x=maxint, in which case log(x)/log(2) + module procedure i1log2_8 + module procedure i1log2_4 + module procedure i1log2_2 + module procedure i1log2_1 + end interface i1log2 + +contains + + !> This function returns log(x+1)/log(2) unless x=maxint, in + !> which case log(x)/log(2) for 8 bit integer numbers. + !> @param[in] ival 8 bit integer numbers. + !> @return value for log(x+1)/log(2) + !> @author Stephen Gilbert @date 2000-06-21 + function i1log2_8(ival) + implicit none + integer(kind=8), value :: ival + integer(kind=8)::i1log2_8 + integer(kind=8), parameter :: one=1 + if(ival+one This function returns log(x+1)/log(2) unless x=maxint, in + !> which case log(x)/log(2) for 4 bit integer numbers. + !> @param[in] ival 4 bit integer numbers. + !> @return value for log(x+1)/log(2) + !> @author Stephen Gilbert @date 2000-06-21 + function i1log2_4(ival) + implicit none + integer(kind=4), value :: ival + integer(kind=4)::i1log2_4 + integer(kind=4), parameter :: one=1 + if(ival+one This function returns log(x+1)/log(2) unless x=maxint, in + !> which case log(x)/log(2) for 2 bit integer numbers. + !> @param[in] ival 2 bit integer numbers. + !> @return value for log(x+1)/log(2) + !> @author Stephen Gilbert @date 2000-06-21 + function i1log2_2(ival) + implicit none + integer(kind=2), value :: ival + integer(kind=2)::i1log2_2 + integer(kind=2), parameter :: one = 1_2 + if(ival+one This function returns log(x+1)/log(2) unless x=maxint, in + !> which case log(x)/log(2) for 1 bit integer numbers. + !> @param[in] ival 1 bit integer numbers. + !> @return value for log(x+1)/log(2) + !> @author Stephen Gilbert @date 2000-06-21 + function i1log2_1(ival) + implicit none + integer(kind=1), value :: ival + integer(kind=1)::i1log2_1 + integer(kind=1), parameter :: one = 1_1 + if(ival+one This function returns log(x)/log(2) for 8 bit integer numbers. + !> @param[in] i_in 8 bit integer numbers. + !> @return value for log(x)/log(2) + !> @author Stephen Gilbert @date 2000-06-21 + function ilog2_8(i_in) + implicit none + integer(kind=8), value :: i_in + integer(kind=8)::ilog2_8,i + ilog2_8=0 + i=i_in + if(i<=0) return + if(iand(i,i-1)/=0) then + !write(0,*) 'iand i-1' + ilog2_8=1 + endif + if(iand(i,Z'FFFFFFFF00000000')/=0) then + ilog2_8=ilog2_8+32 + i=ishft(i,-32) + !write(0,*) 'iand ffffffff',i,ilog2_8 + endif + if(iand(i,Z'00000000FFFF0000')/=0) then + ilog2_8=ilog2_8+16 + i=ishft(i,-16) + !write(0,*) 'iand ffff' ,i,ilog2_8 + endif + if(iand(i,Z'000000000000FF00')/=0) then + ilog2_8=ilog2_8+8 + i=ishft(i,-8) + !write(0,*) 'iand ff',i,ilog2_8 + endif + if(iand(i,Z'00000000000000F0')/=0) then + ilog2_8=ilog2_8+4 + i=ishft(i,-4) + !write(0,*) 'iand f',i,ilog2_8 + endif + if(iand(i,Z'000000000000000C')/=0) then + ilog2_8=ilog2_8+2 + i=ishft(i,-2) + !write(0,*) 'iand c',i,ilog2_8 + endif + if(iand(i,Z'0000000000000002')/=0) then + ilog2_8=ilog2_8+1 + i=ishft(i,-1) + !write(0,*) 'iand 2',i,ilog2_8 + endif + end function ilog2_8 + + !> This function returns log(x)/log(2) for 4 bit integer numbers. + !> @param[in] i_in 4 bit integer numbers. + !> @return value for log(x)/log(2) + !> @author Stephen Gilbert @date 2000-06-21 + function ilog2_4(i_in) + implicit none + integer(kind=4), value :: i_in + integer(kind=4)::ilog2_4,i + ilog2_4=0 + i=i_in + if(i<=0) return + if(iand(i,i-1)/=0) then + !write(0,*) 'iand i-1' + ilog2_4=1 + endif + if(iand(i,Z'FFFF0000')/=0) then + ilog2_4=ilog2_4+16 + i=ishft(i,-16) + !write(0,*) 'iand ffff' ,i,ilog2_4 + endif + if(iand(i,Z'0000FF00')/=0) then + ilog2_4=ilog2_4+8 + i=ishft(i,-8) + !write(0,*) 'iand ff',i,ilog2_4 + endif + if(iand(i,Z'000000F0')/=0) then + ilog2_4=ilog2_4+4 + i=ishft(i,-4) + !write(0,*) 'iand f',i,ilog2_4 + endif + if(iand(i,Z'0000000C')/=0) then + ilog2_4=ilog2_4+2 + i=ishft(i,-2) + !write(0,*) 'iand c',i,ilog2_4 + endif + if(iand(i,Z'00000002')/=0) then + ilog2_4=ilog2_4+1 + i=ishft(i,-1) + !write(0,*) 'iand 2',i,ilog2_4 + endif + end function ilog2_4 + + !> This function returns log(x)/log(2) for 2 bit integer numbers. + !> @param[in] i_in 2 bit integer numbers. + !> @return value for log(x)/log(2) + !> @author Stephen Gilbert @date 2000-06-21 + function ilog2_2(i_in) + implicit none + integer(kind=2), value :: i_in + integer(kind=2)::ilog2_2,i + ilog2_2 = 0_2 + i=i_in + if(i<=0) return + if(iand(i,int(i-1,kind=2))/=0) then + !write(0,*) 'iand i-1' + ilog2_2 = 1_2 + endif + if(iand(i,Z'FF00')/=0) then + ilog2_2 = ilog2_2 + 8_2 + i=ishft(i,-8) + !write(0,*) 'iand ff',i,ilog2_2 + endif + if(iand(i,Z'00F0')/=0) then + ilog2_2 = ilog2_2 + 4_2 + i=ishft(i,-4) + !write(0,*) 'iand f',i,ilog2_2 + endif + if(iand(i,Z'000C')/=0) then + ilog2_2 = ilog2_2 + 2_2 + i=ishft(i,-2) + !write(0,*) 'iand c',i,ilog2_2 + endif + if(iand(i,Z'0002')/=0) then + ilog2_2 = ilog2_2 + 1_2 + i=ishft(i,-1) + !write(0,*) 'iand 2',i,ilog2_2 + endif + end function ilog2_2 + + !> This function returns log(x)/log(2) for 1 bit integer numbers. + !> @param[in] i_in 1 bit integer numbers. + !> @return value for log(x)/log(2) + !> @author Stephen Gilbert @date 2000-06-21 + function ilog2_1(i_in) + implicit none + integer(kind=1), value :: i_in + integer(kind=1)::ilog2_1,i + ilog2_1 = 0_1 + i=i_in + if(i<=0) return + if(iand(i,int(i-1,kind=1))/=0) then + !write(0,*) 'iand i-1' + ilog2_1 = 1_1 + endif + if(iand(i,Z'F0')/=0) then + ilog2_1 = ilog2_1 + 4_1 + i=ishft(i,-4) + !write(0,*) 'iand f',i,ilog2_1 + endif + if(iand(i,Z'0C')/=0) then + ilog2_1 = ilog2_1 + 2_1 + i=ishft(i,-2) + !write(0,*) 'iand c',i,ilog2_1 + endif + if(iand(i,Z'02')/=0) then + ilog2_1 = ilog2_1 + 1_1 + i=ishft(i,-1) + !write(0,*) 'iand 2',i,ilog2_1 + endif + end function ilog2_1 +end module intmath diff --git a/src/intmath.f b/src/intmath.f deleted file mode 100644 index c0fe9da3..00000000 --- a/src/intmath.f +++ /dev/null @@ -1,301 +0,0 @@ -!> @file -!> @brief Define math functions used by compack(), -!> simpack(), and misspack(). -!> @author Stephen Gilbert @date 2000-06-21 - -!> @brief Define math functions used by compack(), -!> simpack(), and misspack(). -!> -!> This module includes the following functions: -!> - ilog2 Calculate log(x)/log(2). -!> - ilog2_8 for 8 bit integer numbers. -!> - ilog2_4 for 4 bit integer numbers. -!> - ilog2_2 for 2 bit integer numbers. -!> - ilog2_1 for 1 bit integer numbers. -!> - i1log2 Calculate log(x+1)/log(2) unless x=maxint, in which case log(x)/log(2). -!> - i1log2_8 for 8 bit integer numbers. -!> - i1log2_4 for 4 bit integer numbers. -!> - i1log2_2 for 2 bit integer numbers. -!> - i1log2_1 for 1 bit integer numbers. -!> -!> @author Stephen Gilbert @date 2000-06-21 - module intmath - implicit none - - interface ilog2 - ! log(x)/log(2) - module procedure ilog2_8 - module procedure ilog2_4 - module procedure ilog2_2 - module procedure ilog2_1 - end interface ilog2 - - interface i1log2 - ! log(x+1)/log(2) unless x=maxint, in which case log(x)/log(2) - module procedure i1log2_8 - module procedure i1log2_4 - module procedure i1log2_2 - module procedure i1log2_1 - end interface i1log2 - - contains - -!> This function returns log(x+1)/log(2) unless x=maxint, in -!> which case log(x)/log(2) for 8 bit integer numbers. -!> @param[in] ival 8 bit integer numbers. -!> @return value for log(x+1)/log(2) -!> @author Stephen Gilbert @date 2000-06-21 - function i1log2_8(ival) - implicit none - integer(kind=8), value :: ival - integer(kind=8)::i1log2_8 - integer(kind=8), parameter :: one=1 - if(ival+one This function returns log(x+1)/log(2) unless x=maxint, in -!> which case log(x)/log(2) for 4 bit integer numbers. -!> @param[in] ival 4 bit integer numbers. -!> @return value for log(x+1)/log(2) -!> @author Stephen Gilbert @date 2000-06-21 - function i1log2_4(ival) - implicit none - integer(kind=4), value :: ival - integer(kind=4)::i1log2_4 - integer(kind=4), parameter :: one=1 - if(ival+one This function returns log(x+1)/log(2) unless x=maxint, in -!> which case log(x)/log(2) for 2 bit integer numbers. -!> @param[in] ival 2 bit integer numbers. -!> @return value for log(x+1)/log(2) -!> @author Stephen Gilbert @date 2000-06-21 - function i1log2_2(ival) - implicit none - integer(kind=2), value :: ival - integer(kind=2)::i1log2_2 - integer(kind=2), parameter :: one = 1_2 - if(ival+one This function returns log(x+1)/log(2) unless x=maxint, in -!> which case log(x)/log(2) for 1 bit integer numbers. -!> @param[in] ival 1 bit integer numbers. -!> @return value for log(x+1)/log(2) -!> @author Stephen Gilbert @date 2000-06-21 - function i1log2_1(ival) - implicit none - integer(kind=1), value :: ival - integer(kind=1)::i1log2_1 - integer(kind=1), parameter :: one = 1_1 - if(ival+one This function returns log(x)/log(2) for 8 bit integer numbers. -!> @param[in] i_in 8 bit integer numbers. -!> @return value for log(x)/log(2) -!> @author Stephen Gilbert @date 2000-06-21 - function ilog2_8(i_in) - implicit none - integer(kind=8), value :: i_in - integer(kind=8)::ilog2_8,i - ilog2_8=0 - i=i_in - if(i<=0) return - if(iand(i,i-1)/=0) then - !write(0,*) 'iand i-1' - ilog2_8=1 - endif - if(iand(i,Z'FFFFFFFF00000000')/=0) then - ilog2_8=ilog2_8+32 - i=ishft(i,-32) - !write(0,*) 'iand ffffffff',i,ilog2_8 - endif - if(iand(i,Z'00000000FFFF0000')/=0) then - ilog2_8=ilog2_8+16 - i=ishft(i,-16) - !write(0,*) 'iand ffff' ,i,ilog2_8 - endif - if(iand(i,Z'000000000000FF00')/=0) then - ilog2_8=ilog2_8+8 - i=ishft(i,-8) - !write(0,*) 'iand ff',i,ilog2_8 - endif - if(iand(i,Z'00000000000000F0')/=0) then - ilog2_8=ilog2_8+4 - i=ishft(i,-4) - !write(0,*) 'iand f',i,ilog2_8 - endif - if(iand(i,Z'000000000000000C')/=0) then - ilog2_8=ilog2_8+2 - i=ishft(i,-2) - !write(0,*) 'iand c',i,ilog2_8 - endif - if(iand(i,Z'0000000000000002')/=0) then - ilog2_8=ilog2_8+1 - i=ishft(i,-1) - !write(0,*) 'iand 2',i,ilog2_8 - endif - end function ilog2_8 - -!> This function returns log(x)/log(2) for 4 bit integer numbers. -!> @param[in] i_in 4 bit integer numbers. -!> @return value for log(x)/log(2) -!> @author Stephen Gilbert @date 2000-06-21 - function ilog2_4(i_in) - implicit none - integer(kind=4), value :: i_in - integer(kind=4)::ilog2_4,i - ilog2_4=0 - i=i_in - if(i<=0) return - if(iand(i,i-1)/=0) then - !write(0,*) 'iand i-1' - ilog2_4=1 - endif - if(iand(i,Z'FFFF0000')/=0) then - ilog2_4=ilog2_4+16 - i=ishft(i,-16) - !write(0,*) 'iand ffff' ,i,ilog2_4 - endif - if(iand(i,Z'0000FF00')/=0) then - ilog2_4=ilog2_4+8 - i=ishft(i,-8) - !write(0,*) 'iand ff',i,ilog2_4 - endif - if(iand(i,Z'000000F0')/=0) then - ilog2_4=ilog2_4+4 - i=ishft(i,-4) - !write(0,*) 'iand f',i,ilog2_4 - endif - if(iand(i,Z'0000000C')/=0) then - ilog2_4=ilog2_4+2 - i=ishft(i,-2) - !write(0,*) 'iand c',i,ilog2_4 - endif - if(iand(i,Z'00000002')/=0) then - ilog2_4=ilog2_4+1 - i=ishft(i,-1) - !write(0,*) 'iand 2',i,ilog2_4 - endif - end function ilog2_4 - -!> This function returns log(x)/log(2) for 2 bit integer numbers. -!> @param[in] i_in 2 bit integer numbers. -!> @return value for log(x)/log(2) -!> @author Stephen Gilbert @date 2000-06-21 - function ilog2_2(i_in) - implicit none - integer(kind=2), value :: i_in - integer(kind=2)::ilog2_2,i - ilog2_2 = 0_2 - i=i_in - if(i<=0) return - if(iand(i,int(i-1,kind=2))/=0) then - !write(0,*) 'iand i-1' - ilog2_2 = 1_2 - endif - if(iand(i,Z'FF00')/=0) then - ilog2_2 = ilog2_2 + 8_2 - i=ishft(i,-8) - !write(0,*) 'iand ff',i,ilog2_2 - endif - if(iand(i,Z'00F0')/=0) then - ilog2_2 = ilog2_2 + 4_2 - i=ishft(i,-4) - !write(0,*) 'iand f',i,ilog2_2 - endif - if(iand(i,Z'000C')/=0) then - ilog2_2 = ilog2_2 + 2_2 - i=ishft(i,-2) - !write(0,*) 'iand c',i,ilog2_2 - endif - if(iand(i,Z'0002')/=0) then - ilog2_2 = ilog2_2 + 1_2 - i=ishft(i,-1) - !write(0,*) 'iand 2',i,ilog2_2 - endif - end function ilog2_2 - -!> This function returns log(x)/log(2) for 1 bit integer numbers. -!> @param[in] i_in 1 bit integer numbers. -!> @return value for log(x)/log(2) -!> @author Stephen Gilbert @date 2000-06-21 - function ilog2_1(i_in) - implicit none - integer(kind=1), value :: i_in - integer(kind=1)::ilog2_1,i - ilog2_1 = 0_1 - i=i_in - if(i<=0) return - if(iand(i,int(i-1,kind=1))/=0) then - !write(0,*) 'iand i-1' - ilog2_1 = 1_1 - endif - if(iand(i,Z'F0')/=0) then - ilog2_1 = ilog2_1 + 4_1 - i=ishft(i,-4) - !write(0,*) 'iand f',i,ilog2_1 - endif - if(iand(i,Z'0C')/=0) then - ilog2_1 = ilog2_1 + 2_1 - i=ishft(i,-2) - !write(0,*) 'iand c',i,ilog2_1 - endif - if(iand(i,Z'02')/=0) then - ilog2_1 = ilog2_1 + 1_1 - i=ishft(i,-1) - !write(0,*) 'iand 2',i,ilog2_1 - endif - end function ilog2_1 - end module intmath -!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! - -c$$$ program test_intmath -c$$$ use intmath -c$$$ implicit none -c$$$ real(kind=16) :: temp -c$$$ real(kind=16), parameter :: alog2=log(2.0_16) -c$$$ integer(kind=8), parameter :: & -c$$$ & one=1,big=Z'7FFFFFFFFFFFFFFF',small=-2000000_8, & -c$$$ & check=Z'1FFFFFFF' -c$$$ integer(kind=8) :: ival, iret -c$$$ !$OMP PARALLEL DO PRIVATE(ival,temp,iret) -c$$$ do ival=small,big -c$$$ 10 format(Z16,' -- MISMATCH: ',I0,'=>',I0,' (',I0,' = ',F0.10,')') -c$$$ 20 format(Z16,' -- OKAY: ',I0,'=>',I0,' (',I0,' = ',F0.10,')') -c$$$ if(ival+one @file +!> @brief Pack up a data field using a GRIB2 algorithm with missing +!> value management. +!> @author Stephen Gilbert @date 2000-06-21 + +!> Pack up a data field using a GRIB2 algorithm with missing value +!> management. +!> +!> This subroutine packs up a data field using a complex packing +!> algorithm as defined in the GRIB2 documention. It supports GRIB2 +!> complex packing templates with or without spatial differences +!> (i.e. Data Representation Tables [5.2] +!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table5-2.shtml) +!> and +!> [5.3](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table5-3.shtml)). +!> +!> This subroutine also fills in GRIB2 Data Representation Template +!> 5.2 or 5.3 with the appropriate values. It assumes that Missing +!> Value Management is being used and that 1 or 2 missing values +!> appear in the data. +!> +!> @param[in] fld Contains the data values to pack. +!> @param[in] ndpts The number of data values in array fld. +!> @param[in] idrsnum Data Representation Template number. Must +!> equal 2 or 3. +!> @param[inout] idrstmpl Contains the array of values for Data +!> Representation Template 5.2 or 5.3 +!> - idrstmpl(1) Reference value - ignored on input set by compack routine. +!> - idrstmpl(2) Binary Scale Factor. +!> - idrstmpl(3) Decimal Scale Factor. +!> - idrstmpl(4) number of bits for each data value - ignored on input. +!> - idrstmpl(5) Original field type, currently ignored on input, set = 0 on +!> !output, Data values assumed to be reals. +!> - idrstmpl(6) = 0 use lossless compression or = 1 use lossy compression. +!> - idrstmpl(7) Missing value management. +!> - idrstmpl(8) Primary missing value. +!> - idrstmpl(9) Secondary missing value. +!> - idrstmpl(17) Order of Spatial Differencing (1 or 2). +!> @param[out] cpack The packed data field (character*1 array). +!> @param[out] lcpack length of packed field cpack. -1 is returned if +!> idrstmpl(7) is not set correctly. +!> +!> @author Stephen Gilbert @date 2000-06-21 +subroutine misspack(fld,ndpts,idrsnum,idrstmpl,cpack,lcpack) + + use intmath + integer,intent(in) :: ndpts,idrsnum + real,intent(in) :: fld(ndpts) + character(len=1),intent(out) :: cpack(*) + integer,intent(inout) :: idrstmpl(*) + integer,intent(out) :: lcpack + + real(4) :: ref, rmin4 + integer(4) :: iref + integer,allocatable :: ifld(:),ifldmiss(:),jfld(:) + integer,allocatable :: jmin(:),jmax(:),lbit(:) + integer,parameter :: zero=0 + integer,allocatable :: gref(:),gwidth(:),glen(:) + integer :: glength,grpwidth + logical :: simple_alg + logical :: have_rmin + + simple_alg = .false. + have_rmin = .false. + bscale=2.0**real(-idrstmpl(2)) + dscale=10.0**real(idrstmpl(3)) + missopt=idrstmpl(7) + if ( missopt.ne.1 .AND. missopt.ne.2 ) then + print *,'misspack: Unrecognized option.' + lcpack=-1 + return + else ! Get missing values + call rdieee(idrstmpl(8),rmissp,1) + if (missopt.eq.2) call rdieee(idrstmpl(9),rmisss,1) + endif + + ! Find min value of non-missing values in the data, + ! AND set up missing value mapping of the field. + allocate(ifldmiss(ndpts)) + + if ( missopt .eq. 1 ) then ! Primary missing value only + do j=1,ndpts + if (fld(j).eq.rmissp) then + ifldmiss(j)=1 + else + ifldmiss(j)=0 + if(have_rmin) then + if (fld(j).lt.rmin) rmin=fld(j) + else + rmin=fld(j) + have_rmin=.true. + endif + endif + enddo + if(.not.have_rmin) rmin=rmissp + endif + if ( missopt .eq. 2 ) then ! Primary and secondary missing values + do j=1,ndpts + if (fld(j).eq.rmissp) then + ifldmiss(j)=1 + elseif (fld(j).eq.rmisss) then + ifldmiss(j)=2 + else + ifldmiss(j)=0 + if(have_rmin) then + if (fld(j).lt.rmin) rmin=fld(j) + else + rmin=fld(j) + have_rmin=.true. + endif + endif + if(.not.have_rmin) rmin=rmissp + enddo + endif + + ! Allocate work arrays: + ! Note: -ifldmiss(j),j=1,ndpts is a map of original field indicating + ! which of the original data values + ! are primary missing (1), sencondary missing (2) or non-missing (0). + ! -jfld(j),j=1,nonmiss is a subarray of just the non-missing values from + ! the original field. + iofst=0 + allocate(ifld(ndpts)) + allocate(jfld(ndpts)) + allocate(gref(ndpts)) + allocate(gwidth(ndpts)) + allocate(glen(ndpts)) + ! + ! Scale original data + ! + nonmiss=0 + if (idrstmpl(2).eq.0) then ! No binary scaling + imin=nint(rmin*dscale) + !imax=nint(rmax*dscale) + rmin=real(imin) + do j=1,ndpts + if (ifldmiss(j).eq.0) then + nonmiss=nonmiss+1 + jfld(nonmiss)=max(0,nint(fld(j)*dscale)-imin) + endif + enddo + else ! Use binary scaling factor + rmin=rmin*dscale + do j=1,ndpts + if (ifldmiss(j).eq.0) then + nonmiss=nonmiss+1 + jfld(nonmiss)=max(0,nint(((fld(j)*dscale)-rmin)*bscale)) + endif + enddo + endif + ! + ! Calculate Spatial differences, if using DRS Template 5.3 + ! + if (idrsnum.eq.3) then ! spatial differences + if (idrstmpl(17).ne.1.and.idrstmpl(17).ne.2) idrstmpl(17)=2 + if (idrstmpl(17).eq.1) then ! first order + if(nonmiss<1) then + ival1=1.0 + else + ival1=jfld(1) + endif + do j=nonmiss,2,-1 + jfld(j)=jfld(j)-jfld(j-1) + enddo + if(nonmiss>0) jfld(1)=0 + elseif (idrstmpl(17).eq.2) then ! second order + if(nonmiss==1) then + ival1=jfld(1) + ival2=jfld(1) + elseif(nonmiss<1) then + ival1=1.0 + ival2=1.0 + else + ival1=jfld(1) + ival2=jfld(2) + endif + do j=nonmiss,3,-1 + jfld(j)=jfld(j)-(2*jfld(j-1))+jfld(j-2) + enddo + if(nonmiss>=1) jfld(1)=0 + if(nonmiss>=2) jfld(2)=0 + endif + + ! subtract min value from spatial diff field + isd=idrstmpl(17)+1 + minsd=minval(jfld(isd:nonmiss)) + do j=isd,nonmiss + jfld(j)=jfld(j)-minsd + enddo + + ! find num of bits need to store minsd and add 1 extra bit + ! to indicate sign + temp=i1log2(abs(minsd)) + nbitsd=ceiling(temp)+1 + + ! find num of bits need to store ifld(1) ( and ifld(2) + ! if using 2nd order differencing ) + maxorig=ival1 + if (idrstmpl(17).eq.2.and.ival2.gt.ival1) maxorig=ival2 + temp=i1log2(maxorig) + nbitorig=ceiling(temp)+1 + if (nbitorig.gt.nbitsd) nbitsd=nbitorig + ! increase number of bits to even multiple of 8 ( octet ) + if (mod(nbitsd,8).ne.0) nbitsd=nbitsd+(8-mod(nbitsd,8)) + + ! Store extra spatial differencing info into the packed + ! data section. + if (nbitsd.ne.0) then + ! pack first original value + if (ival1.ge.0) then + call g2_sbytec(cpack,ival1,iofst,nbitsd) + iofst=iofst+nbitsd + else + call g2_sbytec(cpack,1,iofst,1) + iofst=iofst+1 + call g2_sbytec(cpack,iabs(ival1),iofst,nbitsd-1) + iofst=iofst+nbitsd-1 + endif + if (idrstmpl(17).eq.2) then + ! pack second original value + if (ival2.ge.0) then + call g2_sbytec(cpack,ival2,iofst,nbitsd) + iofst=iofst+nbitsd + else + call g2_sbytec(cpack,1,iofst,1) + iofst=iofst+1 + call g2_sbytec(cpack,iabs(ival2),iofst,nbitsd-1) + iofst=iofst+nbitsd-1 + endif + endif + ! pack overall min of spatial differences + if (minsd.ge.0) then + call g2_sbytec(cpack,minsd,iofst,nbitsd) + iofst=iofst+nbitsd + else + call g2_sbytec(cpack,1,iofst,1) + iofst=iofst+1 + call g2_sbytec(cpack,iabs(minsd),iofst,nbitsd-1) + iofst=iofst+nbitsd-1 + endif + endif + !print *,'SDp ',ival1,ival2,minsd,nbitsd + endif ! end of spatial diff section + + ! Expand non-missing data values to original grid. + miss1=minval(jfld(1:nonmiss))-1 + miss2=miss1-1 + n=0 + do j=1,ndpts + if ( ifldmiss(j).eq.0 ) then + n=n+1 + ifld(j)=jfld(n) + elseif ( ifldmiss(j).eq.1 ) then + ifld(j)=miss1 + elseif ( ifldmiss(j).eq.2 ) then + ifld(j)=miss2 + endif + enddo + if(ndpts<2) simple_alg=.true. + + ! Determine Groups to be used. + if ( simple_alg ) then + ! set group length to 10 : calculate number of groups + ! and length of last group + ngroups=ndpts/10 + glen(1:ngroups)=10 + itemp=mod(ndpts,10) + if (itemp.ne.0) then + ngroups=ngroups+1 + glen(ngroups)=itemp + endif + else + ! Use Dr. Glahn's algorithm for determining grouping. + ! + kfildo=6 + minpk=10 + inc=1 + maxgrps=(ndpts/minpk)+1 + allocate(jmin(maxgrps)) + allocate(jmax(maxgrps)) + allocate(lbit(maxgrps)) + call pack_gp(kfildo,ifld,ndpts,missopt,minpk,inc,miss1,miss2, & + jmin,jmax,lbit,glen,maxgrps,ngroups,ibit,jbit, & + kbit,novref,lbitref,ier) + !print *,'SAGier = ',ier,ibit,jbit,kbit,novref,lbitref + do ng=1,ngroups + glen(ng)=glen(ng)+novref + enddo + deallocate(jmin) + deallocate(jmax) + deallocate(lbit) + endif + + ! For each group, find the group's reference value (min) + ! and the number of bits needed to hold the remaining values + n=1 + do ng=1,ngroups + ! how many of each type? + num0=count(ifldmiss(n:n+glen(ng)-1) .EQ. 0) + num1=count(ifldmiss(n:n+glen(ng)-1) .EQ. 1) + num2=count(ifldmiss(n:n+glen(ng)-1) .EQ. 2) + if ( num0.eq.0 ) then ! all missing values + if ( num1.eq.0 ) then ! all secondary missing + gref(ng)=-2 + gwidth(ng)=0 + elseif ( num2.eq.0 ) then ! all primary missing + gref(ng)=-1 + gwidth(ng)=0 + else ! both primary and secondary + gref(ng)=0 + gwidth(ng)=1 + endif + else ! contains some non-missing data + ! find max and min values of group + gref(ng)=huge(n) + imax=-1*huge(n) + j=n + do lg=1,glen(ng) + if ( ifldmiss(j).eq.0 ) then + if (ifld(j).lt.gref(ng)) gref(ng)=ifld(j) + if (ifld(j).gt.imax) imax=ifld(j) + endif + j=j+1 + enddo + if (missopt.eq.1) imax=imax+1 + if (missopt.eq.2) imax=imax+2 + ! calc num of bits needed to hold data + if ( gref(ng).ne.imax ) then + temp=i1log2(imax-gref(ng)) + gwidth(ng)=ceiling(temp) + else + gwidth(ng)=0 + endif + endif + ! Subtract min from data + j=n + mtemp=2**gwidth(ng) + do lg=1,glen(ng) + if (ifldmiss(j).eq.0) then ! non-missing + ifld(j)=ifld(j)-gref(ng) + elseif (ifldmiss(j).eq.1) then ! primary missing + ifld(j)=mtemp-1 + elseif (ifldmiss(j).eq.2) then ! secondary missing + ifld(j)=mtemp-2 + endif + j=j+1 + enddo + ! increment fld array counter + n=n+glen(ng) + enddo + + ! Find max of the group references and calc num of bits needed + ! to pack each groups reference value, then + ! pack up group reference values + !write(77,*)'GREFS: ',(gref(j),j=1,ngroups) + igmax=maxval(gref(1:ngroups)) + if (missopt.eq.1) igmax=igmax+1 + if (missopt.eq.2) igmax=igmax+2 + if (igmax.ne.0) then + temp=i1log2(igmax) + nbitsgref=ceiling(temp) + ! restet the ref values of any "missing only" groups. + mtemp=2**nbitsgref + do j=1,ngroups + if (gref(j).eq.-1) gref(j)=mtemp-1 + if (gref(j).eq.-2) gref(j)=mtemp-2 + enddo + call g2_sbytesc(cpack,gref,iofst,nbitsgref,0,ngroups) + itemp=nbitsgref*ngroups + iofst=iofst+itemp + ! Pad last octet with Zeros, if necessary, + if (mod(itemp,8).ne.0) then + left=8-mod(itemp,8) + call g2_sbytec(cpack,zero,iofst,left) + iofst=iofst+left + endif + else + nbitsgref=0 + endif + + ! Find max/min of the group widths and calc num of bits needed + ! to pack each groups width value, then + ! pack up group width values + !write(77,*)'GWIDTHS: ',(gwidth(j),j=1,ngroups) + iwmax=maxval(gwidth(1:ngroups)) + ngwidthref=minval(gwidth(1:ngroups)) + if (iwmax.ne.ngwidthref) then + temp=i1log2(iwmax-ngwidthref) + nbitsgwidth=ceiling(temp) + do i=1,ngroups + gwidth(i)=gwidth(i)-ngwidthref + enddo + call g2_sbytesc(cpack,gwidth,iofst,nbitsgwidth,0,ngroups) + itemp=nbitsgwidth*ngroups + iofst=iofst+itemp + ! Pad last octet with Zeros, if necessary, + if (mod(itemp,8).ne.0) then + left=8-mod(itemp,8) + call g2_sbytec(cpack,zero,iofst,left) + iofst=iofst+left + endif + else + nbitsgwidth=0 + gwidth(1:ngroups)=0 + endif + + ! Find max/min of the group lengths and calc num of bits needed + ! to pack each groups length value, then + ! pack up group length values + !write(77,*)'GLENS: ',(glen(j),j=1,ngroups) + ilmax=maxval(glen(1:ngroups-1)) + nglenref=minval(glen(1:ngroups-1)) + if(ngroups>0) then + nglenlast=glen(ngroups) + else + nglenlast=0 + endif + if (ilmax.ne.nglenref) then + temp=i1log2(ilmax-nglenref) + nbitsglen=ceiling(temp) + do i=1,ngroups-1 + glen(i)=glen(i)-nglenref + enddo + call g2_sbytesc(cpack,glen,iofst,nbitsglen,0,ngroups) + itemp=nbitsglen*ngroups + iofst=iofst+itemp + ! Pad last octet with Zeros, if necessary, + if (mod(itemp,8).ne.0) then + left=8-mod(itemp,8) + call g2_sbytec(cpack,zero,iofst,left) + iofst=iofst+left + endif + else + nbitsglen=0 + glen(1:ngroups)=0 + endif + + ! For each group, pack data values + !write(77,*)'IFLDS: ',(ifld(j),j=1,ndpts) + n=1 + ij=0 + do ng=1,ngroups + glength=glen(ng)+nglenref + if (ng.eq.ngroups ) glength=nglenlast + grpwidth=gwidth(ng)+ngwidthref + !write(77,*)'NGP ',ng,grpwidth,glength,gref(ng) + if ( grpwidth.ne.0 ) then + call g2_sbytesc(cpack,ifld(n),iofst,grpwidth,0,glength) + iofst=iofst+(grpwidth*glength) + endif + do kk=1,glength + ij=ij+1 + !write(77,*)'SAG ',ij,fld(ij),ifld(ij),gref(ng),bscale,rmin,dscale + enddo + n=n+glength + enddo + ! Pad last octet with Zeros, if necessary, + if (mod(iofst,8).ne.0) then + left=8-mod(iofst,8) + call g2_sbytec(cpack,zero,iofst,left) + iofst=iofst+left + endif + lcpack=iofst/8 + ! + if ( allocated(ifld) ) deallocate(ifld) + if ( allocated(jfld) ) deallocate(jfld) + if ( allocated(ifldmiss) ) deallocate(ifldmiss) + if ( allocated(gref) ) deallocate(gref) + if ( allocated(gwidth) ) deallocate(gwidth) + if ( allocated(glen) ) deallocate(glen) + + ! Fill in ref value and number of bits in Template 5.2 + rmin4 = real(rmin, 4) + call mkieee(rmin4,ref,1) ! ensure reference value is IEEE format + iref=transfer(ref,iref) + idrstmpl(1)=iref + idrstmpl(4)=nbitsgref + idrstmpl(5)=0 ! original data were reals + idrstmpl(6)=1 ! general group splitting + idrstmpl(10)=ngroups ! Number of groups + idrstmpl(11)=ngwidthref ! reference for group widths + idrstmpl(12)=nbitsgwidth ! num bits used for group widths + idrstmpl(13)=nglenref ! Reference for group lengths + idrstmpl(14)=1 ! length increment for group lengths + idrstmpl(15)=nglenlast ! True length of last group + idrstmpl(16)=nbitsglen ! num bits used for group lengths + if (idrsnum.eq.3) then + idrstmpl(18)=nbitsd/8 ! num bits used for extra spatial + ! differencing values + endif + +end subroutine misspack diff --git a/src/misspack.f b/src/misspack.f deleted file mode 100644 index c8b06de5..00000000 --- a/src/misspack.f +++ /dev/null @@ -1,502 +0,0 @@ -!> @file -!> @brief Pack up a data field using a GRIB2 algorithm with missing -!> value management. -!> @author Stephen Gilbert @date 2000-06-21 - -!> Pack up a data field using a GRIB2 algorithm with missing value -!> management. -!> -!> This subroutine packs up a data field using a complex packing -!> algorithm as defined in the GRIB2 documention. It supports GRIB2 -!> complex packing templates with or without spatial differences -!> (i.e. Data Representation Tables [5.2] -!> (https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table5-2.shtml) -!> and -!> [5.3](https://www.nco.ncep.noaa.gov/pmb/docs/grib2/grib2_doc/grib2_table5-3.shtml)). -!> -!> This subroutine also fills in GRIB2 Data Representation Template -!> 5.2 or 5.3 with the appropriate values. It assumes that Missing -!> Value Management is being used and that 1 or 2 missing values -!> appear in the data. -!> -!> @param[in] fld Contains the data values to pack. -!> @param[in] ndpts The number of data values in array fld. -!> @param[in] idrsnum Data Representation Template number. Must -!> equal 2 or 3. -!> @param[inout] idrstmpl Contains the array of values for Data -!> Representation Template 5.2 or 5.3 -!> - idrstmpl(1) Reference value - ignored on input set by compack routine. -!> - idrstmpl(2) Binary Scale Factor. -!> - idrstmpl(3) Decimal Scale Factor. -!> - idrstmpl(4) number of bits for each data value - ignored on input. -!> - idrstmpl(5) Original field type, currently ignored on input, set = 0 on -!> !output, Data values assumed to be reals. -!> - idrstmpl(6) = 0 use lossless compression or = 1 use lossy compression. -!> - idrstmpl(7) Missing value management. -!> - idrstmpl(8) Primary missing value. -!> - idrstmpl(9) Secondary missing value. -!> - idrstmpl(17) Order of Spatial Differencing (1 or 2). -!> @param[out] cpack The packed data field (character*1 array). -!> @param[out] lcpack length of packed field cpack. -1 is returned if -!> idrstmpl(7) is not set correctly. -!> -!> @author Stephen Gilbert @date 2000-06-21 - subroutine misspack(fld,ndpts,idrsnum,idrstmpl,cpack,lcpack) - - use intmath - integer,intent(in) :: ndpts,idrsnum - real,intent(in) :: fld(ndpts) - character(len=1),intent(out) :: cpack(*) - integer,intent(inout) :: idrstmpl(*) - integer,intent(out) :: lcpack - - real(4) :: ref, rmin4 - integer(4) :: iref - integer,allocatable :: ifld(:),ifldmiss(:),jfld(:) - integer,allocatable :: jmin(:),jmax(:),lbit(:) - integer,parameter :: zero=0 - integer,allocatable :: gref(:),gwidth(:),glen(:) - integer :: glength,grpwidth - logical :: simple_alg - logical :: have_rmin - - simple_alg = .false. - have_rmin = .false. - bscale=2.0**real(-idrstmpl(2)) - dscale=10.0**real(idrstmpl(3)) - missopt=idrstmpl(7) - if ( missopt.ne.1 .AND. missopt.ne.2 ) then - print *,'misspack: Unrecognized option.' - lcpack=-1 - return - else ! Get missing values - call rdieee(idrstmpl(8),rmissp,1) - if (missopt.eq.2) call rdieee(idrstmpl(9),rmisss,1) - endif - -! Find min value of non-missing values in the data, -! AND set up missing value mapping of the field. - allocate(ifldmiss(ndpts)) -c rmin=huge(rmin) - - if ( missopt .eq. 1 ) then ! Primary missing value only - do j=1,ndpts - if (fld(j).eq.rmissp) then - ifldmiss(j)=1 - else - ifldmiss(j)=0 - if(have_rmin) then - if (fld(j).lt.rmin) rmin=fld(j) - else - rmin=fld(j) - have_rmin=.true. - endif - endif - enddo - if(.not.have_rmin) rmin=rmissp - endif - if ( missopt .eq. 2 ) then ! Primary and secondary missing values - do j=1,ndpts - if (fld(j).eq.rmissp) then - ifldmiss(j)=1 - elseif (fld(j).eq.rmisss) then - ifldmiss(j)=2 - else - ifldmiss(j)=0 - if(have_rmin) then - if (fld(j).lt.rmin) rmin=fld(j) - else - rmin=fld(j) - have_rmin=.true. - endif - endif - if(.not.have_rmin) rmin=rmissp - enddo - endif - -! Allocate work arrays: -! Note: -ifldmiss(j),j=1,ndpts is a map of original field indicating -! which of the original data values -! are primary missing (1), sencondary missing (2) or non-missing (0). -! -jfld(j),j=1,nonmiss is a subarray of just the non-missing values from -! the original field. - !if (rmin.ne.rmax) then - iofst=0 - allocate(ifld(ndpts)) - allocate(jfld(ndpts)) - allocate(gref(ndpts)) - allocate(gwidth(ndpts)) - allocate(glen(ndpts)) - ! - ! Scale original data - ! - nonmiss=0 - if (idrstmpl(2).eq.0) then ! No binary scaling - imin=nint(rmin*dscale) - !imax=nint(rmax*dscale) - rmin=real(imin) - do j=1,ndpts - if (ifldmiss(j).eq.0) then - nonmiss=nonmiss+1 - jfld(nonmiss)=max(0,nint(fld(j)*dscale)-imin) - endif - enddo - else ! Use binary scaling factor - rmin=rmin*dscale - !rmax=rmax*dscale - do j=1,ndpts - if (ifldmiss(j).eq.0) then - nonmiss=nonmiss+1 - jfld(nonmiss)=max(0,nint(((fld(j)*dscale)-rmin)*bscale)) - endif - enddo - endif - ! - ! Calculate Spatial differences, if using DRS Template 5.3 - ! - if (idrsnum.eq.3) then ! spatial differences - if (idrstmpl(17).ne.1.and.idrstmpl(17).ne.2) idrstmpl(17)=2 - if (idrstmpl(17).eq.1) then ! first order - if(nonmiss<1) then - ival1=1.0 - else - ival1=jfld(1) - endif - do j=nonmiss,2,-1 - jfld(j)=jfld(j)-jfld(j-1) - enddo - if(nonmiss>0) jfld(1)=0 - elseif (idrstmpl(17).eq.2) then ! second order - if(nonmiss==1) then - ival1=jfld(1) - ival2=jfld(1) - elseif(nonmiss<1) then - ival1=1.0 - ival2=1.0 - else - ival1=jfld(1) - ival2=jfld(2) - endif - do j=nonmiss,3,-1 - jfld(j)=jfld(j)-(2*jfld(j-1))+jfld(j-2) - enddo - if(nonmiss>=1) jfld(1)=0 - if(nonmiss>=2) jfld(2)=0 - endif - - ! subtract min value from spatial diff field - isd=idrstmpl(17)+1 - minsd=minval(jfld(isd:nonmiss)) - do j=isd,nonmiss - jfld(j)=jfld(j)-minsd - enddo - - ! find num of bits need to store minsd and add 1 extra bit - ! to indicate sign - temp=i1log2(abs(minsd)) - nbitsd=ceiling(temp)+1 - - ! find num of bits need to store ifld(1) ( and ifld(2) - ! if using 2nd order differencing ) - maxorig=ival1 - if (idrstmpl(17).eq.2.and.ival2.gt.ival1) maxorig=ival2 - temp=i1log2(maxorig) - nbitorig=ceiling(temp)+1 - if (nbitorig.gt.nbitsd) nbitsd=nbitorig - ! increase number of bits to even multiple of 8 ( octet ) - if (mod(nbitsd,8).ne.0) nbitsd=nbitsd+(8-mod(nbitsd,8)) - - ! Store extra spatial differencing info into the packed - ! data section. - if (nbitsd.ne.0) then - ! pack first original value - if (ival1.ge.0) then - call g2_sbytec(cpack,ival1,iofst,nbitsd) - iofst=iofst+nbitsd - else - call g2_sbytec(cpack,1,iofst,1) - iofst=iofst+1 - call g2_sbytec(cpack,iabs(ival1),iofst,nbitsd-1) - iofst=iofst+nbitsd-1 - endif - if (idrstmpl(17).eq.2) then - ! pack second original value - if (ival2.ge.0) then - call g2_sbytec(cpack,ival2,iofst,nbitsd) - iofst=iofst+nbitsd - else - call g2_sbytec(cpack,1,iofst,1) - iofst=iofst+1 - call g2_sbytec(cpack,iabs(ival2),iofst,nbitsd-1) - iofst=iofst+nbitsd-1 - endif - endif - ! pack overall min of spatial differences - if (minsd.ge.0) then - call g2_sbytec(cpack,minsd,iofst,nbitsd) - iofst=iofst+nbitsd - else - call g2_sbytec(cpack,1,iofst,1) - iofst=iofst+1 - call g2_sbytec(cpack,iabs(minsd),iofst,nbitsd-1) - iofst=iofst+nbitsd-1 - endif - endif - !print *,'SDp ',ival1,ival2,minsd,nbitsd - endif ! end of spatial diff section - - ! Expand non-missing data values to original grid. - miss1=minval(jfld(1:nonmiss))-1 - miss2=miss1-1 - n=0 - do j=1,ndpts - if ( ifldmiss(j).eq.0 ) then - n=n+1 - ifld(j)=jfld(n) - elseif ( ifldmiss(j).eq.1 ) then - ifld(j)=miss1 - elseif ( ifldmiss(j).eq.2 ) then - ifld(j)=miss2 - endif - enddo - if(ndpts<2) simple_alg=.true. - - ! Determine Groups to be used. - if ( simple_alg ) then - ! set group length to 10 : calculate number of groups - ! and length of last group - ngroups=ndpts/10 - glen(1:ngroups)=10 - itemp=mod(ndpts,10) - if (itemp.ne.0) then - ngroups=ngroups+1 - glen(ngroups)=itemp - endif - else - ! Use Dr. Glahn's algorithm for determining grouping. - ! - kfildo=6 - minpk=10 - inc=1 - maxgrps=(ndpts/minpk)+1 - allocate(jmin(maxgrps)) - allocate(jmax(maxgrps)) - allocate(lbit(maxgrps)) - call pack_gp(kfildo,ifld,ndpts,missopt,minpk,inc,miss1,miss2, - & jmin,jmax,lbit,glen,maxgrps,ngroups,ibit,jbit, - & kbit,novref,lbitref,ier) - !print *,'SAGier = ',ier,ibit,jbit,kbit,novref,lbitref - do ng=1,ngroups - glen(ng)=glen(ng)+novref - enddo - deallocate(jmin) - deallocate(jmax) - deallocate(lbit) - endif - - ! For each group, find the group's reference value (min) - ! and the number of bits needed to hold the remaining values - n=1 - do ng=1,ngroups - ! how many of each type? - num0=count(ifldmiss(n:n+glen(ng)-1) .EQ. 0) - num1=count(ifldmiss(n:n+glen(ng)-1) .EQ. 1) - num2=count(ifldmiss(n:n+glen(ng)-1) .EQ. 2) - if ( num0.eq.0 ) then ! all missing values - if ( num1.eq.0 ) then ! all secondary missing - gref(ng)=-2 - gwidth(ng)=0 - elseif ( num2.eq.0 ) then ! all primary missing - gref(ng)=-1 - gwidth(ng)=0 - else ! both primary and secondary - gref(ng)=0 - gwidth(ng)=1 - endif - else ! contains some non-missing data - ! find max and min values of group - gref(ng)=huge(n) - imax=-1*huge(n) - j=n - do lg=1,glen(ng) - if ( ifldmiss(j).eq.0 ) then - if (ifld(j).lt.gref(ng)) gref(ng)=ifld(j) - if (ifld(j).gt.imax) imax=ifld(j) - endif - j=j+1 - enddo - if (missopt.eq.1) imax=imax+1 - if (missopt.eq.2) imax=imax+2 - ! calc num of bits needed to hold data - if ( gref(ng).ne.imax ) then - temp=i1log2(imax-gref(ng)) - gwidth(ng)=ceiling(temp) - else - gwidth(ng)=0 - endif - endif - ! Subtract min from data - j=n - mtemp=2**gwidth(ng) - do lg=1,glen(ng) - if (ifldmiss(j).eq.0) then ! non-missing - ifld(j)=ifld(j)-gref(ng) - elseif (ifldmiss(j).eq.1) then ! primary missing - ifld(j)=mtemp-1 - elseif (ifldmiss(j).eq.2) then ! secondary missing - ifld(j)=mtemp-2 - endif - j=j+1 - enddo - ! increment fld array counter - n=n+glen(ng) - enddo - - ! Find max of the group references and calc num of bits needed - ! to pack each groups reference value, then - ! pack up group reference values - !write(77,*)'GREFS: ',(gref(j),j=1,ngroups) - igmax=maxval(gref(1:ngroups)) - if (missopt.eq.1) igmax=igmax+1 - if (missopt.eq.2) igmax=igmax+2 - if (igmax.ne.0) then - temp=i1log2(igmax) - nbitsgref=ceiling(temp) - ! restet the ref values of any "missing only" groups. - mtemp=2**nbitsgref - do j=1,ngroups - if (gref(j).eq.-1) gref(j)=mtemp-1 - if (gref(j).eq.-2) gref(j)=mtemp-2 - enddo - call g2_sbytesc(cpack,gref,iofst,nbitsgref,0,ngroups) - itemp=nbitsgref*ngroups - iofst=iofst+itemp - ! Pad last octet with Zeros, if necessary, - if (mod(itemp,8).ne.0) then - left=8-mod(itemp,8) - call g2_sbytec(cpack,zero,iofst,left) - iofst=iofst+left - endif - else - nbitsgref=0 - endif - - ! Find max/min of the group widths and calc num of bits needed - ! to pack each groups width value, then - ! pack up group width values - !write(77,*)'GWIDTHS: ',(gwidth(j),j=1,ngroups) - iwmax=maxval(gwidth(1:ngroups)) - ngwidthref=minval(gwidth(1:ngroups)) - if (iwmax.ne.ngwidthref) then - temp=i1log2(iwmax-ngwidthref) - nbitsgwidth=ceiling(temp) - do i=1,ngroups - gwidth(i)=gwidth(i)-ngwidthref - enddo - call g2_sbytesc(cpack,gwidth,iofst,nbitsgwidth,0,ngroups) - itemp=nbitsgwidth*ngroups - iofst=iofst+itemp - ! Pad last octet with Zeros, if necessary, - if (mod(itemp,8).ne.0) then - left=8-mod(itemp,8) - call g2_sbytec(cpack,zero,iofst,left) - iofst=iofst+left - endif - else - nbitsgwidth=0 - gwidth(1:ngroups)=0 - endif - - ! Find max/min of the group lengths and calc num of bits needed - ! to pack each groups length value, then - ! pack up group length values - !write(77,*)'GLENS: ',(glen(j),j=1,ngroups) - ilmax=maxval(glen(1:ngroups-1)) - nglenref=minval(glen(1:ngroups-1)) - if(ngroups>0) then - nglenlast=glen(ngroups) - else - nglenlast=0 - endif - if (ilmax.ne.nglenref) then - temp=i1log2(ilmax-nglenref) - nbitsglen=ceiling(temp) - do i=1,ngroups-1 - glen(i)=glen(i)-nglenref - enddo - call g2_sbytesc(cpack,glen,iofst,nbitsglen,0,ngroups) - itemp=nbitsglen*ngroups - iofst=iofst+itemp - ! Pad last octet with Zeros, if necessary, - if (mod(itemp,8).ne.0) then - left=8-mod(itemp,8) - call g2_sbytec(cpack,zero,iofst,left) - iofst=iofst+left - endif - else - nbitsglen=0 - glen(1:ngroups)=0 - endif - - ! For each group, pack data values - !write(77,*)'IFLDS: ',(ifld(j),j=1,ndpts) - n=1 - ij=0 - do ng=1,ngroups - glength=glen(ng)+nglenref - if (ng.eq.ngroups ) glength=nglenlast - grpwidth=gwidth(ng)+ngwidthref - !write(77,*)'NGP ',ng,grpwidth,glength,gref(ng) - if ( grpwidth.ne.0 ) then - call g2_sbytesc(cpack,ifld(n),iofst,grpwidth,0,glength) - iofst=iofst+(grpwidth*glength) - endif - do kk=1,glength - ij=ij+1 - !write(77,*)'SAG ',ij,fld(ij),ifld(ij),gref(ng),bscale,rmin,dscale - enddo - n=n+glength - enddo - ! Pad last octet with Zeros, if necessary, - if (mod(iofst,8).ne.0) then - left=8-mod(iofst,8) - call g2_sbytec(cpack,zero,iofst,left) - iofst=iofst+left - endif - lcpack=iofst/8 - ! - if ( allocated(ifld) ) deallocate(ifld) - if ( allocated(jfld) ) deallocate(jfld) - if ( allocated(ifldmiss) ) deallocate(ifldmiss) - if ( allocated(gref) ) deallocate(gref) - if ( allocated(gwidth) ) deallocate(gwidth) - if ( allocated(glen) ) deallocate(glen) - !else ! Constant field ( max = min ) - ! nbits=0 - ! lcpack=0 - ! nbitsgref=0 - ! ngroups=0 - !endif - -! Fill in ref value and number of bits in Template 5.2 - 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)=nbitsgref - idrstmpl(5)=0 ! original data were reals - idrstmpl(6)=1 ! general group splitting - idrstmpl(10)=ngroups ! Number of groups - idrstmpl(11)=ngwidthref ! reference for group widths - idrstmpl(12)=nbitsgwidth ! num bits used for group widths - idrstmpl(13)=nglenref ! Reference for group lengths - idrstmpl(14)=1 ! length increment for group lengths - idrstmpl(15)=nglenlast ! True length of last group - idrstmpl(16)=nbitsglen ! num bits used for group lengths - if (idrsnum.eq.3) then - idrstmpl(18)=nbitsd/8 ! num bits used for extra spatial - ! differencing values - endif - - end diff --git a/src/realloc.F90 b/src/realloc.F90 new file mode 100644 index 00000000..335eae5d --- /dev/null +++ b/src/realloc.F90 @@ -0,0 +1,158 @@ +!> @file +!> @brief Reallocate memory, preserving contents. +!> @author Stephen Gilbert @date 2000-10-01 + +!> @brief Reallocate memory, preserving contents. +!> +!> This module contains two subroutines to reallocate integer, and +!> character arrays in memory, preseving the existing contents of the +!> array. +!> +!> @author Stephen Gilbert @date 2000-10-01 +module re_alloc + + interface realloc + module procedure realloc_c1 + module procedure realloc_r + module procedure realloc_i + end interface realloc + +contains + + !> This subroutine reallocates a character array, preserving its + !> contents. + !> + !> @param[inout] c pointer for data in memory. + !> @param[in] n dimension for data in memory. This is how much data + !> is currently stored in the array. + !> @param[in] m dimension for allocatable array. This is the size + !> the array will be after the allocation of memory. + !> @param[out] istat scalar INTEGER variable for allocate. + !> - 0 No error. + !> - 10 Incorrect dimension. + !> - other Allocation error. + !> + !> @author Stephen Gilbert @date 2000-10-01 + subroutine realloc_c1(c,n,m,istat) + character(len=1),pointer,dimension(:) :: c + integer,intent(in) :: n,m + integer,intent(out) :: istat + integer :: num + character(len=1),pointer,dimension(:) :: tmp + + istat=0 + if ( (n<0) .OR. (m<=0) ) then + istat=10 + return + endif + + if ( .not. associated(c) ) then + allocate(c(m),stat=istat) ! allocate new memory + return + endif + + tmp=>c ! save pointer to original mem + nullify(c) + allocate(c(m),stat=istat) ! allocate new memory + if ( istat /= 0 ) then + c=>tmp + return + endif + if ( n /= 0 ) then + num=min(n,m) + c(1:num)=tmp(1:num) ! copy data from orig mem to new loc. + endif + deallocate(tmp) ! deallocate original memory + return + end subroutine realloc_c1 + + !> This subroutine reallocates an integer array, preserving its + !> contents. + !> + !> @param[inout] c pointer for data in memory. + !> @param[in] n dimension for data in memory. This is how much data + !> is currently stored in the array. + !> @param[in] m dimension for allocatable array. This is the size + !> the array will be after the allocation of memory. + !> @param[out] istat scalar INTEGER variable for allocate. + !> - 0 No error. + !> - 10 Incorrect dimension. + !> - other Allocation error. + !> + !> @author Stephen Gilbert @date 2000-10-01 + subroutine realloc_r(c,n,m,istat) + real,pointer,dimension(:) :: c + integer,intent(in) :: n,m + integer,intent(out) :: istat + integer :: num + real,pointer,dimension(:) :: tmp + + istat=0 + if ( (n<0) .OR. (m<=0) ) then + istat=10 + return + endif + + if ( .not. associated(c) ) then + allocate(c(m),stat=istat) ! allocate new memory + return + endif + + tmp=>c ! save pointer to original mem + nullify(c) + allocate(c(m),stat=istat) ! allocate new memory + if ( istat /= 0 ) then + c=>tmp + return + endif + if ( n /= 0 ) then + num=min(n,m) + c(1:num)=tmp(1:num) ! copy data from orig mem to new loc. + endif + deallocate(tmp) ! deallocate original memory + return + end subroutine realloc_r + + !> This subroutine reorganize integer type data in memory into + !> one one dimensional array. + !> @param[inout] c pointer for data in memory. + !> @param[in] n dimension for data in memory. + !> @param[in] m dimension for allocatable array. + !> @param[out] istat scalar INTEGER variable for allocate. + !> @author Stephen Gilbert @date 2000-10-01 + !> + + subroutine realloc_i(c,n,m,istat) + integer,pointer,dimension(:) :: c + integer,intent(in) :: n,m + integer,intent(out) :: istat + integer :: num + integer,pointer,dimension(:) :: tmp + + istat=0 + if ( (n<0) .OR. (m<=0) ) then + istat=10 + return + endif + + if ( .not. associated(c) ) then + allocate(c(m),stat=istat) ! allocate new memory + return + endif + + tmp=>c ! save pointer to original mem + nullify(c) + allocate(c(m),stat=istat) ! allocate new memory + if ( istat /= 0 ) then + c=>tmp + return + endif + if ( n /= 0 ) then + num=min(n,m) + c(1:num)=tmp(1:num) ! copy data from orig mem to new loc. + endif + deallocate(tmp) ! deallocate original memory + return + end subroutine realloc_i + +end module re_alloc diff --git a/src/realloc.f b/src/realloc.f deleted file mode 100644 index 5090abbb..00000000 --- a/src/realloc.f +++ /dev/null @@ -1,158 +0,0 @@ -!> @file -!> @brief Reallocate memory, preserving contents. -!> @author Stephen Gilbert @date 2000-10-01 - -!> @brief Reallocate memory, preserving contents. -!> -!> This module contains two subroutines to reallocate integer, and -!> character arrays in memory, preseving the existing contents of the -!> array. -!> -!> @author Stephen Gilbert @date 2000-10-01 - module re_alloc - - interface realloc - module procedure realloc_c1 - module procedure realloc_r - module procedure realloc_i - end interface - - contains - - !> This subroutine reallocates a character array, preserving its - !> contents. - !> - !> @param[inout] c pointer for data in memory. - !> @param[in] n dimension for data in memory. This is how much data - !> is currently stored in the array. - !> @param[in] m dimension for allocatable array. This is the size - !> the array will be after the allocation of memory. - !> @param[out] istat scalar INTEGER variable for allocate. - !> - 0 No error. - !> - 10 Incorrect dimension. - !> - other Allocation error. - !> - !> @author Stephen Gilbert @date 2000-10-01 - subroutine realloc_c1(c,n,m,istat) - character(len=1),pointer,dimension(:) :: c - integer,intent(in) :: n,m - integer,intent(out) :: istat - integer :: num - character(len=1),pointer,dimension(:) :: tmp - - istat=0 - if ( (n<0) .OR. (m<=0) ) then - istat=10 - return - endif - - if ( .not. associated(c) ) then - allocate(c(m),stat=istat) ! allocate new memory - return - endif - - tmp=>c ! save pointer to original mem - nullify(c) - allocate(c(m),stat=istat) ! allocate new memory - if ( istat /= 0 ) then - c=>tmp - return - endif - if ( n /= 0 ) then - num=min(n,m) - c(1:num)=tmp(1:num) ! copy data from orig mem to new loc. - endif - deallocate(tmp) ! deallocate original memory - return - end subroutine - - !> This subroutine reallocates an integer array, preserving its - !> contents. - !> - !> @param[inout] c pointer for data in memory. - !> @param[in] n dimension for data in memory. This is how much data - !> is currently stored in the array. - !> @param[in] m dimension for allocatable array. This is the size - !> the array will be after the allocation of memory. - !> @param[out] istat scalar INTEGER variable for allocate. - !> - 0 No error. - !> - 10 Incorrect dimension. - !> - other Allocation error. - !> - !> @author Stephen Gilbert @date 2000-10-01 - subroutine realloc_r(c,n,m,istat) - real,pointer,dimension(:) :: c - integer,intent(in) :: n,m - integer,intent(out) :: istat - integer :: num - real,pointer,dimension(:) :: tmp - - istat=0 - if ( (n<0) .OR. (m<=0) ) then - istat=10 - return - endif - - if ( .not. associated(c) ) then - allocate(c(m),stat=istat) ! allocate new memory - return - endif - - tmp=>c ! save pointer to original mem - nullify(c) - allocate(c(m),stat=istat) ! allocate new memory - if ( istat /= 0 ) then - c=>tmp - return - endif - if ( n /= 0 ) then - num=min(n,m) - c(1:num)=tmp(1:num) ! copy data from orig mem to new loc. - endif - deallocate(tmp) ! deallocate original memory - return - end subroutine - -!> This subroutine reorganize integer type data in memory into -!> one one dimensional array. -!> @param[inout] c pointer for data in memory. -!> @param[in] n dimension for data in memory. -!> @param[in] m dimension for allocatable array. -!> @param[out] istat scalar INTEGER variable for allocate. -!> @author Stephen Gilbert @date 2000-10-01 -!> - - subroutine realloc_i(c,n,m,istat) - integer,pointer,dimension(:) :: c - integer,intent(in) :: n,m - integer,intent(out) :: istat - integer :: num - integer,pointer,dimension(:) :: tmp - - istat=0 - if ( (n<0) .OR. (m<=0) ) then - istat=10 - return - endif - - if ( .not. associated(c) ) then - allocate(c(m),stat=istat) ! allocate new memory - return - endif - - tmp=>c ! save pointer to original mem - nullify(c) - allocate(c(m),stat=istat) ! allocate new memory - if ( istat /= 0 ) then - c=>tmp - return - endif - if ( n /= 0 ) then - num=min(n,m) - c(1:num)=tmp(1:num) ! copy data from orig mem to new loc. - endif - deallocate(tmp) ! deallocate original memory - return - end subroutine - - end module re_alloc diff --git a/src/simpack.F90 b/src/simpack.F90 new file mode 100644 index 00000000..f09610a8 --- /dev/null +++ b/src/simpack.F90 @@ -0,0 +1,154 @@ +!> @file +!> @brief Pack a data field using simple packing algorithm. +!> @author Stephen Gilbert @date 2000-06-21 + +!> Pack a data field using a simple packing algorithm. +!> +!> This subroutine also fills in GRIB2 Data Representation Template 5.0 +!> with the appropriate values. +!> +!> @param[in] fld Contains the data values to pack. +!> @param[in] ndpts The number of data values in array fld. +!> @param[inout] idrstmpl Contains the array of values for Data +!> Representation Template 5.2 or 5.3. +!> - (1) Reference value - ignored on input +!> - (2) Binary Scale Factor +!> - (3) Decimal Scale Factor +!> - (4) Number of bits used to pack data, if value is > 0 and <= +!> 31. If this input value is 0 or outside above range then the num +!> of bits is calculated based on given data and scale factors. +!> - (5) Original field type - currently ignored on input Data values +!> assumed to be reals. +!> @param[out] cpack The packed data field (character*1 array). +!> @param[out] lcpack length of packed field cpack. +!> +!> @author Stephen Gilbert @date 2000-06-21 +subroutine simpack(fld,ndpts,idrstmpl,cpack,lcpack) + use intmath + implicit none + + integer,intent(in) :: ndpts + real,intent(in) :: fld(ndpts) + character(len=1),intent(out) :: cpack(*) + integer,intent(inout) :: idrstmpl(*) + integer,intent(out) :: lcpack + + real(4) :: ref,rmin4 + + integer(4) :: iref + integer :: ifld(ndpts) + integer,parameter :: zero=0 + integer :: nbittot, nbits, maxnum, maxdif, left + integer :: imax, imin, j + real :: rmax, rmin, temp, dscale, bscale + + bscale=2.0**real(-idrstmpl(2)) + dscale=10.0**real(idrstmpl(3)) + if (idrstmpl(4).le.0.OR.idrstmpl(4).gt.31) then + nbits=0 + else + nbits=idrstmpl(4) + endif + + ! Find max and min values in the data + if(ndpts<1) then + rmin=0 + rmax=0 + else + rmax=fld(1) + rmin=fld(1) + do j=2,ndpts + if (fld(j).gt.rmax) rmax=fld(j) + if (fld(j).lt.rmin) rmin=fld(j) + enddo + 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) then + + ! Determine which algorithm to use based on user-supplied + ! binary scale factor and number of bits. + if (nbits.eq.0.AND.idrstmpl(2).eq.0) then + + ! No binary scaling and calculate minumum number of + ! bits in which the data will fit. + imin=nint(rmin*dscale) + imax=nint(rmax*dscale) + maxdif=imax-imin + temp=i1log2(maxdif+1) + nbits=ceiling(temp) + rmin=real(imin) + ! scale data + do j=1,ndpts + ifld(j)=nint(fld(j)*dscale)-imin + enddo + elseif (nbits.ne.0.AND.idrstmpl(2).eq.0) then + + ! Use minimum number of bits specified by user and + ! adjust binary scaling factor to accomodate data. + rmin=rmin*dscale + rmax=rmax*dscale + maxnum=(2**nbits)-1 + temp=ilog2(nint(real(maxnum)/(rmax-rmin))) + idrstmpl(2)=ceiling(-1.0*temp) + bscale=2.0**real(-idrstmpl(2)) + ! scale data + do j=1,ndpts + ifld(j)=max(0,nint(((fld(j)*dscale)-rmin)*bscale)) + enddo + elseif (nbits.eq.0.AND.idrstmpl(2).ne.0) then + + ! Use binary scaling factor and calculate minumum number of + ! bits in which the data will fit. + rmin=rmin*dscale + rmax=rmax*dscale + maxdif=nint((rmax-rmin)*bscale) + temp=i1log2(maxdif) + nbits=ceiling(temp) + ! scale data + do j=1,ndpts + ifld(j)=max(0,nint(((fld(j)*dscale)-rmin)*bscale)) + enddo + elseif (nbits.ne.0.AND.idrstmpl(2).ne.0) then + + ! Use binary scaling factor and use minumum number of + ! bits specified by user. Dangerous - may loose + ! information if binary scale factor and nbits not set + ! properly by user. + rmin=rmin*dscale + ! scale data + do j=1,ndpts + ifld(j)=max(0,nint(((fld(j)*dscale)-rmin)*bscale)) + enddo + endif + + ! Pack data, Pad last octet with Zeros, if necessary, + ! and calculate the length of the packed data in bytes + call g2_sbytesc(cpack,ifld,0,nbits,0,ndpts) + nbittot=nbits*ndpts + left=8-mod(nbittot,8) + if (left.ne.8) then + call g2_sbytec(cpack,zero,nbittot,left) ! Pad with zeros to fill Octet + nbittot=nbittot+left + endif + lcpack=nbittot/8 + + else + !print *,'nbits 0' + nbits=0 + lcpack=0 + endif + + + ! 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 + iref=transfer(ref,iref) + idrstmpl(1)=iref + idrstmpl(4)=nbits + idrstmpl(5)=0 ! original data were reals + +end subroutine simpack diff --git a/src/simpack.f b/src/simpack.f deleted file mode 100644 index 3d575980..00000000 --- a/src/simpack.f +++ /dev/null @@ -1,154 +0,0 @@ -!> @file -!> @brief Pack a data field using simple packing algorithm. -!> @author Stephen Gilbert @date 2000-06-21 - -!> Pack a data field using a simple packing algorithm. -!> -!> This subroutine also fills in GRIB2 Data Representation Template 5.0 -!> with the appropriate values. -!> -!> @param[in] fld Contains the data values to pack. -!> @param[in] ndpts The number of data values in array fld. -!> @param[inout] idrstmpl Contains the array of values for Data -!> Representation Template 5.2 or 5.3. -!> - (1) Reference value - ignored on input -!> - (2) Binary Scale Factor -!> - (3) Decimal Scale Factor -!> - (4) Number of bits used to pack data, if value is > 0 and <= -!> 31. If this input value is 0 or outside above range then the num -!> of bits is calculated based on given data and scale factors. -!> - (5) Original field type - currently ignored on input Data values -!> assumed to be reals. -!> @param[out] cpack The packed data field (character*1 array). -!> @param[out] lcpack length of packed field cpack. -!> -!> @author Stephen Gilbert @date 2000-06-21 - subroutine simpack(fld,ndpts,idrstmpl,cpack,lcpack) - use intmath - implicit none - - integer,intent(in) :: ndpts - real,intent(in) :: fld(ndpts) - character(len=1),intent(out) :: cpack(*) - integer,intent(inout) :: idrstmpl(*) - integer,intent(out) :: lcpack - - real(4) :: ref,rmin4 - - integer(4) :: iref - integer :: ifld(ndpts) - integer,parameter :: zero=0 - integer :: nbittot, nbits, maxnum, maxdif, left - integer :: imax, imin, j - real :: rmax, rmin, temp, dscale, bscale - - bscale=2.0**real(-idrstmpl(2)) - dscale=10.0**real(idrstmpl(3)) - if (idrstmpl(4).le.0.OR.idrstmpl(4).gt.31) then - nbits=0 - else - nbits=idrstmpl(4) - endif - -! Find max and min values in the data - if(ndpts<1) then - rmin=0 - rmax=0 - else - rmax=fld(1) - rmin=fld(1) - do j=2,ndpts - if (fld(j).gt.rmax) rmax=fld(j) - if (fld(j).lt.rmin) rmin=fld(j) - enddo - 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) then - - ! Determine which algorithm to use based on user-supplied - ! binary scale factor and number of bits. - if (nbits.eq.0.AND.idrstmpl(2).eq.0) then - - ! No binary scaling and calculate minumum number of - ! bits in which the data will fit. - imin=nint(rmin*dscale) - imax=nint(rmax*dscale) - maxdif=imax-imin - temp=i1log2(maxdif+1) - nbits=ceiling(temp) - rmin=real(imin) - ! scale data - do j=1,ndpts - ifld(j)=nint(fld(j)*dscale)-imin - enddo - elseif (nbits.ne.0.AND.idrstmpl(2).eq.0) then - - ! Use minimum number of bits specified by user and - ! adjust binary scaling factor to accomodate data. - rmin=rmin*dscale - rmax=rmax*dscale - maxnum=(2**nbits)-1 - temp=ilog2(nint(real(maxnum)/(rmax-rmin))) - idrstmpl(2)=ceiling(-1.0*temp) - bscale=2.0**real(-idrstmpl(2)) - ! scale data - do j=1,ndpts - ifld(j)=max(0,nint(((fld(j)*dscale)-rmin)*bscale)) - enddo - elseif (nbits.eq.0.AND.idrstmpl(2).ne.0) then - - ! Use binary scaling factor and calculate minumum number of - ! bits in which the data will fit. - rmin=rmin*dscale - rmax=rmax*dscale - maxdif=nint((rmax-rmin)*bscale) - temp=i1log2(maxdif) - nbits=ceiling(temp) - ! scale data - do j=1,ndpts - ifld(j)=max(0,nint(((fld(j)*dscale)-rmin)*bscale)) - enddo - elseif (nbits.ne.0.AND.idrstmpl(2).ne.0) then - - ! Use binary scaling factor and use minumum number of - ! bits specified by user. Dangerous - may loose - ! information if binary scale factor and nbits not set - ! properly by user. - rmin=rmin*dscale - ! scale data - do j=1,ndpts - ifld(j)=max(0,nint(((fld(j)*dscale)-rmin)*bscale)) - enddo - endif - - ! Pack data, Pad last octet with Zeros, if necessary, - ! and calculate the length of the packed data in bytes - call g2_sbytesc(cpack,ifld,0,nbits,0,ndpts) - nbittot=nbits*ndpts - left=8-mod(nbittot,8) - if (left.ne.8) then - call g2_sbytec(cpack,zero,nbittot,left) ! Pad with zeros to fill Octet - nbittot=nbittot+left - endif - lcpack=nbittot/8 - - else - !print *,'nbits 0' - nbits=0 - lcpack=0 - endif - - -! 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 - iref=transfer(ref,iref) - idrstmpl(1)=iref - idrstmpl(4)=nbits - idrstmpl(5)=0 ! original data were reals - - end diff --git a/tests/test_intmath.F90 b/tests/test_intmath.F90 index 04ef599a..9bb2874a 100644 --- a/tests/test_intmath.F90 +++ b/tests/test_intmath.F90 @@ -35,5 +35,39 @@ program test_intmath print *, 'ok' print *, 'SUCCESS!' + + ! Here is some legacy test code that could be recovered... +! program test_intmath +! use intmath +! implicit none +! real(kind=16) :: temp +! real(kind=16), parameter :: alog2=log(2.0_16) +! integer(kind=8), parameter :: & +! & one=1,big=Z'7FFFFFFFFFFFFFFF',small=-2000000_8, & +! & check=Z'1FFFFFFF' +! integer(kind=8) :: ival, iret +! !$OMP PARALLEL DO PRIVATE(ival,temp,iret) +! do ival=small,big +! 10 format(Z16,' -- MISMATCH: ',I0,'=>',I0,' (',I0,' = ',F0.10,')') +! 20 format(Z16,' -- OKAY: ',I0,'=>',I0,' (',I0,' = ',F0.10,')') +! if(ival+one