gribinfo.f

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

      subroutine gribinfo(cgrib,lcgrib,listsec0,listsec1,
     &                    numlocal,numfields,maxvals,ierr)
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .                                       .
! SUBPROGRAM:    gribinfo 
!   PRGMMR: Gilbert         ORG: W/NP11    DATE: 2000-05-25
!
! ABSTRACT: This subroutine searches through a GRIB2 message and
!   returns the number of Local Use Sections and number of gridded
!   fields found in the message.  It also performs various checks 
!   to see if the message is a valid GRIB2 message.
!   Last, a list of safe array dimensions is returned for use in 
!   allocating return arrays from routines getlocal, gettemplates, and 
!   getfields.  (See maxvals and REMARKS)
!
! PROGRAM HISTORY LOG:
! 2000-05-25  Gilbert
!
! USAGE:    CALL gribinfo(cgrib,lcgrib,listsec0,listsec1,
!     &                    numlocal,numfields,ierr)
!   INPUT ARGUMENT LIST:
!     cgrib    - Character array that contains the GRIB2 message
!     lcgrib   - Length (in bytes) of GRIB message in array cgrib.
!
!   OUTPUT ARGUMENT LIST:      
!     listsec0 - Contains information decoded from GRIB Indicator Section 0.
!                Must be dimensioned >= 2.
!                listsec0(1)=Discipline-GRIB Master Table Number
!                            (see Code Table 0.0)
!                listsec0(2)=GRIB Edition Number (currently 2)
!                listsec0(3)=Length of GRIB message
!     listsec1 - Contains information read from GRIB Identification Section 1.
!                Must be dimensioned >= 13.
!                listsec1(1)=Id of orginating centre (Common Code Table C-1)
!                listsec1(2)=Id of orginating sub-centre (local table)
!                listsec1(3)=GRIB Master Tables Version Number (Code Table 1.0)
!                listsec1(4)=GRIB Local Tables Version Number 
!                listsec1(5)=Significance of Reference Time (Code Table 1.1)
!                listsec1(6)=Reference Time - Year (4 digits)
!                listsec1(7)=Reference Time - Month
!                listsec1(8)=Reference Time - Day
!                listsec1(9)=Reference Time - Hour
!                listsec1(10)=Reference Time - Minute
!                listsec1(11)=Reference Time - Second
!                listsec1(12)=Production status of data (Code Table 1.2)
!                listsec1(13)=Type of processed data (Code Table 1.3)
!     numlocal - The number of Local Use Sections ( Section 2 ) found in 
!                the GRIB message.
!     numfields- The number of gridded fieldse found in the GRIB message.
!     maxvals()- The maximum number of elements that could be returned
!                in various arrays from this GRIB2 message. (see REMARKS)
!                maxvals(1)=max length of local section 2 (for getlocal)
!                maxvals(2)=max length of GDS Template (for gettemplates 
!                                                       and getfield)
!                maxvals(3)=max length of GDS Optional list (for getfield)
!                maxvals(4)=max length of PDS Template (for gettemplates 
!                                                       and getfield)
!                maxvals(5)=max length of PDS Optional list (for getfield)
!                maxvals(6)=max length of DRS Template (for gettemplates 
!                                                       and getfield)
!                maxvals(7)=max number of gridpoints (for getfield)
!     ierr     - Error return code.
!                0 = no error
!                1 = Beginning characters "GRIB" not found.
!                2 = GRIB message is not Edition 2.
!                3 = Could not find Section 1, where expected.
!                4 = End string "7777" found, but not where expected.
!                5 = End string "7777" not found at end of message.
!
! REMARKS: Array maxvals contains the maximum possible 
!          number of values that will be returned in argument arrays
!          for routines getlocal, gettemplates, and getfields.  
!          Users can use this info to determine if their arrays are 
!          dimensioned large enough for the data that may be returned
!          from the above routines, or to dynamically allocate arrays
!          with a reasonable size.
!          NOTE that the actual number of values in these arrays is returned
!          from the routines and will likely be less than the values 
!          calculated by this routine.
!
! ATTRIBUTES:
!   LANGUAGE: Fortran 90
!   MACHINE:  IBM SP
!
!$$$

      character(len=1),intent(in) :: cgrib(lcgrib)
      integer,intent(in) :: lcgrib
      integer,intent(out) :: listsec0(3),listsec1(13),maxvals(7)
      integer,intent(out) :: numlocal,numfields,ierr
      
      character(len=4),parameter :: grib='GRIB',c7777='7777'
      character(len=4) :: ctemp
      integer,parameter :: zero=0,one=1
      integer,parameter :: mapsec1len=13
      integer,parameter :: 
     &        mapsec1(mapsec1len)=(/ 2,2,1,1,1,2,1,1,1,1,1,1,1 /)
      integer iofst,ibeg,istart

      ierr=0
      numlocal=0
      numfields=0
      maxsec2len=1
      maxgdstmpl=1
      maxdeflist=1
      maxpdstmpl=1
      maxcoordlist=1
      maxdrstmpl=1
      maxgridpts=0
