misspack.f

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

      subroutine misspack(fld,ndpts,idrsnum,idrstmpl,cpack,lcpack)
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .                                       .
! SUBPROGRAM:    misspack
!   PRGMMR: Gilbert          ORG: W/NP11    DATE: 2000-06-21
!
! ABSTRACT: 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. DRTs 5.2 and 5.3).
!   It also fills in GRIB2 Data Representation Template 5.2 or 5.3 
!   with the appropriate values.
!   This version assumes that Missing Value Management is being used and that
!   1 or 2 missing values appear in the data.
!
! PROGRAM HISTORY LOG:
! 2000-06-21  Gilbert
! 2004-12-29  Gilbert  -  Corrected bug when encoding secondary missing values.
!
! USAGE:    CALL misspack(fld,ndpts,idrsnum,idrstmpl,cpack,lcpack)
!   INPUT ARGUMENT LIST:
!     fld()    - Contains the data values to pack
!     ndpts    - The number of data values in array fld()
!     idrsnum  - Data Representation Template number 5.N
!                Must equal 2 or 3.
!     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
!                    .
!                    .
!                (7) = Missing value management
!                (8) = Primary missing value
!                (9) = Secondary missing value
!                    .
!                    .
!               (17) = Order of Spatial Differencing  ( 1 or 2 )
!                    .
!                    .
!
!   OUTPUT ARGUMENT LIST: 
!     idrstmpl - Contains the array of values for Data Representation
!                Template 5.3
!                (1) = Reference value - set by misspack routine.
!                (2) = Binary Scale Factor - unchanged from input
!                (3) = Decimal Scale Factor - unchanged from input
!                    .
!                    .
!     cpack    - The packed data field (character*1 array)
!     lcpack   - length of packed field cpack().
!
! REMARKS: None
!
! ATTRIBUTES:
!   LANGUAGE: XL Fortran 90
!   MACHINE:  IBM SP
!
!$$$

      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
      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 = .false.
      
      alog2=alog(2.0)
      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))
      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 (fld(j).lt.rmin) rmin=fld(j)
           endif
         enddo
      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 (fld(j).lt.rmin) rmin=fld(j)
           endif
         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)=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)=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
              ival1=jfld(1)
              do j=nonmiss,2,-1
                 jfld(j)=jfld(j)-jfld(j-1)
              enddo
              jfld(1)=0
           elseif (idrstmpl(17).eq.2) then      ! second order
              ival1=jfld(1)
              ival2=jfld(2)
              do j=nonmiss,3,-1
                 jfld(j)=jfld(j)-(2*jfld(j-1))+jfld(j-2)
              enddo
              jfld(1)=0
              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=alog(real(abs(minsd)+1))/alog2
           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=alog(real(maxorig+1))/alog2
           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 sbyte(cpack,ival1,iofst,nbitsd)
                 iofst=iofst+nbitsd
              else
                 call sbyte(cpack,1,iofst,1)
                 iofst=iofst+1
                 call sbyte(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 sbyte(cpack,ival2,iofst,nbitsd)
                    iofst=iofst+nbitsd
                 else
                    call sbyte(cpack,1,iofst,1)
                    iofst=iofst+1
                    call sbyte(cpack,iabs(ival2),iofst,nbitsd-1)
                    iofst=iofst+nbitsd-1
                 endif
              endif
              !  pack overall min of spatial differences
              if (minsd.ge.0) then
                 call sbyte(cpack,minsd,iofst,nbitsd)
                 iofst=iofst+nbitsd
              else
                 call sbyte(cpack,1,iofst,1)
                 iofst=iofst+1
                 call sbyte(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