Skip to content
Snippets Groups Projects
sfx_sheba.f90 18 KiB
Newer Older
  • Learn to ignore specific revisions
  • #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
    
    数学の武士's avatar
    .  
    数学の武士 committed
    
    
    数学の武士's avatar
    数学の武士 committed
    #if defined(INCLUDE_CXX)
    
    数学の武士's avatar
    数学の武士 committed
        use iso_c_binding, only: C_LOC, C_PTR, C_INT, C_FLOAT
    
    数学の武士's avatar
    .  
    数学の武士 committed
        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
        ! --------------------------------------------------------------------------------
    
    
    数学の武士's avatar
    数学の武士 committed
    #if defined(INCLUDE_CXX)
        type, BIND(C), public :: sfx_sheba_param_C
            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_numericsType_C 
            integer(C_INT) :: maxiters_charnock 
        end type
    
        INTERFACE
    
            SUBROUTINE c_sheba_compute_flux(sfx, meteo, model_param, surface_param, numerics, constants, grid_size) BIND(C, & 
                name="c_sheba_compute_flux")
    
    数学の武士's avatar
    数学の武士 committed
                use sfx_data
                use, intrinsic :: ISO_C_BINDING, ONLY: C_INT, C_PTR
                Import :: sfx_sheba_param_C, sfx_sheba_numericsType_C
                implicit none
                INTEGER(C_INT) :: grid_size
                type(C_PTR), value :: sfx
                type(C_PTR), value :: meteo
                type(sfx_sheba_param_C) :: model_param
                type(sfx_surface_param) :: surface_param
                type(sfx_sheba_numericsType_C) :: numerics
                type(sfx_phys_constants) :: constants
    
            END SUBROUTINE c_sheba_compute_flux
    
    数学の武士's avatar
    数学の武士 committed
        END INTERFACE
    #endif 
    
    
    数学の武士's avatar
    数学の武士 committed
    #if defined(INCLUDE_CXX)
        subroutine set_c_struct_sfx_sheba_param_values(sfx_model_param)
    
    数学の武士's avatar
    数学の武士 committed
            type (sfx_sheba_param_C), intent(inout) :: sfx_model_param
    
    数学の武士's avatar
    数学の武士 committed
            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
            ! ----------------------------------------------------------------------------
    
    数学の武士's avatar
    数学の武士 committed
    #if defined(INCLUDE_CXX)
            type (meteoDataVecTypeC), pointer :: meteo_c         !< meteorological data (input)
            type (sfxDataVecTypeC), pointer :: sfx_c             !< surface flux data (output)
    
    数学の武士's avatar
    数学の武士 committed
            type(C_PTR) :: meteo_c_ptr, sfx_c_ptr
    
    数学の武士's avatar
    数学の武士 committed
            type (sfx_sheba_param_C) :: model_param
    
    数学の武士's avatar
    数学の武士 committed
            type (sfx_surface_param) :: surface_param
    
    数学の武士's avatar
    数学の武士 committed
            type (sfx_sheba_numericsType_C) :: numerics_c
    
    数学の武士's avatar
    数学の武士 committed
            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)
    
    
    数学の武士's avatar
    数学の武士 committed
            meteo_c_ptr = C_LOC(meteo_c)
            sfx_c_ptr   = C_LOC(sfx_c)
    
            call c_sheba_compute_flux(sfx_c_ptr, meteo_c_ptr, model_param, surface_param, numerics_c, phys_constants, n)
    
    数学の武士's avatar
    .  
    数学の武士 committed
    #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
    
    数学の武士's avatar
    .  
    数学の武士 committed
    #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