scale_bulkflux.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
#include "scalelib.h"
module scale_bulkflux
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_io
  use scale_prof
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: bulkflux_setup

  abstract interface
     subroutine bc( &
          Ustar,   & ! (out)
          tstar,   & ! (out)
          qstar,   & ! (out)
          uabs,    & ! (out)
          ra,      & ! (out)
          fracu10, & ! (out)
          fract2,  & ! (out)
          fracq2,  & ! (out)
          t1,      & ! (in)
          t0,      & ! (in)
          p1,      & ! (in)
          p0,      & ! (in)
          q1,      & ! (in)
          q0,      & ! (in)
          u1,      & ! (in)
          v1,      & ! (in)
          z1,      & ! (in)
          pbl,     & ! (in)
          z0m,     & ! (in)
          z0h,     & ! (in)
          z0e      ) ! (in)
       use scale_precision
       implicit none

       real(RP), intent(out) :: ustar   ! friction velocity [m/s]
       real(RP), intent(out) :: tstar   ! friction temperature [K]
       real(RP), intent(out) :: qstar   ! friction mixing rate [kg/kg]
       real(RP), intent(out) :: uabs    ! modified absolute velocity [m/s]
       real(RP), intent(out) :: ra      ! Aerodynamic resistance (=1/Ce)
       real(RP), intent(out) :: fracu10 ! calculation parameter for U10 [-]
       real(RP), intent(out) :: fract2  ! calculation parameter for T2 [-]
       real(RP), intent(out) :: fracq2  ! calculation parameter for Q2 [-]

       real(RP), intent(in) :: t1  ! tempearature at the lowest atmospheric layer [K]
       real(RP), intent(in) :: t0  ! skin temperature [K]
       real(RP), intent(in) :: p1  ! pressure at the lowest atmospheric layer [Pa]
       real(RP), intent(in) :: p0  ! surface pressure [Pa]
       real(RP), intent(in) :: q1  ! mixing ratio at the lowest atmospheric layer [kg/kg]
       real(RP), intent(in) :: q0  ! surface mixing ratio [kg/kg]
       real(RP), intent(in) :: u1  ! zonal wind at the lowest atmospheric layer [m/s]
       real(RP), intent(in) :: v1  ! meridional wind at the lowest atmospheric layer [m/s]
       real(RP), intent(in) :: z1  ! height at the lowest atmospheric layer [m]
       real(RP), intent(in) :: pbl ! the top of atmospheric mixing layer [m]
       real(RP), intent(in) :: z0m ! roughness length of momentum [m]
       real(RP), intent(in) :: z0h ! roughness length of heat [m]
       real(RP), intent(in) :: z0e ! roughness length of moisture [m]
     end subroutine bc
  end interface

  procedure(bc), pointer :: bulkflux => null()
  public :: bulkflux

  !-----------------------------------------------------------------------------
  !
  !++ Public parameters & variables
  !
  !-----------------------------------------------------------------------------
  !
  !++ Private procedure
  !
  private :: bulkflux_u95
  private :: bulkflux_b91w01
  private :: fm_unstable
  private :: fh_unstable
  private :: fm_stable
  private :: fh_stable

  !-----------------------------------------------------------------------------
  !
  !++ Private parameters & variables
  !
  character(len=H_SHORT), private :: bulkflux_type = 'B91W01' ! 'U95', 'B94', and 'B91W01'

  logical,  private :: bulkflux_use_mean = .true.

  integer,  private :: bulkflux_itr_sa_max = 5  ! maximum iteration number for successive approximation
  integer,  private :: bulkflux_itr_nr_max = 10 ! maximum iteration number for Newton-Raphson method

  real(RP), private :: bulkflux_err_min = 1.0e-3_rp ! minimum value of error

  real(RP), private :: bulkflux_wscf ! empirical scaling factor of Wstar (Beljaars 1994)

  ! limiter
  real(RP), private :: bulkflux_uabs_min  = 1.0e-2_rp ! minimum of Uabs [m/s]
  real(RP), private :: bulkflux_wstar_min = 1.0e-4_rp ! minimum of W* [m/s]

  logical,  private :: flag_w01

