Skip to content
Snippets Groups Projects
sfx_surface.f90 24.7 KiB
Newer Older
  • Learn to ignore specific revisions
  • 数学の武士's avatar
    .  
    数学の武士 committed
    #include "../includeF/sfx_def.fi"
    
    数学の武士's avatar
    数学の武士 committed
        !< @brief surface roughness parameterizations
    
    
        ! modules used
        ! --------------------------------------------------------------------------------
        use sfx_phys_const
    
        use sfx_z0m_all_surface
    
        ! --------------------------------------------------------------------------------
    
        ! directives list
        ! --------------------------------------------------------------------------------
        implicit none
        ! --------------------------------------------------------------------------------
    
    
        ! --------------------------------------------------------------------------------
    
    数学の武士's avatar
    数学の武士 committed
    #if defined(INCLUDE_CXX)
        public :: set_c_struct_sfx_surface_param_values
    #endif
    
        !public :: get_charnock_roughness
        !public :: get_thermal_roughness
    
        public :: get_dynamic_roughness_all
        public :: get_dynamic_roughness_definition
    
        public :: get_thermal_roughness_all
        public :: get_thermal_roughness_definition
    
        ! --------------------------------------------------------------------------------
    
    
    数学の武士's avatar
    数学の武士 committed
        integer, public, parameter :: surface_ocean = 0     !< ocean surface
        integer, public, parameter :: surface_land = 1      !< land surface
        integer, public, parameter :: surface_lake = 2      !< lake surface
    
        integer, public, parameter :: surface_snow = 3      !< snow covered surface
    
        integer, public, parameter :: surface_forest = 4      !< forest covered surface
        integer, public, parameter :: surface_user = 5      !< coast covered surface
        integer, public, parameter :: surface_ice = 6      !< ice covered surface
    
        integer, public, parameter :: surface_mosaic = 7      !< mosaic covered surface
    
    
        character(len = 16), parameter :: surface_ocean_tag = 'ocean'
        character(len = 16), parameter :: surface_land_tag = 'land'
        character(len = 16), parameter :: surface_lake_tag = 'lake'
        character(len = 16), parameter :: surface_snow_tag = 'snow'
    
        character(len = 16), parameter :: surface_forest_tag = 'forest'
        character(len = 16), parameter :: surface_user_tag = 'user'
    
        character(len = 16), parameter :: surface_ice_tag = 'ice'
    
        character(len = 16), parameter :: surface_mosaic_tag = 'mosaic'
    
    
        integer, public, parameter :: z0m_ch = 0 
        integer, public, parameter :: z0m_fe = 1   
        integer, public, parameter :: z0m_ow = 2  
    
        integer, public, parameter :: z0m_map_id = 3      
    
        integer, public, parameter :: z0m_user = 4 
    
        integer, public, parameter :: z0m_and = 5 
    
        integer, public, parameter :: z0m_cst1 = 6
        integer, public, parameter :: z0m_cst2 = 7
        integer, public, parameter :: z0m_cst3 = 8 
    
        integer, public, parameter :: z0t_kl_water = 0 
        integer, public, parameter :: z0t_kl_land = 1   
        integer, public, parameter :: z0t_re = 2      
        integer, public, parameter :: z0t_zi = 3         
        integer, public, parameter :: z0t_ca = 4    
        integer, public, parameter :: z0t_cz = 5      
        integer, public, parameter :: z0t_br = 6      
        integer, public, parameter :: z0t_ot = 7      
        integer, public, parameter :: z0t_du = 8      
        integer, public, parameter :: z0t_zm = 9      
        integer, public, parameter :: z0t_mix = 10      
        integer, public, parameter :: z0t_user = 11   
    
        integer, public, parameter :: z0t_map_id = 12
    
        character(len = 16), parameter :: z0m_tag_ch = 'charnock'
        character(len = 16), parameter :: z0m_tag_fe = 'fetch'
        character(len = 16), parameter :: z0m_tag_ow = 'owen'
        character(len = 16), parameter :: z0m_tag_map = 'map'
        character(len = 16), parameter :: z0m_tag_user = 'user'
    
        character(len = 16), parameter :: z0m_tag_and = 'andreas'
    
        character(len = 16), parameter :: z0m_tag_cst1 = 'coast1'
        character(len = 16), parameter :: z0m_tag_cst2 = 'coast2'
        character(len = 16), parameter :: z0m_tag_cst3 = 'coast3'
    
    
        character(len = 16), parameter :: z0t_tag_kl_water = 'kl_water'
        character(len = 16), parameter :: z0t_tag_kl_land = 'kl_land'
        character(len = 16), parameter :: z0t_tag_re = 're'
        character(len = 16), parameter :: z0t_tag_zi = 'zi'
        character(len = 16), parameter :: z0t_tag_ca = 'ca'
        character(len = 16), parameter :: z0t_tag_cz = 'cz'
        character(len = 16), parameter :: z0t_tag_br = 'br'
        character(len = 16), parameter :: z0t_tag_ot = 'ot'
        character(len = 16), parameter :: z0t_tag_du = 'du'
        character(len = 16), parameter :: z0t_tag_zm = 'zm'
        character(len = 16), parameter :: z0t_tag_mix = 'mix'
        character(len = 16), parameter :: z0t_tag_user = 'zt_user'
    
        character(len = 16), parameter :: z0t_tag_map = 'zt_map'
    
        integer, public, parameter :: ocean_z0m_id = z0m_ch     !< ocean surface
    
        integer, public, parameter :: land_z0m_id = z0m_map_id        !< land surface
    
        integer, public, parameter :: lake_z0m_id = z0m_ch        !< lake surface
    
        integer, public, parameter :: snow_z0m_id = z0m_map_id       !< snow covered surface    
    
        integer, public, parameter :: forest_z0m_id = z0m_map_id     !< forest csurface  
    
        integer, public, parameter :: usersf_z0m_id = z0m_ch      !< user surface  
    
        integer, public, parameter :: ice_z0m_id = z0m_map_id       !< ice surface  
    
        integer, public, parameter :: mosaic_z0m_id = z0m_map_id
    
        integer, public, parameter :: ocean_z0t_id = z0t_kl_water     !< ocean surface
    
        integer, public, parameter :: land_z0t_id = z0t_kl_land       !< land surface
    
        integer, public, parameter :: lake_z0t_id = z0t_kl_water        !< lake surface
    
        integer, public, parameter :: snow_z0t_id = z0t_kl_land         !< snow covered surface    
        integer, public, parameter :: forest_z0t_id = z0t_kl_land      !< forest csurface  
    
        integer, public, parameter :: usersf_z0t_id = z0t_kl_water      !< user surface  
    
        integer, public, parameter :: ice_z0t_id = z0t_kl_land     !< user surface  
    
        integer, public, parameter :: mosaic_z0t_id = z0t_kl_land
    
        ! --------------------------------------------------------------------------------
    
    数学の武士's avatar
    数学の武士 committed
        real, parameter, private :: kappa = 0.40         !< von Karman constant [n/d]
    
        ! --------------------------------------------------------------------------------
    
    
    数学の武士's avatar
    数学の武士 committed
        !< Charnock parameters
        !<     z0 = Re_visc_min * (nu / u_dyn) + gamma_c * (u_dyn^2 / g)
    
        ! --------------------------------------------------------------------------------
    
    数学の武士's avatar
    数学の武士 committed
        
    
        !real, parameter :: gamma_c = 0.0144
        !real, parameter :: Re_visc_min = 0.111
    
        !real, parameter :: h_charnock = 10.0
        !real, parameter :: c1_charnock = log(h_charnock * (g / gamma_c))
        !real, parameter :: c2_charnock = Re_visc_min * nu_air * (g / gamma_c)
    
        ! --------------------------------------------------------------------------------
    
    
    数学の武士's avatar
    数学の武士 committed
        !< Re fully roughness minimum value [n/d]
    
    数学の武士's avatar
    数学の武士 committed
        !< roughness model coeff. [n/d]
        !< --- transitional mode
        !<     B = log(z0_m / z0_t) = B1 * log(B3 * Re) + B2
    
        !real, parameter :: B1_rough = 5.0 / 6.0
        !real, parameter :: B2_rough = 0.45
        !real, parameter :: B3_rough = kappa * Pr_m
    
    数学の武士's avatar
    数学の武士 committed
        !< --- fully rough mode (Re > Re_rough_min)
        !<     B = B4 * Re^(B2)
    
        !real, parameter :: B4_rough =(0.14 * (30.0**B2_rough)) * (Pr_m**0.8)
    
        !real, parameter :: B_max_land = 2.0
        !real, parameter :: B_max_ocean = 8.0
        !real, parameter :: B_max_lake = 8.0
    
        ! surface type definition
        ! --------------------------------------------------------------------------------
        function get_surface_id(tag) result(id)
            implicit none
            character(len=*), intent(in) :: tag
            integer :: id
    
            id = - 1
            if (trim(tag) == trim(surface_ocean_tag)) then
                id = surface_ocean
            else if (trim(tag) == trim(surface_land_tag)) then
                id = surface_land
            else if (trim(tag) == trim(surface_lake_tag)) then
                id = surface_lake
            else if (trim(tag) == trim(surface_snow_tag)) then
                id = surface_snow
    
            else if (trim(tag) == trim(surface_forest_tag)) then
                id = surface_forest
            else if (trim(tag) == trim(surface_user_tag)) then
                id = surface_user
    
            else if (trim(tag) == trim(surface_ice_tag)) then
                id = surface_ice
    
            else if (trim(tag) == trim(surface_mosaic_tag)) then
                id = surface_mosaic
    
            end if
    
        end function
    
        function get_surface_tag(id) result(tag)
            implicit none
            integer :: id
            character(len=:), allocatable :: tag
    
            tag = 'undefined'
            if (id == surface_ocean) then
                tag = surface_ocean_tag
            else if (id == surface_land) then
                tag = surface_land_tag
            else if (id == surface_lake) then
                tag = surface_lake_tag
            else if (id == surface_snow) then
                tag = surface_snow_tag
    
            else if (id == surface_forest) then
                tag = surface_forest_tag
            else if (id == surface_user) then
                tag = surface_user_tag
    
            else if (id == surface_ice) then
                tag = surface_ice_tag
    
            else if (id == surface_mosaic) then
                tag = surface_mosaic_tag
    
            end if 
    
        end function
    
    
    
        function get_surface_z0m_id(tag_z0m) result(z0m_id)
            implicit none
            character(len=*), intent(in) :: tag_z0m
            integer :: z0m_id
    
            z0m_id = - 1
            if (trim(tag_z0m) == trim(z0m_tag_ch)) then
                z0m_id = z0m_ch
            else if (trim(tag_z0m) == trim(z0m_tag_fe)) then
                z0m_id = z0m_fe
            else if (trim(tag_z0m) == trim(z0m_tag_ow)) then
                z0m_id = z0m_ow
            else if (trim(tag_z0m) == trim(z0m_tag_map)) then
    
                z0m_id = z0m_map_id
    
            else if (trim(tag_z0m) == trim(z0m_tag_user)) then
                z0m_id = z0m_user
    
            else if (trim(tag_z0m) == trim(z0m_tag_cst1)) then
                z0m_id = z0m_cst1
            else if (trim(tag_z0m) == trim(z0m_tag_cst2)) then
                z0m_id = z0m_cst2
            else if (trim(tag_z0m) == trim(z0m_tag_cst3)) then
                z0m_id = z0m_cst3
    
            end if
    
        end function
    
        function get_surface_z0m_tag(z0m_id) result(tag_z0m)
            implicit none
            integer :: z0m_id
            character(len=:), allocatable :: tag_z0m
    
            tag_z0m = 'undefined'
            if (z0m_id == z0m_ch) then
                tag_z0m = z0m_tag_ch
            else if (z0m_id == z0m_fe) then
                tag_z0m =  z0m_tag_fe
            else if (z0m_id == z0m_ow) then
                tag_z0m = z0m_tag_ow
    
            else if (z0m_id == z0m_map_id) then
    
                tag_z0m = z0m_tag_map
            else if (z0m_id == z0m_user) then
                tag_z0m = z0m_tag_user
    
            else if (z0m_id == z0m_cst1) then
                tag_z0m = z0m_tag_cst1
            else if (z0m_id == z0m_cst2) then
                tag_z0m = z0m_tag_cst2
            else if (z0m_id == z0m_cst3) then
                tag_z0m = z0m_tag_cst3
    
         function get_surface_z0t_id(tag_z0t) result(z0t_id)
            implicit none
            character(len=*), intent(in) :: tag_z0t
            integer :: z0t_id
    
            z0t_id = - 1
            if (trim(tag_z0t) == trim(z0t_tag_kl_water)) then
                z0t_id = z0t_kl_water
            else if (trim(tag_z0t) == trim(z0t_tag_kl_land)) then
                z0t_id = z0t_kl_land
            else if (trim(tag_z0t) == trim(z0t_tag_re)) then
                z0t_id = z0t_re
            else if (trim(tag_z0t) == trim(z0t_tag_zi)) then
                z0t_id = z0t_zi
            else if (trim(tag_z0t) == trim(z0t_tag_ca)) then
                z0t_id = z0t_ca
            else if (trim(tag_z0t) == trim(z0t_tag_cz)) then
                z0t_id = z0t_cz
            else if (trim(tag_z0t) == trim(z0t_tag_br)) then
                z0t_id = z0t_br
            else if (trim(tag_z0t) == trim(z0t_tag_ot)) then
                z0t_id = z0t_ot
            else if (trim(tag_z0t) == trim(z0t_tag_du)) then
                z0t_id = z0t_du
            else if (trim(tag_z0t) == trim(z0t_tag_zm)) then
                z0t_id = z0t_zm
            else if (trim(tag_z0t) == trim(z0t_tag_mix)) then
                z0t_id = z0t_mix
            else if (trim(tag_z0t) == trim(z0t_tag_user)) then
                z0t_id = z0t_user
    
            else if (trim(tag_z0t) == trim(z0t_tag_map)) then
                z0t_id = z0t_map_id
    
        end function
    
        function get_surface_z0t_tag(z0t_id) result(tag_z0t)
            implicit none
            integer :: z0t_id
            character(len=:), allocatable :: tag_z0t
    
            tag_z0t = 'undefined'
            if (z0t_id == z0t_kl_water) then
                tag_z0t = z0t_tag_kl_water
            else if (z0t_id == z0t_kl_land) then
                tag_z0t =  z0t_tag_kl_land
            else if (z0t_id == z0t_re) then
                tag_z0t = z0t_tag_re
            else if (z0t_id == z0t_zi) then
                tag_z0t = z0t_tag_zi
            else if (z0t_id == z0t_ca) then
                tag_z0t = z0t_tag_ca
            else if (z0t_id == z0t_cz) then
                tag_z0t = z0t_tag_cz
            else if (z0t_id == z0t_br) then
                tag_z0t = z0t_tag_br
            else if (z0t_id == z0t_ot) then
                tag_z0t = z0t_tag_ot
            else if (z0t_id == z0t_du) then
                tag_z0t = z0t_tag_du
            else if (z0t_id == z0t_zm) then
                tag_z0t = z0t_tag_zm
            else if (z0t_id == z0t_mix) then
                tag_z0t = z0t_tag_mix
            else if (z0t_id == z0t_user) then
                tag_z0t = z0t_tag_user
            end if 
         end function
    
    
    
    数学の武士's avatar
    数学の武士 committed
    #if defined(INCLUDE_CXX)
        subroutine set_c_struct_sfx_surface_param_values(surface_param)
            use sfx_data
    
            use sfx_z0m_all_surface
            use sfx_z0t_all_surface
    
    数学の武士's avatar
    数学の武士 committed
            implicit none
            type (sfx_surface_param), intent(inout) :: surface_param
            surface_param%surface_ocean = surface_ocean
            surface_param%surface_land = surface_land
            surface_param%surface_lake = surface_lake
    
            surface_param%kappa = kappa
            surface_param%gamma_c = gamma_c
            surface_param%Re_visc_min = Re_visc_min
            surface_param%h_charnock = h_charnock
            surface_param%c1_charnock = c1_charnock
            surface_param%c2_charnock = c2_charnock
            surface_param%Re_rough_min = Re_rough_min
            surface_param%B1_rough = B1_rough
            surface_param%B2_rough = B2_rough
            surface_param%B3_rough = B3_rough
            surface_param%B4_rough = B4_rough
            surface_param%B_max_lake = B_max_lake
            surface_param%B_max_ocean = B_max_ocean
            surface_param%B_max_land = B_max_land
        end subroutine set_c_struct_sfx_surface_param_values
    #endif
    
        ! dynamic roughness definition
    
        ! --------------------------------------------------------------------------------
    
          subroutine 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)
    
        ! ----------------------------------------------------------------------------
            integer, intent(out) :: z0m_id              
    
            integer, intent(in) :: surface_type             
    
            integer, intent(in) :: ocean_z0m_id            
            integer, intent(in) :: land_z0m_id             
            integer, intent(in) :: lake_z0m_id
            integer, intent(in) :: snow_z0m_id
            integer, intent(in) :: forest_z0m_id
            integer, intent(in) :: usersf_z0m_id
    
            integer, intent(in) :: ice_z0m_id
    
            integer, intent(in) :: mosaic_z0m_id
    
     ! ---------------------------------------------------------------------------
    
        if (surface_type == surface_ocean) then
         z0m_id = ocean_z0m_id
        else if (surface_type == surface_land) then
         z0m_id = land_z0m_id
        else if (surface_type == surface_snow) then
         z0m_id = snow_z0m_id
        else if (surface_type == surface_lake) then
        z0m_id = lake_z0m_id
        else if (surface_type == surface_forest) then
         z0m_id = forest_z0m_id
        else if (surface_type == surface_user) then
         z0m_id  = usersf_z0m_id
    
        else if (surface_type == surface_ice) then
        z0m_id  = ice_z0m_id
    
        else if (surface_type == surface_mosaic) then
        z0m_id  = mosaic_z0m_id
    
        !write (*,*) z0m_id, surface_type
    
     end subroutine
    
        ! ----------------------------------------------------------------------------
     subroutine get_dynamic_roughness_all(z0_m, u_dyn0, U, depth, h,&
        maxiters, z0m_map, z0m_id)
    ! ----------------------------------------------------------------------------
           real, intent(out) :: z0_m           !< aerodynamic roughness [m]
           real, intent(out) :: u_dyn0         !< dynamic velocity in neutral conditions [m/s]
    
           real, intent(in) :: U               !< abs(wind speed) [m/s]
           real, intent(in) :: depth           !< depth [m]
           real, intent(in) :: h               !< constant flux layer height [m]
           real, intent(in) :: z0m_map          !< aerodynamic roughness from map [m]
           integer, intent(in) :: maxiters     !< maximum number of iterations
           
           integer, intent(in) :: z0m_id
    ! ---------------------------------------------------------------------------
    
    
      if (z0m_id == z0m_ch) then
          call get_dynamic_roughness_ch(z0_m, u_dyn0, U, h, maxiters)
       else if (z0m_id == z0m_fe) then
           call get_dynamic_roughness_fetch(z0_m, u_dyn0, U, depth, h, maxiters)
       else if (z0m_id == z0m_ow) then
           call get_dynamic_roughness_ow(z0_m, u_dyn0, U, h, maxiters)
    
       else if (z0m_id == z0m_map_id) then
    
           call get_dynamic_roughness_map(z0_m, u_dyn0, U, h, z0m_map)
    
        else if (z0m_id == z0m_and) then
            call get_dynamic_roughness_and(z0_m, u_dyn0, U, h, z0m_map)
    
       else if (z0m_id == z0m_user) then
           write(*, *) 'z0m_user'
    
        else if (z0m_id == z0m_cst1) then
            call get_dynamic_roughness_coast1(z0_m, u_dyn0, U, h, maxiters)
        else if (z0m_id == z0m_cst2) then
            call get_dynamic_roughness_coast2(z0_m, u_dyn0, U, h, maxiters)
        else if (z0m_id == z0m_cst3) then
            call get_dynamic_roughness_coast3(z0_m, u_dyn0, U, h, maxiters)
    
    end subroutine
    ! --------------------------------------------------------------------------------  
    
    
    ! thermal roughness definition
        ! --------------------------------------------------------------------------------
    subroutine 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)
    
    ! ----------------------------------------------------------------------------
        integer, intent(out) :: z0t_id              
    
        integer, intent(in) :: surface_type             
    
        integer, intent(in) :: ocean_z0t_id            
        integer, intent(in) :: land_z0t_id             
        integer, intent(in) :: lake_z0t_id
        integer, intent(in) :: snow_z0t_id
        integer, intent(in) :: forest_z0t_id
        integer, intent(in) :: usersf_z0t_id
    
        integer, intent(in) :: ice_z0t_id
    
        integer, intent(in) :: mosaic_z0t_id
    
    ! ---------------------------------------------------------------------------
    
        if (surface_type == surface_ocean) then
        z0t_id = ocean_z0t_id
        else if (surface_type == surface_land) then
        z0t_id = land_z0t_id
        else if (surface_type == surface_snow) then
        z0t_id = snow_z0t_id
        else if (surface_type == surface_lake) then
        z0t_id = lake_z0t_id
        else if (surface_type == surface_forest) then
        z0t_id = forest_z0t_id
        else if (surface_type == surface_user) then
        z0t_id  = usersf_z0t_id
    
        else if (surface_type == surface_ice) then
        z0t_id  = ice_z0t_id
    
        else if (surface_type == surface_mosaic) then
        z0t_id  = mosaic_z0t_id
    
    
        end subroutine
    
    !---------------------------------------------------------------
        subroutine get_thermal_roughness_all(z0_t, B, &
    
        z0_m, Re, u_dyn, LAI, z0t_map, z0t_id)
    
    ! ----------------------------------------------------------------------------
        real, intent(out) :: z0_t               !< thermal roughness [m]
        real, intent(out) :: B                  !< = log(z0_m / z0_t) [n/d]
    
        real, intent(in) :: z0_m                !< aerodynamic roughness [m]
        real, intent(in) :: Re                  !< roughness Reynolds number [n/d]
        real, intent(in) :: LAI                 !< leaf-area index
        real, intent(in) :: u_dyn              !< dynamic velocity [m/s]
    
    
        integer, intent(in) :: z0t_id
    ! ---------------------------------------------------------------------------
        real  Czm                 !< proportionality coefficient z0_t =Czm*z0_m
        Czm=1
    
        if (z0t_id == z0t_kl_water) then
        call get_thermal_roughness_kl_water(z0_t, B, z0_m, Re)
        else if (z0t_id == z0t_kl_land) then
        call get_thermal_roughness_kl_land(z0_t, B, z0_m, Re)
        else if (z0t_id == z0t_re) then
        call get_thermal_roughness_re(z0_t, B, z0_m, Re)
        else if (z0t_id == z0t_zi) then
        call get_thermal_roughness_zi(z0_t, B, z0_m, Re)
        else if (z0t_id == z0t_ca) then
        call get_thermal_roughness_ca(z0_t, B, z0_m, Re)
        else if (z0t_id == z0t_cz) then
        call get_thermal_roughness_cz(z0_t, B, z0_m, Re)
        else if (z0t_id == z0t_br) then
        call get_thermal_roughness_br(z0_t, B, z0_m, Re)
        else if (z0t_id == z0t_ot) then
        call get_thermal_roughness_ot(z0_t, B, z0_m, Re)
        else if (z0t_id == z0t_du) then
        call get_thermal_roughness_du(z0_t, B, z0_m, u_dyn, LAI)
        else if (z0t_id == z0t_zm) then
        call get_thermal_roughness_zm(z0_t, B, z0_m, Czm)
        else if (z0t_id == z0t_mix) then
        call get_thermal_roughness_mix(z0_t, B, z0_m, u_dyn, Re)
        else if (z0t_id == z0t_user) then
        write(*, *) 'z0t_user'
    
        else if (z0t_id == z0t_map_id) then
        call get_thermal_roughness_map(z0_t, B, z0_m, z0t_map)
    
    
    
    
          subroutine get_charnock_roughness(z0_m, u_dyn0, U, h, maxiters)
    
            ! ----------------------------------------------------------------------------
    
    数学の武士's avatar
    数学の武士 committed
            real, intent(out) :: z0_m           !< aerodynamic roughness [m]
            real, intent(out) :: u_dyn0         !< dynamic velocity in neutral conditions [m/s]
    
    数学の武士's avatar
    数学の武士 committed
            real, intent(in) :: h               !< constant flux layer height [m]
            real, intent(in) :: U               !< abs(wind speed) [m/s]
            integer, intent(in) :: maxiters     !< maximum number of iterations
    
            ! ----------------------------------------------------------------------------
    
            ! --- local variables
            real :: Uc                  ! wind speed at h_charnock [m/s]
    
            real :: a, b, c, c_min
            real :: f
    
            integer :: i, j
            ! ----------------------------------------------------------------------------
    
            Uc = U
            a = 0.0
            b = 25.0
            c_min = log(h_charnock) / kappa
    
            do i = 1, maxiters
                f = c1_charnock - 2.0 * log(Uc)
                do j = 1, maxiters
                    c = (f + 2.0 * log(b)) / kappa
                    ! looks like the check should use U10 instead of U
                    !   but note that a1 should be set = 0 in cycle beforehand
                    if (U <= 8.0e0) a = log(1.0 + c2_charnock * ((b / Uc)**3)) / kappa
                    c = max(c - a, c_min)
                    b = c
                end do
                z0_m = h_charnock * exp(-c * kappa)
                z0_m = max(z0_m, 0.000015e0)
                Uc = U * log(h_charnock / z0_m) / log(h / z0_m)
            end do
    
            ! --- define dynamic velocity in neutral conditions
            u_dyn0 = Uc / c
    
        end subroutine
        ! --------------------------------------------------------------------------------
    
    
        ! thermal roughness definition
        ! --------------------------------------------------------------------------------
        subroutine get_thermal_roughness(z0_t, B, &
                z0_m, Re, surface_type)
            ! ----------------------------------------------------------------------------
    
    数学の武士's avatar
    数学の武士 committed
            real, intent(out) :: z0_t               !< thermal roughness [m]
            real, intent(out) :: B                  !< = log(z0_m / z0_t) [n/d]
    
    数学の武士's avatar
    数学の武士 committed
            real, intent(in) :: z0_m                !< aerodynamic roughness [m]
            real, intent(in) :: Re                  !< roughness Reynolds number [n/d]
            integer, intent(in) :: surface_type     !< = [ocean] || [land] || [lake]
    
            ! ----------------------------------------------------------------------------
    
            ! --- local variables
            ! ----------------------------------------------------------------------------
    
            ! --- define B = log(z0_m / z0_t)
            if (Re <= Re_rough_min) then
                B = B1_rough * alog(B3_rough * Re) + B2_rough
            else
                ! *: B4 takes into account Re value at z' ~ O(10) z0
                B = B4_rough * (Re**B2_rough)
            end if
    
            !   --- apply max restriction based on surface type
            if (surface_type == surface_ocean) then
                B = min(B, B_max_ocean)
            else if (surface_type == surface_lake) then
    
            else if (surface_type == surface_land) then
                B = min(B, B_max_land)
            end if
    
            ! --- define roughness [thermal]
            z0_t = z0_m / exp(B)
    
        end subroutine
        ! --------------------------------------------------------------------------------
    
    end module sfx_surface