most_cpu.f90

Go to the documentation of this file.

! Copyright (c) 2026, The Neko Authors
! All rights reserved.
!
! Redistribution and use in source and binary forms, with or without
! modification, are permitted provided that the following conditions
! are met:
!
!   * Redistributions of source code must retain the above copyright
!     notice, this list of conditions and the following disclaimer.
!
!   * Redistributions in binary form must reproduce the above
!     copyright notice, this list of conditions and the following
!     disclaimer in the documentation and/or other materials provided
!     with the distribution.
!
!   * Neither the name of the authors nor the names of its
!     contributors may be used to endorse or promote products derived
!     from this software without specific prior written permission.
!
! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
! FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
! COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
! INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
! BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
! LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
! CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
! LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN
! ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
! POSSIBILITY OF SUCH DAMAGE.
!
module most_cpu
  use num_types, only : rp
  use utils, only : neko_error, neko_warning
  use logger, only : log_size, neko_log
  use math, only : glsum, glmin, glmax
  implicit none
  private
 
  public :: most_compute_cpu
 
  abstract interface
     function slaw_m_interface(z, L_ob, z0) result(slaw)
       import rp
       real(kind=rp), intent(in) :: z, l_ob, z0
       real(kind=rp) :: slaw
     end function slaw_m_interface
 
     function slaw_h_interface(z, L_ob, z0h) result(slaw)
       import rp
       real(kind=rp), intent(in) :: z, l_ob, z0h
       real(kind=rp) :: slaw
     end function slaw_h_interface
 
     function corr_m_interface(z, L_ob) result(corr)
       import rp
       real(kind=rp), intent(in) :: z, l_ob
       real(kind=rp) :: corr
     end function corr_m_interface
 
     function corr_h_interface(z, L_ob) result(corr)
       import rp
       real(kind=rp), intent(in) :: z, l_ob
       real(kind=rp) :: corr
     end function corr_h_interface
 
     function f_interface(Ri_b, z, z0, z0h, Pr, L_ob, slaw_m, slaw_h) result(f)
       import rp, slaw_m_interface, slaw_h_interface
       real(kind=rp), intent(in) :: ri_b, z, z0, z0h, l_ob, pr
       real(kind=rp) :: f
       procedure(slaw_m_interface) :: slaw_m
       procedure(slaw_h_interface) :: slaw_h
     end function f_interface
 
     function dfdl_interface(l_upper, l_lower, z, z0, z0h, Pr, L_ob,&
          slaw_m, slaw_h, fd_h) result(dfdl)
       import rp, slaw_m_interface, slaw_h_interface
       real(kind=rp), intent(in) :: l_upper, l_lower, z, z0, z0h, l_ob, fd_h, pr
       real(kind=rp) :: dfdl
       procedure(slaw_m_interface) :: slaw_m
       procedure(slaw_h_interface) :: slaw_h
     end function dfdl_interface
  end interface
 
 
  ! These will point to the correct functions
  ! depending on stability regime and bc_type.
  !! @note These are reassigned per node inside the compute loop, so every
  !! thread needs its own copy. They are always set before being used.
  procedure(slaw_m_interface), pointer :: slaw_m_ptr => null()
  procedure(slaw_h_interface), pointer :: slaw_h_ptr => null()
  procedure(corr_m_interface), pointer :: corr_m_ptr => null()
  procedure(corr_h_interface), pointer :: corr_h_ptr => null()
  procedure(f_interface), pointer :: f_ptr => null()
  procedure(dfdl_interface), pointer :: dfdl_ptr => null()
  !$omp threadprivate(slaw_m_ptr, slaw_h_ptr, corr_m_ptr, corr_h_ptr, &
  !$omp& f_ptr, dfdl_ptr)
 