contains

  !-----------------------------------------------------------------------------
  !
  !-----------------------------------------------------------------------------
  subroutine bulkflux_setup( dx )
    use scale_prc, only: &
       prc_abort
    implicit none
    real(RP), intent(in) :: dx

    namelist / param_bulkflux / &
       bulkflux_type,       &
       bulkflux_use_mean,   &
       bulkflux_itr_sa_max, &
       bulkflux_itr_nr_max, &
       bulkflux_err_min,    &
       bulkflux_wscf,       &
       bulkflux_uabs_min,   &
       bulkflux_wstar_min

    integer :: ierr
    !---------------------------------------------------------------------------

    log_newline
    log_info("BULKFLUX_setup",*) 'Setup'

    ! WSCF = 1.2 for dx > 1 km (Beljaars 1994)
    ! lim_{dx->0} WSCF = 0  for LES (Fig. 6 Kitamura and Ito 2016 BLM)
    bulkflux_wscf = 1.2_rp * min(dx*1.e-3_rp, 1.0_rp)


    !--- read namelist
    rewind(io_fid_conf)
    read(io_fid_conf,nml=param_bulkflux,iostat=ierr)

    if( ierr < 0 ) then !--- missing
       log_info("BULKFLUX_setup",*) 'Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       log_error("BULKFLUX_setup",*) 'Not appropriate names in namelist PARAM_BULKFLUX. Check!'
       call prc_abort
    endif
    log_nml(param_bulkflux)

    log_newline
    log_info("BULKFLUX_setup",*) 'Scheme for surface bulk flux : ', trim(bulkflux_type)
    select case(bulkflux_type)
    case('U95')
       log_info_cont(*) '=> Uno et al.(1995)'
       bulkflux => bulkflux_u95
    case('B91W01')
       log_info_cont(*) '=> Beljaars (1991) and Wilson (2001)'
       flag_w01 = .true.
       bulkflux => bulkflux_b91w01
    case('B94')
       log_info_cont(*) '=> Beljaars (1994)'
       flag_w01 = .false.
       bulkflux => bulkflux_b91w01
    case default
       log_error("BULKFLUX_setup",*) 'Unsupported BULKFLUX_type. STOP'
       call prc_abort
    end select

    return
  end subroutine bulkflux_setup

  !-----------------------------------------------------------------------------
  ! ref. Uno et al. (1995)
  !-----------------------------------------------------------------------------
  subroutine bulkflux_u95( &
      Ustar,   & ! (out)
      tstar,   & ! (out)
      qstar,   & ! (out)
      uabs,    & ! (out)
      ra,      & ! (out)
      fracu10, & ! (out)
      fract2,  & ! (out)
      fracq2,  & ! (out)
      t1,      & ! (in)
      t0,      & ! (in)
      p1,      & ! (in)
      p0,      & ! (in)
      q1,      & ! (in)
      q0,      & ! (in)
      u1,      & ! (in)
      v1,      & ! (in)
      z1,      & ! (in)
      pbl,     & ! (in)
      z0m,     & ! (in)
      z0h,     & ! (in)
      z0e      ) ! (in)
    use scale_const, only: &
      grav   => const_grav,   &
      karman => const_karman, &
      rdry   => const_rdry,   &
      cpdry  => const_cpdry,  &
      pre00  => const_pre00
    implicit none

    ! parameter
    real(RP), parameter :: tPrn = 0.74_rp    ! turbulent Prandtl number (Businger et al. 1971)
    real(RP), parameter :: LFb  = 9.4_rp     ! Louis factor b (Louis 1979)
    real(RP), parameter :: LFbp = 4.7_rp     ! Louis factor b' (Louis 1979)
    real(RP), parameter :: LFdm = 7.4_rp     ! Louis factor d for momemtum (Louis 1979)
    real(RP), parameter :: LFdh = 5.3_rp     ! Louis factor d for heat (Louis 1979)

    ! argument
    real(RP), intent(out) :: Ustar   ! friction velocity [m/s]
    real(RP), intent(out) :: Tstar   ! friction temperature [K]
    real(RP), intent(out) :: Qstar   ! friction mixing rate [kg/kg]
    real(RP), intent(out) :: Uabs    ! modified absolute velocity [m/s]
    real(RP), intent(out) :: Ra      ! Aerodynamic resistance (=1/Ce)
    real(RP), intent(out) :: FracU10 ! calculation parameter for U10 [-]
    real(RP), intent(out) :: FracT2  ! calculation parameter for T2 [-]
    real(RP), intent(out) :: FracQ2  ! calculation parameter for Q2 [-]

    real(RP), intent(in) :: T1  ! tempearature at the lowest atmospheric layer [K]
    real(RP), intent(in) :: T0  ! skin temperature [K]
    real(RP), intent(in) :: P1  ! pressure at the lowest atmospheric layer [Pa]
    real(RP), intent(in) :: P0  ! surface pressure [Pa]
    real(RP), intent(in) :: Q1  ! mixing ratio at the lowest atmospheric layer [kg/kg]
    real(RP), intent(in) :: Q0  ! surface mixing ratio [kg/kg]
    real(RP), intent(in) :: U1  ! zonal wind at the lowest atmospheric layer [m/s]
    real(RP), intent(in) :: V1  ! meridional wind at the lowest atmospheric layer [m/s]
    real(RP), intent(in) :: Z1  ! height at the lowest atmospheric layer [m]
    real(RP), intent(in) :: PBL ! the top of atmospheric mixing layer [m]
    real(RP), intent(in) :: Z0M ! roughness length of momentum [m]
    real(RP), intent(in) :: Z0H ! roughness length of heat [m]
    real(RP), intent(in) :: Z0E ! roughness length of moisture [m]

    ! work
    real(RP) :: RiB0, RiB ! bulk Richardson number [-]
    real(RP) :: C0Z1, C010, C002 ! initial drag coefficient [-]
    real(RP) :: CmZ1, ChZ1, CqZ1, fmZ1, fhZ1, t0thZ1, q0qeZ1
    real(RP) :: Cm10, fm10
    real(RP) :: Cm02, Ch02, Cq02, fm02, fh02, t0th02, q0qe02
    real(RP) :: TH1, TH0
    real(RP) :: logZ1Z0M, log10Z0m, log02Z0m
    real(RP) :: logZ0MZ0E
    real(RP) :: logZ0MZ0H
    !---------------------------------------------------------------------------

    logz1z0m = log( z1      / z0m )
    log10z0m = log( 10.0_rp / z0m )
    log02z0m = log( 2.0_rp  / z0m )

    logz0mz0e = max( log( z0m/z0e ), 1.0_rp )
    logz0mz0h = max( log( z0m/z0h ), 1.0_rp )

    uabs = max( sqrt( u1**2 + v1**2 ), bulkflux_uabs_min )
    th1  = t1 * ( p0 / p1 )**( rdry / cpdry )
    th0  = t0

    ! bulk Richardson number
    rib0 = grav * z1 * ( th1 - th0 ) / ( th1 * uabs**2 )

    c0z1 = ( karman / logz1z0m )**2
    c010 = ( karman / log10z0m )**2
    c002 = ( karman / log02z0m )**2

    rib = rib0

    if( rib0 > 0.0_rp ) then
      ! stable condition
      fmz1 = 1.0_rp / ( 1.0_rp + lfbp * rib )**2
      fhz1 = 1.0_rp / ( 1.0_rp + lfbp * rib )**2
      fm10 = 1.0_rp / ( 1.0_rp + lfbp * rib )**2
      fm02 = 1.0_rp / ( 1.0_rp + lfbp * rib )**2
      fh02 = 1.0_rp / ( 1.0_rp + lfbp * rib )**2
    else
      ! unstable condition
      fmz1 = 1.0_rp - lfb * rib / ( 1.0_rp + lfb * lfdm * c0z1 * sqrt( z1      / z0m ) * sqrt( abs( rib ) ) )
      fhz1 = 1.0_rp - lfb * rib / ( 1.0_rp + lfb * lfdh * c0z1 * sqrt( z1      / z0m ) * sqrt( abs( rib ) ) )
      fm10 = 1.0_rp - lfb * rib / ( 1.0_rp + lfb * lfdm * c010 * sqrt( 10.0_rp / z0m ) * sqrt( abs( rib ) ) )
      fm02 = 1.0_rp - lfb * rib / ( 1.0_rp + lfb * lfdm * c002 * sqrt( 2.0_rp  / z0m ) * sqrt( abs( rib ) ) )
      fh02 = 1.0_rp - lfb * rib / ( 1.0_rp + lfb * lfdh * c002 * sqrt( 2.0_rp  / z0m ) * sqrt( abs( rib ) ) )
    end if

    t0thz1 = 1.0_rp / ( 1.0_rp + logz0mz0h / logz1z0m / sqrt( fmz1 ) * fhz1 )
    q0qez1 = 1.0_rp / ( 1.0_rp + logz0mz0e / logz1z0m / sqrt( fmz1 ) * fhz1 )
    t0th02 = 1.0_rp / ( 1.0_rp + logz0mz0h / log02z0m / sqrt( fm02 ) * fh02 )
    q0qe02 = 1.0_rp / ( 1.0_rp + logz0mz0e / log02z0m / sqrt( fm02 ) * fh02 )

    rib = rib * t0thz1

    if( rib0 > 0.0_rp ) then
      ! stable condition
      fmz1 = 1.0_rp / ( 1.0_rp + lfbp * rib )**2
      fhz1 = 1.0_rp / ( 1.0_rp + lfbp * rib )**2
      fm10 = 1.0_rp / ( 1.0_rp + lfbp * rib )**2
      fm02 = 1.0_rp / ( 1.0_rp + lfbp * rib )**2
      fh02 = 1.0_rp / ( 1.0_rp + lfbp * rib )**2
    else
      ! unstable condition
      fmz1 = 1.0_rp - lfb * rib / ( 1.0_rp + lfb * lfdm * c0z1 * sqrt( z1      / z0m ) * sqrt( abs( rib ) ) )
      fhz1 = 1.0_rp - lfb * rib / ( 1.0_rp + lfb * lfdh * c0z1 * sqrt( z1      / z0m ) * sqrt( abs( rib ) ) )
      fm10 = 1.0_rp - lfb * rib / ( 1.0_rp + lfb * lfdm * c010 * sqrt( 10.0_rp / z0m ) * sqrt( abs( rib ) ) )
      fm02 = 1.0_rp - lfb * rib / ( 1.0_rp + lfb * lfdm * c002 * sqrt( 2.0_rp  / z0m ) * sqrt( abs( rib ) ) )
      fh02 = 1.0_rp - lfb * rib / ( 1.0_rp + lfb * lfdh * c002 * sqrt( 2.0_rp  / z0m ) * sqrt( abs( rib ) ) )
    end if

    t0thz1 = 1.0_rp / ( 1.0_rp + logz0mz0h / logz1z0m / sqrt( fmz1 ) * fhz1 )
    q0qez1 = 1.0_rp / ( 1.0_rp + logz0mz0e / logz1z0m / sqrt( fmz1 ) * fhz1 )
    t0th02 = 1.0_rp / ( 1.0_rp + logz0mz0h / log02z0m / sqrt( fm02 ) * fh02 )
    q0qe02 = 1.0_rp / ( 1.0_rp + logz0mz0e / log02z0m / sqrt( fm02 ) * fh02 )

    cmz1 = c0z1 * fmz1
    chz1 = c0z1 * fhz1 * t0thz1 / tprn
    cqz1 = c0z1 * fhz1 * q0qez1 / tprn
    cm10 = c010 * fm10
    cm02 = c002 * fm02
    ch02 = c002 * fh02 * t0th02 / tprn
    cq02 = c002 * fh02 * q0qe02 / tprn

    ustar = sqrt( cmz1 ) * uabs
    tstar = chz1 * uabs / ustar * ( th1 - th0 )
    qstar = cqz1 * uabs / ustar * ( q1  - q0  )

    fracu10 = sqrt( cmz1 / cm10 )
    fract2  = chz1 / ch02 * sqrt( cm02 / cmz1 )
    fracq2  = cqz1 / cq02 * sqrt( cm02 / cmz1 )

    ra = 1.0_rp / ( cqz1 * uabs )

    return
  end subroutine bulkflux_u95

  !-----------------------------------------------------------------------------
  !
  ! refs. Beljaars (1991) and Wilson (2001)
  !
  ! If you want to run with the original Beljaars scheme (Beljaars and Holtslag 1994),
  ! you should fix the stability functions (fm_unstable, fh_unstable, fm_stable, and fh_stable).
  !
  ! Iteration method: refs. JMA-NHM Description Note II, Mar 2008
  !
  !-----------------------------------------------------------------------------
