!***********************************************************************
SUBROUTINE neutral(u10_mps, ustar_mps, obklen_m, &
        u10n_mps)
  !-----------------------------------------------------------------------
  ! SUBROUTINE to compute u10 neutral wind speed
  ! inputs
  ! u10_mps - wind speed at 10 m (m/s)
  ! ustar_mps - friction velocity (m/s)
  ! obklen_m - monin-obukhov length scale (m)
  ! outputs
  ! u10n_mps - wind speed at 10 m under neutral conditions (m/s)
  ! following code assumes reference height Z is 10m, consistent with use
  ! of u10 and u10_neutral.  If not, code
  ! should be changed so that constants of 50. and 160. in equations
  ! below are changed to -5 * Z and -16 * Z respectively.
  ! Reference:  G. L. Geernaert.  'Bulk parameterizations for the
  ! wind stress and heat fluxes,' in Surface Waves and Fluxes, Vol. I,
  ! Current Theory, Geernaert and W.J. Plant, editors, Kluwer Academic
  ! Publishers, Boston, MA, 1990.
  ! SUBROUTINE written Feb 2001 by eg chapman
  ! adapted to LMD-ZT by E. Cosme 310801
  ! Following Will Shaw (PNL, Seattle) the theory applied for flux
  ! calculation with the scheme of Nightingale et al. (2000) does not
  ! hold anymore when -1<obklen<20. In this case, u10n is set to 0,
  ! so that the transfer velocity  computed in nightingale.F will also
  ! be 0. The flux is then set to 0.
  !----------------------------------------------------------------------
  !
  USE dimphy
  INCLUDE "dimensions.h"
  !
  real :: u10_mps(klon), ustar_mps(klon), obklen_m(klon)
  real :: u10n_mps(klon)
  real :: pi, von_karman
  ! parameter (pi = 3.141592653589793, von_karman = 0.4)
  ! pour etre coherent avec vk de bl_for_dms.F
  parameter (pi = 3.141592653589793, von_karman = 0.35)
  !
  real :: phi, phi_inv, phi_inv_sq, f1, f2, f3, dum1, psi
  integer :: i

  psi = 0.
  do i = 1, klon

    if (u10_mps(i) < 0.) u10_mps(i) = 0.0

    if  (obklen_m(i) < 0.) then
      phi = (1. - 160. / obklen_m(i))**(-0.25)
      phi_inv = 1. / phi
      phi_inv_sq = 1. / phi * 1. / phi
      f1 = (1. + phi_inv) / 2.
      f2 = (1. + phi_inv_sq) / 2.
      ! following to avoid numerical overruns. reCALL tan(90deg)=infinity
      dum1 = min (1.e24, phi_inv)
      f3 = atan(dum1)
      psi = 2. * log(f1) + log(f2) - 2. * f3 + pi / 2.
    else if (obklen_m(i) > 0.) then
      psi = -50. / obklen_m(i)
    end if

    u10n_mps(i) = u10_mps(i) + (ustar_mps(i) * psi / von_karman)
    ! u10n set to 0. if -1 < obklen < 20
    if ((obklen_m(i)>-1.).and.(obklen_m(i)<20.)) then
      u10n_mps(i) = 0.
    endif
    if (u10n_mps(i) < 0.) u10n_mps(i) = 0.0

  enddo

end subroutine neutral
!***********************************************************************