contains
 
  subroutine select_bc_operators(bc_type, bc_value, q, ts, ti, kappa, &
       utau, z0h, hi, Pr)
    character(len=*), intent(in) :: bc_type
    real(kind=rp), intent(in) :: hi, ti, kappa, utau, z0h, bc_value, pr
    real(kind=rp), intent(inout) :: q,ts
    select case (bc_type)
    case ("neumann")
       ! ts not used
       q = bc_value
       f_ptr => f_neumann
       dfdl_ptr => dfdl_neumann
    case ("dirichlet")
       ts = bc_value
       q = (kappa/pr)*utau*(ts - ti)/log(hi/z0h)
       f_ptr => f_dirichlet
       dfdl_ptr => dfdl_dirichlet
    case default
       call neko_error("Invalid specified temperature b.c. type " // &
            "('neumann' or 'dirichlet'?)")
    end select
  end subroutine select_bc_operators
 
  subroutine compute_ri_b(bc_type, g_dot_n, hi, ti, ts, magu, kappa, q, &
       Pr, Ri_b)
    character(len=*), intent(in) :: bc_type
    real(kind=rp), intent(in) :: hi, ti, ts, pr
    real(kind=rp), intent(in) :: magu, kappa, g_dot_n
    real(kind=rp), intent(inout) :: q, ri_b
 
    select case (bc_type)
    case ("neumann")
       ri_b = - g_dot_n*hi / ti*q*pr / (magu**3*kappa**2)
    case ("dirichlet")
       ri_b = g_dot_n*hi/ti*(ti - ts)/magu**2
    case default
       call neko_error("Invalid specified temperature b.c. type " // &
            "('neumann' or 'dirichlet'?)")
    end select
  end subroutine compute_ri_b
 
  subroutine set_stability_regime(Ri_b, Ri_threshold)
    real(kind=rp), intent(in) :: ri_b, ri_threshold
 
    if (ri_b > ri_threshold) then
       slaw_m_ptr => slaw_m_stable
       slaw_h_ptr => slaw_h_stable
       corr_m_ptr => corr_m_stable
       corr_h_ptr => corr_h_stable
    elseif (ri_b < -ri_threshold) then
       slaw_m_ptr => slaw_m_convective
       slaw_h_ptr => slaw_h_convective
       corr_m_ptr => corr_m_convective
       corr_h_ptr => corr_h_convective
    else
       slaw_m_ptr => slaw_m_neutral
       slaw_h_ptr => slaw_h_neutral
    end if
  end subroutine set_stability_regime
 
  subroutine most_compute_cpu(u, v, w, temp, temp_w, &
       n_x, n_y, n_z, h, tau_x, tau_y, tau_z, n_nodes, &
       kappa, mu_w, rho_w, g_vec, Pr, z0, z0h_in, bc_type, bc_value, tstep, &
       Ri_b_diagn, L_ob_diagn, utau_diagn, magu_diagn, ti_diagn, ts_diagn,&
       q_diagn)
    integer, intent(in) :: n_nodes, tstep
    real(kind=rp), dimension(n_nodes), intent(in) :: u, v, w, temp, temp_w
    real(kind=rp), dimension(n_nodes), intent(in) :: n_x, n_y, n_z, h
    real(kind=rp), intent(in) :: kappa, z0, z0h_in, bc_value, pr
    real(kind=rp), dimension(3), intent(in) :: g_vec
    real(kind=rp), dimension(n_nodes), intent(in) :: mu_w, rho_w
    real(kind=rp) :: g_dot_n
    character(len=*), intent(in) :: bc_type
    real(kind=rp), dimension(n_nodes), intent(inout) :: tau_x, tau_y, tau_z
    real(kind=rp) :: ui, vi, wi, hi, rho, mu
    real(kind=rp) :: normu, z0h
    real(kind=rp) :: l_upper, l_lower, l_old
    real(kind=rp) :: f, dfdl, fd_h, l_new, l_sign
    integer :: i, count
    integer, parameter :: max_count = 50
    real(kind=rp), parameter :: tol = 0.001_rp
    real(kind=rp), parameter :: nr_step = 0.001_rp
    real(kind=rp), parameter :: ri_threshold = 0.0001_rp
    character(len=LOG_SIZE) :: log_buf
    real(kind=rp) :: utau, ri_b, l_ob, magu, q, ti, ts
    real(kind=rp), dimension(n_nodes), intent(inout) :: ri_b_diagn, l_ob_diagn
    real(kind=rp), dimension(n_nodes), intent(inout) :: utau_diagn, magu_diagn
    real(kind=rp), dimension(n_nodes), intent(inout) :: ti_diagn, ts_diagn
    real(kind=rp), dimension(n_nodes), intent(inout) :: q_diagn
 
    !$omp parallel do private(i, ui, vi, wi, hi, rho, mu, normu, z0h, &
    !$omp& L_upper, L_lower, L_old, f, dfdl, fd_h, L_new, L_sign, count, &
    !$omp& utau, Ri_b, L_ob, magu, q, ti, ts, g_dot_n)
    do i=1, n_nodes
       ! Sample the variables
       ui = u(i)
       vi = v(i)
       wi = w(i)
       ti = temp(i)
       hi = h(i)
       rho = rho_w(i)
       mu = mu_w(i)
 
       ! Project on horizontal directions
       normu = ui * n_x(i) + vi * n_y(i) + wi * n_z(i)
       ui = ui - normu * n_x(i)
       vi = vi - normu * n_y(i)
       wi = wi - normu * n_z(i)
 
       ! Compute velocity magnitude
       magu = sqrt(ui**2 + vi**2 + wi**2)
       magu = max(magu, 1.0e-6_rp)
       utau = magu*kappa / log(hi/z0)
 
       ! Compute thermal roughness length from Zilitinkevich, 1995
       if (z0h_in < 0) then
          ! z0h_in is interpreted as -C_Zil (Zilitinkevich constant) for z0h
          z0h = z0 * exp(z0h_in*sqrt((utau*z0)/(mu/rho)))
       else
          z0h = z0h_in
       end if
 
       ! Get q, Ri_b, f_ptr, dfdl_ptr based on bc_type
       ! Maybe redundant, but needed to initialise Rib
       call select_bc_operators(bc_type, bc_value, q, ts, ti, kappa, &
            utau, z0h, hi, pr)
 
       ! Compute g along the normal (generalisation for hills and similar)
       g_dot_n = abs(g_vec(1)*n_x(i) + g_vec(2)*n_y(i) + g_vec(3)*n_z(i))
 
       ! Compute Richardson and set stability accordingly
       call compute_ri_b(bc_type, g_dot_n, hi, ti, ts, magu, kappa, q, pr, ri_b)
       call set_stability_regime(ri_b, ri_threshold)
 
       if (abs((ri_b)) <= ri_threshold) then
          ! Neutral (L_ob undefined)
          l_ob = 0.0_rp
       else
          ! Determine target regime sign
          if ((ri_b) > 0.0_rp) then
             l_ob = hi / max(ri_b, ri_threshold) ! Stable guess
             l_sign = 1.0_rp
          else
             l_ob = hi / min(ri_b, -ri_threshold) ! Convective guess
             l_sign = -1.0_rp
          end if
 
          l_old = 1.0e10_rp
          count = 0
 
          ! Find Obukhov length
          do while ((abs(l_old - l_ob) / abs(l_ob) > tol) .and. &
               (count < max_count))
             ! Switch between stable and convective based on bulk
             ! Richardson (Ri_b)
             l_old = l_ob
             count = count + 1
             fd_h = nr_step*l_ob
             l_upper = l_ob + fd_h
             l_lower = l_ob - fd_h
 
             ! Compute L_ob based on stability and bc_type
             f = f_ptr(ri_b, hi, z0, z0h, pr, l_ob, slaw_m_ptr, slaw_h_ptr)
             dfdl = dfdl_ptr(l_upper, l_lower, hi, z0, z0h, pr, l_ob, &
                  slaw_m_ptr, slaw_h_ptr, fd_h)
             if (abs(dfdl) < 1.0e-12_rp) then
                call neko_error("Division by zero in dfdl")
             end if
             l_new = l_ob - f/dfdl
             ! Avoid regime crossing during Newton iter (otherwise crash)
             if (l_new*l_sign <= 0.0_rp) then
                ! "damp update" (stay on same side)
                l_new = 0.5_rp * l_ob
             end if
             ! Bound L_ob
             l_ob = sign(max(abs(l_new), 1.0e-8_rp), l_sign)
             l_ob = sign(min(abs(l_ob), 1.0e8_rp), l_sign)
          end do
 
          if (abs(l_ob) > 5e5_rp .or. abs(l_ob) < 1e-6_rp) then
             count = max_count
             !$omp critical
             call neko_warning("Obukhov length did not converge " // &
                  "(MOST wall model)")
             !$omp end critical
          end if
 
          if (.not. associated(f_ptr) .or. .not. associated(dfdl_ptr)) then
             call neko_error("Unassociated pointer for f or dfdl")
          end if
       end if
 
       ! Based on stability and bc_type, compute utau/q
       select case (bc_type)
       case ("neumann")
          ! Compute u* with the new Obukhov length
          utau = kappa*magu/slaw_m_ptr(hi, l_ob, z0)
       case ("dirichlet")
          ! Compute u* with the new Obukhov length
          utau = kappa*magu/slaw_m_ptr(hi, l_ob, z0)
          ! and compute q from here
          q = kappa/pr*utau*(ts - ti)/slaw_h_ptr(hi, l_ob, z0h)
       case default
          call neko_error("Invalid specified temperature b.c. type " // &
               "('neumann' or 'dirichlet'?)")
       end select
 
       ! Distribute according to the velocity vector and bound magu
       ! to avoid 0 division
       magu = max(magu, 1.0e-6_rp)
       tau_x(i) = -rho * utau**2 * ui / magu
       tau_y(i) = -rho * utau**2 * vi / magu
       tau_z(i) = -rho * utau**2 * wi / magu
 
       ri_b_diagn(i) = ri_b
       l_ob_diagn(i) = l_ob
       utau_diagn(i) = utau
       magu_diagn(i) = magu
       ti_diagn(i) = ti
       ts_diagn(i) = temp_w(i)
       q_diagn(i) = q
    end do
    !$omp end parallel do
  end subroutine most_compute_cpu
 
  function slaw_m_stable(z, L_ob, z0) result(slaw)
    real(kind=rp), intent(in) :: z, l_ob, z0
    real(kind=rp) :: slaw
 
    slaw = log(z/z0)-corr_m_stable(z,l_ob)+corr_m_stable(z0,l_ob)
  end function slaw_m_stable
 
  function slaw_h_stable(z,L_ob,z0h) result(slaw)
    real(kind=rp), intent(in) :: z,l_ob,z0h
    real(kind=rp) :: slaw
 
    slaw = log(z/z0h)-corr_h_stable(z,l_ob)+corr_h_stable(z0h,l_ob)
  end function slaw_h_stable
 
  function corr_m_stable(z, L_ob) result(corr)
    real(kind=rp), intent(in) :: z, l_ob
    real(kind=rp) :: corr
    real(kind=rp) :: a, b, c, d
    real(kind=rp) :: zeta
    zeta = z/l_ob
    ! Coefficients specific to Cheng & Brutsaert (2005)
    a = 1.0_rp
    b = 2.0_rp/3.0_rp
    c = 5.0_rp
    d = 0.35_rp
    corr = - a*zeta - b*(zeta-c/d)*exp(-d*zeta) - b*c/d
  end function corr_m_stable
 
  function corr_h_stable(z, L_ob) result(corr)
    real(kind=rp), intent(in) :: z, l_ob
    real(kind=rp) :: corr
    real(kind=rp) :: a, b, c, d
    real(kind=rp) :: zeta
 
    zeta = z/l_ob
    ! Coefficients specific to Cheng & Brutsaert (2005)
    a = 1.0_rp
    b = 2.0_rp/3.0_rp
    c = 5.0_rp
    d = 0.35_rp
    corr = -b * (zeta - c / d) * exp(-d * zeta) - &
         (1.0_rp + 2.0_rp / 3.0_rp * a * zeta)**1.5_rp - b * c / d + &
         1.0_rp
  end function corr_h_stable
 
  function slaw_m_convective(z, L_ob, z0) result(slaw)
    real(kind=rp), intent(in) :: z, l_ob, z0
    real(kind=rp) :: slaw
 
    slaw = log(z/z0) - corr_m_convective(z, l_ob) + corr_m_convective(z0, l_ob)
  end function slaw_m_convective
 
  function slaw_h_convective(z, L_ob, z0h) result(slaw)
    real(kind=rp), intent(in) :: z, l_ob, z0h
    real(kind=rp) :: slaw
 
    slaw = log(z / z0h) - corr_h_convective(z, l_ob) + &
         corr_h_convective(z0h, l_ob)
  end function slaw_h_convective
 
  function corr_m_convective(z, L_ob) result(corr)
    real(kind=rp), intent(in) :: z, l_ob
    real(kind=rp) :: xi, pi, zeta
    real(kind=rp) :: corr
 
    zeta = z/l_ob
    pi = 4*atan(1.0_rp)
    ! Standard Dyer-Businger coefficient gamma = 16.0
    xi = (1.0_rp - 16.0_rp*zeta)**0.25_rp
    corr = 2*log(0.5_rp*(1 + xi)) + log(0.5_rp*(1 + xi**2)) - 2*atan(xi) + pi/2
  end function corr_m_convective
 
  function corr_h_convective(z, L_ob) result(corr)
    real(kind=rp), intent(in) :: z, l_ob
    real(kind=rp) :: zeta, pi, xi
    real(kind=rp) :: corr
 
    zeta = z/l_ob
    pi = 4*atan(1.0_rp)
    ! Standard Dyer-Businger coefficient gamma = 16.0
    xi = (1.0_rp - 16.0_rp*zeta)**0.25_rp
    corr = 2*log(0.5_rp*(1 + xi**2))
  end function corr_h_convective
 
  function slaw_m_neutral(z, L_ob, z0) result(slaw)
    real(kind=rp), intent(in) :: z, l_ob, z0
    real(kind=rp) :: slaw
 
    slaw = log(z/z0)
  end function slaw_m_neutral
 
  function slaw_h_neutral(z, L_ob, z0h) result(slaw)
    real(kind=rp), intent(in) :: z, l_ob, z0h
    real(kind=rp) :: slaw
 
    slaw = log(z/z0h)
  end function slaw_h_neutral
 
  function f_neumann(Ri_b, z, z0, z0h, Pr, L_ob, slaw_m, slaw_h) result(f)
    real(kind=rp), intent(in) :: ri_b, z, z0, z0h, l_ob, pr
    procedure(slaw_m_interface) :: slaw_m
    procedure(slaw_h_interface) :: slaw_h
    real(kind=rp) :: f
 
    f = (ri_b - pr*z/l_ob/slaw_m(z, l_ob, z0)**3)
  end function f_neumann
 
  function dfdl_neumann(l_upper, l_lower, z, z0, z0h, Pr, L_ob, &
       slaw_m, slaw_h, fd_h) result(dfdl)
    real(kind=rp), intent(in) :: l_upper, l_lower, z, z0, z0h, l_ob, fd_h, pr
    procedure(slaw_m_interface) :: slaw_m
    procedure(slaw_h_interface) :: slaw_h
    real(kind=rp) :: dfdl
 
    dfdl = (-z/l_upper/slaw_m(z, l_upper, z0)**3) ! conv
    dfdl = dfdl + (z/l_lower/slaw_m(z, l_lower, z0)**3) ! conv
    dfdl = pr*dfdl/(2*fd_h)
  end function dfdl_neumann
 
  function f_dirichlet(Ri_b, z, z0, z0h, Pr, L_ob, slaw_m, slaw_h) result(f)
    real(kind=rp), intent(in) :: ri_b, z, z0, z0h, l_ob, pr
    procedure(slaw_m_interface) :: slaw_m
    procedure(slaw_h_interface) :: slaw_h
    real(kind=rp) :: f
 
    f = (ri_b - pr*z/l_ob*slaw_h(z, l_ob, z0h)/slaw_m(z, l_ob, z0)**2) ! conv
  end function f_dirichlet
 
  function dfdl_dirichlet(l_upper, l_lower, z, z0, z0h, Pr, L_ob, &
       slaw_m, slaw_h, fd_h) result(dfdl)
    real(kind=rp), intent(in) :: l_upper, l_lower, z, z0, z0h, l_ob, fd_h, pr
    procedure(slaw_m_interface) :: slaw_m
    procedure(slaw_h_interface) :: slaw_h
    real(kind=rp) :: dfdl
 
    dfdl = (-z / l_upper * slaw_h(z, l_upper, z0h) / &
         slaw_m(z, l_upper, z0)**2) ! conv
    dfdl = dfdl + (z / l_lower * slaw_h(z, l_lower, z0h) / &
         slaw_m(z, l_lower, z0)**2) ! conv
    dfdl = pr*dfdl/(2*fd_h)
  end function dfdl_dirichlet
 
 
end module most_cpu