!OCL SERIAL
  subroutine bulkflux_b91w01( &
      Ustar,   & ! (out)
      tstar,   & ! (out)
      qstar,   & ! (out)
      uabs,    & ! (out)
      ra,      & ! (out)
      fracu10, & ! (out)
      fract2,  & ! (out)
      fracq2,  & ! (out)
      t1,      & ! (in)
      t0,      & ! (in)
      p1,      & ! (in)
      p0,      & ! (in)
      q1,      & ! (in)
      q0,      & ! (in)
      u1,      & ! (in)
      v1,      & ! (in)
      z1,      & ! (in)
      pbl,     & ! (in)
      z0m,     & ! (in)
      z0h,     & ! (in)
      z0e      ) ! (in)
    use scale_const, only: &
      grav    => const_grav,    &
      karman  => const_karman,  &
      rdry    => const_rdry,    &
      cpdry   => const_cpdry,   &
      epstvap => const_epstvap, &
      eps     => const_eps,     &
      pre00   => const_pre00
    implicit none

    ! parameter
    real(DP), parameter :: dIL = 1.0e-6_dp ! delta [1/m]

    ! argument
    real(RP), intent(out) :: Ustar   ! friction velocity [m/s]
    real(RP), intent(out) :: Tstar   ! friction temperature [K]
    real(RP), intent(out) :: Qstar   ! friction mixing rate [kg/kg]
    real(RP), intent(out) :: Uabs    ! modified absolute velocity [m/s]
    real(RP), intent(out) :: Ra      ! Aerodynamic resistance (=1/Ce)
    real(RP), intent(out) :: FracU10 ! calculation parameter for U10 [-]
    real(RP), intent(out) :: FracT2  ! calculation parameter for T2 [-]
    real(RP), intent(out) :: FracQ2  ! calculation parameter for Q2 [-]

    real(RP), intent(in) :: T1  ! tempearature at the lowest atmospheric layer [K]
    real(RP), intent(in) :: T0  ! skin temperature [K]
    real(RP), intent(in) :: P1  ! pressure at the lowest atmospheric layer [Pa]
    real(RP), intent(in) :: P0  ! surface pressure [Pa]
    real(RP), intent(in) :: Q1  ! mixing ratio at the lowest atmospheric layer [kg/kg]
    real(RP), intent(in) :: Q0  ! mixing ratio at surface [kg/kg]
    real(RP), intent(in) :: U1  ! zonal wind at the lowest atmospheric layer [m/s]
    real(RP), intent(in) :: V1  ! meridional wind at the lowest atmospheric layer [m/s]
    real(RP), intent(in) :: Z1  ! height at the lowest atmospheric layer [m]
    real(RP), intent(in) :: PBL ! the top of atmospheric mixing layer [m]
    real(RP), intent(in) :: Z0M ! roughness length of momentum [m]
    real(RP), intent(in) :: Z0H ! roughness length of heat [m]
    real(RP), intent(in) :: Z0E ! roughness length of moisture [m]

    ! work
    integer :: n

    real(DP) :: IL   ! inversed Obukhov length [1/m]
    real(DP) :: res  ! residual
    real(DP) :: dres ! d(residual)/dIL

    real(DP) :: RiB0 ! bulk Richardson number [no unit]
    real(DP) :: Wstar, dWstar ! free convection velocity scale [m/s]

    real(DP) :: UabsUS, UabsS, UabsC
    real(DP) :: dUabsUS, dUabsS

    real(DP) :: UstarUS, UstarS, UstarC
    real(DP) :: TstarUS, TstarS, TstarC
    real(DP) :: QstarUS, QstarS, QstarC

    real(DP) :: dUstarUS, dUstarS, dUstarC
    real(DP) :: dTstarUS, dTstarS, dTstarC
    real(DP) :: dQstarUS, dQstarS, dQstarC

    real(DP) :: FracU10US, FracU10S, FracU10C
    real(DP) :: FracT2US,  FracT2S,  FracT2C
    real(DP) :: FracQ2US,  FracQ2S,  FracQ2C

    real(DP) :: TH1, TH0, THM
    real(DP) :: TV1, TV0, TVM
    real(DP) :: QM
    real(DP) :: sw, tmp

    real(DP) :: BFLX, dBFLX

    real(DP) :: DP_Z1, DP_Z0M, DP_Z0H, DP_Z0E
    real(DP) :: log_Z1ovZ0M, log_Z1ovZ0H, log_Z1ovZ0E
    real(DP) :: log_10ovZ0M, log_02ovZ0H, log_02ovZ0E

    real(DP) :: denoM, denoH, denoE
    real(DP) :: RzM, RzH, RzE
    !---------------------------------------------------------------------------

    ! convert to DP
    if ( bulkflux_use_mean ) then
       dp_z1  = real( z1*2.0_rp,  kind=dp )
    else
       dp_z1  = real( z1,  kind=dp )
    end if
    dp_z0m = real( z0m, kind=dp )
    dp_z0h = real( z0h, kind=dp )
    dp_z0e = real( z0e, kind=dp )

    uabsc = max( sqrt( u1**2 + v1**2 ), bulkflux_uabs_min )

    th1 = t1 * ( p0 / p1 )**( rdry / cpdry )
    th0 = t0
    thm = ( th1 + th0 ) * 0.5_rp
    qm  = ( q1  + q0  ) * 0.5_rp
    tv1 = th1 * ( 1.0_dp + epstvap * q1 )
    tv0 = th0 * ( 1.0_dp + epstvap * q0 )
    tvm = ( tv1 + tv0 ) * 0.5_rp

    rzm = 1.0_dp - dp_z0m / dp_z1
    rzh = 1.0_dp - dp_z0h / dp_z1
    rze = 1.0_dp - dp_z0e / dp_z1

    ! make log constant
    log_z1ovz0m = log( dp_z1 / dp_z0m )
    log_z1ovz0h = log( dp_z1 / dp_z0h )
    log_z1ovz0e = log( dp_z1 / dp_z0e )

    log_10ovz0m = log( 10.0_dp / dp_z0m )
    log_02ovz0h = log(  2.0_dp / dp_z0h )
    log_02ovz0e = log(  2.0_dp / dp_z0e )

    ! bulk Richardson number at initial step
    rib0 = grav * dp_z1 * ( tv1 - tv0 ) / ( tvm * uabsc**2 )

    ! inversed Obukhov length at initial step
    il = rib0 / dp_z1 * log_z1ovz0m**2 / log_z1ovz0h

    ! free convection velocity scale at initial step
    wstar  = bulkflux_wstar_min
    dwstar = bulkflux_wstar_min

    ! Successive approximation
    do n = 1, bulkflux_itr_sa_max
       if ( bulkflux_use_mean ) then
          denom = log_z1ovz0m &
                - fmm_unstable(dp_z1,il) + fmm_unstable(dp_z0m,il) * dp_z0m / dp_z1 &
                + rzm * ( fm_unstable(dp_z0m,il) - 1.0_dp )
          denoh = log_z1ovz0h &
                - fhm_unstable(dp_z1,il) + fhm_unstable(dp_z0h,il) * dp_z0h / dp_z1 &
                + rzh * ( fh_unstable(dp_z0h,il) - 1.0_dp )
          denoe = log_z1ovz0e &
                - fhm_unstable(dp_z1,il) + fhm_unstable(dp_z0e,il) * dp_z0e / dp_z1 &
                + rze * ( fh_unstable(dp_z0e,il) - 1.0_dp )
       else
          denom = log_z1ovz0m - fm_unstable(dp_z1,il) + fm_unstable(dp_z0m,il)
          denoh = log_z1ovz0h - fh_unstable(dp_z1,il) + fh_unstable(dp_z0h,il)
          denoe = log_z1ovz0e - fh_unstable(dp_z1,il) + fh_unstable(dp_z0e,il)
       end if
      ! unstable condition
      uabsus  = max( sqrt( u1**2 + v1**2 + (bulkflux_wscf*wstar)**2 ), real( BULKFLUX_Uabs_min, kind=DP ) )
      ustarus = karman / denom * uabsus
      tstarus = karman / denoh * ( th1 - th0 )
      qstarus = karman / denoe * ( q1  - q0  )
