sfx_data.f90

!> @brief surface flux model module data
module sfx_data

    ! modules used
    ! --------------------------------------------------------------------------------
    ! --------------------------------------------------------------------------------
    use iso_c_binding, only: C_FLOAT, C_INT, C_PTR, C_LOC
    ! directives list
    ! --------------------------------------------------------------------------------
    implicit none
    private
    ! --------------------------------------------------------------------------------

    ! public interface
    ! --------------------------------------------------------------------------------
    public :: allocate_meteo_vec, deallocate_meteo_vec
#if defined(INCLUDE_CXX)
    public :: set_meteo_vec_c
#endif
    public :: allocate_sfx_vec, deallocate_sfx_vec
#if defined(INCLUDE_CXX)
    public :: set_sfx_vec_c
#endif

    public :: push_sfx_data
    ! --------------------------------------------------------------------------------

    ! --------------------------------------------------------------------------------
    !> @brief meteorological input for surface flux calculation
    type, public :: meteoDataType

        real(C_FLOAT) :: h       !< constant flux layer height [m]
        real(C_FLOAT) :: U       !< abs(wind speed) at 'h' [m/s]
        real(C_FLOAT) :: dT      !< difference between potential temperature at 'h' and at surface [K]
        real(C_FLOAT) :: Tsemi   !< semi-sum of potential temperature at 'h' and at surface [K]
        real(C_FLOAT) :: dQ      !< difference between humidity at 'h' and at surface [g/g]
        real(C_FLOAT) :: z0_m    !< surface aerodynamic roughness (should be < 0 for water bodies surface)
    end type

    !> @brief meteorological input for surface flux calculation
    !> &details using arrays as input
    type, public :: meteoDataVecType
#if defined(INCLUDE_CXX)
        real, pointer :: h(:)       !< constant flux layer height [m]
        real, pointer :: U(:)       !< abs(wind speed) at 'h' [m/s]
        real, pointer :: dT(:)      !< difference between potential temperature at 'h' and at surface [K]
        real, pointer :: Tsemi(:)   !< semi-sum of potential temperature at 'h' and at surface [K]
        real, pointer :: dQ(:)      !< difference between humidity at 'h' and at surface [g/g]
        real, pointer :: z0_m(:)    !< surface aerodynamic roughness (should be < 0 for water bodies surface)
#else
        real, allocatable :: h(:)       !< constant flux layer height [m]
        real, allocatable :: U(:)       !< abs(wind speed) at 'h' [m/s]
        real, allocatable :: dT(:)      !< difference between potential temperature at 'h' and at surface [K]
        real, allocatable :: Tsemi(:)   !< semi-sum of potential temperature at 'h' and at surface [K]
        real, allocatable :: dQ(:)      !< difference between humidity at 'h' and at surface [g/g]
        real, allocatable :: z0_m(:)    !< surface aerodynamic roughness (should be < 0 for water bodies surface)
#endif
    end type

#if defined(INCLUDE_CXX)
    type, public :: meteoDataVecTypeC

        type(C_PTR) :: h       !< constant flux layer height [m]
        type(C_PTR) :: U       !< abs(wind speed) at 'h' [m/s]
        type(C_PTR) :: dT      !< difference between potential temperature at 'h' and at surface [K]
        type(C_PTR) :: Tsemi   !< semi-sum of potential temperature at 'h' and at surface [K]
        type(C_PTR) :: dQ      !< difference between humidity at 'h' and at surface [g/g]
        type(C_PTR) :: z0_m    !< surface aerodynamic roughness (should be < 0 for water bodies surface)
    end type
