subroutine boxintegral01(f,n1,n2,bmiss,fout,fw,h,kdim)
C revised 06/05/02 see  4th note
C revised 06/01/01
C revised 12/10/00
C JIANT, 04/12/00
C
C Program computes vertical integral of standard depth (SD) level data.
C Data at each SD level weighted with half thickness of the layer defined
C as depth difference between previous and next SD levels except
C first (n1) and last (n2) depth data have a weight of half distance
C to the corresponding nearest SD level within n1,n2 range.
C
C NOTES:
C(1) If value at 1st depth is missing then output is BMISS
C
C(2) Integral will be evaluated from 1st depth to last depth with data
C    which could be less than h(n2)
C
C(3) If value at 2nd depth is missing then
C    output are fw=0.5*(h(n1+1)-h(n1)) and fout=f(n1)*fw
C
C(4) For n1=n2: (06/05/02) fout=f(n1)*0.5*(h(n1+1)-h(n1-1))
C
      dimension f(kdim),h(kdim)

      nsum=0
      fout=bmiss
      fw=bmiss
      hmax=h(kdim)+(h(kdim)-h(kdim-1))

      fval=f(n1)
      if(fval .gt. bmiss)then
         if ( n1 .lt. kdim ) then
          wh=0.5*(h(n1+1)-h(n1))
         else
          wh=0.5*(hmax-h(n1))
         endif
         sum=wh*fval
         wsum=wh
         nsum=nsum+1
      else
         fout=bmiss
         fw=bmiss
         return
      endif

      if(n1.eq.n2)then
       if(n1.gt.1)then
         if ( n1 .lt. kdim ) then
          fw=0.5*(h(n1+1)-h(n1-1))
         else
          fw=0.5*(hmax-h(n1-1))
         endif
         fout=fval*fw
         return
       else
         fout=sum
         fw=wsum
         return
       endif
      endif

      do k=n1+1,n2-1
         if ( k .le. kdim ) then
          fval=f(k)
         else
          fval=f(kdim)
         endif
         if(fval .gt. bmiss)then
          if ( k .lt. kdim ) then
	   wh= 0.5*(h(k+1)-h(k-1)) 
          else
	   wh= 0.5*(hmax-h(k-1)) 
          endif
	  sum=sum + wh*fval 
          wsum=wsum + wh
          nsum=nsum+1
         elseif(nsum.eq.1)then
          fout=sum
          fw=wsum
          return
         else
          if ( k .lt. kdim ) then
           wh= 0.5*(h(k)-h(k-1))
          else
           wh= 0.5*(hmax-h(k-1))
          endif
          fout=sum - wh*f(k-1)
          fw=wsum - wh
          return
         endif
      enddo

      fval=f(n2)
      wh= 0.5*(h(n2)-h(n2-1))
      if(fval .gt. bmiss)then
        fout=sum + wh*fval
        fw=wsum + wh
      else
        fout=sum - wh*f(n2-1)
        fw=wsum - wh
      endif

      return
      end