misspack.f

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

              ifld(j)=miss1
           elseif ( ifldmiss(j).eq.2 ) then
              ifld(j)=miss2
           endif
        enddo
        !
        !   Determine Groups to be used.
        !
        if ( simple_alg ) then
           !  set group length to 10 :  calculate number of groups
           !  and length of last group
           ngroups=ndpts/10
           glen(1:ngroups)=10
           itemp=mod(ndpts,10)
           if (itemp.ne.0) then
              ngroups=ngroups+1
              glen(ngroups)=itemp
           endif
        else
           ! Use Dr. Glahn's algorithm for determining grouping.
           !
           kfildo=6
           minpk=10 
           inc=1
           maxgrps=(ndpts/minpk)+1
           allocate(jmin(maxgrps))
           allocate(jmax(maxgrps))
           allocate(lbit(maxgrps))
           call pack_gp(kfildo,ifld,ndpts,missopt,minpk,inc,miss1,miss2,
     &                  jmin,jmax,lbit,glen,maxgrps,ngroups,ibit,jbit,
     &                  kbit,novref,lbitref,ier)
           !print *,'SAGier = ',ier,ibit,jbit,kbit,novref,lbitref
           do ng=1,ngroups
              glen(ng)=glen(ng)+novref
           enddo
           deallocate(jmin)
           deallocate(jmax)
           deallocate(lbit)
        endif
        !  
        !  For each group, find the group's reference value (min)
        !  and the number of bits needed to hold the remaining values
        !
        n=1
        do ng=1,ngroups
           !  how many of each type?
           num0=count(ifldmiss(n:n+glen(ng)-1) .EQ. 0)
           num1=count(ifldmiss(n:n+glen(ng)-1) .EQ. 1)
           num2=count(ifldmiss(n:n+glen(ng)-1) .EQ. 2)
           if ( num0.eq.0 ) then      ! all missing values
              if ( num1.eq.0 ) then       ! all secondary missing
                gref(ng)=-2
                gwidth(ng)=0
              elseif ( num2.eq.0 ) then       ! all primary missing
                gref(ng)=-1
                gwidth(ng)=0
              else                           ! both primary and secondary
                gref(ng)=0
                gwidth(ng)=1
              endif
           else                       ! contains some non-missing data
             !    find max and min values of group
             gref(ng)=huge(n)
             imax=-1*huge(n)
             j=n
             do lg=1,glen(ng)
                if ( ifldmiss(j).eq.0 ) then
                  if (ifld(j).lt.gref(ng)) gref(ng)=ifld(j) 
                  if (ifld(j).gt.imax) imax=ifld(j) 
                endif
                j=j+1
             enddo
             if (missopt.eq.1) imax=imax+1
             if (missopt.eq.2) imax=imax+2
             !   calc num of bits needed to hold data
             if ( gref(ng).ne.imax ) then
                temp=alog(real(imax-gref(ng)+1))/alog2
                gwidth(ng)=ceiling(temp)
             else
                gwidth(ng)=0
             endif
           endif
           !   Subtract min from data
           j=n
           mtemp=2**gwidth(ng)
           do lg=1,glen(ng)
              if (ifldmiss(j).eq.0) then       ! non-missing
                 ifld(j)=ifld(j)-gref(ng)
              elseif (ifldmiss(j).eq.1) then    ! primary missing
                 ifld(j)=mtemp-1
              elseif (ifldmiss(j).eq.2) then    ! secondary missing
                 ifld(j)=mtemp-2
              endif
              j=j+1
           enddo
           !   increment fld array counter
           n=n+glen(ng)
        enddo
        !  
        !  Find max of the group references and calc num of bits needed 
        !  to pack each groups reference value, then
        !  pack up group reference values
        !
        !write(77,*)'GREFS: ',(gref(j),j=1,ngroups)
        igmax=maxval(gref(1:ngroups))
        if (missopt.eq.1) igmax=igmax+1
        if (missopt.eq.2) igmax=igmax+2
        if (igmax.ne.0) then
           temp=alog(real(igmax+1))/alog2
           nbitsgref=ceiling(temp)
           ! restet the ref values of any "missing only" groups.
           mtemp=2**nbitsgref
           do j=1,ngroups
               if (gref(j).eq.-1) gref(j)=mtemp-1
               if (gref(j).eq.-2) gref(j)=mtemp-2
           enddo
           call sbytes(cpack,gref,iofst,nbitsgref,0,ngroups)
           itemp=nbitsgref*ngroups
           iofst=iofst+itemp
           !         Pad last octet with Zeros, if necessary,
           if (mod(itemp,8).ne.0) then
              left=8-mod(itemp,8)
              call sbyte(cpack,zero,iofst,left)
              iofst=iofst+left
           endif
        else
           nbitsgref=0
        endif
        !
        !  Find max/min of the group widths and calc num of bits needed
        !  to pack each groups width value, then
        !  pack up group width values
        !
        !write(77,*)'GWIDTHS: ',(gwidth(j),j=1,ngroups)
        iwmax=maxval(gwidth(1:ngroups))
        ngwidthref=minval(gwidth(1:ngroups))
        if (iwmax.ne.ngwidthref) then
           temp=alog(real(iwmax-ngwidthref+1))/alog2
           nbitsgwidth=ceiling(temp)
           do i=1,ngroups
              gwidth(i)=gwidth(i)-ngwidthref
           enddo
           call sbytes(cpack,gwidth,iofst,nbitsgwidth,0,ngroups)
           itemp=nbitsgwidth*ngroups
           iofst=iofst+itemp
           !         Pad last octet with Zeros, if necessary,
           if (mod(itemp,8).ne.0) then
              left=8-mod(itemp,8)
              call sbyte(cpack,zero,iofst,left)
              iofst=iofst+left
           endif
        else
           nbitsgwidth=0
           gwidth(1:ngroups)=0
        endif
        !
        !  Find max/min of the group lengths and calc num of bits needed
        !  to pack each groups length value, then
        !  pack up group length values
        !
        !write(77,*)'GLENS: ',(glen(j),j=1,ngroups)
        ilmax=maxval(glen(1:ngroups-1))
        nglenref=minval(glen(1:ngroups-1))
        nglenlast=glen(ngroups)
        if (ilmax.ne.nglenref) then
           temp=alog(real(ilmax-nglenref+1))/alog2
           nbitsglen=ceiling(temp)
           do i=1,ngroups-1
              glen(i)=glen(i)-nglenref
           enddo
           call sbytes(cpack,glen,iofst,nbitsglen,0,ngroups)
           itemp=nbitsglen*ngroups
           iofst=iofst+itemp
           !         Pad last octet with Zeros, if necessary,
           if (mod(itemp,8).ne.0) then
              left=8-mod(itemp,8)
              call sbyte(cpack,zero,iofst,left)
              iofst=iofst+left
           endif
        else
           nbitsglen=0
           glen(1:ngroups)=0
        endif
        !
        !  For each group, pack data values
        !
        !write(77,*)'IFLDS: ',(ifld(j),j=1,ndpts)
        n=1
        ij=0
        do ng=1,ngroups
           glength=glen(ng)+nglenref
           if (ng.eq.ngroups ) glength=nglenlast
           grpwidth=gwidth(ng)+ngwidthref
       !write(77,*)'NGP ',ng,grpwidth,glength,gref(ng)
           if ( grpwidth.ne.0 ) then
              call sbytes(cpack,ifld(n),iofst,grpwidth,0,glength)
              iofst=iofst+(grpwidth*glength)
           endif
           do kk=1,glength
              ij=ij+1
        !write(77,*)'SAG ',ij,fld(ij),ifld(ij),gref(ng),bscale,rmin,dscale
           enddo
           n=n+glength
        enddo
        !         Pad last octet with Zeros, if necessary,
        if (mod(iofst,8).ne.0) then
           left=8-mod(iofst,8)
           call sbyte(cpack,zero,iofst,left)
           iofst=iofst+left
        endif
        lcpack=iofst/8
        !
        if ( allocated(ifld) ) deallocate(ifld)
        if ( allocated(jfld) ) deallocate(jfld)
        if ( allocated(ifldmiss) ) deallocate(ifldmiss)
        if ( allocated(gref) ) deallocate(gref)
        if ( allocated(gwidth) ) deallocate(gwidth)
        if ( allocated(glen) ) deallocate(glen)
      !else           !   Constant field ( max = min )
      !  nbits=0
      !  lcpack=0
      !  nbitsgref=0
      !  ngroups=0
      !endif

!
!  Fill in ref value and number of bits in Template 5.2
!
      call mkieee(rmin,ref,1)   ! ensure reference value is IEEE format
!      call gbyte(ref,idrstmpl(1),0,32)
      iref=transfer(ref,iref)
      idrstmpl(1)=iref
      idrstmpl(4)=nbitsgref
      idrstmpl(5)=0         ! original data were reals
      idrstmpl(6)=1         ! general group splitting
      idrstmpl(10)=ngroups          ! Number of groups
      idrstmpl(11)=ngwidthref       ! reference for group widths
      idrstmpl(12)=nbitsgwidth      ! num bits used for group widths
      idrstmpl(13)=nglenref         ! Reference for group lengths
      idrstmpl(14)=1                ! length increment for group lengths
      idrstmpl(15)=nglenlast        ! True length of last group
      idrstmpl(16)=nbitsglen        ! num bits used for group lengths
      if (idrsnum.eq.3) then
         idrstmpl(18)=nbitsd/8      ! num bits used for extra spatial
                                    ! differencing values
      endif

      return
      end