addfield.f

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

      !   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 sbyte(cgrib,ipdstmpl(i),iofst,nbits)
        else
          call sbyte(cgrib,one,iofst,1)
          call sbyte(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
        call mkieee(coordlist,coordieee,numcoord)
        call sbytes(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 sbyte(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(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
      else 
         ndpts=ngrdpts;
         pfld=>fld;
      endif
      lcpack=0
      allocate(cpack(ndpts*4),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 gbyte(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
#ifdef USE_JPEG2000
      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
        call jpcpack(pfld,width,height,idrstmpl,cpack,lcpack)
#endif  /* USE_JPEG2000 */
#ifdef USE_PNG
      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
        call pngpack(pfld,width,height,idrstmpl,cpack,lcpack)
#endif   /* USE_PNG */
      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 sbyte(cgrib,five,iofst,8)     ! Store section number ( 5 )
      iofst=iofst+8
      call sbyte(cgrib,ndpts,iofst,32)    ! Store num of actual data points
      iofst=iofst+32
      call sbyte(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 sbyte(cgrib,idrstmpl(i),iofst,nbits)
        else
          call sbyte(cgrib,one,iofst,1)
          call sbyte(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 sbyte(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 sbyte(cgrib,six,iofst,8)     ! Store section number ( 6 )
      iofst=iofst+8
      call sbyte(cgrib,ibmap,iofst,8)    ! Store Bit Map indicator
      iofst=iofst+8
      !
      !  Store bitmap, if supplied
      !
      if (ibmap.eq.0) then
        call sbytes(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 sbyte(cgrib,zero,iofst,left)     ! Pad with zeros to fill Octet
        iofst=iofst+left
      endif
      lensec6=(iofst-ibeg)/8
      call sbyte(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 sbyte(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 sbyte(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 sbyte(cgrib,newlen,96,32)

      return
      end