sfx_esm.f90

#include "../includeF/sfx_def.fi"

module sfx_esm
    !< @brief main Earth System Model surface flux module

    ! modules used
    ! --------------------------------------------------------------------------------
#ifdef SFX_CHECK_NAN
    use sfx_common
#endif
    use sfx_data
    use sfx_surface
    use sfx_esm_param
#if defined(INCLUDE_CXX)
    use iso_c_binding, only: C_LOC, C_PTR, C_INT, C_FLOAT
    use C_FUNC
#endif
    ! --------------------------------------------------------------------------------

    ! directives list
    ! --------------------------------------------------------------------------------
    implicit none
    private
    ! --------------------------------------------------------------------------------

    ! public interface
    ! --------------------------------------------------------------------------------
    public :: get_surface_fluxes
    public :: get_surface_fluxes_vec
    integer z0m_id
    integer z0t_id
    ! --------------------------------------------------------------------------------

    ! --------------------------------------------------------------------------------
    type, public :: numericsType
        integer :: maxiters_convection = 10    !< maximum (actual) number of iterations in convection
        integer :: maxiters_charnock = 10      !< maximum (actual) number of iterations in charnock roughness
    end type
    ! --------------------------------------------------------------------------------

#if defined(INCLUDE_CXX)
    type, BIND(C), public :: sfx_esm_param_C 
        real(C_FLOAT) :: kappa           
        real(C_FLOAT) :: Pr_t_0_inv
        real(C_FLOAT) :: Pr_t_inf_inv

        real(C_FLOAT) :: alpha_m           
        real(C_FLOAT) :: alpha_h
        real(C_FLOAT) :: alpha_h_fix
        real(C_FLOAT) :: beta_m
        real(C_FLOAT) :: beta_h
        real(C_FLOAT) :: Rib_max
    end type

    type, BIND(C), public :: sfx_esm_numericsType_C 
        integer(C_INT) :: maxiters_convection           
        integer(C_INT) :: maxiters_charnock 
    end type

    INTERFACE
        SUBROUTINE c_esm_compute_flux(sfx, meteo, model_param, surface_param, numerics, constants, grid_size) BIND(C, & 
            name="c_esm_compute_flux")
            use sfx_data
            use, intrinsic :: ISO_C_BINDING, ONLY: C_INT, C_PTR
            Import :: sfx_esm_param_C, sfx_esm_numericsType_C
            implicit none
            integer(C_INT) :: grid_size
            type(C_PTR), value :: sfx
            type(C_PTR), value :: meteo
            type(sfx_esm_param_C) :: model_param
            type(sfx_surface_param) :: surface_param
            type(sfx_esm_numericsType_C) :: numerics
            type(sfx_phys_constants) :: constants
        END SUBROUTINE c_esm_compute_flux
    END INTERFACE
#endif 

contains

    ! --------------------------------------------------------------------------------
#if defined(INCLUDE_CXX)
    subroutine set_c_struct_sfx_esm_param_values(sfx_model_param)
        type (sfx_esm_param_C), intent(inout) :: sfx_model_param
        sfx_model_param%kappa = kappa
        sfx_model_param%Pr_t_0_inv = Pr_t_0_inv
        sfx_model_param%Pr_t_inf_inv = Pr_t_inf_inv

        sfx_model_param%alpha_m = alpha_m
        sfx_model_param%alpha_h = alpha_h
        sfx_model_param%alpha_h_fix = alpha_h_fix
        sfx_model_param%beta_m = beta_m
        sfx_model_param%beta_h = beta_h
        sfx_model_param%Rib_max = Rib_max
    end subroutine set_c_struct_sfx_esm_param_values
#endif

    subroutine get_surface_fluxes_vec(sfx, meteo, numerics, n)
        !< @brief surface flux calculation for array data
        !< @details contains C/C++ & CUDA interface
        ! ----------------------------------------------------------------------------
        type (sfxDataVecType), intent(inout) :: sfx

        type (meteoDataVecType), intent(in) :: meteo
        type (numericsType), intent(in) :: numerics
        integer, intent(in) :: n
        ! ----------------------------------------------------------------------------

        ! --- local variables
        type (meteoDataType)  meteo_cell
        type (sfxDataType) sfx_cell
        integer i
        ! ----------------------------------------------------------------------------
