gf_getfld.f

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

            character(len=1),intent(in) :: cgrib(lcgrib)
            integer,intent(in) :: lcgrib,ndpts,idrsnum,igdsnum
            integer,intent(inout) :: iofst
            integer,pointer,dimension(:) :: idrstmpl,igdstmpl
            integer,intent(out) :: ierr
            real,pointer,dimension(:) :: fld
         end subroutine gf_unpack7
      end interface

      have3=.false.
      have4=.false.
      have5=.false.
      have6=.false.
      have7=.false.
      ierr=0
      numfld=0
      gfld%locallen=0
      nullify(gfld%list_opt,gfld%igdtmpl,gfld%ipdtmpl)
      nullify(gfld%coord_list,gfld%idrtmpl,gfld%bmap,gfld%fld)
!
!  Check for valid request number
!  
      if (ifldnum.le.0) then
        print *,'gf_getfld: Request for field number must be positive.'
        ierr=3
        return
      endif
!
!  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 *,'gf_getfld:  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
      lensec0=16
      ipos=istart+lensec0
!
!  Currently handles only GRIB Edition 2.
!  
      if (listsec0(2).ne.2) then
        print *,'gf_getfld: can only decode GRIB edition 2.'
        ierr=2
        return
      endif
!
!  Loop through the remaining sections keeping track of the 
!  length of each.  Also keep the latest Grid Definition Section info.
!  Unpack the requested field number.
!
      do
        !    Check to see if we are at end of GRIB message
        ctemp=cgrib(ipos)//cgrib(ipos+1)//cgrib(ipos+2)//cgrib(ipos+3)
        if (ctemp.eq.c7777 ) then
          ipos=ipos+4
          !    If end of GRIB message not where expected, issue error
          if (ipos.ne.(istart+lengrib)) then
            print *,'gf_getfld: "7777" found, but not where expected.'
            ierr=4
            return
          endif
          exit
        endif
        !     Get length of Section and Section number
        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
        !print *,' lensec= ',lensec,'    secnum= ',isecnum
        !
        !  Check to see if section number is valid
        !
        if ( (isecnum.lt.1).OR.(isecnum.gt.7) ) then
          print *,'gf_getfld: Unrecognized Section Encountered=',isecnum     
          ierr=8
          return
        endif
        !
        !   If found Section 1, decode elements in Identification Section
        !
        if (isecnum.eq.1) then
          iofst=iofst-40       ! reset offset to beginning of section
          call gf_unpack1(cgrib,lcgrib,iofst,gfld%idsect,
     &                    gfld%idsectlen,jerr)
          if (jerr.ne.0) then
            ierr=15
            return
          endif
        endif
        !
        !   If found Section 2, Grab local section
        !   Save in case this is the latest one before the requested field.
        !
        if (isecnum.eq.2) then
          iofst=iofst-40       ! reset offset to beginning of section
          if (associated(gfld%local)) deallocate(gfld%local)
          call gf_unpack2(cgrib,lcgrib,iofst,gfld%locallen,
     &                    gfld%local,jerr)
          if (jerr.ne.0) then
            ierr=16
            return
          endif
        endif
        !
        !   If found Section 3, unpack the GDS info using the 
        !   appropriate template.  Save in case this is the latest
        !   grid before the requested field.
        !
        if (isecnum.eq.3) then
          iofst=iofst-40       ! reset offset to beginning of section
          if (associated(gfld%igdtmpl)) deallocate(gfld%igdtmpl)
          if (associated(gfld%list_opt)) deallocate(gfld%list_opt)
          call gf_unpack3(cgrib,lcgrib,iofst,igds,gfld%igdtmpl,
     &                 gfld%igdtlen,gfld%list_opt,gfld%num_opt,jerr)
          if (jerr.eq.0) then
            have3=.true.
            gfld%griddef=igds(1)
            gfld%ngrdpts=igds(2)
            gfld%numoct_opt=igds(3)
            gfld%interp_opt=igds(4)
            gfld%igdtnum=igds(5)
          else
            ierr=10
            return
          endif
        endif
        !
        !   If found Section 4, check to see if this field is the
        !   one requested.
        !
        if (isecnum.eq.4) then
          numfld=numfld+1
          if (numfld.eq.ifldnum) then
            gfld%discipline=listsec0(1)
            gfld%version=listsec0(2)
            gfld%ifldnum=ifldnum
            gfld%unpacked=unpack
            gfld%expanded=.false.
            iofst=iofst-40       ! reset offset to beginning of section
            call gf_unpack4(cgrib,lcgrib,iofst,gfld%ipdtnum,
     &                      gfld%ipdtmpl,gfld%ipdtlen,gfld%coord_list,
     &                      gfld%num_coord,jerr)
            if (jerr.eq.0) then
              have4=.true.
            else
              ierr=11
              return
            endif
          endif
        endif
        !
        !   If found Section 5, check to see if this field is the
        !   one requested.
        !
        if ((isecnum.eq.5).and.(numfld.eq.ifldnum)) then
          iofst=iofst-40       ! reset offset to beginning of section
          call gf_unpack5(cgrib,lcgrib,iofst,gfld%ndpts,gfld%idrtnum,
     &                    gfld%idrtmpl,gfld%idrtlen,jerr)
          if (jerr.eq.0) then
            have5=.true.
          else
            ierr=12
            return
          endif
        endif
        !
        !   If found Section 6, Unpack bitmap.
        !   Save in case this is the latest
        !   bitmap before the requested field.
        !
        if (isecnum.eq.6) then
          if (unpack) then   ! unpack bitmap
            iofst=iofst-40       ! reset offset to beginning of section
            bmpsave=>gfld%bmap      ! save pointer to previous bitmap
            call gf_unpack6(cgrib,lcgrib,iofst,gfld%ngrdpts,gfld%ibmap,
     &                   gfld%bmap,jerr)
            if (jerr.eq.0) then
              have6=.true.
              if (gfld%ibmap .eq. 254) then    ! use previously specified bitmap
                 if ( associated(bmpsave) ) then
                    gfld%bmap=>bmpsave
                 else
                    print *,'gf_getfld:  Previous bit-map specified,',
     &                       ' but none exists,'
                    ierr=17
                    return
                 endif
              else                             ! get rid of it
                 if ( associated(bmpsave) ) deallocate(bmpsave)
              endif
            else
              ierr=13
              return
            endif
          else    ! do not unpack bitmap
            call gbyte(cgrib,gfld%ibmap,iofst,8)      ! Get BitMap Indicator
            have6=.true.
          endif
        endif
        !
        !   If found Section 7, check to see if this field is the
        !   one requested.
        !
        if ((isecnum.eq.7).and.(numfld.eq.ifldnum).and.unpack) then
          iofst=iofst-40       ! reset offset to beginning of section
          call gf_unpack7(cgrib,lcgrib,iofst,gfld%igdtnum,
     &                    gfld%igdtmpl,gfld%idrtnum,
     &                    gfld%idrtmpl,gfld%ndpts,
     &                    gfld%fld,jerr)
          if (jerr.eq.0) then
            have7=.true.
            !  If bitmap is used with this field, expand data field
            !  to grid, if possible.
            if ( gfld%ibmap .ne. 255 .AND. associated(gfld%bmap) ) then
               if ( expand ) then
                  allocate(newfld(gfld%ngrdpts))
                  !newfld(1:gfld%ngrdpts)=0.0
                  !newfld=unpack(gfld%fld,gfld%bmap,newfld)
                  n=1
                  do j=1,gfld%ngrdpts
                      if ( gfld%bmap(j) ) then
                        newfld(j)=gfld%fld(n)
                        n=n+1
                      else
                        newfld(j)=0.0
                      endif
                  enddo
                  deallocate(gfld%fld);
                  gfld%fld=>newfld;
                  gfld%expanded=.true.
               else 
                  gfld%expanded=.false.
               endif
            else 
               gfld%expanded=.true.
            endif
          else
            print *,'gf_getfld: return from gf_unpack7 = ',jerr
            ierr=14
            return
          endif
        endif
        !
        !   Check to see if we read pass the end of the GRIB
        !   message and missed the terminator string '7777'.
        !
        ipos=ipos+lensec                 ! Update beginning of section pointer
        if (ipos.gt.(istart+lengrib)) then
          print *,'gf_getfld: "7777"  not found at end of GRIB message.'
          ierr=7
          return
        endif
        !
        !  If unpacking requested, return when all sections have been
        !  processed
        !
        if (unpack.and.have3.and.have4.and.have5.and.have6.and.have7)
     &      return
        !
        !  If unpacking is not requested, return when sections 
        !  3 through 6 have been processed
        !
        if ((.NOT.unpack).and.have3.and.have4.and.have5.and.have6)
     &      return
        
      enddo

!
!  If exited from above loop, the end of the GRIB message was reached
!  before the requested field was found.
!
      print *,'gf_getfld: GRIB message contained ',numlocal,
     &        ' different fields.'
      print *,'gf_getfld: The request was for the ',ifldnum,
     &        ' field.'
      ierr=6

      return
      end