Skip to content
Snippets Groups Projects
sfx_log.f90 8.7 KiB
Newer Older
  • Learn to ignore specific revisions
  • 数学の武士's avatar
    .  
    数学の武士 committed
    #include "../includeF/sfx_def.fi"
    
    数学の武士's avatar
    数学の武士 committed
        !< @brief simple log-roughness surface flux module
    
    
        ! modules used
        ! --------------------------------------------------------------------------------
    #ifdef SFX_CHECK_NAN
        use sfx_common
    #endif
        use sfx_data
    
        use sfx_log_param
        ! --------------------------------------------------------------------------------
    
        ! directives list
        ! --------------------------------------------------------------------------------
        implicit none
        private
        ! --------------------------------------------------------------------------------
    
        ! public interface
        ! --------------------------------------------------------------------------------
        public :: get_surface_fluxes
        public :: get_surface_fluxes_vec
    
        ! --------------------------------------------------------------------------------
    
        ! --------------------------------------------------------------------------------
        type, public :: numericsType
    
    数学の武士's avatar
    数学の武士 committed
            integer :: maxiters_charnock = 10      !< maximum (actual) number of iterations in charnock roughness
    
        end type
        ! --------------------------------------------------------------------------------
    
    contains
    
        ! --------------------------------------------------------------------------------
        subroutine get_surface_fluxes_vec(sfx, meteo, numerics, n)
    
    数学の武士's avatar
    数学の武士 committed
            !< @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
            ! ----------------------------------------------------------------------------
    
            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), z0_t = meteo%z0_t(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)
            end do
    
        end subroutine get_surface_fluxes_vec
        ! --------------------------------------------------------------------------------
    
        ! --------------------------------------------------------------------------------
        subroutine get_surface_fluxes(sfx, meteo, numerics)
    
    数学の武士's avatar
    数学の武士 committed
            !< @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
            ! ----------------------------------------------------------------------------
    
    数学の武士's avatar
    数学の武士 committed
            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]
    
    数学の武士's avatar
    数学の武士 committed
            real B                  !< = ln(z0_m / z0_t) [n/d]
            real h0_m, h0_t         !< = h / z0_m, h / z0_h [n/d]
    
    数学の武士's avatar
    数学の武士 committed
            real u_dyn0             !< dynamic velocity in neutral conditions [m/s]
            real Re                 !< roughness Reynolds number = u_dyn0 * z0_m / nu [n/d]
    
    数学の武士's avatar
    数学の武士 committed
            real Rib                !< bulk Richardson number
    
    数学の武士's avatar
    数学の武士 committed
            real psi_m, psi_h       !< universal functions (momentum) & (heat) [n/d]
            real phi_m, phi_h       !< stability functions (momentum) & (heat) [n/d]
    
    数学の武士's avatar
    数学の武士 committed
            real Km                 !< eddy viscosity coeff. at h [m^2/s]
            real Pr_t_inv           !< invese Prandt number [n/d]
    
    数学の武士's avatar
    数学の武士 committed
            real Cm, Ct             !< transfer coeff. for (momentum) & (heat) [n/d]
    
    数学の武士's avatar
    数学の武士 committed
            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_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, mosaic_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, mosaic_z0t_id, z0t_id)
    
            call get_thermal_roughness_all(z0_t, B, z0_m, Re, u_dyn0, lai, z0_t1, z0t_id)
    
            ! --- 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_psi_neutral(psi_m, psi_h, h0_m, h0_t, B)
    
    Evgeny Mortikov's avatar
    Evgeny Mortikov committed
            ! ----------------------------------------------------------------------------
    
    
            phi_m = 1.0
            phi_h = 1.0 / Pr_t_0_inv
    
            ! --- 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 = 0.0, 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
        ! --------------------------------------------------------------------------------
    
        ! universal functions
        ! --------------------------------------------------------------------------------
    
        subroutine get_psi_neutral(psi_m, psi_h, h0_m, h0_t, B)
    
    数学の武士's avatar
    数学の武士 committed
            !< @brief universal functions (momentum) & (heat): neutral case
    
            ! ----------------------------------------------------------------------------
    
    数学の武士's avatar
    数学の武士 committed
            real, intent(out) :: psi_m, psi_h   !< universal functions
    
    数学の武士's avatar
    数学の武士 committed
            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
        ! --------------------------------------------------------------------------------