!      UstarUS = KARMAN / ( log_Z1ovZ0M - fm_unstable(DP_Z1,IL) + fm_unstable(DP_Z0M,IL) ) * UabsUS
!      TstarUS = KARMAN / ( log_Z1ovZ0H - fh_unstable(DP_Z1,IL) + fh_unstable(DP_Z0H,IL) ) * ( TH1 - TH0 )
!      QstarUS = KARMAN / ( log_Z1ovZ0E - fh_unstable(DP_Z1,IL) + fh_unstable(DP_Z0E,IL) ) * ( Q1  - Q0  )

      ! stable condition
       if ( bulkflux_use_mean ) then
          denom = log_z1ovz0m &
                - fmm_stable(dp_z1,il) + fmm_stable(dp_z0m,il) * dp_z0m / dp_z1 &
                + rzm * ( fm_stable(dp_z0m,il) - 1.0_dp )
          denoh = log_z1ovz0h &
                - fhm_stable(dp_z1,il) + fhm_stable(dp_z0h,il) * dp_z0h / dp_z1 &
                + rzh * ( fh_stable(dp_z0h,il) - 1.0_dp )
          denoe = log_z1ovz0e &
                - fhm_stable(dp_z1,il) + fhm_stable(dp_z0e,il) * dp_z0e / dp_z1 &
                + rze * ( fh_stable(dp_z0e,il) - 1.0_dp )
       else
          denom = log_z1ovz0m - fm_stable(dp_z1,il) + fm_stable(dp_z0m,il)
          denoh = log_z1ovz0h - fh_stable(dp_z1,il) + fh_stable(dp_z0h,il)
          denoe = log_z1ovz0e - fh_stable(dp_z1,il) + fh_stable(dp_z0e,il)
       end if
       uabss  = max( sqrt( u1**2 + v1**2 ), bulkflux_uabs_min )
       ustars = karman / denom * uabss
       tstars = karman / denoh * ( th1 - th0 )
       qstars = karman / denoe * ( q1  - q0  )
!      UstarS = KARMAN / ( log_Z1ovZ0M - fm_stable(DP_Z1,IL) + fm_stable(DP_Z0M,IL) ) * UabsS
!      TstarS = KARMAN / ( log_Z1ovZ0H - fh_stable(DP_Z1,IL) + fh_stable(DP_Z0H,IL) ) * ( TH1 - TH0 )
!      QstarS = KARMAN / ( log_Z1ovZ0E - fh_stable(DP_Z1,IL) + fh_stable(DP_Z0E,IL) ) * ( Q1  - Q0  )

      sw = 0.5_dp - sign( 0.5_dp, il ) ! if unstable, sw = 1

      ustarc = ( sw ) * ustarus + ( 1.0_dp-sw ) * ustars
      tstarc = ( sw ) * tstarus + ( 1.0_dp-sw ) * tstars
      qstarc = ( sw ) * qstarus + ( 1.0_dp-sw ) * qstars

      ! estimate buoyancy flux
      bflx = - ustarc * tstarc * ( 1.0_rp + epstvap * qm ) - epstvap * ustarc * qstarc * thm

      ! update free convection velocity scale
      tmp = pbl * grav / t1 * bflx
      sw  = 0.5_dp + sign( 0.5_dp, tmp ) ! if tmp is plus, sw = 1

      wstar = ( tmp * sw )**( 1.0_dp / 3.0_dp )

      ! avoid zero division with UstarC = 0
      sw = 0.5_dp + sign( 0.5_dp, abs( ustarc ) - eps )

      ustarc = ( sw ) * ustarc + ( 1.0_dp-sw ) * eps

      ! estimate the inversed Obukhov length
      il = - karman * grav * bflx / ( ustarc**3 * thm )
    end do


    ! Newton-Raphson method
    do n = 1, bulkflux_itr_nr_max
      ! unstable condition
       if ( bulkflux_use_mean ) then
          denom = log_z1ovz0m &
                - fmm_unstable(dp_z1,il) + fmm_unstable(dp_z0m,il) * dp_z0m / dp_z1 &
                + rzm * ( fm_unstable(dp_z0m,il) - 1.0_dp )
          denoh = log_z1ovz0h &
                - fhm_unstable(dp_z1,il) + fhm_unstable(dp_z0h,il) * dp_z0h / dp_z1 &
                + rzh * ( fh_unstable(dp_z0h,il) - 1.0_dp )
          denoe = log_z1ovz0e &
                - fhm_unstable(dp_z1,il) + fhm_unstable(dp_z0e,il) * dp_z0e / dp_z1 &
                + rze * ( fh_unstable(dp_z0e,il) - 1.0_dp )
       else
          denom = log_z1ovz0m - fm_unstable(dp_z1,il) + fm_unstable(dp_z0m,il)
          denoh = log_z1ovz0h - fh_unstable(dp_z1,il) + fh_unstable(dp_z0h,il)
          denoe = log_z1ovz0e - fh_unstable(dp_z1,il) + fh_unstable(dp_z0e,il)
       end if
       uabsus  = max( sqrt( u1**2 + v1**2 + (bulkflux_wscf*wstar)**2 ), real( BULKFLUX_Uabs_min, kind=DP ) )
       ustarus = karman / denom * uabsus
       tstarus = karman / denoh * ( th1 - th0 )
       qstarus = karman / denoe * ( q1  - q0  )
!      UstarUS = KARMAN / ( log_Z1ovZ0M - fm_unstable(DP_Z1,IL) + fm_unstable(DP_Z0M,IL) ) * UabsUS
!      TstarUS = KARMAN / ( log_Z1ovZ0H - fh_unstable(DP_Z1,IL) + fh_unstable(DP_Z0H,IL) ) * ( TH1 - TH0 )
!      QstarUS = KARMAN / ( log_Z1ovZ0E - fh_unstable(DP_Z1,IL) + fh_unstable(DP_Z0E,IL) ) * ( Q1  - Q0  )

      ! stable condition
       if ( bulkflux_use_mean ) then
          denom = log_z1ovz0m &
                - fmm_stable(dp_z1,il) + fmm_stable(dp_z0m,il) * dp_z0m / dp_z1 &
                + rzm * ( fm_stable(dp_z0m,il) - 1.0_dp )
          denoh = log_z1ovz0h &
                - fhm_stable(dp_z1,il) + fhm_stable(dp_z0h,il) * dp_z0h / dp_z1 &
                + rzh * ( fh_stable(dp_z0h,il) - 1.0_dp )
          denoe = log_z1ovz0e &
                - fhm_stable(dp_z1,il) + fhm_stable(dp_z0e,il) * dp_z0e / dp_z1 &
                + rze * ( fh_stable(dp_z0e,il) - 1.0_dp )
       else
          denom = log_z1ovz0m - fm_stable(dp_z1,il) + fm_stable(dp_z0m,il)
          denoh = log_z1ovz0h - fh_stable(dp_z1,il) + fh_stable(dp_z0h,il)
          denoe = log_z1ovz0e - fh_stable(dp_z1,il) + fh_stable(dp_z0e,il)
       end if
       uabss  = max( sqrt( u1**2 + v1**2 ), bulkflux_uabs_min )
       ustars = karman / denom * uabss
       tstars = karman / denoh * ( th1 - th0 )
       qstars = karman / denoe * ( q1  - q0  )