#if defined(INCLUDE_CXX)
        type (meteoDataVecTypeC), target :: meteo_c         !< meteorological data (input)
        type (sfxDataVecTypeC), target :: sfx_c             !< surface flux data (output)
        type(C_PTR) :: meteo_c_ptr, sfx_c_ptr
        type (sfx_esm_param_C) :: model_param
        type (sfx_surface_param) :: surface_param
        type (sfx_esm_numericsType_C) :: numerics_c
        type (sfx_phys_constants) :: phys_constants

        numerics_c%maxiters_convection = numerics%maxiters_convection
        numerics_c%maxiters_charnock = numerics%maxiters_charnock

        phys_constants%Pr_m = Pr_m;
        phys_constants%nu_air = nu_air;
        phys_constants%g = g;

        call set_c_struct_sfx_esm_param_values(model_param)
        call set_c_struct_sfx_surface_param_values(surface_param)
        call set_meteo_vec_c(meteo, meteo_c)
        call set_sfx_vec_c(sfx, sfx_c)
        meteo_c_ptr = C_LOC(meteo_c)
        sfx_c_ptr   = C_LOC(sfx_c)
        
        call c_esm_compute_flux(sfx_c_ptr, meteo_c_ptr, model_param, surface_param, numerics_c, phys_constants, n)
#else
        do i = 1, n
#ifdef SFX_FORCE_DEPRECATED_ESM_CODE
#else
            meteo_cell = meteoDataType(&
                    h = meteo%h(i), &
                    U = meteo%U(i), dT = meteo%dT(i), Tsemi = meteo%Tsemi(i), dQ = meteo%dQ(i), &
                    z0_m = meteo%z0_m(i), depth=meteo%depth(i), lai=meteo%lai(i), surface_type=meteo%surface_type(i))

            call get_surface_fluxes(sfx_cell, meteo_cell, numerics)

            call push_sfx_data(sfx, sfx_cell, i)
#endif
        end do
#endif

    end subroutine get_surface_fluxes_vec
    ! --------------------------------------------------------------------------------

    ! --------------------------------------------------------------------------------
    subroutine get_surface_fluxes(sfx, meteo, numerics)
        !< @brief surface flux calculation for single cell
        !< @details contains C/C++ interface
        ! ----------------------------------------------------------------------------
#ifdef SFX_CHECK_NAN
        use ieee_arithmetic
#endif

        type (sfxDataType), intent(out) :: sfx

        type (meteoDataType), intent(in) :: meteo
        type (numericsType), intent(in) :: numerics
        ! ----------------------------------------------------------------------------

        ! --- meteo derived datatype name shadowing
        ! ----------------------------------------------------------------------------
        real :: h       !< constant flux layer height [m]
        real :: U       !< abs(wind speed) at 'h' [m/s]
        real :: dT      !< difference between potential temperature at 'h' and at surface [K]
        real :: Tsemi   !< semi-sum of potential temperature at 'h' and at surface [K]
        real :: dQ      !< difference between humidity at 'h' and at surface [g/g]
        real :: z0_m    !< surface aerodynamic roughness (should be < 0 for water bodies surface)
        real :: depth   
        real :: lai
        integer :: surface_type
        ! ----------------------------------------------------------------------------

        ! --- local variables
        ! ----------------------------------------------------------------------------
        real z0_t               !< thermal roughness [m]
        real B                  !< = ln(z0_m / z0_t) [n/d]
        real h0_m, h0_t         !< = h / z0_m, h / z0_h [n/d]

        real u_dyn0             !< dynamic velocity in neutral conditions [m/s]
        real Re                 !< roughness Reynolds number = u_dyn0 * z0_m / nu [n/d]

        real zeta               !< = z/L [n/d]
        real Rib                !< bulk Richardson number

        real zeta_conv_lim      !< z/L critical value for matching free convection limit [n/d]
        real Rib_conv_lim       !< Ri-bulk critical value for matching free convection limit [n/d]

        real f_m_conv_lim       !< stability function (momentum) value in free convection limit [n/d]
        real f_h_conv_lim       !< stability function (heat) value in free convection limit [n/d]

        real psi_m, psi_h       !< universal functions (momentum) & (heat) [n/d]
        real phi_m, phi_h       !< stability functions (momentum) & (heat) [n/d]

        real Km                 !< eddy viscosity coeff. at h [m^2/s]
        real Pr_t_inv           !< invese Prandt number [n/d]

        real Cm, Ct             !< transfer coeff. for (momentum) & (heat) [n/d]


        real fval               !< just a shortcut for partial calculations
        real z0_m1