#endif
    ! --------------------------------------------------------------------------------

    ! --------------------------------------------------------------------------------
    !> @brief surface flux output data
    type, public :: sfxDataType
        
        real(C_FLOAT) :: zeta            !< = z/L [n/d]
        real(C_FLOAT) :: Rib             !< bulk Richardson number [n/d]
        real(C_FLOAT) :: Re              !< Reynolds number [n/d]
        real(C_FLOAT) :: B               !< = log(z0_m / z0_h) [n/d]
        real(C_FLOAT) :: z0_m            !< aerodynamic roughness [m]
        real(C_FLOAT) :: z0_t            !< thermal roughness [m]
        real(C_FLOAT) :: Rib_conv_lim    !< Ri-bulk convection critical value [n/d]
        real(C_FLOAT) :: Cm              !< transfer coefficient for momentum [n/d]
        real(C_FLOAT) :: Ct              !< transfer coefficient for heat [n/d]
        real(C_FLOAT) :: Km              !< eddy viscosity coeff. at h [m^2/s]
        real(C_FLOAT) :: Pr_t_inv        !< inverse turbulent Prandtl number at h [n/d]
    end type

    !> @brief surface flux output data
    !> &details using arrays as output
    type, public :: sfxDataVecType
#if defined(INCLUDE_CXX)
        real, pointer :: zeta(:)            !< = z/L [n/d]
        real, pointer :: Rib(:)             !< bulk Richardson number [n/d]
        real, pointer :: Re(:)              !< Reynolds number [n/d]
        real, pointer :: B(:)               !< = log(z0_m / z0_h) [n/d]
        real, pointer :: z0_m(:)            !< aerodynamic roughness [m]
        real, pointer :: z0_t(:)            !< thermal roughness [m]
        real, pointer :: Rib_conv_lim(:)    !< Ri-bulk convection critical value [n/d]
        real, pointer :: Cm(:)              !< transfer coefficient for momentum [n/d]
        real, pointer :: Ct(:)              !< transfer coefficient for heat [n/d]
        real, pointer :: Km(:)              !< eddy viscosity coeff. at h [m^2/s]
        real, pointer :: Pr_t_inv(:)        !< inverse turbulent Prandtl number at h [n/d]
#else
        real, allocatable :: zeta(:)            !< = z/L [n/d]
        real, allocatable :: Rib(:)             !< bulk Richardson number [n/d]
        real, allocatable :: Re(:)              !< Reynolds number [n/d]
        real, allocatable :: B(:)               !< = log(z0_m / z0_h) [n/d]
        real, allocatable :: z0_m(:)            !< aerodynamic roughness [m]
        real, allocatable :: z0_t(:)            !< thermal roughness [m]
        real, allocatable :: Rib_conv_lim(:)    !< Ri-bulk convection critical value [n/d]
        real, allocatable :: Cm(:)              !< transfer coefficient for momentum [n/d]
        real, allocatable :: Ct(:)              !< transfer coefficient for heat [n/d]
        real, allocatable :: Km(:)              !< eddy viscosity coeff. at h [m^2/s]
        real, allocatable :: Pr_t_inv(:)        !< inverse turbulent Prandtl number at h [n/d]
#endif
    end type