!
!  Check for beginning of GRIB message in the first 100 bytes
!
      istart=0
      do j=1,100
        ctemp=cgrib(j)//cgrib(j+1)//cgrib(j+2)//cgrib(j+3)
        if (ctemp.eq.grib ) then
          istart=j
          exit
        endif
      enddo
      if (istart.eq.0) then
        print *,'gribinfo:  Beginning characters GRIB not found.'
        ierr=1
        return
      endif
!
!  Unpack Section 0 - Indicator Section 
!
      iofst=8*(istart+5)
      call gbyte(cgrib,listsec0(1),iofst,8)     ! Discipline
      iofst=iofst+8
      call gbyte(cgrib,listsec0(2),iofst,8)     ! GRIB edition number
      iofst=iofst+8
      iofst=iofst+32
      call gbyte(cgrib,lengrib,iofst,32)        ! Length of GRIB message
      iofst=iofst+32
      listsec0(3)=lengrib
      lensec0=16
      ipos=istart+lensec0
!
!  Currently handles only GRIB Edition 2.
!  
      if (listsec0(2).ne.2) then
        print *,'gribinfo: can only decode GRIB edition 2.'
        ierr=2
        return
      endif
!
!  Unpack Section 1 - Identification Section
!
      call gbyte(cgrib,lensec1,iofst,32)        ! Length of Section 1
      iofst=iofst+32
      call gbyte(cgrib,isecnum,iofst,8)         ! Section number ( 1 )
      iofst=iofst+8
      if (isecnum.ne.1) then
        print *,'gribinfo: Could not find section 1.'
        ierr=3
        return
      endif
      !
      !   Unpack each input value in array listsec1 into the
      !   the appropriate number of octets, which are specified in
      !   corresponding entries in array mapsec1.
      !
      do i=1,mapsec1len
        nbits=mapsec1(i)*8
        call gbyte(cgrib,listsec1(i),iofst,nbits)
        iofst=iofst+nbits
      enddo
      ipos=ipos+lensec1
!
!  Loop through the remaining sections keeping track of the 
!  length of each.  Also count the number of times Section 2
!  and Section 4 appear.
!
      do
        ctemp=cgrib(ipos)//cgrib(ipos+1)//cgrib(ipos+2)//cgrib(ipos+3)
        if (ctemp.eq.c7777 ) then
          ipos=ipos+4
          if (ipos.ne.(istart+lengrib)) then
            print *,'gribinfo: "7777" found, but not where expected.'
            ierr=4
            return
          endif
          exit
        endif
        iofst=(ipos-1)*8
        call gbyte(cgrib,lensec,iofst,32)        ! Get Length of Section
        iofst=iofst+32
        call gbyte(cgrib,isecnum,iofst,8)         ! Get Section number
        iofst=iofst+8
        ipos=ipos+lensec                 ! Update beginning of section pointer
        if (ipos.gt.(istart+lengrib)) then
          print *,'gribinfo: "7777"  not found at end of GRIB message.'
          ierr=5
          return
        endif
        if (isecnum.eq.2) then     ! Local Section 2
           !   increment counter for total number of local sections found
           !   and determine largest Section 2 in message
           numlocal=numlocal+1
           lenposs=lensec-5
           if ( lenposs.gt.maxsec2len ) maxsec2len=lenposs
        elseif (isecnum.eq.3) then
           iofst=iofst+8                      ! skip source of grid def.
           call gbyte(cgrib,ngdpts,iofst,32)         ! Get Num of Grid Points
           iofst=iofst+32
           call gbyte(cgrib,nbyte,iofst,8)      ! Get Num octets for opt. list
           iofst=iofst+8
           if (ngdpts.gt.maxgridpts) maxgridpts=ngdpts
           lenposs=lensec-14
           if ( lenposs.gt.maxgdstmpl ) maxgdstmpl=lenposs
           if (nbyte.ne.0) then
              lenposs=lenposs/nbyte
              if ( lenposs.gt.maxdeflist ) maxdeflist=lenposs
           endif
        elseif (isecnum.eq.4) then
          numfields=numfields+1
           call gbyte(cgrib,numcoord,iofst,16)      ! Get Num of Coord Values
           iofst=iofst+16
           if (numcoord.ne.0) then
              if (numcoord.gt.maxcoordlist) maxcoordlist=numcoord
           endif
           lenposs=lensec-9
           if ( lenposs.gt.maxpdstmpl ) maxpdstmpl=lenposs
        elseif (isecnum.eq.5) then
           lenposs=lensec-11
           if ( lenposs.gt.maxdrstmpl ) maxdrstmpl=lenposs
        endif
        
      enddo

      maxvals(1)=maxsec2len
      maxvals(2)=maxgdstmpl
      maxvals(3)=maxdeflist
      maxvals(4)=maxpdstmpl
      maxvals(5)=maxcoordlist
      maxvals(6)=maxdrstmpl
      maxvals(7)=maxgridpts

      return
      end