#ifdef SFX_CHECK_NAN
        real NaN
#endif
        ! ----------------------------------------------------------------------------

#ifdef SFX_CHECK_NAN
        ! --- checking if arguments are finite
        if (.not.(is_finite(meteo%U).and.is_finite(meteo%Tsemi).and.is_finite(meteo%dT).and.is_finite(meteo%dQ) &
                .and.is_finite(meteo%z0_m).and.is_finite(meteo%h))) then

            NaN = ieee_value(0.0, ieee_quiet_nan)   ! setting NaN
            sfx = sfxDataType(zeta = NaN, Rib = NaN, &
                    Re = NaN, B = NaN, z0_m = NaN, z0_t = NaN, &
                    Rib_conv_lim = NaN, &
                    Cm = NaN, Ct = NaN, Km = NaN, Pr_t_inv = NaN)
            return
        end if
#endif

        ! --- shadowing names for clarity
        U = meteo%U
        Tsemi = meteo%Tsemi
        dT = meteo%dT
        dQ = meteo%dQ
        h = meteo%h
        z0_m1 = meteo%z0_m
        depth = meteo%depth
        lai = meteo%lai
        surface_type=meteo%surface_type
       
        
        call get_dynamic_roughness_definition(surface_type, ocean_z0m_id, land_z0m_id, lake_z0m_id, snow_z0m_id, &
        forest_z0m_id, usersf_z0m_id, ice_z0m_id, z0m_id)
        
        
        call get_dynamic_roughness_all(z0_m, u_dyn0, U, depth, h, numerics%maxiters_charnock, z0_m1, z0m_id)
        
        call get_thermal_roughness_definition(surface_type, ocean_z0t_id, land_z0t_id, lake_z0t_id, snow_z0t_id, &
        forest_z0t_id, usersf_z0t_id, ice_z0t_id, z0t_id)
        
        
        Re = u_dyn0 * z0_m / nu_air
        
        call get_thermal_roughness_all(z0_t, B, z0_m, Re, u_dyn0, lai, z0t_id)
       
       ! --- define relative height
            h0_m = h / z0_m
        ! --- define relative height [thermal]
        h0_t = h / z0_t

        ! --- define Ri-bulk
        Rib = (g / Tsemi) * h * (dT + 0.61e0 * Tsemi * dQ) / U**2

        ! --- define free convection transition zeta = z/L value
        call get_convection_lim(zeta_conv_lim, Rib_conv_lim, f_m_conv_lim, f_h_conv_lim, &
                h0_m, h0_t, B)

        ! --- get the fluxes
        ! ----------------------------------------------------------------------------
        if (Rib > 0.0) then
            ! --- stable stratification block

            !   --- restrict bulk Ri value
            !   *: note that this value is written to output
            Rib = min(Rib, Rib_max)
            call get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B)

            fval = beta_m * zeta
            phi_m = 1.0 + fval
            phi_h = 1.0/Pr_t_0_inv + fval

        else if (Rib < Rib_conv_lim) then
            ! --- strong instability block

            call get_psi_convection(psi_m, psi_h, zeta, Rib, &
                    zeta_conv_lim, f_m_conv_lim, f_h_conv_lim, h0_m, h0_t, B, numerics%maxiters_convection)

            fval = (zeta_conv_lim / zeta)**(1.0/3.0)
            phi_m = fval / f_m_conv_lim
            phi_h = fval / (Pr_t_0_inv * f_h_conv_lim)

        else if (Rib > -0.001) then
            ! --- nearly neutral [-0.001, 0] block

            call get_psi_neutral(psi_m, psi_h, h0_m, h0_t, B)

            zeta = 0.0
            phi_m = 1.0
            phi_h = 1.0 / Pr_t_0_inv
        else
            ! --- weak & semistrong instability block

            call get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, numerics%maxiters_convection)

            phi_m = (1.0 - alpha_m * zeta)**(-0.25)
            phi_h = 1.0 / (Pr_t_0_inv * sqrt(1.0 - alpha_h_fix * zeta))
        end if
        ! ----------------------------------------------------------------------------

        ! --- define transfer coeff. (momentum) & (heat)
        Cm = kappa / psi_m
        Ct = kappa / psi_h

        ! --- define eddy viscosity & inverse Prandtl number
        Km = kappa * Cm * U * h / phi_m
        Pr_t_inv = phi_m / phi_h

        ! --- setting output
        sfx = sfxDataType(zeta = zeta, Rib = Rib, &
            Re = Re, B = B, z0_m = z0_m, z0_t = z0_t, &
            Rib_conv_lim = Rib_conv_lim, &
            Cm = Cm, Ct = Ct, Km = Km, Pr_t_inv = Pr_t_inv)

    end subroutine get_surface_fluxes
    ! --------------------------------------------------------------------------------

    ! convection universal functions shortcuts
    ! --------------------------------------------------------------------------------
    function f_m_conv(zeta)
        ! ----------------------------------------------------------------------------
        real :: f_m_conv
        real, intent(in) :: zeta
        ! ----------------------------------------------------------------------------

        f_m_conv = (1.0 - alpha_m * zeta)**0.25

    end function f_m_conv

    function f_h_conv(zeta)
        ! ----------------------------------------------------------------------------
        real :: f_h_conv
        real, intent(in) :: zeta
        ! ----------------------------------------------------------------------------

        f_h_conv = (1.0 - alpha_h * zeta)**0.5

    end function f_h_conv
    ! --------------------------------------------------------------------------------

    ! universal functions
    ! --------------------------------------------------------------------------------
    subroutine get_psi_neutral(psi_m, psi_h, h0_m, h0_t, B)
        !< @brief universal functions (momentum) & (heat): neutral case
        ! ----------------------------------------------------------------------------
        real, intent(out) :: psi_m, psi_h   !< universal functions

        real, intent(in) :: h0_m, h0_t      !< = z/z0_m, z/z0_h
        real, intent(in) :: B               !< = log(z0_m / z0_h)
        ! ----------------------------------------------------------------------------

        psi_m = log(h0_m)
        psi_h = log(h0_t) / Pr_t_0_inv
        !*: this looks redundant z0_t = z0_m in case |B| ~ 0
        if (abs(B) < 1.0e-10) psi_h = psi_m / Pr_t_0_inv

    end subroutine

    subroutine get_psi_stable(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B)
        !< @brief universal functions (momentum) & (heat): stable case
        ! ----------------------------------------------------------------------------
        real, intent(out) :: psi_m, psi_h   !< universal functions [n/d]
        real, intent(out) :: zeta           !< = z/L [n/d]

        real, intent(in) :: Rib             !< bulk Richardson number [n/d]
        real, intent(in) :: h0_m, h0_t      !< = z/z0_m, z/z0_h [n/d]
        real, intent(in) :: B               !< = log(z0_m / z0_h) [n/d]
        ! ----------------------------------------------------------------------------

        ! --- local variables
        real :: Rib_coeff
        real :: psi0_m, psi0_h
        real :: phi
        real :: c
        ! ----------------------------------------------------------------------------

        psi0_m = alog(h0_m)
        psi0_h = B / psi0_m

        Rib_coeff = beta_m * Rib
        c = (psi0_h + 1.0) / Pr_t_0_inv - 2.0 * Rib_coeff
        zeta = psi0_m * (sqrt(c**2 + 4.0 * Rib_coeff * (1.0 - Rib_coeff)) - c) / &
                (2.0 * beta_m * (1.0 - Rib_coeff))

        phi = beta_m * zeta

        psi_m = psi0_m + phi
        psi_h = (psi0_m + B) / Pr_t_0_inv + phi

    end subroutine

    subroutine get_psi_semi_convection(psi_m, psi_h, zeta, Rib, h0_m, h0_t, B, maxiters)
        !< @brief universal functions (momentum) & (heat): semi-strong convection case
        ! ----------------------------------------------------------------------------
        real, intent(out) :: psi_m, psi_h   !< universal functions [n/d]
        real, intent(out) :: zeta           !< = z/L [n/d]

        real, intent(in) :: Rib             !< bulk Richardson number [n/d]
        real, intent(in) :: h0_m, h0_t      !< = z/z0_m, z/z0_h [n/d]
        real, intent(in) :: B               !< = log(z0_m / z0_h) [n/d]
        integer, intent(in) :: maxiters     !< maximum number of iterations

        ! --- local variables
        real :: zeta0_m, zeta0_h
        real :: f0_m, f0_h
        real :: f_m, f_h

        integer :: i
        ! ----------------------------------------------------------------------------

        psi_m = log(h0_m)
        psi_h = log(h0_t)
        if (abs(B) < 1.0e-10) psi_h = psi_m

        zeta = Rib * Pr_t_0_inv * psi_m**2 / psi_h

        do i = 1, maxiters + 1
            zeta0_m = zeta / h0_m
            zeta0_h = zeta / h0_t
            if (abs(B) < 1.0e-10) zeta0_h = zeta0_m

            f_m = (1.0 - alpha_m * zeta)**0.25e0
            f_h = sqrt(1.0 - alpha_h_fix * zeta)

            f0_m = (1.0 - alpha_m * zeta0_m)**0.25e0
            f0_h = sqrt(1.0 - alpha_h_fix * zeta0_h)

            f0_m = max(f0_m, 1.000001e0)
            f0_h = max(f0_h, 1.000001e0)

            psi_m = log((f_m - 1.0e0)*(f0_m + 1.0e0)/((f_m + 1.0e0)*(f0_m - 1.0e0))) + 2.0e0*(atan(f_m) - atan(f0_m))
            psi_h = log((f_h - 1.0e0)*(f0_h + 1.0e0)/((f_h + 1.0e0)*(f0_h - 1.0e0))) / Pr_t_0_inv

            ! *: don't update zeta = z/L at last iteration
            if (i == maxiters + 1) exit

            zeta = Rib * psi_m**2 / psi_h
        end do

    end subroutine

    subroutine get_psi_convection(psi_m, psi_h, zeta, Rib, &
            zeta_conv_lim, f_m_conv_lim, f_h_conv_lim, &
            h0_m, h0_t, B, maxiters)
        !< @brief universal functions (momentum) & (heat): fully convective case
        ! ----------------------------------------------------------------------------
        real, intent(out) :: psi_m, psi_h               !< universal functions [n/d]
        real, intent(out) :: zeta                       !< = z/L [n/d]

        real, intent(in) :: Rib                         !< bulk Richardson number [n/d]
        real, intent(in) :: h0_m, h0_t                  !< = z/z0_m, z/z0_h [n/d]
        real, intent(in) :: B                           !< = log(z0_m / z0_h) [n/d]
        integer, intent(in) :: maxiters                 !< maximum number of iterations

        real, intent(in) :: zeta_conv_lim               !< convective limit zeta
        real, intent(in) :: f_m_conv_lim, f_h_conv_lim  !< universal function shortcuts in limiting case
        ! ----------------------------------------------------------------------------

        ! --- local variables
        real :: zeta0_m, zeta0_h
        real :: f0_m, f0_h
        real :: p_m, p_h
        real :: a_m, a_h
        real :: c_lim, f

        integer :: i
        ! ----------------------------------------------------------------------------

        p_m = 2.0 * atan(f_m_conv_lim) + log((f_m_conv_lim - 1.0) / (f_m_conv_lim + 1.0))
        p_h = log((f_h_conv_lim - 1.0) / (f_h_conv_lim + 1.0))

        zeta = zeta_conv_lim
        do i = 1, maxiters + 1
            zeta0_m = zeta / h0_m
            zeta0_h = zeta / h0_t
            if (abs(B) < 1.0e-10) zeta0_h = zeta0_m

            f0_m = (1.0 - alpha_m * zeta0_m)**0.25
            f0_h = sqrt(1.0 - alpha_h_fix * zeta0_h)

            a_m = -2.0*atan(f0_m) + log((f0_m + 1.0)/(f0_m - 1.0))
            a_h = log((f0_h + 1.0)/(f0_h - 1.0))

            c_lim = (zeta_conv_lim / zeta)**(1.0 / 3.0)
            f = 3.0 * (1.0 - c_lim)

            psi_m = f / f_m_conv_lim + p_m + a_m
            psi_h = (f / f_h_conv_lim + p_h + a_h) / Pr_t_0_inv

            ! *: don't update zeta = z/L at last iteration
            if (i == maxiters + 1) exit

            zeta = Rib * psi_m**2 / psi_h
        end do

    end subroutine
    ! --------------------------------------------------------------------------------

    ! convection limit definition
    ! --------------------------------------------------------------------------------
    subroutine get_convection_lim(zeta_lim, Rib_lim, f_m_lim, f_h_lim, &
            h0_m, h0_t, B)
        ! ----------------------------------------------------------------------------
        real, intent(out) :: zeta_lim           !< limiting value of z/L
        real, intent(out) :: Rib_lim            !< limiting value of Ri-bulk
        real, intent(out) :: f_m_lim, f_h_lim   !< limiting values of universal functions shortcuts

        real, intent(in) :: h0_m, h0_t          !< = z/z0_m, z/z0_h [n/d]
        real, intent(in) :: B                   !< = log(z0_m / z0_h) [n/d]
        ! ----------------------------------------------------------------------------

        ! --- local variables
        real :: psi_m, psi_h
        real :: f_m, f_h
        real :: c
        ! ----------------------------------------------------------------------------

        ! --- define limiting value of zeta = z / L
        c = (Pr_t_inf_inv / Pr_t_0_inv)**4
        zeta_lim = (2.0 * alpha_h - c * alpha_m - &
                sqrt((c * alpha_m)**2 + 4.0 * c * alpha_h * (alpha_h - alpha_m))) / (2.0 * alpha_h**2)

        f_m_lim = f_m_conv(zeta_lim)
        f_h_lim = f_h_conv(zeta_lim)

        ! --- universal functions
        f_m = zeta_lim / h0_m
        f_h = zeta_lim / h0_t
        if (abs(B) < 1.0e-10) f_h = f_m

        f_m = (1.0 - alpha_m * f_m)**0.25
        f_h = sqrt(1.0 - alpha_h_fix * f_h)

        psi_m = 2.0 * (atan(f_m_lim) - atan(f_m)) + alog((f_m_lim - 1.0) * (f_m + 1.0)/((f_m_lim + 1.0) * (f_m - 1.0)))
        psi_h = alog((f_h_lim - 1.0) * (f_h + 1.0)/((f_h_lim + 1.0) * (f_h - 1.0))) / Pr_t_0_inv

        ! --- bulk Richardson number
        Rib_lim = zeta_lim * psi_h / (psi_m * psi_m)

    end subroutine
    ! --------------------------------------------------------------------------------

end module sfx_esm