!      UstarS = KARMAN / ( log_Z1ovZ0M - fm_stable(DP_Z1,IL) + fm_stable(DP_Z0M,IL) ) * UabsS
!      TstarS = KARMAN / ( log_Z1ovZ0H - fh_stable(DP_Z1,IL) + fh_stable(DP_Z0H,IL) ) * ( TH1 - TH0 )
!      QstarS = KARMAN / ( log_Z1ovZ0E - fh_stable(DP_Z1,IL) + fh_stable(DP_Z0E,IL) ) * ( Q1  - Q0  )

      sw = 0.5_dp - sign( 0.5_dp, il ) ! if unstable, sw = 1

      ustarc = ( sw ) * ustarus + ( 1.0_dp-sw ) * ustars
      tstarc = ( sw ) * tstarus + ( 1.0_dp-sw ) * tstars
      qstarc = ( sw ) * qstarus + ( 1.0_dp-sw ) * qstars

      ! estimate buoyancy flux
      bflx = - ustarc * tstarc * ( 1.0_rp + epstvap * qm ) - epstvap * ustarc * qstarc * thm

      ! update free convection velocity scale
      tmp = pbl * grav / t1 * bflx
      sw  = 0.5_dp + sign( 0.5_dp, tmp ) ! if tmp is plus, sw = 1

      wstar = ( tmp * sw )**( 1.0_dp / 3.0_dp )

      ! avoid zero division with UstarC = 0
      sw = 0.5_dp + sign( 0.5_dp, abs( ustarc ) - eps )

      ustarc = ( sw ) * ustarc + ( 1.0_dp-sw ) * eps

      ! calculate residual
      res = il + karman * grav * bflx / ( ustarc**3 * thm )

      ! unstable condition
       if ( bulkflux_use_mean ) then
          denom = log_z1ovz0m &
                - fmm_unstable(dp_z1,il+dil) + fmm_unstable(dp_z0m,il+dil) * dp_z0m / dp_z1 &
                + rzm * ( fm_unstable(dp_z0m,il+dil) - 1.0_dp )
          denoh = log_z1ovz0h &
                - fhm_unstable(dp_z1,il+dil) + fhm_unstable(dp_z0h,il+dil) * dp_z0h / dp_z1 &
                + rzh * ( fh_unstable(dp_z0h,il+dil) - 1.0_dp )
          denoe = log_z1ovz0e &
                - fhm_unstable(dp_z1,il+dil) + fhm_unstable(dp_z0e,il+dil) * dp_z0e / dp_z1 &
                + rze * ( fh_unstable(dp_z0e,il+dil) - 1.0_dp )
       else
          denom = log_z1ovz0m - fm_unstable(dp_z1,il+dil) + fm_unstable(dp_z0m,il+dil)
          denoh = log_z1ovz0h - fh_unstable(dp_z1,il+dil) + fh_unstable(dp_z0h,il+dil)
          denoe = log_z1ovz0e - fh_unstable(dp_z1,il+dil) + fh_unstable(dp_z0e,il+dil)
       end if
       duabsus  = max( sqrt( u1**2 + v1**2 + (bulkflux_wscf*wstar)**2 ), real( BULKFLUX_Uabs_min, kind=DP ) )
       dustarus = karman / denom * uabsus
       dtstarus = karman / denoh * ( th1 - th0 )
       dqstarus = karman / denoe * ( q1  - q0  )
!      dUstarUS = KARMAN / ( log_Z1ovZ0M - fm_unstable(DP_Z1,IL+dIL) + fm_unstable(DP_Z0M,IL+dIL) ) * dUabsUS
!      dTstarUS = KARMAN / ( log_Z1ovZ0H - fh_unstable(DP_Z1,IL+dIL) + fh_unstable(DP_Z0H,IL+dIL) ) * ( TH1 - TH0 )
!      dQstarUS = KARMAN / ( log_Z1ovZ0E - fh_unstable(DP_Z1,IL+dIL) + fh_unstable(DP_Z0E,IL+dIL) ) * ( Q1  - Q0  )

      ! stable condition
       if ( bulkflux_use_mean ) then
          denom = log_z1ovz0m &
                - fmm_stable(dp_z1,il+dil) + fmm_stable(dp_z0m,il+dil) * dp_z0m / dp_z1 &
                + rzm * ( fm_stable(dp_z0m,il+dil) - 1.0_dp )
          denoh = log_z1ovz0h &
                - fhm_stable(dp_z1,il+dil) + fhm_stable(dp_z0h,il+dil) * dp_z0h / dp_z1 &
                + rzh * ( fh_stable(dp_z0h,il+dil) - 1.0_dp )
          denoe = log_z1ovz0e &
                - fhm_stable(dp_z1,il+dil) + fhm_stable(dp_z0e,il+dil) * dp_z0e / dp_z1 &
                + rze * ( fh_stable(dp_z0e,il+dil) - 1.0_dp )
       else
          denom = log_z1ovz0m - fm_stable(dp_z1,il+dil) + fm_stable(dp_z0m,il+dil)
          denoh = log_z1ovz0h - fh_stable(dp_z1,il+dil) + fh_stable(dp_z0h,il+dil)
          denoe = log_z1ovz0e - fh_stable(dp_z1,il+dil) + fh_stable(dp_z0e,il+dil)
       end if
       duabss  = max( sqrt( u1**2 + v1**2 ), bulkflux_uabs_min )
       dustars = karman / denom * uabss
       dtstars = karman / denoh * ( th1 - th0 )
       dqstars = karman / denoe * ( q1  - q0  )
!      dUstarS = KARMAN / ( log_Z1ovZ0M - fm_stable(DP_Z1,IL+dIL) + fm_stable(DP_Z0M,IL+dIL) ) * dUabsS
!      dTstarS = KARMAN / ( log_Z1ovZ0H - fh_stable(DP_Z1,IL+dIL) + fh_stable(DP_Z0H,IL+dIL) ) * ( TH1 - TH0 )
!      dQstarS = KARMAN / ( log_Z1ovZ0E - fh_stable(DP_Z1,IL+dIL) + fh_stable(DP_Z0E,IL+dIL) ) * ( Q1  - Q0  )

      sw = 0.5_dp - sign( 0.5_dp, il+dil ) ! if unstable, sw = 1

      dustarc = ( sw ) * dustarus + ( 1.0_dp-sw ) * dustars
      dtstarc = ( sw ) * dtstarus + ( 1.0_dp-sw ) * dtstars
      dqstarc = ( sw ) * dqstarus + ( 1.0_dp-sw ) * dqstars

      ! estimate buoyancy flux
      dbflx = - dustarc * dtstarc * ( 1.0_rp + epstvap * qm ) - epstvap * dustarc * dqstarc * thm

      ! update d(free convection velocity scale)
      tmp = pbl * grav / t1 * dbflx
      sw  = 0.5_dp + sign( 0.5_dp, tmp ) ! if tmp is plus, sw = 1

      dwstar = ( tmp * sw )**( 1.0_dp / 3.0_dp )

      ! avoid zero division with dUstarC = 0
      sw = 0.5_dp + sign( 0.5_dp, abs( dustarc ) - eps )

      dustarc = ( sw ) * dustarc + ( 1.0_dp-sw ) * eps

      ! calculate d(residual)/dIL
      dres = 1.0_dp + karman * grav / ( thm * dil ) * ( dbflx / dustarc**3 - bflx / ustarc**3 )

      ! stop iteration to prevent numerical error
      if( abs( dres ) < eps ) exit

      ! convergence test with error levels
      if( abs( res/dres ) < bulkflux_err_min ) exit

      ! avoid sign changing
      if( il * ( il - res / dres ) < 0.0_rp ) exit

      ! update the inversed Obukhov length
      il = il - res / dres
    end do

    ! Successive approximation after Newton-Raphson method
    if( .NOT. abs( res/dres ) < bulkflux_err_min ) then
      ! unstable condition
       if ( bulkflux_use_mean ) then
          denom = log_z1ovz0m &
                - fmm_unstable(dp_z1,il) + fmm_unstable(dp_z0m,il) * dp_z0m / dp_z1 &
                + rzm * ( fm_unstable(dp_z0m,il) - 1.0_dp )
          denoh = log_z1ovz0h &
                - fhm_unstable(dp_z1,il) + fhm_unstable(dp_z0h,il) * dp_z0h / dp_z1 &
                + rzh * ( fh_unstable(dp_z0h,il) - 1.0_dp )
          denoe = log_z1ovz0e &
                - fhm_unstable(dp_z1,il) + fhm_unstable(dp_z0e,il) * dp_z0e / dp_z1 &
                + rze * ( fh_unstable(dp_z0e,il) - 1.0_dp )
       else
          denom = log_z1ovz0m - fm_unstable(dp_z1,il) + fm_unstable(dp_z0m,il)
          denoh = log_z1ovz0h - fh_unstable(dp_z1,il) + fh_unstable(dp_z0h,il)
          denoe = log_z1ovz0e - fh_unstable(dp_z1,il) + fh_unstable(dp_z0e,il)
       end if
       uabsus  = max( sqrt( u1**2 + v1**2 + (bulkflux_wscf*wstar)**2 ), real( BULKFLUX_Uabs_min, kind=DP ) )
       ustarus = karman / denom * uabsus
       tstarus = karman / denoh * ( th1 - th0 )
       qstarus = karman / denoe * ( q1  - q0  )
