getfield.f

计算线性趋势 回归系数 主要用于气象站点值的线性趋势计算

      real,intent(out) :: coordlist(*)
      integer,intent(out) :: ipdsnum,ipdstmpl(*)
      integer,intent(out) :: ierr,numcoord

      real(4),allocatable :: coordieee(:)
      integer,allocatable :: mappds(:)
      integer :: mappdslen
      logical needext

      ierr=0

      call gbyte(cgrib,lensec,iofst,32)        ! Get Length of Section
      iofst=iofst+32
      iofst=iofst+8     ! skip section number
      allocate(mappds(lensec))

      call gbyte(cgrib,numcoord,iofst,16)    ! Get num of coordinate values
      iofst=iofst+16
      call gbyte(cgrib,ipdsnum,iofst,16)    ! Get 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
      !
      !   Unpack each value into array ipdstmpl from 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 ) then
          call gbyte(cgrib,ipdstmpl(i),iofst,nbits)
        else
          call gbyte(cgrib,isign,iofst,1)
          call gbyte(cgrib,ipdstmpl(i),iofst+1,nbits-1)
          if (isign.eq.1) ipdstmpl(i)=-ipdstmpl(i)
        endif
        iofst=iofst+nbits
      enddo
      !
      !   Check to see if the Product Definition Template needs to be
      !   extended.
      !   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,newmappdslen,mappds)
        !   Unpack the rest of the Product Definition Template
        do i=mappdslen+1,newmappdslen
          nbits=iabs(mappds(i))*8
          if ( mappds(i).ge.0 ) then
            call gbyte(cgrib,ipdstmpl(i),iofst,nbits)
          else
            call gbyte(cgrib,isign,iofst,1)
            call gbyte(cgrib,ipdstmpl(i),iofst+1,nbits-1)
            if (isign.eq.1) ipdstmpl(i)=-ipdstmpl(i)
          endif
          iofst=iofst+nbits
        enddo
        mappdslen=newmappdslen
      endif
      !
      !   Get Optional list of vertical coordinate values
      !   after the Product Definition Template, if necessary.
      !
      if ( numcoord .ne. 0 ) then
        allocate (coordieee(numcoord))
        call gbytes(cgrib,coordieee,iofst,32,0,numcoord)
        call rdieee(coordieee,coordlist,numcoord)
        deallocate (coordieee)
        iofst=iofst+(32*numcoord)
      endif
      if( allocated(mappds) ) deallocate(mappds)
      return    ! End of Section 4 processing
      end


      subroutine unpack5(cgrib,lcgrib,iofst,ndpts,idrsnum,idrstmpl,
     &                   mapdrslen,ierr)
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .                                       .
! SUBPROGRAM:    unpack5 
!   PRGMMR: Gilbert         ORG: W/NP11    DATE: 2000-05-26
!
! ABSTRACT: This subroutine unpacks Section 5 (Data Representation Section)
!   starting at octet 6 of that Section.  
!
! PROGRAM HISTORY LOG:
! 2000-05-26  Gilbert
!
! USAGE:    CALL unpack5(cgrib,lcgrib,iofst,ndpts,idrsnum,idrstmpl,
!                        mapdrslen,ierr)
!   INPUT ARGUMENT LIST:
!     cgrib    - Character array that contains the GRIB2 message
!     lcgrib   - Length (in bytes) of GRIB message array cgrib.
!     iofst    - Bit offset of the beginning of Section 5.
!
!   OUTPUT ARGUMENT LIST:      
!     iofst    - Bit offset at the end of Section 5, returned.
!     ndpts    - Number of data points unpacked and returned.
!     idrsnum  - Data Representation Template Number ( see Code Table 5.0)
!     idrstmpl - Contains the data values for the specified Data Representation
!                Template ( N=idrsnum ).  Each element of this integer
!                array contains an entry (in the order specified) of Data
!                Representation Template 5.N
!     mapdrslen- Number of elements in idrstmpl().  i.e. number of entries
!                in Data Representation Template 5.N  ( N=idrsnum ).
!     ierr     - Error return code.
!                0 = no error
!                7 = "GRIB" message contains an undefined Data
!                    Representation Template.
!
! REMARKS: None
!
! ATTRIBUTES:
!   LANGUAGE: Fortran 90
!   MACHINE:  IBM SP
!
!$$$

      use drstemplates

      character(len=1),intent(in) :: cgrib(lcgrib)
      integer,intent(in) :: lcgrib
      integer,intent(inout) :: iofst
      integer,intent(out) :: ndpts,idrsnum,idrstmpl(*)
      integer,intent(out) :: ierr

