subroutine dispthicknessmo(output)
! output is an array of length (ib:ie)) containing displacement thickness values
use modglobal, only : ib,ie,kb,ke,dzf,xf,Uinf,numol,grav,prandtlmoli
use modsurfdata, only : thls
implicit none
real, dimension(ib:ie), intent(out) :: output !< dispacement thickness
! real, dimension(kb:ke) :: dthick
real :: dispm
real :: disp2m
real :: xfdispm
real :: ustar,tstar,blth
real :: B = 5.0 ! Wake parameter
real :: C = 0.5 ! Coles parameter
real :: kappa = 0.41 ! Von k�rm�n constant
real :: cmo = 0.702 ! constant in MO theory (0.135*5.2)
real :: lam ! = Uinf/ustar
real :: func,dfunc,utaunu,lmo
integer :: i,n
blth = di ! initial value
do i=ib,ie
ustar = sqrt(abs(2.*numol* Utav(i,kb)/dzf(kb))) ! average streamwise friction at x-location
tstar = numol*prandtlmoli* 2.*(Ttav(i,kb)-thls)/(dzf(kb)*ustar) ! average shear temp. at x-location
lmo = (thls*ustar**2)/(kappa*grav*tstar) ! obukhov length at this x-location
if ((lmo >= 10000.) .or. (lmo <= 0.01)) then
lmo = 1000.
end if
! lmo = 0.3 !! TEMPORARY
utaunu = ustar / numol
lam = Uinf / ustar
do n=1,10 ! Newton Raphson method to find BL height
! write(6,*) 'blth,ustar,tstar,Lmo =',blth,ustar,tstar,lmo
func = log(blth) + (cmo*blth/lmo) + log(utaunu) - kappa*(lam-B) +2.*C
! func = log(blth) + log(utaunu) - kappa*(lam-B) +2.*C
dfunc = 1./blth + cmo/lmo
blth = blth - (func / dfunc)
if (blth <= 0.) then
blth = di
end if
end do
output(i) = ((1. + C + 0.5*cmo*blth/lmo) / (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 dispthicknessmo