#if defined(INCLUDE_CXX)
    type, public :: sfxDataVecTypeC
        type(C_PTR) :: zeta            !< = z/L [n/d]
        type(C_PTR) :: Rib             !< bulk Richardson number [n/d]
        type(C_PTR) :: Re              !< Reynolds number [n/d]
        type(C_PTR) :: B               !< = log(z0_m / z0_h) [n/d]
        type(C_PTR) :: z0_m            !< aerodynamic roughness [m]
        type(C_PTR) :: z0_t            !< thermal roughness [m]
        type(C_PTR) :: Rib_conv_lim    !< Ri-bulk convection critical value [n/d]
        type(C_PTR) :: Cm              !< transfer coefficient for momentum [n/d]
        type(C_PTR) :: Ct              !< transfer coefficient for heat [n/d]
        type(C_PTR) :: Km              !< eddy viscosity coeff. at h [m^2/s]
        type(C_PTR) :: Pr_t_inv        !< inverse turbulent Prandtl number at h [n/d]
    end type

    type, BIND(C), public :: sfx_surface_param 
        integer(C_INT) :: surface_ocean           
        integer(C_INT) :: surface_land
        integer(C_INT) :: surface_lake

        real(C_FLOAT)  :: gamma_c;
        real(C_FLOAT)  :: Re_visc_min;
        real(C_FLOAT)  :: h_charnock;
        real(C_FLOAT)  :: c1_charnock;
        real(C_FLOAT)  :: c2_charnock;
        real(C_FLOAT)  :: Re_rough_min;
        real(C_FLOAT)  :: B1_rough;
        real(C_FLOAT)  :: B2_rough;
        real(C_FLOAT)  :: B3_rough;
        real(C_FLOAT)  :: B4_rough;
        real(C_FLOAT)  :: B_max_lake;
        real(C_FLOAT)  :: B_max_ocean;
        real(C_FLOAT)  :: B_max_land;
    end type

    type, BIND(C), public :: sfx_phys_constants 
        real(C_FLOAT)  :: Pr_m;
        real(C_FLOAT)  :: g;
        real(C_FLOAT)  :: nu_air;
    end type

    ! ESM strcuctures
    type, BIND(C), public :: sfx_esm_param 
        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_numericsTypeC 
        integer(C_INT) :: maxiters_convection           
        integer(C_INT) :: maxiters_charnock 
    end type

    ! SHEBA strcuctures
    type, BIND(C), public :: sfx_sheba_param 
        real(C_FLOAT) :: kappa           
        real(C_FLOAT) :: Pr_t_0_inv

        real(C_FLOAT) :: alpha_m           
        real(C_FLOAT) :: alpha_h
        real(C_FLOAT) :: a_m
        real(C_FLOAT) :: b_m
        real(C_FLOAT) :: a_h
        real(C_FLOAT) :: b_h
        real(C_FLOAT) :: c_h
    end type

    type, BIND(C), public :: sfx_sheba_numericsTypeC 
        integer(C_INT) :: maxiters_charnock 
    end type
#endif
    ! --------------------------------------------------------------------------------

contains

    ! --------------------------------------------------------------------------------
    subroutine allocate_meteo_vec(meteo, n)
        !> @brief allocate meteo data vector
        ! ----------------------------------------------------------------------------
        type (meteoDataVecType), intent(inout) :: meteo

        integer, intent(in) :: n
        ! ----------------------------------------------------------------------------

        allocate(meteo%h(n))
        allocate(meteo%U(n))
        allocate(meteo%dT(n))
        allocate(meteo%Tsemi(n))
        allocate(meteo%dQ(n))
        allocate(meteo%z0_m(n))

    end subroutine allocate_meteo_vec

#if defined(INCLUDE_CXX)
    subroutine set_meteo_vec_c(meteo, meteo_C)
        !> @brief allocate meteo data vector
        ! ----------------------------------------------------------------------------
        type (meteoDataVecType)  , intent(in) :: meteo
        type (meteoDataVecTypeC), pointer, intent(inout) :: meteo_C

        allocate(meteo_C)

        meteo_C%h = c_loc(meteo%h)
        meteo_C%U = c_loc(meteo%U)
        meteo_C%dT = c_loc(meteo%dT)
        meteo_C%Tsemi = c_loc(meteo%Tsemi)
        meteo_C%dQ = c_loc(meteo%dQ)
        meteo_C%z0_m = c_loc(meteo%z0_m)
    end subroutine set_meteo_vec_c
