Changeset 4630


Ignore:
Timestamp:
Jul 30, 2020 2:54:34 PM (4 years ago)
Author:
suehring
Message:

Land-surface model: bugfix in level 3 initialization of root-area-density; Avoid double classifiation of vertical walls (at surfaces that are alo covered by buildings); Land/urban surface: bugfix in resistance calculation - avoid potential divisions by zero; init_grid: in case of ASCII topography flag grid points as terrain and building to allow application of land/urban-surface model

Location:
palm/trunk/SOURCE
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/init_grid.f90

    r4601 r4630  
    2525! -----------------
    2626! $Id$
     27! In case of ASCII topography input flag grid points as terrain and building.
     28!
     29! 4601 2020-07-14 12:06:09Z suehring
    2730! Minor formatting adjustments
    2831!
     
    13471350!--             Flag topography for all grid points which are below
    13481351!--             the local topography height.
    1349 !--             Note, each topography is flagged as building.
     1352!--             Note, each topography is flagged as building (bit 2) as well as
     1353!--             terrain (bit 1) in order to employ urban-surface as well as
     1354!--             land-surface model.
    13501355                IF ( zu(k) - ocean_offset <= buildings_f%var_2d(j,i) )  THEN
    13511356                   topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 0 )
    1352                    topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 ) !indicates building
     1357                   topo_3d(k,j,i) = IBCLR( topo_3d(k,j,i), 1 )
     1358                   topo_3d(k,j,i) = IBSET( topo_3d(k,j,i), 2 )
    13531359                ENDIF
    13541360             ENDDO
  • palm/trunk/SOURCE/land_surface_model_mod.f90

    r4602 r4630  
    2525! -----------------
    2626! $Id$
     27! - Bugfix in level 3 initialization of root-area-density
     28! - Bugfix in resistance calculation - avoid potential divisions by zero
     29! - Avoid double classifiation of vertical walls (at surfaces that are also
     30!   covered by buildings)
     31! - Minor formatting adjustment to increase readability
     32!
     33! 4602 2020-07-14 14:49:45Z suehring
    2734! - Bugfix in level 3 initialization of pavements - wrongly assumed existence of
    2835!   pavement_subsurface_pars
     
    15791586    ddz_soil(nzb_soil:nzt_soil) = 1.0_wp / dz_soil(nzb_soil:nzt_soil)
    15801587
    1581 
    1582 
    15831588 END SUBROUTINE lsm_check_parameters
    15841589
     
    16071612    LOGICAL      ::  horizontal !< Flag indicating horizontal or vertical surfaces
    16081613
    1609     REAL(wp) :: c_surface_tmp,& !< temporary variable for storing the volumetric heat capacity of the surface
    1610                 f1,          & !< resistance correction term 1
    1611                 f2,          & !< resistance correction term 2
    1612                 f3,          & !< resistance correction term 3
    1613                 m_min,       & !< minimum soil moisture
    1614                 e,           & !< water vapour pressure
    1615                 e_s,         & !< water vapour saturation pressure
    1616                 e_s_dt,      & !< derivate of e_s with respect to T
    1617                 tend,        & !< tendency
    1618                 dq_s_dt,     & !< derivate of q_s with respect to T
    1619                 coef_1,      & !< coef. for prognostic equation
    1620                 coef_2,      & !< coef. for prognostic equation
    1621                 f_qsws,      & !< factor for qsws
    1622                 f_qsws_veg,  & !< factor for qsws_veg
    1623                 f_qsws_soil, & !< factor for qsws_soil
    1624                 f_qsws_liq,  & !< factor for qsws_liq
    1625                 f_shf,       & !< factor for shf
    1626                 lambda_soil, & !< Thermal conductivity of the uppermost soil layer (W/m2/K)
     1614    REAL(wp) :: c_surface_tmp,  & !< temporary variable for storing the volumetric heat capacity of the surface
     1615                f1,             & !< resistance correction term 1
     1616                f2,             & !< resistance correction term 2
     1617                f3,             & !< resistance correction term 3
     1618                m_min,          & !< minimum soil moisture
     1619                e,              & !< water vapour pressure
     1620                e_s,            & !< water vapour saturation pressure
     1621                e_s_dt,         & !< derivate of e_s with respect to T
     1622                tend,           & !< tendency
     1623                dq_s_dt,        & !< derivate of q_s with respect to T
     1624                coef_1,         & !< coef. for prognostic equation
     1625                coef_2,         & !< coef. for prognostic equation
     1626                f_qsws,         & !< factor for qsws
     1627                f_qsws_veg,     & !< factor for qsws_veg
     1628                f_qsws_soil,    & !< factor for qsws_soil
     1629                f_qsws_liq,     & !< factor for qsws_liq
     1630                f_shf,          & !< factor for shf
     1631                lambda_soil,    & !< Thermal conductivity of the uppermost soil layer (W/m2/K)
    16271632                lambda_surface, & !< Current value of lambda_surface (W/m2/K)
    1628                 m_liq_max      !< maxmimum value of the liq. water reservoir
     1633                m_liq_max,      & !< maxmimum value of the liq. water reservoir
     1634                ueff              !< limited near-surface wind speed - used for calculation of resistance
    16291635
    16301636    TYPE(surf_type_lsm), POINTER ::  surf_t_surface
     
    17951801!--    to positive values.
    17961802       IF ( horizontal  .OR.  .NOT. aero_resist_kray )  THEN
    1797           surf%r_a(m) = ABS( ( surf%pt1(m) - surf%pt_surface(m) ) /            &
     1803          surf%r_a(m) = ABS( ( surf%pt1(m) - surf%pt_surface(m) ) /                               &
    17981804                             ( surf%ts(m) * surf%us(m) + 1.0E-20_wp ) )
    17991805       ELSE
    1800           surf%r_a(m) = rho_cp / ( surf%z0(m) * 1000.0_wp                      &
    1801                         * ( 11.8_wp + 4.2_wp *                                 &
    1802                         SQRT( MAX( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 + &
    1803                                    ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 + &
    1804                                    ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2,  &
    1805                               0.01_wp ) )                                      &
    1806                            )  - 4.0_wp  )
     1806!
     1807!--       Limit wind velocity in order to avoid division by zero.
     1808!--       The nominator can become <= 0.0 for values z0 < 3*10E-4.
     1809          ueff        = MAX( SQRT( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 +                    &
     1810                                    ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 +                   &
     1811                                    ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2 ),                  &
     1812                             1.0_wp / 4.2_wp * ( 4.0_wp / ( surf%z0(m) * 1000.0_wp ) - 11.8_wp ), &
     1813                             0.1_wp                                                               &
     1814                           )
     1815          surf%r_a(m) = rho_cp / ( surf%z0(m) * 1000.0_wp                                         &
     1816                                              * ( 11.8_wp + 4.2_wp * ueff )  - 4.0_wp  )
    18071817       ENDIF
    18081818!
     
    20202030       ENDIF
    20212031
    2022        surf%pt_surface(m)          = surf_t_surface_p%var_1d(m) / exner(nzb)
     2032       surf%pt_surface(m) = surf_t_surface_p%var_1d(m) / exner(nzb)
    20232033
    20242034!
     
    20342044
    20352045       surf%ghf(m) = lambda_surface * ( surf_t_surface_p%var_1d(m)             &
    2036                                              - surf_t_soil%var_2d(nzb_soil,m) )
     2046                                      - surf_t_soil%var_2d(nzb_soil,m) )
    20372047
    20382048       surf%shf(m) = - f_shf * ( surf%pt1(m) - surf%pt_surface(m) ) / c_p
    2039 
    20402049!
    20412050! update the 3d field of rad_lw_out array to have consistent output
     
    21982207       CALL calc_q_surface
    21992208    ENDIF
    2200 
    22012209!
    22022210!-- Calculate new roughness lengths (for water surfaces only)
     
    27212729                           building_type_f%var(j,i) /= building_type_f%fill )  &
    27222730                      THEN
    2723                          vegetation_type_f%var(j,i)                 = 1 ! bare soil
    2724                          soil_type_f%var_2d(j,i)                    = 1
     2731                         vegetation_type_f%var(j,i) = 1 ! bare soil
     2732                         soil_type_f%var_2d(j,i)    = 1
     2733
     2734                         water_type_f%var(j,i)      = water_type_f%fill
     2735                         pavement_type_f%var(j,i)   = pavement_type_f%fill
    27252736!
    27262737!--                      If surface_fraction is provided in static input,
     
    27382749!--                Normally proceed with setting surface types.
    27392750                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
    2740                                             surf_lsm_v(l)%building_covered(m) )
     2751                                                   surf_lsm_v(l)%building_covered(m) )
    27412752                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
    2742                                             surf_lsm_v(l)%building_covered(m) )
     2753                                                   surf_lsm_v(l)%building_covered(m) )
    27432754                   IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill ) &
    27442755                      surf_lsm_v(l)%vegetation_surface(m) = .TRUE.
     
    38213832             nzt_pavement = pavement_depth_level
    38223833          ENDIF
    3823 
    38243834       ENDIF
    38253835!
     
    49044914          DO  m = 1, surf_lsm_h%ns
    49054915             IF ( surf_lsm_h%vegetation_surface(m) )  THEN
    4906                 i = surf_lsm_h%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,            &
    4907                                              surf_lsm_v(l)%building_covered(m) )
    4908                 j = surf_lsm_h%j(m) + MERGE( 0, surf_lsm_v(l)%joff,            &
    4909                                              surf_lsm_v(l)%building_covered(m) )
     4916                i = surf_lsm_h%i(m)
     4917                j = surf_lsm_h%j(m)
    49104918                DO  k = nzb_soil, nzt_soil
    49114919                   surf_lsm_h%root_fr(k,m) = root_area_density_lsm_f%var(k,j,i)
     
    54535461!--          true pavement depth is considered
    54545462             DO  k = nzb_soil, nzt_soil-1
    5455                    surf%lambda_h(k,m) = ( lambda_temp(k+1) + lambda_temp(k) )  &
    5456                                         * 0.5_wp
     5463                surf%lambda_h(k,m) = ( lambda_temp(k+1) + lambda_temp(k) ) * 0.5_wp
    54575464             ENDDO
    54585465             surf%lambda_h(nzt_soil,m) = lambda_temp(nzt_soil)
     
    54625469             tend(:) = 0.0_wp
    54635470
    5464              tend(nzb_soil) = ( 1.0_wp / surf%rho_c_total(nzb_soil,m) ) *            &
    5465                     ( surf%lambda_h(nzb_soil,m) * ( surf_t_soil%var_2d(nzb_soil+1,m) &
    5466                       - surf_t_soil%var_2d(nzb_soil,m) ) * ddz_soil_center(nzb_soil) &
    5467                       + surf%ghf(m) ) * ddz_soil(nzb_soil)
     5471             tend(nzb_soil) = ( 1.0_wp / surf%rho_c_total(nzb_soil,m) ) *                          &
     5472                    ( surf%lambda_h(nzb_soil,m) *                                                  &
     5473                      ( surf_t_soil%var_2d(nzb_soil+1,m) - surf_t_soil%var_2d(nzb_soil,m) )        &
     5474                      * ddz_soil_center(nzb_soil)                                                  &
     5475                    + surf%ghf(m) ) * ddz_soil(nzb_soil)
    54685476
    54695477             DO  k = nzb_soil+1, nzt_soil
     
    54725480                     * ( surf_t_soil%var_2d(k+1,m) - surf_t_soil%var_2d(k,m) ) &
    54735481                     * ddz_soil_center(k)                                      &
    5474                      - surf%lambda_h(k-1,m)                                    &
     5482                              - surf%lambda_h(k-1,m)                           &
    54755483                     * ( surf_t_soil%var_2d(k,m) - surf_t_soil%var_2d(k-1,m) ) &
    54765484                     * ddz_soil_center(k-1)                                    &
  • palm/trunk/SOURCE/urban_surface_mod.f90

    r4602 r4630  
    2727! -----------------
    2828! $Id$
     29! - Bugfix in resistance calculation - avoid potential divisions by zero
     30! - Minor formatting adjustment
     31!
     32! 4602 2020-07-14 14:49:45Z suehring
    2933! Add missing initialization of albedo type with values given from static input file
    3034!
     
    463467    INTEGER(iwp) ::  wall_category = 2               !< default category for wall surface over pedestrian zone
    464468
     469    REAL(wp)     ::  d_roughness_concrete            !< inverse roughness length of average concrete surface
    465470    REAL(wp)     ::  roughness_concrete = 0.001_wp   !< roughness length of average concrete surface
    466471
     
    35833588    ENDIF
    35843589!
     3590!-- Calculate constant values
     3591    d_roughness_concrete = 1.0_wp / roughness_concrete
     3592!
    35853593!-- Flag surface elements belonging to the ground floor level. Therefore, use terrain height array
    35863594!-- from file, if available. This flag is later used to control initialization of surface attributes.
     
    75397547                  q_s,                     &  !< saturation specific humidity
    75407548                  rho_lv,                  &  !< frequently used parameter for green layers
    7541                   tend                      !< tendency
     7549                  tend,                    &  !< tendency
     7550                  ueff                        !< limited near-surface wind speed - used for calculation of resistance
    75427551
    75437552
     
    80468055!--        normal component is obtained by simple linear interpolation. (An alternative would be an
    80478056!--        logarithmic interpolation.) Parameter roughness_concrete (default value = 0.001) is used
    8048 !--        to calculation of roughness relative to concrete
    8049            surf_usm_v(l)%r_a(m) = rho_cp / ( surf_usm_v(l)%z0(m) / roughness_concrete              &
    8050                                   * ( 11.8_wp + 4.2_wp                                             &
    8051                                       * SQRT( MAX( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 +     &
    8052                                           ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 +              &
    8053                                           ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2, 0.01_wp ) )   &
    8054                                     )  - 4.0_wp  )
     8057!--        to calculation of roughness relative to concrete. Note, wind velocity is limited
     8058!--        to avoid division by zero. The nominator can become <= 0.0 for values z0 < 3*10E-4.
     8059           ueff        = MAX ( SQRT( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 +                   &
     8060                                     ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 +                   &
     8061                                     ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2 ),                  &
     8062                               1.0_wp / 4.2_wp                                                     &
     8063                               * ( 4.0_wp / ( surf_usm_v(l)%z0(m) * d_roughness_concrete )         &
     8064                               - 11.8_wp ),                                                        &
     8065                               0.1_wp                                                              &
     8066                              )
     8067
     8068           surf_usm_v(l)%r_a(m) = rho_cp / ( surf_usm_v(l)%z0(m) * d_roughness_concrete            &
     8069                                  * ( 11.8_wp + 4.2_wp * ueff )  - 4.0_wp  )
    80558070!
    80568071!--        Limit aerodynamic resistance
Note: See TracChangeset for help on using the changeset viewer.