!      UstarUS = KARMAN / ( log_Z1ovZ0M - fm_unstable(DP_Z1,IL) + fm_unstable(DP_Z0M,IL) ) * UabsUS
!      TstarUS = KARMAN / ( log_Z1ovZ0H - fh_unstable(DP_Z1,IL) + fh_unstable(DP_Z0H,IL) ) * ( TH1 - TH0 )
!      QstarUS = KARMAN / ( log_Z1ovZ0E - fh_unstable(DP_Z1,IL) + fh_unstable(DP_Z0E,IL) ) * ( Q1  - Q0  )

      ! stable condition
       if ( bulkflux_use_mean ) then
          denom = log_z1ovz0m &
                - fmm_stable(dp_z1,il) + fmm_stable(dp_z0m,il) * dp_z0m / dp_z1 &
                + rzm * ( fm_stable(dp_z0m,il) - 1.0_dp )
          denoh = log_z1ovz0h &
                - fhm_stable(dp_z1,il) + fhm_stable(dp_z0h,il) * dp_z0h / dp_z1 &
                + rzh * ( fh_stable(dp_z0h,il) - 1.0_dp )
          denoe = log_z1ovz0e &
                - fhm_stable(dp_z1,il) + fhm_stable(dp_z0e,il) * dp_z0e / dp_z1 &
                + rze * ( fh_stable(dp_z0e,il) - 1.0_dp )
       else
          denom = log_z1ovz0m - fm_stable(dp_z1,il) + fm_stable(dp_z0m,il)
          denoh = log_z1ovz0h - fh_stable(dp_z1,il) + fh_stable(dp_z0h,il)
          denoe = log_z1ovz0e - fh_stable(dp_z1,il) + fh_stable(dp_z0e,il)
       end if
       uabss  = max( sqrt( u1**2 + v1**2 ), bulkflux_uabs_min )
       ustars = karman / denom * uabss
       tstars = karman / denoh * ( th1 - th0 )
       qstars = karman / denoe * ( q1  - q0  )
!      UstarS = KARMAN / ( log_Z1ovZ0M - fm_stable(DP_Z1,IL) + fm_stable(DP_Z0M,IL) ) * UabsS
!      TstarS = KARMAN / ( log_Z1ovZ0H - fh_stable(DP_Z1,IL) + fh_stable(DP_Z0H,IL) ) * ( TH1 - TH0 )
!      QstarS = KARMAN / ( log_Z1ovZ0E - fh_stable(DP_Z1,IL) + fh_stable(DP_Z0E,IL) ) * ( Q1  - Q0  )

      sw = 0.5_dp - sign( 0.5_dp, il ) ! if unstable, sw = 1

      ustarc = ( sw ) * ustarus + ( 1.0_dp-sw ) * ustars
      tstarc = ( sw ) * tstarus + ( 1.0_dp-sw ) * tstars
      qstarc = ( sw ) * qstarus + ( 1.0_dp-sw ) * qstars

      ! estimate buoyancy flux
      bflx = - ustarc * tstarc * ( 1.0_rp + epstvap * qm ) - epstvap * ustarc * qstarc * thm

      ! update free convection velocity scale
      tmp = pbl * grav / t1 * bflx
      sw  = 0.5_dp + sign( 0.5_dp, tmp ) ! if tmp is plus, sw = 1

      wstar = ( tmp * sw )**( 1.0_dp / 3.0_dp )

      ! avoid zero division with UstarC = 0
      sw = 0.5_dp + sign( 0.5_dp, abs( ustarc ) - eps )

      ustarc = ( sw ) * ustarc + ( 1.0_dp-sw ) * eps

      ! estimate the inversed Obukhov length
      il = - karman * grav * bflx / ( ustarc**3 * thm )
    end if


    ! calculate Ustar, Tstar, and Qstar based on IL

    ! unstable condition
    if ( bulkflux_use_mean ) then
       denom = log_z1ovz0m &
             - fmm_unstable(dp_z1,il) + fmm_unstable(dp_z0m,il) * dp_z0m / dp_z1 &
             + rzm * ( fm_unstable(dp_z0m,il) - 1.0_dp )
       denoh = log_z1ovz0h &
             - fhm_unstable(dp_z1,il) + fhm_unstable(dp_z0h,il) * dp_z0h / dp_z1 &
             + rzh * ( fh_unstable(dp_z0h,il) - 1.0_dp )
       denoe = log_z1ovz0e &
             - fhm_unstable(dp_z1,il) + fhm_unstable(dp_z0e,il) * dp_z0e / dp_z1 &
             + rze * ( fh_unstable(dp_z0e,il) - 1.0_dp )
    else
       denom = log_z1ovz0m - fm_unstable(dp_z1,il) + fm_unstable(dp_z0m,il)
       denoh = log_z1ovz0h - fh_unstable(dp_z1,il) + fh_unstable(dp_z0h,il)
       denoe = log_z1ovz0e - fh_unstable(dp_z1,il) + fh_unstable(dp_z0e,il)
    end if
    uabsus  = max( sqrt( u1**2 + v1**2 + (bulkflux_wscf*wstar)**2 ), real( BULKFLUX_Uabs_min, kind=DP ) )
    ustarus = karman / denom * uabsus
    tstarus = karman / denoh * ( th1 - th0 )
    qstarus = karman / denoe * ( q1  - q0  )
