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

module sfx_sheba
    !< @brief SHEBA surface flux module

    ! modules used
    ! --------------------------------------------------------------------------------
#ifdef SFX_CHECK_NAN
    use sfx_common
#endif
    use sfx_data
    use sfx_surface
    use sfx_sheba_param

#if defined(INCLUDE_CXX)
    use iso_c_binding, only: C_LOC, C_PTR
    use C_FUNC
#endif
    ! --------------------------------------------------------------------------------

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

    ! public interface
    ! --------------------------------------------------------------------------------
    public :: get_surface_fluxes
    public :: get_surface_fluxes_vec
    ! --------------------------------------------------------------------------------

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

contains

#if defined(INCLUDE_CXX)
    subroutine set_c_struct_sfx_sheba_param_values(sfx_model_param)
        type (sfx_sheba_param), intent(inout) :: sfx_model_param
        sfx_model_param%kappa = kappa
        sfx_model_param%Pr_t_0_inv = Pr_t_0_inv

        sfx_model_param%alpha_m = alpha_m
        sfx_model_param%alpha_h = alpha_h
        sfx_model_param%a_m = a_m
        sfx_model_param%b_m = b_m
        sfx_model_param%a_h = a_h
        sfx_model_param%b_h = b_h
        sfx_model_param%c_h = c_h
    end subroutine set_c_struct_sfx_sheba_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), pointer :: meteo_c         !< meteorological data (input)
        type (sfxDataVecTypeC), pointer :: sfx_c             !< surface flux data (output)
        type(C_PTR) :: meteo_c_ptr, sfx_c_ptr
        type (sfx_sheba_param) :: model_param
        type (sfx_surface_param) :: surface_param
        type (sfx_sheba_numericsTypeC) :: numerics_c
        type (sfx_phys_constants) :: phys_constants

        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_sheba_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 get_surface_fluxes_sheba(sfx_c_ptr, meteo_c_ptr, model_param, surface_param, numerics_c, phys_constants, n)
#else
        do i = 1, n
            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))

            call get_surface_fluxes(sfx_cell, meteo_cell, numerics)

            call push_sfx_data(sfx, sfx_cell, i)
        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)
        ! ----------------------------------------------------------------------------

        ! --- 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 Udyn, Tdyn, Qdyn   !< dynamic scales

        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]

        integer surface_type    !< surface type = (ocean || land)


