dispthickness Subroutine

  subroutine dispthickness(output)
! output is an array of length (ib:ie)) containing displacement thickness values
    use modglobal, only : ib,ie,kb,ke,dzf,xf,Uinf,numol
    implicit none

       real, dimension(ib:ie),       intent(out) :: output  !< dispacement thickness
!       real, dimension(kb:ke)                    :: dthick
       real    :: dispm
       real    :: disp2m
       real    :: xfdispm
       real    :: ustar,blth
       real    :: B     = 5.0       ! Wake parameter
       real    :: C     = 0.5       ! Coles parameter
       real    :: kappa = 0.41      ! Von k�r�n constant
       real    :: lam               ! = Uinf/ustar
       integer :: i
   

       do i=ib,ie
         ustar     =  sqrt(abs(2*numol* Utav(i,kb)/dzf(kb)))              ! average streamwise friction
         lam       = Uinf / ustar
         blth      = (lam*numol/Uinf) * exp( kappa * (lam - B) - 2.*C)      ! See App. Lund et al.
         output(i) = ((1. + C) / (kappa*lam) ) * blth                 
       end do
       dispm  = sum(output(ib:ie))/(ie-ib+1)      ! mean(displacement)
       disp2m = sum(output(ib:ie)**2.)/(ie-ib+1)   ! mean(displacement^2)
       xfdispm  = sum(xf(ib:ie)*output(ib:ie))/(ie-ib+1)     ! mean(xf*displ)

       ddispdx = (xfdispm - (xfm*dispm)) / (xf2m -xfm**2.)    ! this is d/dx(delta*)


  end subroutine dispthickness