!    UstarUS = KARMAN / ( log_Z1ovZ0M - fm_unstable(DP_Z1,IL) + fm_unstable(DP_Z0M,IL) ) * UabsUS
!    TstarUS = KARMAN / ( log_Z1ovZ0H - fh_unstable(DP_Z1,IL) + fh_unstable(DP_Z0H,IL) ) * ( TH1 - TH0 )
!    QstarUS = KARMAN / ( log_Z1ovZ0E - fh_unstable(DP_Z1,IL) + fh_unstable(DP_Z0E,IL) ) * ( Q1  - Q0  )

    fracu10us = ( log_10ovz0m - fm_unstable(10.0_dp,il) + fm_unstable(dp_z0m,il) ) &
              / ( log_z1ovz0m - fm_unstable(  dp_z1,il) + fm_unstable(dp_z0m,il) )
    fract2us  = ( log_02ovz0h - fm_unstable( 2.0_dp,il) + fm_unstable(dp_z0h,il) ) &
              / ( log_z1ovz0h - fm_unstable(  dp_z1,il) + fm_unstable(dp_z0h,il) )
    fracq2us  = ( log_02ovz0e - fm_unstable( 2.0_dp,il) + fm_unstable(dp_z0e,il) ) &
              / ( log_z1ovz0e - fm_unstable(  dp_z1,il) + fm_unstable(dp_z0e,il) )

    ! stable condition
    if ( bulkflux_use_mean ) then
       denom = log_z1ovz0m &
             - fmm_stable(dp_z1,il) + fmm_stable(dp_z0m,il) * dp_z0m / dp_z1 &
             + rzm * ( fm_stable(dp_z0m,il) - 1.0_dp )
       denoh = log_z1ovz0h &
             - fhm_stable(dp_z1,il) + fhm_stable(dp_z0h,il) * dp_z0h / dp_z1 &
             + rzh * ( fh_stable(dp_z0h,il) - 1.0_dp )
       denoe = log_z1ovz0e &
             - fhm_stable(dp_z1,il) + fhm_stable(dp_z0e,il) * dp_z0e / dp_z1 &
             + rze * ( fh_stable(dp_z0e,il) - 1.0_dp )
    else
       denom = log_z1ovz0m - fm_stable(dp_z1,il) + fm_stable(dp_z0m,il)
       denoh = log_z1ovz0h - fh_stable(dp_z1,il) + fh_stable(dp_z0h,il)
       denoe = log_z1ovz0e - fh_stable(dp_z1,il) + fh_stable(dp_z0e,il)
    end if
    uabss  = max( sqrt( u1**2 + v1**2 ), bulkflux_uabs_min )
    ustars = karman / denom * uabss
    tstars = karman / denoh * ( th1 - th0 )
    qstars = karman / denoe * ( q1  - q0  )
!    UstarS = KARMAN / ( log_Z1ovZ0M - fm_stable(DP_Z1,IL) + fm_stable(DP_Z0M,IL) ) * UabsS
!    TstarS = KARMAN / ( log_Z1ovZ0H - fh_stable(DP_Z1,IL) + fh_stable(DP_Z0H,IL) ) * ( TH1 - TH0 )
!    QstarS = KARMAN / ( log_Z1ovZ0E - fh_stable(DP_Z1,IL) + fh_stable(DP_Z0E,IL) ) * ( Q1  - Q0  )

    fracu10s = ( log_10ovz0m - fm_stable(10.0_dp,il) + fm_stable(dp_z0m,il) ) &
             / ( log_z1ovz0m - fm_stable(  dp_z1,il) + fm_stable(dp_z0m,il) )
    fract2s  = ( log_02ovz0h - fm_stable( 2.0_dp,il) + fm_stable(dp_z0h,il) ) &
             / ( log_z1ovz0h - fm_stable(  dp_z1,il) + fm_stable(dp_z0h,il) )
    fracq2s  = ( log_02ovz0e - fm_stable( 2.0_dp,il) + fm_stable(dp_z0e,il) ) &
             / ( log_z1ovz0e - fm_stable(  dp_z1,il) + fm_stable(dp_z0e,il) )

    sw = 0.5_dp - sign( 0.5_dp, il ) ! if unstable, sw = 1

    ustarc   = ( sw ) * ustarus   + ( 1.0_dp-sw ) * ustars
    tstarc   = ( sw ) * tstarus   + ( 1.0_dp-sw ) * tstars
    qstarc   = ( sw ) * qstarus   + ( 1.0_dp-sw ) * qstars
    uabsc    = ( sw ) * uabsus    + ( 1.0_dp-sw ) * uabss
    fracu10c = ( sw ) * fracu10us + ( 1.0_dp-sw ) * fracu10s
    fract2c  = ( sw ) * fract2us  + ( 1.0_dp-sw ) * fract2s
    fracq2c  = ( sw ) * fracq2us  + ( 1.0_dp-sw ) * fracq2s

    ! revert to RP
    ustar   = real( ustarc,   kind=rp )
    tstar   = real( tstarc,   kind=rp )
    qstar   = real( qstarc,   kind=rp )
    uabs    = real( uabsc,    kind=rp )
    fracu10 = real( fracu10c, kind=rp )
    fract2  = real( fract2c,  kind=rp )
    fracq2  = real( fracq2c,  kind=rp )

    ra = max( ( q1 - q0 ) / real(UstarC * QstarC + EPS,kind=RP), eps )


    return
  end subroutine bulkflux_b91w01

  !-----------------------------------------------------------------------------
  ! stability function for momemtum in unstable condition
  function fm_unstable( Z, IL )
    use scale_const, only: &
         pi => const_pi
    implicit none

    ! argument
    real(DP), intent(in) :: Z
    real(DP), intent(in) :: IL

    ! function
    real(DP) :: fm_unstable

    ! Wilson (2001)
    real(DP), parameter :: gamma = 3.6_dp

    ! works
    real(DP) :: R
    real(DP) :: r4R
    !---------------------------------------------------------------------------

    r = min( z * il, 0.0_dp )

    if ( flag_w01 ) then
       ! Wilson (2001)
       fm_unstable = 3.0_dp * log( ( 1.0_dp + sqrt( 1.0_dp + gamma * (-r)**(2.0_dp/3.0_dp) ) ) * 0.5_dp )
    else
       ! Beljaars and Holtslag (1994), originally Paulson (1974) and Dyer (1974)
       r4r = sqrt( sqrt( 1.0_dp - 16.0_dp * r ) )
       fm_unstable = log( ( 1.0_dp + r4r )**2 * ( 1.0_dp + r4r * r4r ) * 0.125_dp ) - 2.0_dp * atan( r4r ) + pi * 0.5_dp
    end if

    return
  end function fm_unstable
  function fmm_unstable( Z, IL )
    use scale_const, only: &
      pi  => const_pi, &
      eps => const_eps
    implicit none

    ! argument
    real(DP), intent(in) :: Z
    real(DP), intent(in) :: IL

    ! function
    real(DP) :: fmm_unstable

    ! Wilson (2001)
    real(DP), parameter :: gamma = 3.6_dp

    ! works
    real(DP) :: R
    real(DP) :: f, r3
    real(DP) :: r4R, r2R
    !---------------------------------------------------------------------------

    r = min( z * il, 0.0_dp )

    if ( flag_w01 ) then
       ! Wilson (2001)
       r3 = (-r)**(1.0_dp/3.0_dp)
       if ( r > -eps ) then
          fmm_unstable = 9.0_dp / 20.0_dp * gamma * r3**2
       else
          f = sqrt( 1.0_dp + gamma * r3**2 )
          fmm_unstable = 3.0_dp * log( ( 1.0_dp + f ) * 0.5_dp ) &
                       + 1.5_dp / ( sqrt(gamma)**3 * r) * asinh( sqrt(gamma) * r3 ) &
                       + 1.5_dp * f / ( gamma * r3**2 ) &
                       - 1.0_dp
       end if
    else
       ! Beljaars and Holtslag (1994), originally Paulson (1974) and Dyer (1974)
       if ( r > -eps ) then ! |R|<EPS, now R < 0
          fmm_unstable = - 2.0_dp * r
       else
          r2r = sqrt( 1.0_dp - 16.0_dp * r )
          r4r = sqrt( r2r )
          fmm_unstable = log( ( 1.0_dp + r4r )**2 * ( 1.0_dp + r2r ) * 0.125_dp ) &
                        - 2.0_dp * atan( r4r ) &
                        + ( 1.0_dp - r4r*r2r ) / ( 12.0_dp * r ) &
                        + pi * 0.5_dp - 1.0_dp
       end if
    end if

    return
  end function fmm_unstable

  !-----------------------------------------------------------------------------
  ! stability function for heat/vapor in unstable condition
  function fh_unstable( Z, IL )
    implicit none

    ! argument
    real(DP), intent(in) :: Z
    real(DP), intent(in) :: IL

    ! function
    real(DP) :: fh_unstable

    ! Wilson (2001)
    real(DP), parameter :: Pt = 0.95_dp ! turbulent Prandtl number
    real(DP), parameter :: gamma = 7.9_dp

    ! works
    real(DP) :: R
    !---------------------------------------------------------------------------

    r = min( z * il, 0.0_dp )

    if ( flag_w01 ) then
       ! Wilson (2001)
       fh_unstable = 3.0_dp * log( ( 1.0_dp + sqrt( 1.0_dp + gamma * (-r)**(2.0_dp/3.0_dp) ) ) * 0.5_dp ) * pt &
                   + log( z ) * ( 1.0_dp - pt )
    else
       ! Beljaars and Holtslag (1994), originally Paulson (1974); Dyer (1974)
       fh_unstable = 2.0_dp * log( ( 1.0_dp + sqrt( 1.0_dp - 16.0_dp * r ) ) * 0.5_dp )
    end if

    return
  end function fh_unstable
  function fhm_unstable( Z, IL )
    use scale_const, only: &
      eps => const_eps
    implicit none

    ! argument
    real(DP), intent(in) :: Z
    real(DP), intent(in) :: IL

    ! function
    real(DP) :: fhm_unstable

    ! Wilson (2001)
    real(DP), parameter :: Pt = 0.95_dp ! turbulent Prandtl number
    real(DP), parameter :: gamma = 7.9_dp

    ! works
    real(DP) :: R
    real(DP) :: f, r3
    real(DP) :: r2R
    !---------------------------------------------------------------------------

    r = min( z * il, 0.0_dp )

    if ( flag_w01 ) then
       ! Wilson (2001)
       r3 = (-r)**(1.0_dp/3.0_dp)
       if ( r > -eps ) then
          fhm_unstable = 9.0_dp / 20.0_dp * gamma * r3**2
       else
          f = sqrt( 1.0_dp + gamma * r3**2 )
          fhm_unstable = 3.0_dp * log( ( 1.0_dp + f ) * 0.5_dp ) &
                       + 1.5_dp / ( sqrt(gamma)**3 * r) * asinh( sqrt(gamma) * r3 ) &
                       + 1.5_dp * f / ( gamma * r3**2 ) &
                       - 1.0_dp
          fhm_unstable = fhm_unstable * pt + ( log( z ) - 1.0_dp ) * ( 1.0_dp - pt )
       end if
    else
       ! Beljaars and Holtslag (1994), originally Paulson (1974); Dyer (1974)
       if ( r > -eps ) then ! |R| < EPS, now R < 0
          fhm_unstable = - 4.0_dp * r
       else
          r2r = sqrt( 1.0_dp - 16.0_dp * r )
          fhm_unstable = 2.0_dp * log( ( 1.0_dp + r2r ) * 0.5_dp ) &
                       + ( 1.0_dp - r2r ) / ( 8.0_dp * r ) &
                       - 1.0_dp
       end if
    end if

    return
  end function fhm_unstable

  !-----------------------------------------------------------------------------
  ! stability function for momemtum in stable condition
  function fm_stable( Z, IL )
    implicit none

    ! argument
    real(DP), intent(in) :: Z
    real(DP), intent(in) :: IL

    ! function
    real(DP) :: fm_stable

    ! parameters of stability functions (Beljaars and Holtslag 1991)
    real(DP), parameter :: a = 1.0_dp
    real(DP), parameter :: b = 0.667_dp
    real(DP), parameter :: c = 5.0_dp
    real(DP), parameter :: d = 0.35_dp

    ! works
    real(DP) :: R
    !---------------------------------------------------------------------------

    r = max( z * il, 0.0_dp )

    ! Holtslag and DeBruin (1988)