#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_m = meteo%z0_m

        ! --- define surface type
        if (z0_m < 0.0) then
            surface_type = surface_ocean
        else
            surface_type = surface_land
        end if

        if (surface_type == surface_ocean) then
            ! --- define surface roughness [momentum] & dynamic velocity in neutral conditions
            call get_charnock_roughness(z0_m, u_dyn0, U, h, numerics%maxiters_charnock)
            ! --- define relative height
            h0_m = h / z0_m
        endif
        if (surface_type == surface_land) then
            ! --- define relative height
            h0_m = h / z0_m
            ! --- define dynamic velocity in neutral conditions
            u_dyn0 = U * kappa / log(h0_m)
        end if

        ! --- define thermal roughness & B = log(z0_m / z0_h)
        Re = u_dyn0 * z0_m / nu_air
        call get_thermal_roughness(z0_t, B, z0_m, Re, surface_type)

        ! --- define relative height [thermal]
        h0_t = h / z0_t

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

        ! --- get the fluxes
        ! ----------------------------------------------------------------------------
        call get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta, &
                U, Tsemi, dT, dQ, h, z0_m, z0_t, (g / Tsemi), 10)
        ! ----------------------------------------------------------------------------

        call get_phi(phi_m, phi_h, zeta)
        ! ----------------------------------------------------------------------------

        ! --- define transfer coeff. (momentum) & (heat)
        Cm = 0.0
        if (U > 0.0) then
            Cm = Udyn / U
        end if
        Ct = 0.0
        if (abs(dT) > 0.0) then
            Ct = Tdyn / dT
        end if

        ! --- 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 = 0.0, &
                Cm = Cm, Ct = Ct, Km = Km, Pr_t_inv = Pr_t_inv)

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

    !< @brief get dynamic scales
    ! --------------------------------------------------------------------------------
    subroutine get_dynamic_scales(Udyn, Tdyn, Qdyn, zeta, &
            U, Tsemi, dT, dQ, z, z0_m, z0_t, beta, maxiters)
        ! ----------------------------------------------------------------------------
        real, intent(out) :: Udyn, Tdyn, Qdyn   !< dynamic scales
        real, intent(out) :: zeta               !< = z/L

        real, intent(in) :: U                   !< abs(wind speed) at z
        real, intent(in) :: Tsemi               !< semi-sum of temperature at z and at surface
        real, intent(in) :: dT, dQ              !< temperature & humidity difference between z and at surface
        real, intent(in) :: z                   !< constant flux layer height
        real, intent(in) :: z0_m, z0_t          !< roughness parameters
        real, intent(in) :: beta                !< buoyancy parameter

        integer, intent(in) :: maxiters         !< maximum number of iterations
        ! ----------------------------------------------------------------------------

        ! --- local variables
        real, parameter :: gamma = 0.61

        real :: psi_m, psi_h
        real :: psi0_m, psi0_h
        real :: Linv
        integer :: i
        ! ----------------------------------------------------------------------------


        Udyn = kappa * U / log(z / z0_m)
        Tdyn = kappa * dT * Pr_t_0_inv / log(z / z0_t)
        Qdyn = kappa * dQ * Pr_t_0_inv / log(z / z0_t)
        zeta = 0.0

        ! --- no wind
        if (Udyn < 1e-5) return

        Linv = kappa * beta * (Tdyn + gamma * Qdyn * Tsemi) / (Udyn * Udyn)
        zeta = z * Linv

        ! --- near neutral case
        if (Linv < 1e-5) return

        do i = 1, maxiters

            call get_psi(psi_m, psi_h, zeta)
            call get_psi_mh(psi0_m, psi0_h, z0_m * Linv, z0_t * Linv)

            Udyn = kappa * U / (log(z / z0_m) - (psi_m - psi0_m))
            Tdyn = kappa * dT * Pr_t_0_inv / (log(z / z0_t) - (psi_h - psi0_h))
            Qdyn = kappa * dQ * Pr_t_0_inv / (log(z / z0_t) - (psi_h - psi0_h))

            if (Udyn < 1e-5) exit

            Linv = kappa * beta * (Tdyn + gamma * Qdyn * Tsemi) / (Udyn * Udyn)
            zeta = z * Linv
        end do

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

    ! stability functions
    ! --------------------------------------------------------------------------------
    subroutine get_phi(phi_m, phi_h, zeta)
        !< @brief stability functions (momentum) & (heat): neutral case
        ! ----------------------------------------------------------------------------
        real, intent(out) :: phi_m, phi_h   !< stability functions

        real, intent(in) :: zeta            !< = z/L
        ! ----------------------------------------------------------------------------


        if (zeta >= 0.0) then
            phi_m = 1.0 + (a_m * zeta * (1.0 + zeta)**(1.0 / 3.0)) / (1.0 + b_m * zeta)
            phi_h = 1.0 + (a_h * zeta + b_h * zeta * zeta) / (1.0 + c_h * zeta + zeta * zeta)
        else
            phi_m = (1.0 - alpha_m * zeta)**(-0.25)
            phi_h = (1.0 - alpha_h * zeta)**(-0.5)
        end if

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

    ! universal functions
    ! --------------------------------------------------------------------------------
    subroutine get_psi(psi_m, psi_h, zeta)
        !< @brief universal functions (momentum) & (heat): neutral case
        ! ----------------------------------------------------------------------------
        real, intent(out) :: psi_m, psi_h   !< universal functions

        real, intent(in) :: zeta            !< = z/L
        ! ----------------------------------------------------------------------------

        ! --- local variables
        real :: x_m, x_h
        real :: q_m, q_h
        ! ----------------------------------------------------------------------------


        if (zeta >= 0.0) then

            q_m = ((1.0 - b_m) / b_m)**(1.0 / 3.0)
            q_h = sqrt(c_h * c_h - 4.0)

            x_m = (1.0 + zeta)**(1.0 / 3.0)
            x_h = zeta

            psi_m = -3.0 * (a_m / b_m) * (x_m - 1.0) + 0.5 * (a_m / b_m) * q_m * (&
                    2.0 * log((x_m + q_m) / (1.0 + q_m)) - &
                            log((x_m * x_m - x_m * q_m + q_m * q_m) / (1.0 - q_m + q_m * q_m)) + &
                            2.0 * sqrt(3.0) * (&
                                    atan((2.0 * x_m - q_m) / (sqrt(3.0) * q_m)) - &
                                            atan((2.0 - q_m) / (sqrt(3.0) * q_m))))

            psi_h = -0.5 * b_h * log(1.0 + c_h * x_h + x_h * x_h) + &
                    ((-a_h / q_h) + ((b_h * c_h) / (2.0 * q_h))) * (&
                            log((2.0 * x_h + c_h - q_h) / (2.0 * x_h + c_h + q_h)) - &
                                    log((c_h - q_h) / (c_h + q_h)))
        else
            x_m = (1.0 - alpha_m * zeta)**(0.25)
            x_h = (1.0 - alpha_h * zeta)**(0.25)

            psi_m = (4.0 * atan(1.0) / 2.0) + 2.0 * log(0.5 * (1.0 + x_m)) + log(0.5 * (1.0 + x_m * x_m)) - 2.0 * atan(x_m)
            psi_h = 2.0 * log(0.5 * (1.0 + x_h * x_h))
        end if

    end subroutine


    subroutine get_psi_mh(psi_m, psi_h, zeta_m, zeta_h)
        !< @brief universal functions (momentum) & (heat): neutral case
        ! ----------------------------------------------------------------------------
        real, intent(out) :: psi_m, psi_h   !< universal functions

        real, intent(in) :: zeta_m, zeta_h  !< = z/L
        ! ----------------------------------------------------------------------------

        ! --- local variables
        real :: x_m, x_h
        real :: q_m, q_h
        ! ----------------------------------------------------------------------------


        if (zeta_m >= 0.0) then
            q_m = ((1.0 - b_m) / b_m)**(1.0 / 3.0)
            x_m = (1.0 + zeta_m)**(1.0 / 3.0)

            psi_m = -3.0 * (a_m / b_m) * (x_m - 1.0) + 0.5 * (a_m / b_m) * q_m * (&
                    2.0 * log((x_m + q_m) / (1.0 + q_m)) - &
                            log((x_m * x_m - x_m * q_m + q_m * q_m) / (1.0 - q_m + q_m * q_m)) + &
                            2.0 * sqrt(3.0) * (&
                                    atan((2.0 * x_m - q_m) / (sqrt(3.0) * q_m)) - &
                                            atan((2.0 - q_m) / (sqrt(3.0) * q_m))))
        else
            x_m = (1.0 - alpha_m * zeta_m)**(0.25)
            psi_m = (4.0 * atan(1.0) / 2.0) + 2.0 * log(0.5 * (1.0 + x_m)) + log(0.5 * (1.0 + x_m * x_m)) - 2.0 * atan(x_m)
        end if

        if (zeta_h >= 0.0) then
            q_h = sqrt(c_h * c_h - 4.0)
            x_h = zeta_h

            psi_h = -0.5 * b_h * log(1.0 + c_h * x_h + x_h * x_h) + &
                    ((-a_h / q_h) + ((b_h * c_h) / (2.0 * q_h))) * (&
                            log((2.0 * x_h + c_h - q_h) / (2.0 * x_h + c_h + q_h)) - &
                                    log((c_h - q_h) / (c_h + q_h)))
        else
            x_h = (1.0 - alpha_h * zeta_h)**(0.25)
            psi_h = 2.0 * log(0.5 * (1.0 + x_h * x_h))
        end if

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

end module sfx_sheba