scale_bulkflux.F90

Go to the documentation of this file.

!-------------------------------------------------------------------------------
!-------------------------------------------------------------------------------
module scale_bulkflux
  !-----------------------------------------------------------------------------
  !
  !++ used modules
  !
  use scale_precision
  use scale_stdio
  !-----------------------------------------------------------------------------
  implicit none
  private
  !-----------------------------------------------------------------------------
  !
  !++ Public procedure
  !
  public :: bulkflux_setup

  abstract interface
     subroutine bc( &
          Ustar,   & ! (out)
          tstar,   & ! (out)
          qstar,   & ! (out)
          uabs,    & ! (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) :: 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'

  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]

contains

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

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

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

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '++++++ Module[BULKFLUX] / Categ[COUPLER] / Origin[SCALElib]'

    ! 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
       if( io_l ) write(io_fid_log,*) '*** Not found namelist. Default used.'
    elseif( ierr > 0 ) then !--- fatal error
       write(*,*) 'xxx Not appropriate names in namelist PARAM_BULKFLUX. Check!'
       call prc_mpistop
    endif
    if( io_nml ) write(io_fid_nml,nml=param_bulkflux)

    if( io_l ) write(io_fid_log,*)
    if( io_l ) write(io_fid_log,*) '*** Scheme for surface bulk flux : ', trim(bulkflux_type)
    select case(bulkflux_type)
    case('U95')
       if( io_l ) write(io_fid_log,*) '*** => Uno et al.(1995)'
       bulkflux => bulkflux_u95
    case('B91W01')
       if( io_l ) write(io_fid_log,*) '*** => Beljaars (1991) and Wilson (2001)'
       bulkflux => bulkflux_b91w01
    case default
       write(*,*) 'xxx Unsupported BULKFLUX_type. STOP'
       call prc_mpistop
    end select

    return
  end subroutine bulkflux_setup

  !-----------------------------------------------------------------------------
  ! ref. Uno et al. (1995)
  !-----------------------------------------------------------------------------
  subroutine bulkflux_u95( &
      Ustar,   & ! (out)
      tstar,   & ! (out)
      qstar,   & ! (out)
      uabs,    & ! (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) :: 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 * ( pre00 / p1 )**( rdry / cpdry )
    th0  = t0 * ( pre00 / p0 )**( rdry / cpdry )

    ! 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 )

    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
  !
  !-----------------------------------------------------------------------------
  subroutine bulkflux_b91w01( &
      Ustar,   & ! (out)
      tstar,   & ! (out)
      qstar,   & ! (out)
      uabs,    & ! (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) :: 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
    !---------------------------------------------------------------------------

    ! convert to DP
    dp_z1  = real( z1,  kind=dp )
    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 * ( pre00 / p1 )**( rdry / cpdry )
    th0 = t0 * ( pre00 / p0 )**( rdry / cpdry )
    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

    ! 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
      ! unstable condition
      uabsus  = max( sqrt( u1**2 + v1**2 + (bulkflux_wscf*wstar)**2 ), real( BULKFLUX_Uabs_min, kind=DP ) )
      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
      uabss  = max( sqrt( u1**2 + v1**2 ), bulkflux_uabs_min )
      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
      uabsus  = max( sqrt( u1**2 + v1**2 + (bulkflux_wscf*wstar)**2 ), real( BULKFLUX_Uabs_min, kind=DP ) )
      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
      uabss  = max( sqrt( u1**2 + v1**2 ), bulkflux_uabs_min )
      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
      duabsus  = max( sqrt( u1**2 + v1**2 + (bulkflux_wscf*dwstar)**2 ), real( BULKFLUX_Uabs_min, kind=DP ) )
      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
      duabss  = max( sqrt( u1**2 + v1**2 ), bulkflux_uabs_min )
      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
      uabsus  = max( sqrt( u1**2 + v1**2 + (bulkflux_wscf*wstar)**2 ), real( BULKFLUX_Uabs_min, kind=DP ) )
      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
      uabss  = max( sqrt( u1**2 + v1**2 ), bulkflux_uabs_min )
      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
    uabsus  = max( sqrt( u1**2 + v1**2 + (bulkflux_wscf*wstar)**2 ), real( BULKFLUX_Uabs_min, kind=DP ) )
    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
    uabss  = max( sqrt( u1**2 + v1**2 ), bulkflux_uabs_min )
    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 )

    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

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

    r = min( z * il, 0.0_dp )

    ! Wilson (2001)
    fm_unstable = 3.0_dp * log( ( 1.0_dp + sqrt( 1.0_dp + 3.6_dp * (-r)**(2.0_dp/3.0_dp) ) ) * 0.5_dp )

    ! If you want to run with the original Beljaars scheme (Beljaars and Holtslag 1994),
    ! you should comment out the above line (Wilson 2001) and uncomment the below lines (Paulson 1974; Dyer 1974).
    !
    !! Paulson (1974); Dyer (1974)
    !r4R = ( 1.0_DP - 16.0_DP * R )**0.25_DP
    !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

    return
  end function fm_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

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

    r = min( z * il, 0.0_dp )

    ! Wilson (2001)
    fh_unstable = 3.0_dp * log( ( 1.0_dp + sqrt( 1.0_dp + 7.9_dp * (-r)**(2.0_dp/3.0_dp) ) ) * 0.5_dp )
    fh_unstable = pt * fh_unstable  + ( 1.0_dp - pt ) * log( z )

    ! If you want to run with the original Beljaars scheme (Beljaars and Holtslag 1994),
    ! you should comment out the above line (Wilson 2001) and uncomment the below lines (Paulson 1974; Dyer 1974).
    !
    !! Paulson (1974); Dyer (1974)
    !fh_unstable = 2.0_DP * log( ( 1.0_DP + sqrt( 1.0_DP - 16.0_DP * R ) ) * 0.5_DP )

    return
  end function fh_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

  !-----------------------------------------------------------------------------
  ! 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

end module scale_bulkflux