#if defined(PGI) || defined(SX)
    fm_stable = - a*r - b*( r - c/d )*exp( -min( d*r, 1.e+3_rp ) ) - b*c/d ! apply exp limiter
#else
    fm_stable = - a*r - b*( r - c/d )*exp( -d*r ) - b*c/d
#endif

    return
  end function fm_stable
  function fmm_stable( Z, IL )
    use scale_const, only: &
       eps => const_eps
    implicit none

    ! argument
    real(DP), intent(in) :: Z
    real(DP), intent(in) :: IL

    ! function
    real(DP) :: fmm_stable

    ! parameters of stability functions (Beljaars and Holtslag 1991)
    real(DP), parameter :: a = 1.0_dp
    real(DP), parameter :: b = 0.667_dp
    real(DP), parameter :: c = 5.0_dp
    real(DP), parameter :: d = 0.35_dp

    ! works
    real(DP) :: R
    !---------------------------------------------------------------------------

    r = max( z * il, 0.0_dp )

    ! Holtslag and DeBruin (1988)
    if ( r < eps ) then
       fmm_stable = - 0.5_dp * ( a + b * c + d ) * r
    else
       fmm_stable = b * ( d*r - c + 1.0_dp ) / ( d**2 * r ) &
#if defined(__PGI) || defined(__ES2)
    ! apply exp limiter
                    * exp( -min( d*r, 1.e+3_dp ) ) &
#else
                    * exp( -d*r ) &
#endif
                  - a * r * 0.5_dp &
                  - b * ( c*d*r - c + 1.0_dp ) / ( d**2 * r )
    end if

    return
  end function fmm_stable

  !-----------------------------------------------------------------------------
  ! stability function for heat/vapor in stable condition
  function fh_stable( Z, IL )
    implicit none

    ! argument
    real(DP), intent(in) :: Z
    real(DP), intent(in) :: IL

    ! function
    real(DP) :: fh_stable

    ! parameters of stability functions (Beljaars and Holtslag 1991)
    real(DP), parameter :: a = 1.0_dp
    real(DP), parameter :: b = 0.667_dp
    real(DP), parameter :: c = 5.0_dp
    real(DP), parameter :: d = 0.35_dp

    ! works
    real(DP) :: R
    !---------------------------------------------------------------------------

    r = max( z * il, 0.0_dp )

    ! Beljaars and Holtslag (1991)
#if defined(PGI) || defined(SX)
    fh_stable = 1.0_dp - ( 1.0_dp + 2.0_dp/3.0_dp * a*r )**1.5_dp - b*( r - c/d )*exp( -min( d*r, 1.e+3_rp ) ) - b*c/d ! apply exp limiter
#else
    fh_stable = 1.0_dp - ( 1.0_dp + 2.0_dp/3.0_dp * a*r )**1.5_dp - b*( r - c/d )*exp( -d*r ) - b*c/d
#endif

    return
  end function fh_stable
  function fhm_stable( Z, IL )
    use scale_const, only: &
       eps => const_eps
    implicit none

    ! argument
    real(DP), intent(in) :: Z
    real(DP), intent(in) :: IL

    ! function
    real(DP) :: fhm_stable

    ! parameters of stability functions (Beljaars and Holtslag 1991)
    real(DP), parameter :: a = 1.0_dp
    real(DP), parameter :: b = 0.667_dp
    real(DP), parameter :: c = 5.0_dp
    real(DP), parameter :: d = 0.35_dp

    ! works
    real(DP) :: R
    !---------------------------------------------------------------------------

    r = max( z * il, 0.0_dp )

    ! Beljaars and Holtslag (1991)
    if ( r < eps ) then
       fhm_stable = - 0.5_dp * ( a + b*c + b ) * r
    else
       fhm_stable = b * ( d*r - c + 1.0_dp ) / ( d**2 * r ) &
#if defined(__PGI) || defined(__ES2)
                    * exp( -min( d*r, 1.e+3_dp) ) &
#else
                    * exp( -d*r ) &
#endif
                  - 3.0_dp * sqrt( 1.0_dp + 2.0_dp*a*r/3.0_dp )**5 / ( 5.0_dp * a * r ) &
                  - b * ( c*d*r - c + 1.0_dp ) / ( d**2 * r ) &
                  + 0.6_rp / ( a * r ) + 1.0_dp
    end if

    return
  end function fhm_stable

end module scale_bulkflux