subroutine boxintegralSE(f,n1,n2,bmiss,fout,fw)
C
C 06/18/07 updated for n1=n2 > 1 (to avoid negative FOUT)
C revised 06/05/02 see 4th note
C based on sub_depth_integral.f 
C JIANT, 06/04/01
C
C Computes Root_Mean_Square errors of vertically integrated data. 
C
C  Definition of vertical integration:
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) see sub_depth_integral.f
C
      dimension f(33)
      double precision dpp,wh,sum,wsum,h(34)
      data h/0.d0,10.d0,20.d0,30.d0,50.d0,75.d0,100.d0,
     *      125.d0,150.d0,200.d0,250.d0,300.d0,400.d0,500.d0,
     *      600.d0,700.d0,800.d0,
     *      900.d0,1000.d0,1100.d0,1200.d0,1300.d0,1400.d0,
     *      1500.d0,1750.d0,2000.d0,2500.d0,3000.d0,
     *      3500.d0,4000.d0,4500.d0,5000.d0,5500.d0,6000.d0/

      nsum=0
      fout=bmiss
      fw=bmiss

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

      if(n1.eq.n2)then
       if(n1.gt.1)then
         fw= 0.5*(h(n1+1)-h(n1-1))  
         fout= fval*fw
C line below added on 06/18/07
         fout=sqrt(fout*fout)
         return
       else
         fout=real( dsqrt(sum) )
         fw=real( dsqrt(wsum) )
         return
       endif
      endif

      do k=n1+1,n2-1
         fval=f(k)
         if(fval .gt. bmiss)then
	  wh= 0.5d0*(h(k+1)-h(k-1)) 
	  sum=sum + wh*wh*dble(fval*fval) 
          wsum=wsum + wh*wh
          nsum=nsum+1
         elseif(nsum.eq.1)then
          fout=real( dsqrt(sum) )
          fw=real( dsqrt(wsum) )
          return
         else
          wh= 0.5d0*(h(k)-h(k-1))
          dpp=  sum - wh*wh*dble(f(k-1)*f(k-1))
          fout= real( dsqrt(dpp) )
          fw= real( dsqrt(wsum - wh*wh) )
          return
         endif
      enddo

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

      return
      end