#endif
    ! --------------------------------------------------------------------------------

    ! --------------------------------------------------------------------------------
    subroutine deallocate_meteo_vec(meteo)
        !> @brief deallocate meteo data vector
        ! ----------------------------------------------------------------------------
        type (meteoDataVecType), intent(inout) :: meteo
        ! ----------------------------------------------------------------------------

        deallocate(meteo%h)
        deallocate(meteo%U)
        deallocate(meteo%dT)
        deallocate(meteo%Tsemi)
        deallocate(meteo%dQ)
        deallocate(meteo%z0_m)

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

    ! --------------------------------------------------------------------------------
    subroutine allocate_sfx_vec(sfx, n)
        !> @brief allocate surface fluxes data vector
        ! ----------------------------------------------------------------------------
        type (sfxDataVecType), intent(inout) :: sfx

        integer, intent(in) :: n
        ! ----------------------------------------------------------------------------

        allocate(sfx%zeta(n))
        allocate(sfx%Rib(n))
        allocate(sfx%Re(n))
        allocate(sfx%B(n))
        allocate(sfx%z0_m(n))
        allocate(sfx%z0_t(n))
        allocate(sfx%Rib_conv_lim(n))
        allocate(sfx%Cm(n))
        allocate(sfx%Ct(n))
        allocate(sfx%Km(n))
        allocate(sfx%Pr_t_inv(n))

    end subroutine allocate_sfx_vec

#if defined(INCLUDE_CXX)
    subroutine set_sfx_vec_c(sfx, sfx_C)
        !> @brief allocate surface fluxes data vector
        ! ----------------------------------------------------------------------------
        type (sfxDataVecType), intent(inout) :: sfx
        type (sfxDataVecTypeC), pointer, intent(inout) :: sfx_C

        allocate(sfx_C)

        sfx_C%zeta = c_loc(sfx%zeta)
        sfx_C%Rib = c_loc(sfx%Rib)
        sfx_C%Re = c_loc(sfx%Re)
        sfx_C%B = c_loc(sfx%B)
        sfx_C%z0_m = c_loc(sfx%z0_m)
        sfx_C%z0_t = c_loc(sfx%z0_t)
        sfx_C%Rib_conv_lim = c_loc(sfx%Rib_conv_lim)
        sfx_C%Cm = c_loc(sfx%Cm)
        sfx_C%Ct = c_loc(sfx%Ct)
        sfx_C%Km = c_loc(sfx%Km)
        sfx_C%Pr_t_inv = c_loc(sfx%Pr_t_inv)
    end subroutine set_sfx_vec_c
#endif
    ! --------------------------------------------------------------------------------

    ! --------------------------------------------------------------------------------
    !> @brief deallocate surface fluxes data vector
    subroutine deallocate_sfx_vec(sfx)
        ! ----------------------------------------------------------------------------
        type (sfxDataVecType), intent(inout) :: sfx
        ! ----------------------------------------------------------------------------

        deallocate(sfx%zeta)
        deallocate(sfx%Rib)
        deallocate(sfx%Re)
        deallocate(sfx%B)
        deallocate(sfx%z0_m)
        deallocate(sfx%z0_t)
        deallocate(sfx%Rib_conv_lim)
        deallocate(sfx%Cm)
        deallocate(sfx%Ct)
        deallocate(sfx%Km)
        deallocate(sfx%Pr_t_inv)

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

    ! --------------------------------------------------------------------------------
    !> @brief helper subroutine for copying data in sfxDataVecType
    subroutine push_sfx_data(sfx, sfx_cell, idx)
        ! ----------------------------------------------------------------------------
        type (sfxDataVecType), intent(inout) :: sfx

        type (sfxDataType), intent(in) :: sfx_cell
        integer, intent(in) :: idx
        ! ----------------------------------------------------------------------------

        sfx%zeta(idx) = sfx_cell%zeta
        sfx%Rib(idx) = sfx_cell%Rib
        sfx%Re(idx) = sfx_cell%Re
        sfx%B(idx) = sfx_cell%B
        sfx%z0_m(idx) = sfx_cell%z0_m
        sfx%z0_t(idx) = sfx_cell%z0_t
        sfx%Rib_conv_lim(idx) = sfx_cell%Rib_conv_lim
        sfx%Cm(idx) = sfx_cell%Cm
        sfx%Ct(idx) = sfx_cell%Ct
        sfx%Km(idx) = sfx_cell%Km
        sfx%Pr_t_inv(idx) = sfx_cell%Pr_t_inv

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

end module sfx_data