C      integer,allocatable :: mapdrs(:)
      integer,allocatable :: mapdrs(:)
      integer :: mapdrslen
      logical needext

      ierr=0

      call gbyte(cgrib,lensec,iofst,32)        ! Get Length of Section
      iofst=iofst+32
      iofst=iofst+8     ! skip section number
      allocate(mapdrs(lensec))

      call gbyte(cgrib,ndpts,iofst,32)    ! Get num of data points
      iofst=iofst+32
      call gbyte(cgrib,idrsnum,iofst,16)     ! Get Data Rep Template Num.
      iofst=iofst+16
      !   Gen Data Representation Template
      call getdrstemplate(idrsnum,mapdrslen,mapdrs,needext,iret)
      if (iret.ne.0) then
        ierr=7
        return
      endif
      !
      !   Unpack each value into array ipdstmpl from the
      !   the appropriate number of octets, which are specified in
      !   corresponding entries in array mappds.
      !
      do i=1,mapdrslen
        nbits=iabs(mapdrs(i))*8
        if ( mapdrs(i).ge.0 ) then
          call gbyte(cgrib,idrstmpl(i),iofst,nbits)
        else
          call gbyte(cgrib,isign,iofst,1)
          call gbyte(cgrib,idrstmpl(i),iofst+1,nbits-1)
          if (isign.eq.1) idrstmpl(i)=-idrstmpl(i)
        endif
        iofst=iofst+nbits
      enddo
      !
      !   Check to see if the Data Representation Template needs to be
      !   extended.
      !   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 extdrstemplate(idrsnum,idrstmpl,newmapdrslen,mapdrs)
        !   Unpack the rest of the Data Representation Template
        do i=mapdrslen+1,newmapdrslen
          nbits=iabs(mapdrs(i))*8
          if ( mapdrs(i).ge.0 ) then
            call gbyte(cgrib,idrstmpl(i),iofst,nbits)
          else
            call gbyte(cgrib,isign,iofst,1)
            call gbyte(cgrib,idrstmpl(i),iofst+1,nbits-1)
            if (isign.eq.1) idrstmpl(i)=-idrstmpl(i)
          endif
          iofst=iofst+nbits
        enddo
        mapdrslen=newmapdrslen
      endif
      if( allocated(mapdrs) ) deallocate(mapdrs)
      return    ! End of Section 5 processing
      end


      subroutine unpack6(cgrib,lcgrib,iofst,ngpts,ibmap,bmap,ierr)
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .                                       .
! SUBPROGRAM:    unpack6 
!   PRGMMR: Gilbert         ORG: W/NP11    DATE: 2000-05-26
!
! ABSTRACT: This subroutine unpacks Section 6 (Bit-Map Section)
!   starting at octet 6 of that Section.  
!
! PROGRAM HISTORY LOG:
! 2000-05-26  Gilbert
!
! USAGE:    CALL unpack6(cgrib,lcgrib,iofst,ngpts,ibmap,bmap,ierr)
!   INPUT ARGUMENT LIST:
!     cgrib    - Character array that contains the GRIB2 message
!     lcgrib   - Length (in bytes) of GRIB message array cgrib.
!     iofst    - Bit offset of the beginning of Section 6.
!     ngpts    - Number of grid points specified in the bit-map
!
!   OUTPUT ARGUMENT LIST:      
!     iofst    - Bit offset at the end of Section 6, returned.
!     ibmap    - Bitmap indicator ( see Code Table 6.0 )
!                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.
!     bmap()   - Logical*1 array containing decoded bitmap. ( if ibmap=0 )
!     ierr     - Error return code.
!                0 = no error
!                4 = Unrecognized pre-defined bit-map.
!
! REMARKS: None
!
! ATTRIBUTES:
!   LANGUAGE: Fortran 90
!   MACHINE:  IBM SP
!
!$$$

      character(len=1),intent(in) :: cgrib(lcgrib)
      integer,intent(in) :: lcgrib,ngpts
      integer,intent(inout) :: iofst
      integer,intent(out) :: ibmap
      integer,intent(out) :: ierr
      logical*1,intent(out) :: bmap(ngpts)

      integer :: intbmap(ngpts)

      ierr=0

      iofst=iofst+32    ! skip Length of Section
      iofst=iofst+8     ! skip section number

      call gbyte(cgrib,ibmap,iofst,8)    ! Get bit-map indicator
      iofst=iofst+8

      if (ibmap.eq.0) then               ! Unpack bitmap
        call gbytes(cgrib,intbmap,iofst,1,0,ngpts)
        iofst=iofst+ngpts
        do j=1,ngpts
          bmap(j)=.true.
          if (intbmap(j).eq.0) bmap(j)=.false.
        enddo
      elseif (ibmap.eq.254) then               ! Use previous bitmap
        return
      elseif (ibmap.eq.255) then               ! No bitmap in message
        bmap(1:ngpts)=.true.
      else
        print *,'unpack6: Predefined bitmap ',ibmap,' not recognized.'
        ierr=4
      endif
      
      return    ! End of Section 6 processing
      end