blthicknessmo Subroutine

  subroutine blthicknessmo(output,ustar,lmo)
! This routine compute the BL thicknes for a buoyancy affected boundary layer:
! Newton-Raphson method is used 
    use modglobal, only : numol,Uinf
    implicit none

       real, intent(in)                   :: ustar      ! friction velocity
       real, intent(in)                   :: lmo        ! Obukhov length
       real, intent(inout)                :: output     !< BL thickness based on law of the wake
!       real, dimension(kb:ke)             :: mthick
!       real                               :: ucrit
!       real                               :: umax
!       integer :: k
       real :: B     = 5.0       ! Wake parameter
       real :: C     = 0.5       ! Coles parameter
       real :: kappa = 0.41      ! Von k�r�n constant
       real :: cmo   = 0.702     ! Constant in MO theory (0.135*5.2)
       real :: lam               ! = Uinf/ustar
       real :: func,dfunc,utaunu
       integer :: n
       
       utaunu = ustar / numol 
       lam = Uinf / ustar
!       write(6,*) 'Initial delta, Lmo =', output,lmo
       do n=1,10
         func   = log(output) + (cmo*output/lmo) + log(utaunu) - kappa*(lam-B) +2.*C
!         func   = log(output) + log(utaunu) - kappa*(lam-B) +2.*C
         dfunc  = 1./output + cmo/lmo
         output = output - (func / dfunc)
         if (output <= 0.) then
           output = di
         end if
       end do
!       write(6,*) 'Computed delta, Lmo =', output,lmo

   end subroutine blthicknessmo