Ignore:
Timestamp:
Jun 9, 2017 12:18:47 PM (4 years ago)
Author:
maronga
Message:

major revisions in land surface model

File:
1 edited

Legend:

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

    r2249 r2270  
    2020! Current revisions:
    2121! -----------------
    22 !
    2322!
    2423!
     
    2625! -----------------
    2726! $Id$
     27! Revised parameterization of heat conductivity between skin layer and soil.
     28! Temperature and moisture are now defined at the center of the layers.
     29! Renamed veg_type to vegetation_type and pave_type to pavement_type_name
     30! Renamed and reduced the number of look-up tables (vegetation_pars, soil_pars)
     31! Revised land surface model initialization
     32! Removed output of shf_eb and qsws_eb and removed _eb throughout code
     33! Removed Clapp & Hornberger parameterization
     34!
     35! 2249 2017-06-06 13:58:01Z sward
    2836!
    2937! 2248 2017-06-06 13:52:54Z sward $
     
    5058!
    5159! 2149 2017-02-09 16:57:03Z scharf
    52 ! Land surface parameters II corrected for veg_type 18 and 19
     60! Land surface parameters II corrected for vegetation_type 18 and 19
    5361!
    5462! 2031 2016-10-21 15:11:58Z knoop
     
    177185! ------------
    178186!> Land surface model, consisting of a solver for the energy balance at the
    179 !> surface and a four layer soil scheme. The scheme is similar to the TESSEL
     187!> surface and a multi layer soil scheme. The scheme is similar to the TESSEL
    180188!> scheme implemented in the ECMWF IFS model, with modifications according to
    181189!> H-TESSEL. The implementation is based on the formulation implemented in the
     
    195203!>       with considerable precipitation.
    196204!> @todo Fix crashes with radiation_scheme == 'constant'
    197 !>
     205!> @todo Revise calculation of f2 when wilting point is non-constant in the
     206!>       soil
     207!> @todo Solve problems with soil heat flux paramterization
     208!> @todo Make pavement_depth variable for each surface element
     209!> @todo Allow for zero soil moisture (currently, it is set to wilting point)
    198210!> @note No time step criterion is required as long as the soil layers do not
    199211!>       become too thin.
     
    211223               initializing_actions, intermediate_timestep_count_max,          &
    212224               land_surface, max_masks, precipitation, pt_surface,             &
    213                rho_surface, roughness_length, surface_pressure,                &
    214                timestep_scheme, tsc, z0h_factor, time_since_reference_point
     225               rho_surface, surface_pressure, timestep_scheme, tsc,            &
     226               time_since_reference_point
    215227
    216228    USE indices,                                                               &
     
    223235    USE radiation_model_mod,                                                   &
    224236        ONLY:  force_radiation_call, rad_net, rad_sw_in, rad_lw_out,           &
    225                rad_lw_out_change_0, radiation_scheme, unscheduled_radiation_calls
     237               rad_lw_out_change_0, radiation_scheme,                          &
     238               unscheduled_radiation_calls
    226239       
    227240    USE statistics,                                                            &
    228241        ONLY:  hom, statistic_regions
    229242
    230     USE surface_mod,                                                        &
     243    USE surface_mod,                                                           &
    231244        ONLY :  surf_lsm_h, surf_lsm_v, surf_type
    232245
     
    243256    REAL(wp), PARAMETER  ::                    &
    244257              b_ch               = 6.04_wp,    & ! Clapp & Hornberger exponent
    245               lambda_h_dry       = 0.19_wp,    & ! heat conductivity for dry soil  
    246               lambda_h_sm        = 3.44_wp,    & ! heat conductivity of the soil matrix
    247               lambda_h_water     = 0.57_wp,    & ! heat conductivity of water
     258              lambda_h_dry       = 0.19_wp,    & ! heat conductivity for dry soil (W/m/K) 
     259              lambda_h_sm        = 3.44_wp,    & ! heat conductivity of the soil matrix (W/m/K)
     260              lambda_h_water     = 0.57_wp,    & ! heat conductivity of water (W/m/K)
    248261              psi_sat            = -0.388_wp,  & ! soil matrix potential at saturation
    249               rho_c_soil         = 2.19E6_wp,  & ! volumetric heat capacity of soil
    250               rho_c_water        = 4.20E6_wp,  & ! volumetric heat capacity of water
     262              rho_c_soil         = 2.19E6_wp,  & ! volumetric heat capacity of soil (J/m3/K)
     263              rho_c_water        = 4.20E6_wp,  & ! volumetric heat capacity of water (J/m3/K)
    251264              m_max_depth        = 0.0002_wp     ! Maximum capacity of the water reservoir (m)
    252265
     
    260273!
    261274!-- LSM variables
    262     INTEGER(iwp) :: nzb_soil = 0,  & !< bottom of the soil model (Earth's surface)
    263                     nzt_soil = 0,  & !< top of the soil model
    264                     nzs = 0,       & !< number of soil layers
    265                     veg_type  = 2, & !< default NAMELIST veg_type_2d
    266                     soil_type = 3    !< default NAMELIST soil_type_2d
     275    CHARACTER(10) :: surface_type = 'netcdf'  !< general classification. Allowed are:
     276                                             !< 'vegetation', 'pavement', ('building'),
     277                                             !< 'water', and 'netcdf'
     278
     279
     280
     281    INTEGER(iwp) :: nzb_soil = 0,        & !< bottom of the soil model (Earth's surface)
     282                    nzt_soil = 7,        & !< top of the soil model
     283                    nzs = 8,             & !< number of soil layers
     284                    pavement_type = 1,   & !< default NAMELIST pavement_type                 
     285                    soil_type = 3,       & !< default NAMELIST soil_type
     286                    vegetation_type = 2, & !< default NAMELIST vegetation_type
     287                    water_type = 1         !< default NAMELISt water_type
     288                   
    267289   
    268290       
    269291    LOGICAL :: conserve_water_content = .TRUE.,  & !< open or closed bottom surface for the soil model
     292               constant_roughness = .FALSE.,     & !< use fixed/dynamic roughness lengths for water surfaces
    270293               force_radiation_call_l = .FALSE., & !< flag to force calling of radiation routine
    271294               aero_resist_kray = .TRUE.           !< flag to control parametrization of aerodynamic resistance at vertical surface elements
     
    275298    REAL(wp) :: alpha_vangenuchten = 9999999.9_wp,      & !< NAMELIST alpha_vg
    276299                canopy_resistance_coefficient = 9999999.9_wp, & !< NAMELIST g_d
    277                 c_surface   = 20000.0_wp,               & !< Surface (skin) heat capacity
     300                c_surface = 9999999.9_wp,               & !< Surface (skin) heat capacity (J/m2/K)
    278301                drho_l_lv,                              & !< (rho_l * l_v)**-1
    279302                exn,                                    & !< value of the Exner function
     
    283306                hydraulic_conductivity = 9999999.9_wp,  & !< NAMELIST gamma_w_sat
    284307                ke = 0.0_wp,                            & !< Kersten number
    285                 lambda_h_sat = 0.0_wp,                  & !< heat conductivity for saturated soil
    286                 lambda_surface_stable = 9999999.9_wp,   & !< NAMELIST lambda_surface_s
    287                 lambda_surface_unstable = 9999999.9_wp, & !< NAMELIST lambda_surface_u
     308                lambda_h_sat = 0.0_wp,                  & !< heat conductivity for saturated soil (W/m/K)
     309                lambda_surface_stable = 9999999.9_wp,   & !< NAMELIST lambda_surface_s (W/m2/K)
     310                lambda_surface_unstable = 9999999.9_wp, & !< NAMELIST lambda_surface_u (W/m2/K)
    288311                leaf_area_index = 9999999.9_wp,         & !< NAMELIST lai
    289312                l_vangenuchten = 9999999.9_wp,          & !< NAMELIST l_vg
     
    292315                m_total = 0.0_wp,                       & !< weighted total water content of the soil (m3/m3)
    293316                n_vangenuchten = 9999999.9_wp,          & !< NAMELIST n_vg
    294                 pave_depth = 9999999.9_wp,              & !< depth of the pavement
    295                 pave_heat_capacity = 1.94E6_wp,         & !< volumetric heat capacity of pavement (e.g. roads)
    296                 pave_heat_conductivity = 1.00_wp,       & !< heat conductivity for pavements (e.g. roads)
     317                pavement_depth = 9999999.9_wp,          & !< depth of the pavement
     318                pavement_heat_capacity = 9999999.9_wp,  & !< volumetric heat capacity of pavement (e.g. roads) (J/m3/K)
     319                pavement_heat_conduct  = 9999999.9_wp,  & !< heat conductivity for pavements (e.g. roads) (W/m/K)
    297320                q_s = 0.0_wp,                           & !< saturation specific humidity
    298321                residual_moisture = 9999999.9_wp,       & !< NAMELIST m_res
     
    303326                skip_time_do_lsm = 0.0_wp,              & !< LSM is not called before this time
    304327                vegetation_coverage = 9999999.9_wp,     & !< NAMELIST c_veg
     328                water_temperature = 9999999.9_wp,       & !< water temperature
    305329                wilting_point = 9999999.9_wp,           & !< NAMELIST m_wilt
    306                 z0_eb  = 9999999.9_wp,                  & !< NAMELIST z0 (lsm_par)
    307                 z0h_eb = 9999999.9_wp,                  & !< NAMELIST z0h (lsm_par)
    308                 z0q_eb = 9999999.9_wp                     !< NAMELIST z0q (lsm_par)
    309 
     330                z0_vegetation  = 9999999.9_wp,          & !< NAMELIST z0 (lsm_par)
     331                z0h_vegetation = 9999999.9_wp,          & !< NAMELIST z0h (lsm_par)
     332                z0q_vegetation = 9999999.9_wp,          & !< NAMELIST z0q (lsm_par)
     333                z0_pavement    = 9999999.9_wp,          & !< NAMELIST z0 (lsm_par)
     334                z0h_pavement   = 9999999.9_wp,          & !< NAMELIST z0h (lsm_par)
     335                z0q_pavement   = 9999999.9_wp,          & !< NAMELIST z0q (lsm_par)
     336                z0_water       = 9999999.9_wp,          & !< NAMELIST z0 (lsm_par)
     337                z0h_water      = 9999999.9_wp,          & !< NAMELIST z0h (lsm_par)
     338                z0q_water      = 9999999.9_wp             !< NAMELIST z0q (lsm_par) 
     339               
     340               
    310341    REAL(wp), DIMENSION(:), ALLOCATABLE  :: ddz_soil,         & !< 1/dz_soil
    311                                             ddz_soil_stag,    & !< 1/dz_soil_stag
    312                                             dz_soil,          & !< soil grid spacing (edge-edge)
    313                                             dz_soil_stag,     & !< soil grid spacing (center-center)
     342                                            ddz_soil_layer,   & !< 1/dz_soil_layer
     343                                            dz_soil,          & !< soil grid spacing (center-center)
     344                                            dz_soil_layer,    & !< soil grid spacing (edge-edge)
    314345                                            root_extr           !< root extraction
    315346
    316347
    317348                                           
    318     REAL(wp), DIMENSION(0:20)  ::  root_fraction = 9999999.9_wp,               & !< distribution of root surface area to the individual soil layers
     349    REAL(wp), DIMENSION(0:20)  ::  root_fraction = 9999999.9_wp,               & !< (NAMELIST) distribution of root surface area to the individual soil layers
    319350                                   soil_moisture = 0.0_wp,                     & !< NAMELIST soil moisture content (m3/m3)
    320351                                   soil_temperature = 300.0_wp,                & !< NAMELIST soil temperature (K) +1
    321                                    zs = 9999999.9_wp                             !< soil layer depths
    322                                    
     352                                   zs = 9999999.9_wp,                          & !< (NAMELIST) soil layer depths (center)
     353                                   zs_layer = 9999999.9_wp                       !< soil layer depths (edge)
     354                                 
    323355#if defined( __nopointer )
    324356    TYPE(surf_type_lsm), TARGET  ::  t_soil_h,    & !< Soil temperature (K), horizontal surface elements
     
    359391    TYPE(surf_type_lsm), TARGET   ::  t_surface_h,    & !< surface temperature (K), horizontal surface elements
    360392                                      t_surface_h_p,  & !< progn. surface temperature (K), horizontal surface elements
    361                                       m_liq_eb_h,     & !< liquid water reservoir (m), horizontal surface elements
    362                                       m_liq_eb_h_p      !< progn. liquid water reservoir (m), horizontal surface elements
     393                                      m_liq_h,     & !< liquid water reservoir (m), horizontal surface elements
     394                                      m_liq_h_p      !< progn. liquid water reservoir (m), horizontal surface elements
    363395
    364396    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET   ::  &
    365397                                      t_surface_v,    & !< surface temperature (K), vertical surface elements
    366398                                      t_surface_v_p,  & !< progn. surface temperature (K), vertical surface elements
    367                                       m_liq_eb_v,     & !< liquid water reservoir (m), vertical surface elements
    368                                       m_liq_eb_v_p      !< progn. liquid water reservoir (m), vertical surface elements
     399                                      m_liq_v,     & !< liquid water reservoir (m), vertical surface elements
     400                                      m_liq_v_p      !< progn. liquid water reservoir (m), vertical surface elements
    369401#else
    370402    TYPE(surf_type_lsm), POINTER  ::  t_surface_h,    & !< surface temperature (K), horizontal surface elements
    371403                                      t_surface_h_p,  & !< progn. surface temperature (K), horizontal surface elements
    372                                       m_liq_eb_h,     & !< liquid water reservoir (m), horizontal surface elements
    373                                       m_liq_eb_h_p      !< progn. liquid water reservoir (m), horizontal surface elements
     404                                      m_liq_h,     & !< liquid water reservoir (m), horizontal surface elements
     405                                      m_liq_h_p      !< progn. liquid water reservoir (m), horizontal surface elements
    374406
    375407    TYPE(surf_type_lsm), TARGET   ::  t_surface_h_1,  & !<
    376408                                      t_surface_h_2,  & !<
    377                                       m_liq_eb_h_1,   & !<
    378                                       m_liq_eb_h_2      !<
     409                                      m_liq_h_1,   & !<
     410                                      m_liq_h_2      !<
    379411
    380412    TYPE(surf_type_lsm), DIMENSION(:), POINTER  ::    &
    381413                                      t_surface_v,    & !< surface temperature (K), vertical surface elements
    382414                                      t_surface_v_p,  & !< progn. surface temperature (K), vertical surface elements
    383                                       m_liq_eb_v,     & !< liquid water reservoir (m), vertical surface elements
    384                                       m_liq_eb_v_p      !< progn. liquid water reservoir (m), vertical surface elements
     415                                      m_liq_v,     & !< liquid water reservoir (m), vertical surface elements
     416                                      m_liq_v_p      !< progn. liquid water reservoir (m), vertical surface elements
    385417
    386418    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET   ::  &
    387419                                      t_surface_v_1,  & !<
    388420                                      t_surface_v_2,  & !<
    389                                       m_liq_eb_v_1,   & !<
    390                                       m_liq_eb_v_2      !<
     421                                      m_liq_v_1,   & !<
     422                                      m_liq_v_2      !<
    391423#endif
    392424
    393425#if defined( __nopointer )
    394     REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: m_liq_eb_av
     426    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: m_liq_av
    395427#else
    396     REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: m_liq_eb_av
     428    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: m_liq_av
    397429#endif
    398430
     
    405437#endif
    406438
    407     TYPE(surf_type_lsm), TARGET ::  tm_liq_eb_h_m   !< liquid water reservoir tendency (m), horizontal surface elements
     439    TYPE(surf_type_lsm), TARGET ::  tm_liq_h_m   !< liquid water reservoir tendency (m), horizontal surface elements
    408440    TYPE(surf_type_lsm), TARGET ::  tt_surface_h_m  !< surface temperature tendency (K), horizontal surface elements
    409441    TYPE(surf_type_lsm), TARGET ::  tt_soil_h_m     !< t_soil storage array, horizontal surface elements
    410442    TYPE(surf_type_lsm), TARGET ::  tm_soil_h_m     !< m_soil storage array, horizontal surface elements
    411443
    412     TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tm_liq_eb_v_m   !< liquid water reservoir tendency (m), vertical surface elements
     444    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tm_liq_v_m   !< liquid water reservoir tendency (m), vertical surface elements
    413445    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tt_surface_v_m  !< surface temperature tendency (K), vertical surface elements
    414446    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tt_soil_v_m     !< t_soil storage array, vertical surface elements
     
    421453              c_soil_av,        & !< average of c_soil
    422454              c_veg_av,         & !< average of c_veg
    423               ghf_eb_av,        & !< average of ghf_eb
     455              ghf_av,        & !< average of ghf
    424456              lai_av,           & !< average of lai
    425               qsws_eb_av,       & !< average of qsws_eb
    426               qsws_liq_eb_av,   & !< average of qsws_liq_eb
    427               qsws_soil_eb_av,  & !< average of qsws_soil_eb
    428               qsws_veg_eb_av,   & !< average of qsws_veg_eb
     457              qsws_liq_av,   & !< average of qsws_liq
     458              qsws_soil_av,  & !< average of qsws_soil
     459              qsws_veg_av,   & !< average of qsws_veg
    429460              r_a_av,           & !< average of r_a
    430               r_s_av,           & !< average of r_s
    431               shf_eb_av           !< average of shf_eb
    432 
    433 
    434 !
    435 !-- Predefined Land surface classes (veg_type)
    436     CHARACTER(26), DIMENSION(0:20), PARAMETER :: veg_type_name = (/ &
    437                                    'user defined              ',    & ! 0
    438                                    'crops, mixed farming      ',    & !  1
    439                                    'short grass               ',    & !  2
    440                                    'evergreen needleleaf trees',    & !  3
    441                                    'deciduous needleleaf trees',    & !  4
    442                                    'evergreen broadleaf trees ',    & !  5
    443                                    'deciduous broadleaf trees ',    & !  6
    444                                    'tall grass                ',    & !  7
    445                                    'desert                    ',    & !  8
    446                                    'tundra                    ',    & !  9
    447                                    'irrigated crops           ',    & ! 10
    448                                    'semidesert                ',    & ! 11
    449                                    'ice caps and glaciers     ',    & ! 12
    450                                    'bogs and marshes          ',    & ! 13
    451                                    'inland water              ',    & ! 14
    452                                    'ocean                     ',    & ! 15
    453                                    'evergreen shrubs          ',    & ! 16
    454                                    'deciduous shrubs          ',    & ! 17
    455                                    'mixed forest/woodland     ',    & ! 18
    456                                    'interrupted forest        ',    & ! 19
    457                                    'pavements/roads           '     & ! 20
     461              r_s_av              !< average of r_s
     462                   
     463
     464!
     465!-- Predefined Land surface classes (vegetation_type)
     466    CHARACTER(26), DIMENSION(0:18), PARAMETER :: vegetation_type_name = (/ &
     467                                   'user defined              ',           & !  0
     468                                   'bare soil                 ',           & !  1                           
     469                                   'crops, mixed farming      ',           & !  2
     470                                   'short grass               ',           & !  3
     471                                   'evergreen needleleaf trees',           & !  4
     472                                   'deciduous needleleaf trees',           & !  5
     473                                   'evergreen broadleaf trees ',           & !  6
     474                                   'deciduous broadleaf trees ',           & !  7
     475                                   'tall grass                ',           & !  8
     476                                   'desert                    ',           & !  9
     477                                   'tundra                    ',           & ! 10
     478                                   'irrigated crops           ',           & ! 11
     479                                   'semidesert                ',           & ! 12
     480                                   'ice caps and glaciers     ',           & ! 13
     481                                   'bogs and marshes          ',           & ! 14
     482                                   'evergreen shrubs          ',           & ! 15
     483                                   'deciduous shrubs          ',           & ! 16
     484                                   'mixed forest/woodland     ',           & ! 17
     485                                   'interrupted forest        '            & ! 18
    458486                                                                 /)
    459487
    460488!
    461489!-- Soil model classes (soil_type)
    462     CHARACTER(12), DIMENSION(0:7), PARAMETER :: soil_type_name = (/ &
     490    CHARACTER(12), DIMENSION(0:6), PARAMETER :: soil_type_name = (/ &
    463491                                   'user defined',                  & ! 0
    464492                                   'coarse      ',                  & ! 1
     
    467495                                   'fine        ',                  & ! 4
    468496                                   'very fine   ',                  & ! 5
    469                                    'organic     ',                  & ! 6
    470                                    'loamy (CH)  '                   & ! 7
     497                                   'organic     '                   & ! 6
    471498                                                                 /)
    472 !
    473 !-- Land surface parameters according to the respective classes (veg_type)
    474 
    475 !
    476 !-- Land surface parameters I
    477 !--                          r_canopy_min,     lai,   c_veg,     g_d
    478     REAL(wp), DIMENSION(0:3,1:20), PARAMETER :: veg_pars = RESHAPE( (/ &
    479                                  180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,  & !  1
    480                                  110.0_wp, 2.00_wp, 0.85_wp, 0.00_wp,  & !  2
    481                                  500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  3
    482                                  500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  4
    483                                  175.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & !  5
    484                                  240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp,  & !  6
    485                                  100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp,  & !  7
    486                                  250.0_wp, 0.05_wp, 0.00_wp, 0.00_wp,  & !  8
    487                                   80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp,  & !  9
    488                                  180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,  & ! 10
    489                                  150.0_wp, 0.50_wp, 0.10_wp, 0.00_wp,  & ! 11
    490                                    0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 12
    491                                  240.0_wp, 4.00_wp, 0.60_wp, 0.00_wp,  & ! 13
    492                                    0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 14
    493                                    0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp,  & ! 15
    494                                  225.0_wp, 3.00_wp, 0.50_wp, 0.00_wp,  & ! 16
    495                                  225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp,  & ! 17
    496                                  250.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,  & ! 18
    497                                  175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp,  & ! 19
    498                                    0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp   & ! 20
    499                                  /), (/ 4, 20 /) )
    500 
    501 !
    502 !-- Land surface parameters II          z0,         z0h,         z0q
    503     REAL(wp), DIMENSION(0:2,1:20), PARAMETER :: roughness_par = RESHAPE( (/ &
    504                                    0.25_wp,  0.25E-2_wp,  0.25E-2_wp,       & !  1
    505                                    0.20_wp,  0.20E-2_wp,  0.20E-2_wp,       & !  2
    506                                    2.00_wp,     2.00_wp,     2.00_wp,       & !  3
    507                                    2.00_wp,     2.00_wp,     2.00_wp,       & !  4
    508                                    2.00_wp,     2.00_wp,     2.00_wp,       & !  5
    509                                    2.00_wp,     2.00_wp,     2.00_wp,       & !  6
    510                                    0.47_wp,  0.47E-2_wp,  0.47E-2_wp,       & !  7
    511                                   0.013_wp, 0.013E-2_wp, 0.013E-2_wp,       & !  8
    512                                   0.034_wp, 0.034E-2_wp, 0.034E-2_wp,       & !  9
    513                                     0.5_wp,  0.50E-2_wp,  0.50E-2_wp,       & ! 10
    514                                    0.17_wp,  0.17E-2_wp,  0.17E-2_wp,       & ! 11
    515                                  1.3E-3_wp,   1.3E-4_wp,   1.3E-4_wp,       & ! 12
    516                                    0.83_wp,  0.83E-2_wp,  0.83E-2_wp,       & ! 13
    517                                    0.00_wp,     0.00_wp,     0.00_wp,       & ! 14
    518                                    0.00_wp,     0.00_wp,     0.00_wp,       & ! 15
    519                                    0.10_wp,  0.10E-2_wp,  0.10E-2_wp,       & ! 16
    520                                    0.25_wp,  0.25E-2_wp,  0.25E-2_wp,       & ! 17
    521                                    2.00_wp,     2.00_wp,     2.00_wp,       & ! 18
    522                                    1.10_wp,     1.10_wp,     1.10_wp,       & ! 19
    523                                  1.0E-4_wp,   1.0E-5_wp,   1.0E-5_wp        & ! 20
    524                                  /), (/ 3, 20 /) )
    525 
    526 !
    527 !-- Land surface parameters III lambda_surface_s, lambda_surface_u, f_sw_in
    528     REAL(wp), DIMENSION(0:2,1:20), PARAMETER :: surface_pars = RESHAPE( (/ &
    529                                       10.0_wp,       10.0_wp, 0.05_wp,     & !  1
    530                                       10.0_wp,       10.0_wp, 0.05_wp,     & !  2
    531                                       20.0_wp,       15.0_wp, 0.03_wp,     & !  3
    532                                       20.0_wp,       15.0_wp, 0.03_wp,     & !  4
    533                                       20.0_wp,       15.0_wp, 0.03_wp,     & !  5
    534                                       20.0_wp,       15.0_wp, 0.03_wp,     & !  6
    535                                       10.0_wp,       10.0_wp, 0.05_wp,     & !  7
    536                                       15.0_wp,       15.0_wp, 0.00_wp,     & !  8
    537                                       10.0_wp,       10.0_wp, 0.05_wp,     & !  9
    538                                       10.0_wp,       10.0_wp, 0.05_wp,     & ! 10
    539                                       10.0_wp,       10.0_wp, 0.05_wp,     & ! 11
    540                                       58.0_wp,       58.0_wp, 0.00_wp,     & ! 12
    541                                       10.0_wp,       10.0_wp, 0.05_wp,     & ! 13
    542                                     1.0E10_wp,     1.0E10_wp, 0.00_wp,     & ! 14
    543                                     1.0E10_wp,     1.0E10_wp, 0.00_wp,     & ! 15
    544                                       10.0_wp,       10.0_wp, 0.05_wp,     & ! 16
    545                                       10.0_wp,       10.0_wp, 0.05_wp,     & ! 17
    546                                       20.0_wp,       15.0_wp, 0.03_wp,     & ! 18
    547                                       20.0_wp,       15.0_wp, 0.03_wp,     & ! 19
    548                                        0.0_wp,        0.0_wp, 0.00_wp      & ! 20
    549                                       /), (/ 3, 20 /) )
    550 
    551 !
    552 !-- Root distribution (sum = 1)  level 1-8
    553 !--                               
    554     REAL(wp), DIMENSION(0:7,1:20), PARAMETER :: root_distribution = RESHAPE( (/&
     499
     500!
     501!-- Pavement classes
     502    CHARACTER(20), DIMENSION(0:7), PARAMETER :: pavement_type_name = (/ &
     503                                   'user defined        ',          & ! 0
     504                                   'asphalt             ',          & ! 1
     505                                   'concrete            ',          & ! 2
     506                                   'asphalt/concrete mix',          & ! 3
     507                                   'brick pavers        ',          & ! 4
     508                                   'cobblestone pavers  ',          & ! 5
     509                                   'sett pavers         ',          & ! 6
     510                                   'gravel pavers       '           & ! 7
     511                                                                 /)                                                                 
     512                                                                 
     513!
     514!-- Water classes
     515    CHARACTER(12), DIMENSION(0:4), PARAMETER :: water_type_name = (/ &
     516                                   'user defined',                   & ! 0
     517                                   'lake        ',                   & ! 1
     518                                   'river       ',                   & ! 2
     519                                   'ocean       ',                   & ! 3
     520                                   'pond        '                    & ! 4
     521                                                                  /)                                                                 
     522                                                                 
     523                                                                 !
     524!-- Land surface parameters according to the respective classes (vegetation_type)
     525
     526!
     527!-- Land surface parameters
     528!-- r_canopy_min,     lai,   c_veg,     g_d         z0,         z0h, lambda_s_s, lambda_s_u, f_sw_in,  c_surface
     529    REAL(wp), DIMENSION(0:9,1:18), PARAMETER :: vegetation_pars = RESHAPE( (/ &
     530          0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, 0.5E-2_wp,   0.5E-4_wp,     0.0_wp,     0.0_wp, 0.00_wp, 0.00_wp, & !  1
     531        180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,   0.25_wp,  0.25E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & !  2
     532        110.0_wp, 2.00_wp, 0.85_wp, 0.00_wp,   0.20_wp,  0.20E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & !  3
     533        500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,   2.00_wp,     2.00_wp,   280.0_wp,   196.0_wp, 0.03_wp, 0.00_wp, & !  4
     534        500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,   2.00_wp,     2.00_wp,   280.0_wp,   196.0_wp, 0.03_wp, 0.00_wp, & !  5
     535        175.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,   2.00_wp,     2.00_wp,   280.0_wp,   196.0_wp, 0.03_wp, 0.00_wp, & !  6
     536        240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp,   2.00_wp,     2.00_wp,   280.0_wp,   196.0_wp, 0.03_wp, 0.00_wp, & !  7
     537        100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp,   0.47_wp,  0.47E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & !  8
     538        250.0_wp, 0.05_wp, 0.00_wp, 0.00_wp,  0.013_wp, 0.013E-2_wp,   196.0_wp,   196.0_wp, 0.00_wp, 0.00_wp, & !  9
     539         80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp,  0.034_wp, 0.034E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & ! 10
     540        180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,    0.5_wp,  0.50E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & ! 11
     541        150.0_wp, 0.50_wp, 0.10_wp, 0.00_wp,   0.17_wp,  0.17E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & ! 12
     542          0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, 1.3E-3_wp,   1.3E-4_wp,   812.0_wp,   812.0_wp, 0.00_wp, 0.00_wp, & ! 13
     543        240.0_wp, 4.00_wp, 0.60_wp, 0.00_wp,   0.83_wp,  0.83E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & ! 14
     544        225.0_wp, 3.00_wp, 0.50_wp, 0.00_wp,   0.10_wp,  0.10E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & ! 15
     545        225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp,   0.25_wp,  0.25E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & ! 16
     546        250.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,   2.00_wp,     2.00_wp,   280.0_wp,   196.0_wp, 0.03_wp, 0.00_wp, & ! 17
     547        175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp,   1.10_wp,     1.10_wp,   280.0_wp,   196.0_wp, 0.03_wp, 0.00_wp  & ! 18
     548                                                               /), (/ 10, 18 /) )
     549
     550                                   
     551!
     552!-- Root distribution for default soil layer configuration (sum = 1)
     553!--                                level 1 - level 8
     554    REAL(wp), DIMENSION(0:7,1:18), PARAMETER :: root_distribution = RESHAPE( (/&
     555                                   1.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     &
     556                                   0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     & ! 1
    555557                                   0.035_wp, 0.069_wp, 0.069_wp, 0.108_wp,     &
    556                                    0.195_wp, 0.214_wp, 0.284_wp, 0.026_wp,     & ! 1
     558                                   0.195_wp, 0.214_wp, 0.284_wp, 0.026_wp,     & ! 2
    557559                                   0.050_wp, 0.100_wp, 0.100_wp, 0.136_wp,     &
    558                                    0.181_wp, 0.192_wp, 0.215_wp, 0.026_wp,     & ! 2
     560                                   0.181_wp, 0.192_wp, 0.215_wp, 0.026_wp,     & ! 3
    559561                                   0.038_wp, 0.075_wp, 0.075_wp, 0.111_wp,     &
    560                                    0.185_wp, 0.203_wp, 0.273_wp, 0.040_wp,     & ! 3
     562                                   0.185_wp, 0.203_wp, 0.273_wp, 0.040_wp,     & ! 4
    561563                                   0.038_wp, 0.075_wp, 0.075_wp, 0.110_wp,     &
    562                                    0.180_wp, 0.199_wp, 0.277_wp, 0.046_wp,     & ! 4
     564                                   0.180_wp, 0.199_wp, 0.277_wp, 0.046_wp,     & ! 5
    563565                                   0.035_wp, 0.069_wp, 0.069_wp, 0.105_wp,     &
    564                                    0.180_wp, 0.201_wp, 0.295_wp, 0.046_wp,     & ! 5
     566                                   0.180_wp, 0.201_wp, 0.295_wp, 0.046_wp,     & ! 6
    565567                                   0.035_wp, 0.072_wp, 0.072_wp, 0.105_wp,     &
    566                                    0.161_wp, 0.180_wp, 0.282_wp, 0.093_wp,     & ! 6
     568                                   0.161_wp, 0.180_wp, 0.282_wp, 0.093_wp,     & ! 7
    567569                                   0.040_wp, 0.077_wp, 0.077_wp, 0.112_wp,     &
    568                                    0.176_wp, 0.192_wp, 0.266_wp, 0.060_wp,     & ! 7
     570                                   0.176_wp, 0.192_wp, 0.266_wp, 0.060_wp,     & ! 8
    569571                                   0.142_wp, 0.286_wp, 0.286_wp, 0.286_wp,     &
    570                                    0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     & ! 8
     572                                   0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     & ! 9
    571573                                   0.068_wp, 0.134_wp, 0.134_wp, 0.177_wp,     &
    572                                    0.214_wp, 0.203_wp, 0.070_wp, 0.000_wp,     & ! 9
     574                                   0.214_wp, 0.203_wp, 0.070_wp, 0.000_wp,     & ! 10
    573575                                   0.035_wp, 0.068_wp, 0.068_wp, 0.108_wp,     &
    574                                    0.195_wp, 0.215_wp, 0.285_wp, 0.026_wp,     & ! 10
     576                                   0.195_wp, 0.215_wp, 0.285_wp, 0.026_wp,     & ! 11
    575577                                   0.025_wp, 0.048_wp, 0.048_wp, 0.078_wp,     &
    576                                    0.147_wp, 0.175_wp, 0.353_wp, 0.126_wp,     & ! 11
     578                                   0.147_wp, 0.175_wp, 0.353_wp, 0.126_wp,     & ! 12
    577579                                   0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     &
    578                                    0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     & ! 12
     580                                   0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     & ! 13
    579581                                   0.036_wp, 0.072_wp, 0.072_wp, 0.103_wp,     &
    580                                    0.163_wp, 0.180_wp, 0.273_wp, 0.074_wp,     & ! 13
    581                                    0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     &
    582                                    0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     & ! 14
    583                                    0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     &
    584                                    0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     & ! 15
     582                                   0.163_wp, 0.180_wp, 0.273_wp, 0.074_wp,     & ! 14
     583                                   0.032_wp, 0.066_wp, 0.066_wp, 0.100_wp,     &
     584                                   0.172_wp, 0.192_wp, 0.299_wp, 0.073_wp,     & ! 15
    585585                                   0.032_wp, 0.066_wp, 0.066_wp, 0.100_wp,     &
    586586                                   0.172_wp, 0.192_wp, 0.299_wp, 0.073_wp,     & ! 16
    587                                    0.032_wp, 0.066_wp, 0.066_wp, 0.100_wp,     &
    588                                    0.172_wp, 0.192_wp, 0.299_wp, 0.073_wp,     & ! 17
    589587                                   0.028_wp, 0.055_wp, 0.055_wp, 0.087_wp,     &
    590                                    0.166_wp, 0.195_wp, 0.348_wp, 0.066_wp,     & ! 18
     588                                   0.166_wp, 0.195_wp, 0.348_wp, 0.066_wp,     & ! 17
    591589                                   0.028_wp, 0.055_wp, 0.055_wp, 0.087_wp,     &
    592                                    0.166_wp, 0.195_wp, 0.348_wp, 0.066_wp,     & ! 19
    593                                    0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp,     &
    594                                    0.000_wp, 0.000_wp, 0.000_wp, 0.000_wp      & ! 20
    595                                  /), (/ 8, 20 /) )
     590                                   0.166_wp, 0.195_wp, 0.348_wp, 0.066_wp      & ! 18
     591                                 /), (/ 8, 18 /) )
    596592
    597593
     
    600596
    601597!
    602 !-- Soil parameters I           alpha_vg,      l_vg,    n_vg, gamma_w_sat
    603     REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: soil_pars = RESHAPE( (/     &
    604                                  3.83_wp,  1.250_wp, 1.38_wp,  6.94E-6_wp, & ! 1
    605                                  3.14_wp, -2.342_wp, 1.28_wp,  1.16E-6_wp, & ! 2
    606                                  0.83_wp, -0.588_wp, 1.25_wp,  0.26E-6_wp, & ! 3
    607                                  3.67_wp, -1.977_wp, 1.10_wp,  2.87E-6_wp, & ! 4
    608                                  2.65_wp,  2.500_wp, 1.10_wp,  1.74E-6_wp, & ! 5
    609                                  1.30_wp,  0.400_wp, 1.20_wp,  0.93E-6_wp, & ! 6
    610                                  0.00_wp,  0.00_wp,  0.00_wp,  0.57E-6_wp  & ! 7
    611                                  /), (/ 4, 7 /) )
    612 
    613 !
    614 !-- Soil parameters II              m_sat,     m_fc,   m_wilt,    m_res 
    615     REAL(wp), DIMENSION(0:3,1:7), PARAMETER :: m_soil_pars = RESHAPE( (/ &
    616                                  0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp, & ! 1
    617                                  0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp, & ! 2
    618                                  0.430_wp, 0.383_wp, 0.133_wp, 0.010_wp, & ! 3
    619                                  0.520_wp, 0.448_wp, 0.279_wp, 0.010_wp, & ! 4
    620                                  0.614_wp, 0.541_wp, 0.335_wp, 0.010_wp, & ! 5
    621                                  0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp, & ! 6
    622                                  0.472_wp, 0.323_wp, 0.171_wp, 0.000_wp  & ! 7
    623                                  /), (/ 4, 7 /) )
    624 
    625 
     598!-- Soil parameters             alpha_vg,      l_vg,    n_vg, gamma_w_sat,    m_sat,     m_fc,   m_wilt,    m_res
     599    REAL(wp), DIMENSION(0:7,1:6), PARAMETER :: soil_pars = RESHAPE( (/     &
     600                                 3.83_wp,  1.250_wp, 1.38_wp,  6.94E-6_wp, 0.403_wp, 0.244_wp, 0.059_wp, 0.025_wp,& ! 1
     601                                 3.14_wp, -2.342_wp, 1.28_wp,  1.16E-6_wp, 0.439_wp, 0.347_wp, 0.151_wp, 0.010_wp,& ! 2
     602                                 0.83_wp, -0.588_wp, 1.25_wp,  0.26E-6_wp, 0.430_wp, 0.383_wp, 0.133_wp, 0.010_wp,& ! 3
     603                                 3.67_wp, -1.977_wp, 1.10_wp,  2.87E-6_wp, 0.520_wp, 0.448_wp, 0.279_wp, 0.010_wp,& ! 4
     604                                 2.65_wp,  2.500_wp, 1.10_wp,  1.74E-6_wp, 0.614_wp, 0.541_wp, 0.335_wp, 0.010_wp,& ! 5
     605                                 1.30_wp,  0.400_wp, 1.20_wp,  0.93E-6_wp, 0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp & ! 6
     606                                 /), (/ 8, 6 /) )
     607
     608!
     609!-- TO BE FILLED
     610!-- Pavement parameters  depth,        z0,       z0h, lambda_h,      rho_c 
     611    REAL(wp), DIMENSION(0:4,1:7), PARAMETER :: pavement_pars = RESHAPE( (/ &
     612                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, & ! 1
     613                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, & ! 2
     614                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, & ! 3                                 
     615                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, & ! 4
     616                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, & ! 5
     617                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, & ! 6
     618                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp  & ! 7
     619                                                                 /), (/ 5, 7 /) )                             
     620
     621!
     622!-- TO BE FILLED
     623!-- Water parameters                    temperature,     z0,      z0h
     624    REAL(wp), DIMENSION(0:2,1:4), PARAMETER :: water_pars = RESHAPE( (/ &
     625                                           283.0_wp, 0.01_wp, 0.001_wp, & ! 1
     626                                           283.0_wp, 0.01_wp, 0.001_wp, & ! 2
     627                                           283.0_wp, 0.01_wp, 0.001_wp, & ! 3
     628                                           283.0_wp, 10.01_wp, 0.001_wp & ! 4
     629                                                                /), (/ 3, 4 /) )                                                                     
     630                                                                                                                                     
    626631    SAVE
    627632
     
    637642           lsm_header, lsm_init, lsm_init_arrays, lsm_parin, lsm_soil_model,   &
    638643           lsm_swap_timelevel, lsm_read_restart_data, lsm_last_actions
    639 !
     644! !vegetat
    640645!-- Public parameters, constants and initial values
    641646    PUBLIC aero_resist_kray, skip_time_do_lsm
     
    754759          unit = 'K'   
    755760             
    756        CASE ( 'lai*', 'c_liq*', 'c_soil*', 'c_veg*', 'ghf_eb*', 'm_liq_eb*',&
    757               'qsws_eb*', 'qsws_liq_eb*', 'qsws_soil_eb*', 'qsws_veg_eb*',  &
    758               'r_a*', 'r_s*', 'shf_eb*' )
     761       CASE ( 'lai*', 'c_liq*', 'c_soil*', 'c_veg*', 'ghf*', 'm_liq*',&
     762              'qsws_liq*', 'qsws_soil*', 'qsws_veg*',  &
     763              'r_a*', 'r_s*' )
    759764          IF ( k == 0  .OR.  data_output(i)(ilen-2:ilen) /= '_xy' )  THEN
    760765             message_string = 'illegal value for data_output: "' //         &
     
    783788             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    784789          ENDIF
    785           IF ( TRIM( var ) == 'ghf_eb*'  .AND.  .NOT.  land_surface )  THEN
     790          IF ( TRIM( var ) == 'ghf*'  .AND.  .NOT.  land_surface )  THEN
    786791             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    787792                              'res land_surface = .TRUE.'
    788793             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    789794          ENDIF
    790           IF ( TRIM( var ) == 'm_liq_eb*'  .AND.  .NOT.  land_surface )  THEN
     795          IF ( TRIM( var ) == 'm_liq*'  .AND.  .NOT.  land_surface )  THEN
    791796             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    792797                              'res land_surface = .TRUE.'
    793798             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    794799          ENDIF
    795           IF ( TRIM( var ) == 'qsws_eb*'  .AND.  .NOT.  land_surface )  THEN
    796              message_string = 'output of "' // TRIM( var ) // '" requi' //  &
    797                               'res land_surface = .TRUE.'
    798              CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    799           ENDIF
    800           IF ( TRIM( var ) == 'qsws_liq_eb*'  .AND.  .NOT. land_surface )   &
     800          IF ( TRIM( var ) == 'qsws_liq*'  .AND.  .NOT. land_surface )   &
    801801          THEN
    802802             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     
    804804             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    805805          ENDIF
    806           IF ( TRIM( var ) == 'qsws_soil_eb*'  .AND.  .NOT.  land_surface ) &
     806          IF ( TRIM( var ) == 'qsws_soil*'  .AND.  .NOT.  land_surface ) &
    807807          THEN
    808808             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     
    810810             CALL message( 'check_parameters', 'PA0404', 1, 2, 0, 6, 0 )
    811811          ENDIF
    812           IF ( TRIM( var ) == 'qsws_veg_eb*'  .AND.  .NOT. land_surface )   &
     812          IF ( TRIM( var ) == 'qsws_veg*'  .AND.  .NOT. land_surface )   &
    813813          THEN
    814814             message_string = 'output of "' // TRIM( var ) // '" requi' //  &
     
    833833          IF ( TRIM( var ) == 'c_soil*')  unit = 'none'
    834834          IF ( TRIM( var ) == 'c_veg*' )  unit = 'none'
    835           IF ( TRIM( var ) == 'ghf_eb*')  unit = 'W/m2'
    836           IF ( TRIM( var ) == 'm_liq_eb*'     )  unit = 'm'
    837           IF ( TRIM( var ) == 'qsws_eb*'      ) unit = 'W/m2'
    838           IF ( TRIM( var ) == 'qsws_liq_eb*'  ) unit = 'W/m2'
    839           IF ( TRIM( var ) == 'qsws_soil_eb*' ) unit = 'W/m2'
    840           IF ( TRIM( var ) == 'qsws_veg_eb*'  ) unit = 'W/m2'
     835          IF ( TRIM( var ) == 'ghf*')  unit = 'W/m2'
     836          IF ( TRIM( var ) == 'm_liq*'     )  unit = 'm'
     837          IF ( TRIM( var ) == 'qsws_liq*'  ) unit = 'W/m2'
     838          IF ( TRIM( var ) == 'qsws_soil*' ) unit = 'W/m2'
     839          IF ( TRIM( var ) == 'qsws_veg*'  ) unit = 'W/m2'
    841840          IF ( TRIM( var ) == 'r_a*')     unit = 's/m'     
    842841          IF ( TRIM( var ) == 'r_s*')     unit = 's/m'
    843           IF ( TRIM( var ) == 'shf_eb*')  unit = 'W/m2'
    844842             
    845843       CASE DEFAULT
     
    888886             dopr_index(var_count) = 89
    889887             dopr_unit     = 'K'
    890              hom(0:nzs-1,2,89,:)  = SPREAD( - zs, 2, statistic_regions+1 )
     888             hom(0:nzs-1,2,89,:)  = SPREAD( - zs(nzb_soil:nzt_soil), 2, statistic_regions+1 )
    891889             IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
    892890                dopr_initial_index(var_count) = 90
    893                 hom(0:nzs-1,2,90,:)   = SPREAD( - zs, 2, statistic_regions+1 )
     891                hom(0:nzs-1,2,90,:)   = SPREAD( - zs(nzb_soil:nzt_soil), 2, statistic_regions+1 )
    894892                data_output_pr(var_count)     = data_output_pr(var_count)(2:)
    895893             ENDIF
     
    906904             dopr_index(var_count) = 91
    907905             dopr_unit     = 'm3/m3'
    908              hom(0:nzs-1,2,91,:)  = SPREAD( - zs, 2, statistic_regions+1 )
     906             hom(0:nzs-1,2,91,:)  = SPREAD( - zs(nzb_soil:nzt_soil), 2, statistic_regions+1 )
    909907             IF ( data_output_pr(var_count)(1:1) == '#' )  THEN
    910908                dopr_initial_index(var_count) = 92
    911                 hom(0:nzs-1,2,92,:)   = SPREAD( - zs, 2, statistic_regions+1 )
     909                hom(0:nzs-1,2,92,:)   = SPREAD( - zs(nzb_soil:nzt_soil), 2, statistic_regions+1 )
    912910                data_output_pr(var_count)     = data_output_pr(var_count)(2:)
    913911             ENDIF
     
    943941
    944942    INTEGER(iwp) ::  k        !< running index, z-dimension
     943
     944!
     945!-- Check for a valid setting of surface_type. The default value is 'netcdf'.
     946!-- In that case, the surface types are read from NetCDF file
     947    IF ( TRIM( surface_type ) /= 'vegetation'  .AND.                           &
     948         TRIM( surface_type ) /= 'pavement'    .AND.                           &
     949         TRIM( surface_type ) /= 'water'       .AND.                           &
     950         TRIM( surface_type ) /= 'netcdf' )  THEN 
     951          message_string = 'unknown surface type surface_type = "' //     &
     952                        TRIM( surface_type ) // '"'
     953                CALL message( 'check_parameters', 'PA0019', 1, 2, 0, 6, 0 )
     954    ENDIF
     955
    945956!
    946957!-- Dirichlet boundary conditions are required as the surface fluxes are
     
    948959!-- model
    949960    IF ( bc_pt_b == 'neumann'  .OR.  bc_q_b == 'neumann' )  THEN
    950        message_string = 'lsm requires setting of'//                         &
    951                         'bc_pt_b = "dirichlet" and '//                      &
     961       message_string = 'lsm requires setting of'//                            &
     962                        'bc_pt_b = "dirichlet" and '//                         &
    952963                        'bc_q_b  = "dirichlet"'
    953964       CALL message( 'check_parameters', 'PA0399', 1, 2, 0, 6, 0 )
     
    955966
    956967    IF (  .NOT.  constant_flux_layer )  THEN
    957        message_string = 'lsm requires '//                                   &
     968       message_string = 'lsm requires '//                                      &
    958969                        'constant_flux_layer = .T.'
    959970       CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
    960971    ENDIF
    961972
    962     IF ( ( veg_type == 14  .OR.  veg_type == 15 ) .AND.                     &
    963            most_method == 'lookup' )  THEN
    964         WRITE( message_string, * ) 'veg_type = ', veg_type, ' is not ',     &
    965                                    'allowed in combination with ',          &
    966                                    'most_method = ', most_method
    967        CALL message( 'check_parameters', 'PA0417', 1, 2, 0, 6, 0 )
     973    IF ( TRIM( surface_type ) == 'vegetation' )  THEN
     974   
     975       IF ( vegetation_type == 0 )  THEN
     976          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
     977             message_string = 'vegetation_type = 0 (user defined)'//              &
     978                              'requires setting of min_canopy_resistance'//       &
     979                              '/= 9999999.9'
     980             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
     981          ENDIF
     982
     983          IF ( leaf_area_index == 9999999.9_wp )  THEN
     984             message_string = 'vegetation_type = 0 (user_defined)'//              &
     985                              'requires setting of leaf_area_index'//             &
     986                              '/= 9999999.9'
     987             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
     988          ENDIF
     989
     990          IF ( vegetation_coverage == 9999999.9_wp )  THEN
     991             message_string = 'vegetation_type = 0 (user_defined)'//              &
     992                              'requires setting of vegetation_coverage'//         &
     993                              '/= 9999999.9'
     994                CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
     995          ENDIF
     996
     997          IF ( canopy_resistance_coefficient == 9999999.9_wp)  THEN
     998             message_string = 'vegetation_type = 0 (user_defined)'//              &
     999                              'requires setting of'//                             &
     1000                              'canopy_resistance_coefficient /= 9999999.9'
     1001             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
     1002          ENDIF
     1003
     1004          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
     1005             message_string = 'vegetation_type = 0 (user_defined)'//              &
     1006                              'requires setting of lambda_surface_stable'//       &
     1007                              '/= 9999999.9'
     1008             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
     1009          ENDIF
     1010
     1011          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
     1012             message_string = 'vegetation_type = 0 (user_defined)'//              &
     1013                              'requires setting of lambda_surface_unstable'//     &
     1014                              '/= 9999999.9'
     1015             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
     1016          ENDIF
     1017
     1018          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
     1019             message_string = 'vegetation_type = 0 (user_defined)'//              &
     1020                              'requires setting of f_shortwave_incoming'//        &
     1021                              '/= 9999999.9'
     1022             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
     1023          ENDIF
     1024
     1025          IF ( z0_vegetation == 9999999.9_wp )  THEN
     1026             message_string = 'vegetation_type = 0 (user_defined)'//              &
     1027                              'requires setting of z0_vegetation'//               &
     1028                              '/= 9999999.9'
     1029             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
     1030          ENDIF
     1031
     1032          IF ( z0h_vegetation == 9999999.9_wp )  THEN
     1033             message_string = 'vegetation_type = 0 (user_defined)'//              &
     1034                              'requires setting of z0h_vegetation'//              &
     1035                              '/= 9999999.9'
     1036             CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
     1037          ENDIF
     1038       ENDIF
     1039       
    9681040    ENDIF
    969 
    970     IF ( veg_type == 0 )  THEN
    971        IF ( min_canopy_resistance == 9999999.9_wp )  THEN
    972           message_string = 'veg_type = 0 (user defined)'//                  &
    973                            'requires setting of min_canopy_resistance'//    &
    974                            '/= 9999999.9'
    975           CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    976        ENDIF
    977 
    978        IF ( leaf_area_index == 9999999.9_wp )  THEN
    979           message_string = 'veg_type = 0 (user_defined)'//                  &
    980                            'requires setting of leaf_area_index'//          &
    981                            '/= 9999999.9'
    982           CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    983        ENDIF
    984 
    985        IF ( vegetation_coverage == 9999999.9_wp )  THEN
    986           message_string = 'veg_type = 0 (user_defined)'//                  &
    987                            'requires setting of vegetation_coverage'//      &
    988                            '/= 9999999.9'
    989              CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    990        ENDIF
    991 
    992        IF ( canopy_resistance_coefficient == 9999999.9_wp)  THEN
    993           message_string = 'veg_type = 0 (user_defined)'//                  &
    994                            'requires setting of'//                          &
    995                            'canopy_resistance_coefficient /= 9999999.9'
    996           CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    997        ENDIF
    998 
    999        IF ( lambda_surface_stable == 9999999.9_wp )  THEN
    1000           message_string = 'veg_type = 0 (user_defined)'//                  &
    1001                            'requires setting of lambda_surface_stable'//    &
    1002                            '/= 9999999.9'
    1003           CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    1004        ENDIF
    1005 
    1006        IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
    1007           message_string = 'veg_type = 0 (user_defined)'//                  &
    1008                            'requires setting of lambda_surface_unstable'//  &
    1009                            '/= 9999999.9'
    1010           CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    1011        ENDIF
    1012 
    1013        IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
    1014           message_string = 'veg_type = 0 (user_defined)'//                  &
    1015                            'requires setting of f_shortwave_incoming'//     &
    1016                            '/= 9999999.9'
    1017           CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    1018        ENDIF
    1019 
    1020        IF ( z0_eb == 9999999.9_wp )  THEN
    1021           message_string = 'veg_type = 0 (user_defined)'//                  &
    1022                            'requires setting of z0_eb'//                    &
    1023                            '/= 9999999.9'
    1024           CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    1025        ENDIF
    1026 
    1027        IF ( z0h_eb == 9999999.9_wp )  THEN
    1028           message_string = 'veg_type = 0 (user_defined)'//                  &
    1029                            'requires setting of z0h_eb'//                   &
    1030                            '/= 9999999.9'
    1031           CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    1032        ENDIF
    1033 
    1034 
     1041   
     1042    IF ( TRIM( surface_type ) == 'water' )  THEN
     1043
     1044       IF ( TRIM( most_method ) == 'lookup' )  THEN   
     1045          WRITE( message_string, * ) 'surface_type = ', surface_type,          &
     1046                                     ' is not allowed in combination with ',   &
     1047                                     'most_method = ', most_method
     1048          CALL message( 'check_parameters', 'PA0417', 1, 2, 0, 6, 0 )
     1049       ENDIF
     1050
     1051       IF ( water_type == 0 )  THEN 
     1052       
     1053          IF ( z0_water == 9999999.9_wp )  THEN
     1054             message_string = 'water_type = 0 (user_defined)'//                &
     1055                              'requires setting of z0_water'//                 &
     1056                              '/= 9999999.9'
     1057             CALL message( 'check_parameters', 'PA0415', 1, 2, 0, 6, 0 )
     1058          ENDIF
     1059
     1060          IF ( z0h_water == 9999999.9_wp )  THEN
     1061             message_string = 'water_type = 0 (user_defined)'//                &
     1062                              'requires setting of z0h_water'//                &
     1063                              '/= 9999999.9'
     1064             CALL message( 'check_parameters', 'PA0392', 1, 2, 0, 6, 0 )
     1065          ENDIF
     1066         
     1067          IF ( water_temperature == 9999999.9_wp )  THEN
     1068             message_string = 'water_type = 0 (user_defined)'//                &
     1069                              'requires setting of water_temperature'//        &
     1070                              '/= 9999999.9'
     1071             CALL message( 'check_parameters', 'PA0379', 1, 2, 0, 6, 0 )
     1072          ENDIF       
     1073         
     1074       ENDIF
     1075       
    10351076    ENDIF
    1036 
     1077   
     1078   
     1079     
     1080    IF ( TRIM( surface_type ) == 'pavement' )  THEN
     1081
     1082       IF ( pavement_type == 0 )  THEN 
     1083       
     1084          IF ( z0_pavement == 9999999.9_wp )  THEN
     1085             message_string = 'pavement_type = 0 (user_defined)'//             &
     1086                              'requires setting of z0_pavement'//              &
     1087                              '/= 9999999.9'
     1088             CALL message( 'check_parameters', 'PA0352', 1, 2, 0, 6, 0 )
     1089          ENDIF
     1090
     1091          IF ( z0h_pavement == 9999999.9_wp )  THEN
     1092             message_string = 'pavement_type = 0 (user_defined)'//             &
     1093                              'requires setting of z0h_pavement'//             &
     1094                              '/= 9999999.9'
     1095             CALL message( 'check_parameters', 'PA0353', 1, 2, 0, 6, 0 )
     1096          ENDIF
     1097
     1098           IF ( pavement_depth == 9999999.9_wp )  THEN
     1099             message_string = 'pavement_type = 0 (user_defined)'//             &
     1100                              'requires setting of pavement_depth'//           &
     1101                              '/= 9999999.9'
     1102             CALL message( 'check_parameters', 'PA0341', 1, 2, 0, 6, 0 )
     1103          ENDIF         
     1104         
     1105           IF ( pavement_heat_conduct == 9999999.9_wp )  THEN
     1106             message_string = 'pavement_type = 0 (user_defined)'//             &
     1107                              'requires setting of pavement_heat_conduct'//    &
     1108                              '/= 9999999.9'
     1109             CALL message( 'check_parameters', 'PA0342', 1, 2, 0, 6, 0 )
     1110          ENDIF
     1111         
     1112           IF ( pavement_heat_capacity == 9999999.9_wp )  THEN
     1113             message_string = 'pavement_type = 0 (user_defined)'//             &
     1114                              'requires setting of pavement_heat_capacity'//   &
     1115                              '/= 9999999.9'
     1116             CALL message( 'check_parameters', 'PA0139', 1, 2, 0, 6, 0 )
     1117          ENDIF
     1118
     1119       ENDIF
     1120   
     1121    ENDIF
     1122   
    10371123    IF ( soil_type == 0 )  THEN
    10381124
     
    10941180
    10951181    ENDIF
    1096 !
    1097 !-- Check for proper setting of soil moisture, must not be larger than its
    1098 !-- saturation value.
    1099     DO  k = nzb_soil, nzt_soil
    1100        IF ( soil_moisture(k) > m_soil_pars(0,soil_type) )  THEN
    1101           message_string = 'soil_moisture must not exceed its saturation' // &
    1102                            ' value'
    1103           CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
    1104        ENDIF
    1105     ENDDO
    1106 
    1107     IF (  .NOT.  radiation )  THEN
    1108        message_string = 'lsm requires '//                                   &
    1109                         'radiation = .T.'
    1110        CALL message( 'check_parameters', 'PA0400', 1, 2, 0, 6, 0 )
    1111     ENDIF
    1112        
    1113 
    11141182
    11151183!
     
    11171185!-- root fraction is prescribed
    11181186    nzb_soil = 0
    1119     nzt_soil = 0
     1187    nzt_soil = -1
    11201188    IF ( ALL( zs == 9999999.9_wp ) )  THEN
    11211189       nzt_soil = 7
     
    11361204    ENDIF
    11371205
    1138     IF ( veg_type == 0 )  THEN
     1206    IF ( vegetation_type == 0 )  THEN
    11391207       IF ( SUM( root_fraction(nzb_soil:nzt_soil) ) /= 1.0_wp )  THEN
    1140           message_string = 'veg_type = 0 (user_defined)'//                  &
    1141                            'requires setting of root_fraction'//            &
     1208          message_string = 'vegetation_type = 0 (user_defined)'//              &
     1209                           'requires setting of root_fraction'//               &
    11421210                           '/= 9999999.9 and SUM(root_fraction) = 1'
    11431211          CALL message( 'check_parameters', 'PA0401', 1, 2, 0, 6, 0 )
    11441212       ENDIF
    1145     ENDIF
     1213    ENDIF   
    11461214   
     1215   
     1216!
     1217!-- Check for proper setting of soil moisture, must not be larger than its
     1218!-- saturation value.
     1219    DO  k = nzb_soil, nzt_soil
     1220       IF ( soil_moisture(k) > soil_pars(4,soil_type) )  THEN
     1221          message_string = 'soil_moisture must not exceed its saturation' // &
     1222                           ' value'
     1223          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
     1224       ENDIF
     1225    ENDDO
    11471226 
    11481227 END SUBROUTINE lsm_check_parameters
     
    11551234 SUBROUTINE lsm_energy_balance( horizontal, l )
    11561235
     1236 
     1237    USE pegrid
    11571238
    11581239    IMPLICIT NONE
     
    11831264                coef_1,      & !< coef. for prognostic equation
    11841265                coef_2,      & !< coef. for prognostic equation
    1185                 f_qsws,      & !< factor for qsws_eb
    1186                 f_qsws_veg,  & !< factor for qsws_veg_eb
    1187                 f_qsws_soil, & !< factor for qsws_soil_eb
    1188                 f_qsws_liq,  & !< factor for qsws_liq_eb
    1189                 f_shf,       & !< factor for shf_eb
    1190                 lambda_surface, & !< Current value of lambda_surface
    1191                 m_liq_eb_max,   & !< maxmimum value of the liq. water reservoir
     1266                f_qsws,      & !< factor for qsws
     1267                f_qsws_veg,  & !< factor for qsws_veg
     1268                f_qsws_soil, & !< factor for qsws_soil
     1269                f_qsws_liq,  & !< factor for qsws_liq
     1270                f_shf,       & !< factor for shf
     1271                lambda_surface, & !< Current value of lambda_surface (W/m2/K)
     1272                m_liq_max,   & !< maxmimum value of the liq. water reservoir
    11921273                pt1,         & !< potential temperature at first grid level
    11931274                qv1            !< specific humidity at first grid level
     
    11961277    TYPE(surf_type_lsm), POINTER ::  surf_t_surface_p
    11971278    TYPE(surf_type_lsm), POINTER ::  surf_tt_surface_m
    1198     TYPE(surf_type_lsm), POINTER ::  surf_m_liq_eb
    1199     TYPE(surf_type_lsm), POINTER ::  surf_m_liq_eb_p
    1200     TYPE(surf_type_lsm), POINTER ::  surf_tm_liq_eb_m
     1279    TYPE(surf_type_lsm), POINTER ::  surf_m_liq
     1280    TYPE(surf_type_lsm), POINTER ::  surf_m_liq_p
     1281    TYPE(surf_type_lsm), POINTER ::  surf_tm_liq_m
    12011282
    12021283    TYPE(surf_type_lsm), POINTER ::  surf_m_soil
     
    12111292       surf_t_surface_p  => t_surface_h_p
    12121293       surf_tt_surface_m => tt_surface_h_m
    1213        surf_m_liq_eb     => m_liq_eb_h
    1214        surf_m_liq_eb_p   => m_liq_eb_h_p
    1215        surf_tm_liq_eb_m  => tm_liq_eb_h_m
     1294       surf_m_liq     => m_liq_h
     1295       surf_m_liq_p   => m_liq_h_p
     1296       surf_tm_liq_m  => tm_liq_h_m
    12161297       surf_m_soil       => m_soil_h
    12171298       surf_t_soil       => t_soil_h
     
    12261307       surf_t_surface_p  => t_surface_v_p(l)
    12271308       surf_tt_surface_m => tt_surface_v_m(l)
    1228        surf_m_liq_eb     => m_liq_eb_v(l)
    1229        surf_m_liq_eb_p   => m_liq_eb_v_p(l)
    1230        surf_tm_liq_eb_m  => tm_liq_eb_v_m(l)
     1309       surf_m_liq     => m_liq_v(l)
     1310       surf_m_liq_p   => m_liq_v_p(l)
     1311       surf_tm_liq_m  => tm_liq_v_m(l)
    12311312       surf_m_soil       => m_soil_v(l)
    12321313       surf_t_soil       => t_soil_v(l)
     
    12631344!--    radiation can be obtained at respective height level
    12641345       k_rad = MERGE( 0, k + k_off, radiation_scheme /= 'rrtmg' )
    1265 !
    1266 !--    Set lambda_surface according to stratification between skin layer and soil
    1267        IF (  .NOT.  surf%pave_surface(m) )  THEN
    1268 
    1269           c_surface_tmp = c_surface
    1270 
    1271           IF ( surf_t_surface%var_1d(m) >= surf_t_soil%var_2d(nzb_soil,m))  THEN
    1272              lambda_surface = surf%lambda_surface_s(m)
    1273           ELSE
    1274              lambda_surface = surf%lambda_surface_u(m)
    1275           ENDIF
     1346
     1347!
     1348!--    Define heat conductivity between surface and soil depending on surface
     1349!--    type. For vegetation, a skin layer parameterization is used. For bare
     1350!--    soil and pavements, no skin layer is applied. In these cases, the
     1351!--    temperature is assumed to be constant between the surface and the first
     1352!--    soil layer. The heat conductivity is then derived from the soil/
     1353!--    pavement properties. For water surfaces, the conductivity is already set
     1354!--    to 1E10.
     1355!--    Moreover, the heat capacity is set. For bare soil the heat capacity is
     1356!--    the capacity of the uppermost soil layer
     1357       IF ( surf%lambda_surface_s(m) == 0.0_wp )  THEN
     1358       
     1359          lambda_h_sat = lambda_h_sm**(1.0_wp - surf%m_sat(nzb_soil,m)) *      &
     1360                         lambda_h_water ** surf_m_soil%var_2d(nzb_soil,m)
     1361                         
     1362          ke = 1.0_wp + LOG10( MAX( 0.1_wp, surf_m_soil%var_2d(nzb_soil,m) /   &
     1363                                                     surf%m_sat(nzb_soil,m) ) )                   
     1364                         
     1365          lambda_surface = (ke * (lambda_h_sat - lambda_h_dry) + lambda_h_dry )&
     1366                           * ddz_soil_layer(nzb_soil) * 2.0_wp
     1367         
     1368          c_surface_tmp = (rho_c_soil * ( 1.0_wp - surf%m_sat(nzb_soil,m) )    &
     1369                           + rho_c_water * surf_m_soil%var_2d(nzb_soil,m) )    &
     1370                           * dz_soil_layer(nzb_soil) * 0.25_wp               
     1371
     1372       ELSEIF ( surf_t_surface%var_1d(m) >= surf_t_soil%var_2d(nzb_soil,m)) &
     1373       THEN
     1374
     1375          lambda_surface = surf%lambda_surface_s(m)
     1376          c_surface_tmp = surf%c_surface(m)
     1377         
    12761378       ELSE
    12771379
    1278           c_surface_tmp = pave_heat_capacity * dz_soil(nzb_soil) * 0.5_wp
    1279           lambda_surface = pave_heat_conductivity * ddz_soil(nzb_soil)
    1280 
     1380          lambda_surface = surf%lambda_surface_u(m)
     1381          c_surface_tmp = surf%c_surface(m)
     1382         
    12811383       ENDIF
    12821384
     
    13331435       DO  ks = nzb_soil, nzt_soil
    13341436           m_total = m_total + surf%root_fr(ks,m)                              &
    1335                      * MAX( surf_m_soil%var_2d(ks,m), surf%m_wilt(m) )
     1437                     * MAX( surf_m_soil%var_2d(ks,m), surf%m_wilt(ks,m) )
    13361438       ENDDO
    13371439
    1338        IF ( m_total > surf%m_wilt(m)  .AND.                                    &
    1339             m_total < surf%m_fc(m) )  THEN
    1340           f2 = ( m_total - surf%m_wilt(m) ) /                                  &
    1341                ( surf%m_fc(m) - surf%m_wilt(m) )
    1342        ELSEIF ( m_total >= surf%m_fc(m) )  THEN
     1440!
     1441!--    The calculation of f2 is based on only one wilting point value for all
     1442!--    soil layers. The value at k=nzb_soil is used here as a proxy but might
     1443!--    need refinement in the future.
     1444       IF ( m_total > surf%m_wilt(nzb_soil,m)  .AND.                           &
     1445            m_total < surf%m_fc(nzb_soil,m) )  THEN
     1446          f2 = ( m_total - surf%m_wilt(nzb_soil,m) ) /                         &
     1447               ( surf%m_fc(nzb_soil,m) - surf%m_wilt(nzb_soil,m) )
     1448       ELSEIF ( m_total >= surf%m_fc(nzb_soil,m) )  THEN
    13431449          f2 = 1.0_wp
    13441450       ELSE
     
    13661472                              ( surf%lai(m) * f1 * f2 * f3 + 1.0E-20_wp )
    13671473!
    1368 !--    Third step: calculate bare soil resistance r_soil. The Clapp &
    1369 !--    Hornberger parametrization does not consider c_veg.
    1370        IF ( surf%soil_type_2d(m) /= 7 )  THEN
    1371           m_min = surf%c_veg(m) * surf%m_wilt(m) +                             &
    1372                          ( 1.0_wp - surf%c_veg(m) ) * surf%m_res(m)
    1373        ELSE
    1374           m_min = surf%m_wilt(m)
    1375        ENDIF
    1376 
    1377        f2 = ( surf_m_soil%var_2d(nzb_soil,m) - m_min ) / ( surf%m_fc(m) - m_min )
     1474!--    Third step: calculate bare soil resistance r_soil.
     1475       m_min = surf%c_veg(m) * surf%m_wilt(nzb_soil,m) +                       &
     1476                         ( 1.0_wp - surf%c_veg(m) ) * surf%m_res(nzb_soil,m)
     1477
     1478
     1479       f2 = ( surf_m_soil%var_2d(nzb_soil,m) - m_min ) / ( surf%m_fc(nzb_soil,m) - m_min )
    13781480       f2 = MAX( f2, 1.0E-20_wp )
    13791481       f2 = MIN( f2, 1.0_wp     )
    13801482
    13811483       surf%r_soil(m) = surf%r_soil_min(m) / f2
    1382 
     1484       
    13831485!
    13841486!--    Calculate the maximum possible liquid water amount on plants and
     
    13881490!--    Noilhan & Planton (1989), while the ECMWF formulation is used for
    13891491!--    vegetated surfaces and bare soils.
    1390        IF ( surf%pave_surface(m) )  THEN
    1391           m_liq_eb_max = m_max_depth * 5.0_wp
    1392           surf%c_liq(m) = MIN( 1.0_wp, ( surf_m_liq_eb%var_1d(m) / m_liq_eb_max)**0.67 )
     1492       IF ( surf%pavement_surface(m) )  THEN
     1493          m_liq_max = m_max_depth * 5.0_wp
     1494          surf%c_liq(m) = MIN( 1.0_wp, ( surf_m_liq%var_1d(m) / m_liq_max)**0.67 )
    13931495       ELSE
    1394           m_liq_eb_max = m_max_depth * ( surf%c_veg(m) * surf%lai(m)&
    1395                             + ( 1.0_wp - surf%c_veg(m) ) )
    1396           surf%c_liq(m) = MIN( 1.0_wp, surf_m_liq_eb%var_1d(m) / m_liq_eb_max )
     1496          m_liq_max = m_max_depth * ( surf%c_veg(m) * surf%lai(m)              &
     1497                      + ( 1.0_wp - surf%c_veg(m) ) )
     1498          surf%c_liq(m) = MIN( 1.0_wp, surf_m_liq%var_1d(m) / m_liq_max )
    13971499       ENDIF
    13981500!
     
    14151517          f_qsws_soil = 0.0_wp
    14161518          f_qsws_liq  = rho_lv / surf%r_a(m)
    1417        ELSEIF ( surf%pave_surface (m) )  THEN
     1519       ELSEIF ( surf%pavement_surface (m) )  THEN
    14181520          f_qsws_veg  = 0.0_wp
    14191521          f_qsws_soil = 0.0_wp
     
    14281530                              surf%r_a(m)
    14291531       ENDIF
    1430 !
    1431 !--       If soil moisture is below wilting point, plants do no longer
    1432 !--       transpirate.
    1433 !           IF ( m_soil(k,j,i) < m_wilt(j,i) )  THEN
    1434 !              f_qsws_veg = 0.0_wp
    1435 !           ENDIF
    14361532
    14371533       f_shf  = rho_cp / surf%r_a(m)
     
    15351631                     ( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) )
    15361632
    1537        surf%ghf_eb(m) = lambda_surface * ( surf_t_surface_p%var_1d(m)          &
     1633       surf%ghf(m) = lambda_surface * ( surf_t_surface_p%var_1d(m)          &
    15381634                                             - surf_t_soil%var_2d(nzb_soil,m) )
    15391635
    1540        surf%shf_eb(m) = - f_shf * ( pt1 - surf%pt_surface(m) )
    1541 
    1542        surf%shf(m) = surf%shf_eb(m) / rho_cp
     1636       surf%shf(m) = - f_shf * ( pt1 - surf%pt_surface(m) ) / cp
     1637
    15431638
    15441639       IF ( humidity )  THEN
    1545           surf%qsws_eb(m)  = - f_qsws * ( qv1 - q_s + dq_s_dt                  &
     1640          surf%qsws(m)  = - f_qsws * ( qv1 - q_s + dq_s_dt                  &
    15461641                          * surf_t_surface%var_1d(m) - dq_s_dt *               &
    15471642                            surf_t_surface_p%var_1d(m) )
    15481643
    1549           surf%qsws(m) = surf%qsws_eb(m) / rho_lv
    1550 
    1551           surf%qsws_veg_eb(m)  = - f_qsws_veg  * ( qv1 - q_s                   &
     1644          surf%qsws_veg(m)  = - f_qsws_veg  * ( qv1 - q_s                   &
    15521645                              + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt   &
    15531646                              * surf_t_surface_p%var_1d(m) )
    15541647
    1555           surf%qsws_soil_eb(m) = - f_qsws_soil * ( qv1 - q_s                   &
     1648          surf%qsws_soil(m) = - f_qsws_soil * ( qv1 - q_s                   &
    15561649                              + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt   &
    15571650                              * surf_t_surface_p%var_1d(m) )
    15581651
    1559           surf%qsws_liq_eb(m)  = - f_qsws_liq  * ( qv1 - q_s                   &
     1652          surf%qsws_liq(m)  = - f_qsws_liq  * ( qv1 - q_s                   &
    15601653                              + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt   &
    15611654                              * surf_t_surface_p%var_1d(m) )
     
    15641657!
    15651658!--    Calculate the true surface resistance
    1566        IF ( surf%qsws_eb(m) == 0.0_wp )  THEN
     1659       IF ( .NOT.  humidity )  THEN
    15671660          surf%r_s(m) = 1.0E10_wp
    15681661       ELSE
     
    15701663                          * surf_t_surface%var_1d(m) - dq_s_dt *               &
    15711664                            surf_t_surface_p%var_1d(m) ) /                     &
    1572                             surf%qsws_eb(m) - surf%r_a(m)
     1665                            surf%qsws(m) - surf%r_a(m)
    15731666       ENDIF
    15741667
     
    15781671       IF ( humidity )  THEN
    15791672!
    1580 !--       If precipitation is activated, add rain water to qsws_liq_eb
    1581 !--       and qsws_soil_eb according the the vegetation coverage.
     1673!--       If precipitation is activated, add rain water to qsws_liq
     1674!--       and qsws_soil according the the vegetation coverage.
    15821675!--       precipitation_rate is given in mm.
    15831676          IF ( precipitation )  THEN
     
    15871680!--          Otherwise, add the water to soil. In case of
    15881681!--          pavements, the exceeding water amount is implicitely removed
    1589 !--          as runoff as qsws_soil_eb is then not used in the soil model
    1590              IF ( surf_m_liq_eb%var_1d(m) /= m_liq_eb_max )  THEN
    1591                 surf%qsws_liq_eb(m) = surf%qsws_liq_eb(m)                      &
     1682!--          as runoff as qsws_soil is then not used in the soil model
     1683             IF ( surf_m_liq%var_1d(m) /= m_liq_max )  THEN
     1684                surf%qsws_liq(m) = surf%qsws_liq(m)                      &
    15921685                                 + surf%c_veg(m) * prr(k+k_off,j+j_off,i+i_off)&
    15931686                                 * hyrho(k+k_off)                              &
    15941687                                 * 0.001_wp * rho_l * l_v
    15951688             ELSE
    1596                 surf%qsws_soil_eb(m) = surf%qsws_soil_eb(m)                    &
     1689                surf%qsws_soil(m) = surf%qsws_soil(m)                    &
    15971690                                 + surf%c_veg(m) * prr(k+k_off,j+j_off,i+i_off)&
    15981691                                 * hyrho(k+k_off)                              &
     
    16021695!--          Add precipitation to bare soil according to the bare soil
    16031696!--          coverage.
    1604              surf%qsws_soil_eb(m) = surf%qsws_soil_eb(m) + ( 1.0_wp            &
     1697             surf%qsws_soil(m) = surf%qsws_soil(m) + ( 1.0_wp            &
    16051698                               - surf%c_veg(m) ) * prr(k+k_off,j+j_off,i+i_off)&
    16061699                               * hyrho(k+k_off)                                &
     
    16101703!
    16111704!--       If the air is saturated, check the reservoir water level
    1612           IF ( surf%qsws_eb(m) < 0.0_wp )  THEN
    1613 !
    1614 !--          Check if reservoir is full (avoid values > m_liq_eb_max)
    1615 !--          In that case, qsws_liq_eb goes to qsws_soil_eb. In this
    1616 !--          case qsws_veg_eb is zero anyway (because c_liq = 1),       
     1705          IF ( surf%qsws(m) < 0.0_wp )  THEN
     1706!
     1707!--          Check if reservoir is full (avoid values > m_liq_max)
     1708!--          In that case, qsws_liq goes to qsws_soil. In this
     1709!--          case qsws_veg is zero anyway (because c_liq = 1),       
    16171710!--          so that tend is zero and no further check is needed
    1618              IF ( surf_m_liq_eb%var_1d(m) == m_liq_eb_max )  THEN
    1619                 surf%qsws_soil_eb(m) = surf%qsws_soil_eb(m) + surf%qsws_liq_eb(m)
    1620 
    1621                 surf%qsws_liq_eb(m)  = 0.0_wp
     1711             IF ( surf_m_liq%var_1d(m) == m_liq_max )  THEN
     1712                surf%qsws_soil(m) = surf%qsws_soil(m) + surf%qsws_liq(m)
     1713
     1714                surf%qsws_liq(m)  = 0.0_wp
    16221715             ENDIF
    16231716
    16241717!
    1625 !--          In case qsws_veg_eb becomes negative (unphysical behavior),
     1718!--          In case qsws_veg becomes negative (unphysical behavior),
    16261719!--          let the water enter the liquid water reservoir as dew on the
    16271720!--          plant
    1628              IF ( surf%qsws_veg_eb(m) < 0.0_wp )  THEN
    1629                 surf%qsws_liq_eb(m) = surf%qsws_liq_eb(m) + surf%qsws_veg_eb(m)
    1630                 surf%qsws_veg_eb(m) = 0.0_wp
     1721             IF ( surf%qsws_veg(m) < 0.0_wp )  THEN
     1722                surf%qsws_liq(m) = surf%qsws_liq(m) + surf%qsws_veg(m)
     1723                surf%qsws_veg(m) = 0.0_wp
    16311724             ENDIF
    16321725          ENDIF                   
    16331726 
    1634           tend = - surf%qsws_liq_eb(m) * drho_l_lv
    1635           surf_m_liq_eb_p%var_1d(m) = surf_m_liq_eb%var_1d(m) + dt_3d *        &
     1727          surf%qsws(m) = surf%qsws(m) / l_v
     1728 
     1729          tend = - surf%qsws_liq(m) * drho_l_lv
     1730          surf_m_liq_p%var_1d(m) = surf_m_liq%var_1d(m) + dt_3d *        &
    16361731                                        ( tsc(2) * tend +                      &
    1637                                           tsc(3) * surf_tm_liq_eb_m%var_1d(m) )
     1732                                          tsc(3) * surf_tm_liq_m%var_1d(m) )
    16381733!
    16391734!--       Check if reservoir is overfull -> reduce to maximum
    16401735!--       (conservation of water is violated here)
    1641           surf_m_liq_eb_p%var_1d(m) = MIN( surf_m_liq_eb_p%var_1d(m),m_liq_eb_max )
     1736          surf_m_liq_p%var_1d(m) = MIN( surf_m_liq_p%var_1d(m),m_liq_max )
    16421737
    16431738!
    16441739!--       Check if reservoir is empty (avoid values < 0.0)
    16451740!--       (conservation of water is violated here)
    1646           surf_m_liq_eb_p%var_1d(m) = MAX( surf_m_liq_eb_p%var_1d(m), 0.0_wp )
    1647 !
    1648 !--       Calculate m_liq_eb tendencies for the next Runge-Kutta step
     1741          surf_m_liq_p%var_1d(m) = MAX( surf_m_liq_p%var_1d(m), 0.0_wp )
     1742!
     1743!--       Calculate m_liq tendencies for the next Runge-Kutta step
    16491744          IF ( timestep_scheme(1:5) == 'runge' )  THEN
    16501745             IF ( intermediate_timestep_count == 1 )  THEN
    1651                 surf_tm_liq_eb_m%var_1d(m) = tend
     1746                surf_tm_liq_m%var_1d(m) = tend
    16521747             ELSEIF ( intermediate_timestep_count <                            &
    16531748                      intermediate_timestep_count_max )  THEN
    1654                 surf_tm_liq_eb_m%var_1d(m) = -9.5625_wp * tend +               &
    1655                                               5.3125_wp * surf_tm_liq_eb_m%var_1d(m)
     1749                surf_tm_liq_m%var_1d(m) = -9.5625_wp * tend +               &
     1750                                              5.3125_wp * surf_tm_liq_m%var_1d(m)
    16561751             ENDIF
    16571752          ENDIF
     
    16841779!
    16851780!-- Calculate new roughness lengths (for water surfaces only)
    1686     IF ( horizontal )  CALL calc_z0_water_surface
     1781    IF ( horizontal  .AND.  .NOT. constant_roughness )  CALL calc_z0_water_surface
    16871782
    16881783    CONTAINS
     
    17991894       ENDIF
    18001895
    1801        WRITE( io, 4 ) TRIM( veg_type_name(veg_type) ),                         &
     1896       WRITE( io, 4 ) TRIM( vegetation_type_name(vegetation_type) ),           &
    18021897                        TRIM (soil_type_name(soil_type) )
    18031898       WRITE( io, 5 ) TRIM( soil_depth_chr ), TRIM( t_soil_chr ),              &
     
    18291924    SUBROUTINE lsm_init
    18301925   
    1831 
     1926       USE control_parameters,                                                 &
     1927           ONLY:  message_string
     1928   
    18321929       IMPLICIT NONE
    18331930
     
    18641961       tt_soil_h_m%var_2d    = 0.0_wp
    18651962       tm_soil_h_m%var_2d    = 0.0_wp
    1866        tm_liq_eb_h_m%var_1d  = 0.0_wp
     1963       tm_liq_h_m%var_1d  = 0.0_wp
    18671964       surf_lsm_h%c_liq = 0.0_wp
    18681965
    1869        surf_lsm_h%ghf_eb = 0.0_wp
    1870        surf_lsm_h%shf_eb = rho_cp * surf_lsm_h%shf
    1871 
    1872        IF ( humidity )  THEN
    1873           surf_lsm_h%qsws_eb = rho_lv * surf_lsm_h%qsws
    1874        ELSE
    1875           surf_lsm_h%qsws_eb = 0.0_wp
    1876        ENDIF
    1877 
    1878        surf_lsm_h%qsws_liq_eb  = 0.0_wp
    1879        surf_lsm_h%qsws_soil_eb = 0.0_wp
    1880        surf_lsm_h%qsws_veg_eb  = 0.0_wp
     1966       surf_lsm_h%ghf = 0.0_wp
     1967
     1968       surf_lsm_h%qsws_liq  = 0.0_wp
     1969       surf_lsm_h%qsws_soil = 0.0_wp
     1970       surf_lsm_h%qsws_veg  = 0.0_wp
    18811971
    18821972       surf_lsm_h%r_a        = 50.0_wp
     
    18901980          tt_soil_v_m(l)%var_2d    = 0.0_wp
    18911981          tm_soil_v_m(l)%var_2d    = 0.0_wp
    1892           tm_liq_eb_v_m(l)%var_1d  = 0.0_wp
     1982          tm_liq_v_m(l)%var_1d  = 0.0_wp
    18931983          surf_lsm_v(l)%c_liq = 0.0_wp
    18941984
    1895           surf_lsm_v(l)%ghf_eb = 0.0_wp
    1896           surf_lsm_v(l)%shf_eb = rho_cp * surf_lsm_v(l)%shf
    1897 
    1898           IF ( humidity )  THEN
    1899              surf_lsm_v(l)%qsws_eb = rho_lv * surf_lsm_v(l)%qsws
    1900           ELSE
    1901              surf_lsm_v(l)%qsws_eb = 0.0_wp
    1902           ENDIF
    1903 
    1904           surf_lsm_v(l)%qsws_liq_eb  = 0.0_wp
    1905           surf_lsm_v(l)%qsws_soil_eb = 0.0_wp
    1906           surf_lsm_v(l)%qsws_veg_eb  = 0.0_wp
     1985          surf_lsm_v(l)%ghf = 0.0_wp
     1986
     1987          surf_lsm_v(l)%qsws_liq  = 0.0_wp
     1988          surf_lsm_v(l)%qsws_soil = 0.0_wp
     1989          surf_lsm_v(l)%qsws_veg  = 0.0_wp
    19071990
    19081991          surf_lsm_v(l)%r_a        = 50.0_wp
     
    19111994          surf_lsm_v(l)%r_soil     = 0.0_wp
    19121995       ENDDO
    1913 
    1914        
     1996   
    19151997!
    19161998!--    Allocate 3D soil model arrays
    19171999!--    First, for horizontal surfaces
     2000       ALLOCATE ( surf_lsm_h%alpha_vg(nzb_soil:nzt_soil,1:surf_lsm_h%ns)    )
     2001       ALLOCATE ( surf_lsm_h%gamma_w_sat(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
     2002       ALLOCATE ( surf_lsm_h%lambda_h(nzb_soil:nzt_soil,1:surf_lsm_h%ns)    )
     2003       ALLOCATE ( surf_lsm_h%l_vg(nzb_soil:nzt_soil,1:surf_lsm_h%ns)        )
     2004       ALLOCATE ( surf_lsm_h%m_fc(nzb_soil:nzt_soil,1:surf_lsm_h%ns)        )
     2005       ALLOCATE ( surf_lsm_h%m_res(nzb_soil:nzt_soil,1:surf_lsm_h%ns)       )
     2006       ALLOCATE ( surf_lsm_h%m_sat(nzb_soil:nzt_soil,1:surf_lsm_h%ns)       )
     2007       ALLOCATE ( surf_lsm_h%m_wilt(nzb_soil:nzt_soil,1:surf_lsm_h%ns)      )
     2008       ALLOCATE ( surf_lsm_h%n_vg(nzb_soil:nzt_soil,1:surf_lsm_h%ns)        )
     2009       ALLOCATE ( surf_lsm_h%rho_c_total(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
    19182010       ALLOCATE ( surf_lsm_h%root_fr(nzb_soil:nzt_soil,1:surf_lsm_h%ns)     )
    1919        ALLOCATE ( surf_lsm_h%lambda_h(nzb_soil:nzt_soil,1:surf_lsm_h%ns)    )
    1920        ALLOCATE ( surf_lsm_h%rho_c_total(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
    1921 
    1922        surf_lsm_h%lambda_h = 0.0_wp
     2011   
     2012       surf_lsm_h%lambda_h     = 0.0_wp
    19232013!
    19242014!--    If required, allocate humidity-related variables for the soil model
     
    19322022!--    For vertical surfaces
    19332023       DO  l = 0, 3
     2024          ALLOCATE ( surf_lsm_v(l)%alpha_vg(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)    )
     2025          ALLOCATE ( surf_lsm_v(l)%gamma_w_sat(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
     2026          ALLOCATE ( surf_lsm_v(l)%lambda_h(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)    )
     2027          ALLOCATE ( surf_lsm_v(l)%l_vg(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)        )
     2028          ALLOCATE ( surf_lsm_v(l)%m_fc(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)        )
     2029          ALLOCATE ( surf_lsm_v(l)%m_res(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)       )
     2030          ALLOCATE ( surf_lsm_v(l)%m_sat(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)       )
     2031          ALLOCATE ( surf_lsm_v(l)%m_wilt(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)      )
     2032          ALLOCATE ( surf_lsm_v(l)%n_vg(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)        )
     2033          ALLOCATE ( surf_lsm_v(l)%rho_c_total(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) ) 
    19342034          ALLOCATE ( surf_lsm_v(l)%root_fr(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)     )
    1935           ALLOCATE ( surf_lsm_v(l)%lambda_h(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)    )
    1936           ALLOCATE ( surf_lsm_v(l)%rho_c_total(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
    1937 
    1938           surf_lsm_v(l)%lambda_h = 0.0_wp
     2035
     2036          surf_lsm_v(l)%lambda_h     = 0.0_wp
     2037         
    19392038!
    19402039!--       If required, allocate humidity-related variables for the soil model
     
    19492048!
    19502049!--    Calculate grid spacings. Temperature and moisture are defined at
    1951 !--    the edges of the soil layers (_stag), whereas gradients/fluxes are defined
    1952 !--    at the centers
    1953        dz_soil(nzb_soil) = zs(nzb_soil)
    1954 
     2050!--    the center of the soil layers, whereas gradients/fluxes are
     2051!--    defined at the edges (_layer)
     2052       DO  k = nzb_soil, nzt_soil-1
     2053          dz_soil(k) = zs(k+1) - zs(k)
     2054         
     2055          IF ( dz_soil(k) == 0.0_wp )  THEN
     2056             message_string = 'invalid soil layer configuration found ' //     &
     2057                              '(dz_soil(k) = 0.0)'
     2058             CALL message( 'lsm_read_restart_data', 'PA0140', 1, 2, 0, 6, 0 )
     2059          ENDIF
     2060         
     2061       ENDDO
     2062
     2063       dz_soil_layer(nzb_soil) = 2.0_wp * zs(nzb_soil)
     2064       zs_layer(nzb_soil)      = 2.0_wp * zs(nzb_soil)
     2065
     2066
     2067       
     2068       
    19552069       DO  k = nzb_soil+1, nzt_soil
    1956           dz_soil(k) = zs(k) - zs(k-1)
     2070          dz_soil_layer(k) = ( zs(k) - zs_layer(k-1) ) * 2.0_wp
     2071          zs_layer(k) = zs_layer(k-1) + dz_soil_layer(k)
    19572072       ENDDO
    1958        dz_soil(nzt_soil+1) = dz_soil(nzt_soil)
    1959 
    1960        DO  k = nzb_soil, nzt_soil-1
    1961           dz_soil_stag(k) = 0.5_wp * (dz_soil(k+1) + dz_soil(k))
     2073       dz_soil_layer(nzt_soil+1) = 0.5_wp * dz_soil_layer(nzt_soil)
     2074       dz_soil(nzt_soil) = zs_layer(k-1) + dz_soil_layer(k) - zs(nzt_soil)
     2075
     2076       
     2077       ddz_soil       = 1.0_wp / dz_soil
     2078       ddz_soil_layer = 1.0_wp / dz_soil_layer
     2079
     2080
     2081!
     2082!--    Set flag parameter for the prescribed surface type depending on user
     2083!--    input.
     2084       DO  m = 1, surf_lsm_h%ns
     2085
     2086          SELECT CASE ( TRIM( surface_type ) )
     2087         
     2088             CASE ( 'vegetation' )
     2089         
     2090                surf_lsm_h%vegetation_surface(m)    = .TRUE.
     2091   
     2092             CASE ( 'water' )
     2093             
     2094                surf_lsm_h%water_surface(m)      = .TRUE.
     2095
     2096             CASE ( 'pavement' )
     2097             
     2098                surf_lsm_h%pavement_surface(m)   = .TRUE.
     2099   
     2100          END SELECT
     2101 
    19622102       ENDDO
    1963        dz_soil_stag(nzt_soil) = dz_soil(nzt_soil)
    1964 
    1965        ddz_soil      = 1.0_wp / dz_soil
    1966        ddz_soil_stag = 1.0_wp / dz_soil_stag
    1967 
     2103       
     2104       DO  l = 0, 3
     2105          SELECT CASE ( TRIM( surface_type ) )
     2106         
     2107             CASE ( 'vegetation' )
     2108         
     2109                surf_lsm_v(l)%vegetation_surface = .TRUE.
     2110   
     2111             CASE ( 'water' )
     2112
     2113                surf_lsm_v(l)%water_surface   = .TRUE.
     2114
     2115             CASE ( 'pavement' )
     2116             
     2117                surf_lsm_h%pavement_surface   = .TRUE.
     2118
     2119          END SELECT   
     2120
     2121       ENDDO
    19682122!
    19692123!--    Initialize standard soil types. It is possible to overwrite each
     
    19892143
    19902144          IF ( saturation_moisture == 9999999.9_wp )  THEN
    1991              saturation_moisture = m_soil_pars(0,soil_type)           
     2145             saturation_moisture = soil_pars(4,soil_type)           
    19922146          ENDIF
    19932147
    19942148          IF ( field_capacity == 9999999.9_wp )  THEN
    1995              field_capacity = m_soil_pars(1,soil_type)           
     2149             field_capacity = soil_pars(5,soil_type)           
    19962150          ENDIF
    19972151
    19982152          IF ( wilting_point == 9999999.9_wp )  THEN
    1999              wilting_point = m_soil_pars(2,soil_type)           
     2153             wilting_point = soil_pars(6,soil_type)           
    20002154          ENDIF
    20012155
    20022156          IF ( residual_moisture == 9999999.9_wp )  THEN
    2003              residual_moisture = m_soil_pars(3,soil_type)       
     2157             residual_moisture = soil_pars(7,soil_type)       
    20042158          ENDIF
    20052159
    20062160       ENDIF   
    2007 
    2008 !
    2009 !--    Map values to the respective 2D arrays
     2161 
     2162!
     2163!--    Check whether parameters from the lookup tables are to be used
     2164       IF ( vegetation_type /= 0 )  THEN
     2165
     2166          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
     2167             min_canopy_resistance = vegetation_pars(0,vegetation_type)
     2168          ENDIF
     2169
     2170          IF ( leaf_area_index == 9999999.9_wp )  THEN
     2171             leaf_area_index = vegetation_pars(1,vegetation_type)         
     2172          ENDIF
     2173
     2174          IF ( vegetation_coverage == 9999999.9_wp )  THEN
     2175             vegetation_coverage = vegetation_pars(2,vegetation_type)     
     2176          ENDIF
     2177
     2178          IF ( canopy_resistance_coefficient == 9999999.9_wp )  THEN
     2179              canopy_resistance_coefficient= vegetation_pars(3,vegetation_type)     
     2180          ENDIF
     2181
     2182          IF ( z0_vegetation == 9999999.9_wp )  THEN
     2183             z0_vegetation  = vegetation_pars(4,vegetation_type)
     2184          ENDIF
     2185
     2186          IF ( z0h_vegetation == 9999999.9_wp )  THEN
     2187             z0h_vegetation = vegetation_pars(5,vegetation_type)
     2188          ENDIF   
     2189         
     2190!
     2191!--       The heat conductivity between canopy and first soil layer must
     2192!--       depend on the depth of the soil layer. The value in the lookup
     2193!--       table refers to the first soil temperature at zs = -0.005 m. Lambda
     2194!--       is thus interpolated to the actual soil layer depth used.
     2195          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
     2196             lambda_surface_stable = vegetation_pars(6,vegetation_type)        &
     2197                                     * 0.005_wp * ddz_soil_layer(nzb_soil)     &
     2198                                     * 2.0_wp   
     2199          ENDIF
     2200
     2201          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
     2202             lambda_surface_unstable = vegetation_pars(7,vegetation_type)      &
     2203                                       * 0.005_wp  * ddz_soil_layer(nzb_soil)  &
     2204                                       * 2.0_wp             
     2205          ENDIF
     2206
     2207          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
     2208             f_shortwave_incoming = vegetation_pars(8,vegetation_type)       
     2209          ENDIF
     2210
     2211          IF ( c_surface == 9999999.9_wp )  THEN
     2212             c_surface = vegetation_pars(9,vegetation_type)       
     2213          ENDIF
     2214             
     2215       ENDIF
     2216
     2217       IF ( water_type /= 0 )  THEN
     2218
     2219          IF ( water_temperature == 9999999.9_wp )  THEN
     2220             water_temperature = water_pars(0,water_type)       
     2221          ENDIF
     2222
     2223          IF ( z0_water == 9999999.9_wp )  THEN
     2224             z0_water = water_pars(1,water_type)       
     2225          ENDIF       
     2226
     2227          IF ( z0h_water == 9999999.9_wp )  THEN
     2228             z0h_water = water_pars(2,water_type)       
     2229          ENDIF 
     2230
     2231       ENDIF
     2232       
     2233       IF ( pavement_type /= 0 )  THEN 
     2234
     2235          IF ( pavement_depth == 9999999.9_wp )  THEN
     2236             pavement_depth = pavement_pars(0,pavement_type)
     2237          ELSE
     2238             pavement_depth = MAX( zs(nzb_soil), pavement_depth )
     2239          ENDIF
     2240
     2241          IF ( z0_pavement == 9999999.9_wp )  THEN
     2242             z0_pavement  = pavement_pars(1,pavement_type)
     2243          ENDIF
     2244
     2245          IF ( z0h_pavement == 9999999.9_wp )  THEN
     2246             z0h_pavement = pavement_pars(2,pavement_type)
     2247          ENDIF
     2248
     2249          IF ( pavement_heat_conduct == 9999999.9_wp )  THEN
     2250             pavement_heat_conduct = pavement_pars(3,pavement_type)
     2251          ENDIF
     2252
     2253          IF ( pavement_heat_capacity == 9999999.9_wp )  THEN
     2254             pavement_heat_capacity = pavement_pars(4,pavement_type)
     2255          ENDIF   
     2256   
     2257       ENDIF
     2258!
     2259!--    Map values to the respective 2D/3D arrays
    20102260!--    Horizontal surfaces
    20112261       surf_lsm_h%alpha_vg      = alpha_vangenuchten
     
    20322282       ENDDO
    20332283
    2034 
    2035 
    20362284!
    20372285!--    Initial run actions
     
    20392287!
    20402288!--       First, for horizontal surfaces
    2041           t_soil_h%var_2d    = 0.0_wp
    2042           m_soil_h%var_2d    = 0.0_wp
    2043           m_liq_eb_h%var_1d  = 0.0_wp
    2044 
    2045 !
    2046 !--       Map user settings of T and q for each soil layer
    2047 !--       (make sure that the soil moisture does not drop below the permanent
    2048 !--       wilting point) -> problems with devision by zero)
    2049           DO  k = nzb_soil, nzt_soil
    2050              t_soil_h%var_2d(k,:)      = soil_temperature(k)
    2051              m_soil_h%var_2d(k,:)      = MAX(soil_moisture(k),surf_lsm_h%m_wilt(:))
    2052              soil_moisture(k) = MAX(soil_moisture(k),wilting_point)
    2053           ENDDO
    2054           t_soil_h%var_2d(nzt_soil+1,:) = soil_temperature(nzt_soil+1)
    2055 
    2056 !
    2057 !--       Calculate surface temperature
    2058           t_surface_h%var_1d(:)    = pt_surface * exn
    2059           surf_lsm_h%pt_surface(:) = pt_surface
    2060 
    20612289!
    20622290!--       Set artifical values for ts and us so that r_a has its initial value
    20632291!--       for the first time step. Only for interior core domain, not for ghost points
    20642292          DO  m = 1, surf_lsm_h%ns
     2293
     2294             t_soil_h%var_2d(:,m)    = 0.0_wp
     2295             m_soil_h%var_2d(:,m)    = 0.0_wp
     2296             m_liq_h%var_1d(m)       = 0.0_wp
     2297
     2298!--          Map user settings of T and q for each soil layer
     2299!--          (make sure that the soil moisture does not drop below the permanent
     2300!--          wilting point) -> problems with devision by zero)
     2301             IF ( surf_lsm_h%water_surface(m) )  THEN
     2302
     2303                surf_lsm_h%z0(m)          = z0_water
     2304                surf_lsm_h%z0h(m)         = z0h_water
     2305                surf_lsm_h%z0q(m)         = z0h_water
     2306                t_soil_h%var_2d(:,m)      = water_temperature
     2307                surf_lsm_h%lambda_surface_s(m) = 1.0E10_wp
     2308                surf_lsm_h%lambda_surface_u(m) = 1.0E10_wp               
     2309                surf_lsm_h%c_surface(m)        = 0.0_wp
     2310               
     2311             ELSEIF ( surf_lsm_h%pavement_surface(m) )  THEN
     2312
     2313                surf_lsm_h%z0(m)          = z0_pavement
     2314                surf_lsm_h%z0h(m)         = z0h_pavement
     2315                surf_lsm_h%z0q(m)         = z0h_pavement   
     2316                DO  k = nzb_soil, nzt_soil 
     2317                   t_soil_h%var_2d(k,m)   = soil_temperature(k)
     2318                   m_soil_h%var_2d(k,m)   = MAX(soil_moisture(k),surf_lsm_h%m_wilt(k,m))
     2319                   soil_moisture(k)       = MAX(soil_moisture(k),wilting_point)           
     2320                ENDDO
     2321                t_soil_h%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1)
     2322                surf_lsm_h%lambda_surface_s(m) = pavement_heat_conduct         &
     2323                                                 * ddz_soil_layer(nzb_soil)    &
     2324                                                 * 2.0_wp
     2325                surf_lsm_h%lambda_surface_u(m) = surf_lsm_h%lambda_surface_s(m)                                         
     2326                surf_lsm_h%c_surface(m)        = pavement_heat_capacity        &
     2327                                                 * dz_soil_layer(nzb_soil)     &
     2328                                                 * 0.25_wp
     2329                surf_lsm_h%lambda_h_def(m)    = pavement_heat_conduct                                     
     2330                surf_lsm_h%rho_c_def(m)       = pavement_heat_capacity                                               
     2331               
     2332             ELSEIF ( surf_lsm_h%vegetation_surface(m) )  THEN
     2333
     2334                surf_lsm_h%z0(m)          = z0_vegetation
     2335                surf_lsm_h%z0h(m)         = z0h_vegetation
     2336                surf_lsm_h%z0q(m)         = z0h_vegetation   
     2337                DO  k = nzb_soil, nzt_soil                       
     2338                   t_soil_h%var_2d(k,m)   = soil_temperature(k)
     2339                   m_soil_h%var_2d(k,m)   = MAX(soil_moisture(k),surf_lsm_h%m_wilt(k,m))
     2340                   soil_moisture(k)       = MAX(soil_moisture(k),wilting_point)         
     2341                ENDDO
     2342                t_soil_h%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1)
     2343                surf_lsm_h%lambda_surface_s(m) = lambda_surface_stable
     2344                surf_lsm_h%lambda_surface_u(m) = lambda_surface_unstable
     2345                surf_lsm_h%c_surface(m)        = c_surface
     2346
     2347             ENDIF
     2348 
    20652349             i   = surf_lsm_h%i(m)           
    20662350             j   = surf_lsm_h%j(m)
    20672351             k   = surf_lsm_h%k(m)
    2068 
     2352!
     2353!--          Calculate surface temperature. In case of bare soil, the surface
     2354!--          temperature must be reset to the soil temperature in the first soil
     2355!--          layer
     2356             IF ( surf_lsm_h%lambda_surface_s(m) == 0.0_wp )  THEN
     2357                t_surface_h%var_1d(:)    = t_soil_h%var_2d(nzb_soil,m)
     2358                surf_lsm_h%pt_surface(:) = t_soil_h%var_2d(nzb_soil,m) / exn
     2359             ELSE
     2360                t_surface_h%var_1d(:)    = pt_surface * exn
     2361                surf_lsm_h%pt_surface(:) = pt_surface 
     2362             ENDIF
     2363             
    20692364             IF ( cloud_physics )  THEN
    20702365                pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
     
    20792374             surf_lsm_h%us(m)   = 0.1_wp
    20802375             surf_lsm_h%ts(m)   = ( pt1 - surf_lsm_h%pt_surface(m) ) / surf_lsm_h%r_a(m)
    2081              surf_lsm_h%shf(m)  = - surf_lsm_h%us(m) * surf_lsm_h%ts(m)
    2082           ENDDO
    2083 !
    2084 !--       Vertical surfaces
    2085           DO  l = 0, 3
    2086 
    2087              t_soil_v(l)%var_2d    = 0.0_wp
    2088              m_soil_v(l)%var_2d    = 0.0_wp
    2089              m_liq_eb_v(l)%var_1d  = 0.0_wp
    2090 
    2091 !
    2092 !--          Map user settings of T and q for each soil layer
    2093 !--          (make sure that the soil moisture does not drop below the permanent
    2094 !--          wilting point) -> problems with devision by zero)
    2095              DO  k = nzb_soil, nzt_soil
    2096                 t_soil_v(l)%var_2d(k,:)      = soil_temperature(k)
    2097                 m_soil_v(l)%var_2d(k,:)      = MAX(soil_moisture(k),surf_lsm_v(l)%m_wilt(:))
    2098                 soil_moisture(k) = MAX(soil_moisture(k),wilting_point)
    2099              ENDDO
    2100              t_soil_v(l)%var_2d(nzt_soil+1,:) = soil_temperature(nzt_soil+1)
    2101 
    2102 !
    2103 !--          Calculate surface temperature
    2104              t_surface_v(l)%var_1d(:)   = pt_surface * exn
    2105              surf_lsm_h%pt_surface(:)   = pt_surface
    2106 
    2107 !
    2108 !--          Set artifical values for ts and us so that r_a has its initial value
    2109 !--          for the first time step. Only for interior core domain, not for ghost points
     2376             surf_lsm_h%shf(m)  = - surf_lsm_h%us(m) * surf_lsm_h%ts(m) * rho_surface
     2377
     2378         ENDDO
     2379
     2380!
     2381!--      Vertical surfaces
     2382         DO  l = 0, 3
    21102383             DO  m = 1, surf_lsm_v(l)%ns
     2384             
     2385                t_soil_v(l)%var_2d    = 0.0_wp
     2386                m_soil_v(l)%var_2d    = 0.0_wp
     2387                m_liq_v(l)%var_1d  = 0.0_wp
     2388
     2389             
     2390!--             Map user settings of T and q for each soil layer
     2391!--             (make sure that the soil moisture does not drop below the permanent
     2392!--             wilting point) -> problems with devision by zero)
     2393                IF ( surf_lsm_v(l)%water_surface(m) )  THEN
     2394                   surf_lsm_v(l)%z0(m)          = z0_water
     2395                   surf_lsm_v(l)%z0h(m)         = z0h_water
     2396                   surf_lsm_v(l)%z0q(m)         = z0h_water                   
     2397                   t_soil_v(l)%var_2d(:,m)      = water_temperature
     2398                   surf_lsm_v(l)%lambda_surface_s(m) = 1.0E10_wp
     2399                   surf_lsm_v(l)%lambda_surface_u(m) = 1.0E10_wp             
     2400                   surf_lsm_v(l)%c_surface(m)        = 0.0_wp
     2401                   
     2402                ELSEIF ( surf_lsm_v(l)%pavement_surface(m) )  THEN
     2403                   surf_lsm_v(l)%z0(m)          = z0_pavement
     2404                   surf_lsm_v(l)%z0h(m)         = z0h_pavement
     2405                   surf_lsm_v(l)%z0q(m)         = z0h_pavement 
     2406                   DO  k = nzb_soil, nzt_soil 
     2407                      t_soil_v(l)%var_2d(k,m)      = soil_temperature(k)
     2408                      m_soil_v(l)%var_2d(k,m)   = MAX(soil_moisture(k),surf_lsm_v(l)%m_wilt(k,m))
     2409                      soil_moisture(k)          = MAX(soil_moisture(k),wilting_point)   
     2410                   ENDDO
     2411                   t_soil_v(l)%var_2d(nzt_soil+1,m)  = soil_temperature(nzt_soil+1)
     2412                   surf_lsm_v(l)%lambda_surface_s(m) = pavement_heat_conduct      &
     2413                                                       * ddz_soil_layer(nzb_soil) &
     2414                                                       * 2.0_wp
     2415                   surf_lsm_v(l)%lambda_surface_u(m) = surf_lsm_v(l)%lambda_surface_s(m)
     2416                   surf_lsm_v(l)%c_surface(m) = pavement_heat_capacity         &
     2417                                                * dz_soil_layer(nzb_soil)      &
     2418                                                * 0.25_wp
     2419                   surf_lsm_v(l)%lambda_h_def(m)   = pavement_heat_conduct                                     
     2420                   surf_lsm_v(l)%rho_c_def(m)      = pavement_heat_capacity                           
     2421                   
     2422                ELSEIF ( surf_lsm_v(l)%vegetation_surface(m) )  THEN
     2423
     2424                   surf_lsm_v(l)%z0(m)          = z0_vegetation
     2425                   surf_lsm_v(l)%z0h(m)         = z0h_vegetation
     2426                   surf_lsm_v(l)%z0q(m)         = z0h_vegetation   
     2427                   DO  k = nzb_soil, nzt_soil                       
     2428                      t_soil_v(l)%var_2d(k,m)   = soil_temperature(k)
     2429                      m_soil_v(l)%var_2d(k,m)   = MAX(soil_moisture(k),surf_lsm_v(l)%m_wilt(k,m))
     2430                      soil_moisture(k)          = MAX(soil_moisture(k),wilting_point)         
     2431                   ENDDO
     2432                   t_soil_v(l)%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1)
     2433                   surf_lsm_v(l)%lambda_surface_s(m) = lambda_surface_stable
     2434                   surf_lsm_v(l)%lambda_surface_u(m) = lambda_surface_unstable
     2435                   surf_lsm_v(l)%c_surface(m)        = c_surface
     2436                ENDIF       
     2437         
     2438!
     2439!--             Calculate surface temperature. In case of bare soil, the surface
     2440!--             temperature must be reset to the soil temperature in the first soil
     2441!--             layer
     2442                IF ( surf_lsm_v(l)%lambda_surface_s(m) == 0.0_wp )  THEN
     2443                   t_surface_v(l)%var_1d(:)      = t_soil_v(l)%var_2d(nzb_soil,m)
     2444                   surf_lsm_v(l)%pt_surface(:)   = t_soil_v(l)%var_2d(nzb_soil,m) / exn
     2445                ELSE
     2446                   t_surface_v(l)%var_1d(:)      = pt_surface * exn
     2447                   surf_lsm_v(l)%pt_surface(:)   = pt_surface               
     2448                ENDIF
     2449
     2450!
     2451!--             Set artifical values for ts and us so that r_a has its initial value
     2452!--             for the first time step. Only for interior core domain, not for ghost points
     2453
    21112454                i   = surf_lsm_v(l)%i(m)           
    21122455                j   = surf_lsm_v(l)%j(m)
     
    21212464!
    21222465!--             Assure that r_a cannot be zero at model start
    2123                 IF ( pt1 == surf_lsm_h%pt_surface(m) )  pt1 = pt1 + 1.0E-10_wp
     2466                IF ( pt1 == surf_lsm_v(l)%pt_surface(m) )  pt1 = pt1 + 1.0E-10_wp
    21242467
    21252468                surf_lsm_v(l)%us(m)   = 0.1_wp
    2126                 surf_lsm_v(l)%ts(m)   = ( pt1 - surf_lsm_h%pt_surface(m) ) / surf_lsm_v(l)%r_a(m)
    2127                 surf_lsm_v(l)%shf(m)  = - surf_lsm_v(l)%us(m) * surf_lsm_v(l)%ts(m)
     2469                surf_lsm_v(l)%ts(m)   = ( pt1 - surf_lsm_v(l)%pt_surface(m) ) / surf_lsm_v(l)%r_a(m)
     2470                surf_lsm_v(l)%shf(m)  = - surf_lsm_v(l)%us(m) * surf_lsm_v(l)%ts(m) * rho_surface
    21282471             ENDDO
    21292472          ENDDO
     
    21912534       ENDDO
    21922535
    2193        IF ( veg_type /= 0 )  THEN
    2194           IF ( min_canopy_resistance == 9999999.9_wp )  THEN
    2195              min_canopy_resistance = veg_pars(0,veg_type)
    2196           ENDIF
    2197           IF ( leaf_area_index == 9999999.9_wp )  THEN
    2198              leaf_area_index = veg_pars(1,veg_type)         
    2199           ENDIF
    2200           IF ( vegetation_coverage == 9999999.9_wp )  THEN
    2201              vegetation_coverage = veg_pars(2,veg_type)     
    2202           ENDIF
    2203           IF ( canopy_resistance_coefficient == 9999999.9_wp )  THEN
    2204               canopy_resistance_coefficient= veg_pars(3,veg_type)     
    2205           ENDIF
    2206           IF ( lambda_surface_stable == 9999999.9_wp )  THEN
    2207              lambda_surface_stable = surface_pars(0,veg_type)         
    2208           ENDIF
    2209           IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
    2210              lambda_surface_unstable = surface_pars(1,veg_type)       
    2211           ENDIF
    2212           IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
    2213              f_shortwave_incoming = surface_pars(2,veg_type)       
    2214           ENDIF
    2215           IF ( z0_eb == 9999999.9_wp )  THEN
    2216              roughness_length = roughness_par(0,veg_type)
    2217              z0_eb            = roughness_par(0,veg_type)
    2218           ENDIF
    2219           IF ( z0h_eb == 9999999.9_wp )  THEN
    2220              z0h_eb = roughness_par(1,veg_type)
    2221           ENDIF
    2222           IF ( z0q_eb == 9999999.9_wp )  THEN
    2223              z0q_eb = roughness_par(2,veg_type)
    2224           ENDIF
    2225           z0h_factor = z0h_eb / ( z0_eb + 1.0E-20_wp )
    2226 
    2227           IF ( ALL( root_fraction == 9999999.9_wp ) )  THEN
    2228              DO  m = 1, surf_lsm_h%ns
    2229                 i   = surf_lsm_h%i(m)           
    2230                 j   = surf_lsm_h%j(m)
    2231 
    2232 
    2233                 DO  k = nzb_soil, nzt_soil
    2234                    surf_lsm_h%root_fr(k,m) = root_distribution(k,veg_type)
    2235                    root_fraction(k)        = root_distribution(k,veg_type)
     2536
     2537
     2538       IF ( ALL( root_fraction == 9999999.9_wp ) )  THEN
     2539          DO  m = 1, surf_lsm_h%ns
     2540             i   = surf_lsm_h%i(m)           
     2541             j   = surf_lsm_h%j(m)
     2542
     2543
     2544             DO  k = nzb_soil, nzt_soil
     2545                surf_lsm_h%root_fr(k,m) = root_distribution(k,vegetation_type)
     2546                root_fraction(k)        = root_distribution(k,vegetation_type)
     2547             ENDDO
     2548          ENDDO
     2549          DO  l = 0, 3
     2550             DO  m = 1, surf_lsm_v(l)%ns
     2551                i   = surf_lsm_v(l)%i(m)           
     2552                j   = surf_lsm_v(l)%j(m)
     2553
     2554                DO  k = nzb_soil, nzt_soil
     2555                   surf_lsm_v(l)%root_fr(k,m) = root_distribution(k,vegetation_type)
     2556                   root_fraction(k)           = root_distribution(k,vegetation_type)
    22362557                ENDDO
    22372558             ENDDO
    2238              DO  l = 0, 3
    2239                 DO  m = 1, surf_lsm_v(l)%ns
    2240                    i   = surf_lsm_v(l)%i(m)           
    2241                    j   = surf_lsm_v(l)%j(m)
    2242 
    2243                    DO  k = nzb_soil, nzt_soil
    2244                       surf_lsm_v(l)%root_fr(k,m) = root_distribution(k,veg_type)
    2245                       root_fraction(k)    = root_distribution(k,veg_type)
    2246                    ENDDO
    2247                 ENDDO
    2248              ENDDO
    2249           ENDIF
    2250 
    2251        ELSE
    2252 
    2253           IF ( z0_eb == 9999999.9_wp )  THEN
    2254              z0_eb = roughness_length
    2255           ENDIF
    2256           IF ( z0h_eb == 9999999.9_wp )  THEN
    2257              z0h_eb = z0_eb * z0h_factor
    2258           ENDIF
    2259           IF ( z0q_eb == 9999999.9_wp )  THEN
    2260              z0q_eb = z0_eb * z0h_factor
    2261           ENDIF
    2262 
    2263        ENDIF
    2264 
    2265 !
    2266 !--    For surfaces covered with pavement, set depth of the pavement (with dry
    2267 !--    soil below). The depth must be greater than the first soil layer depth
    2268        IF ( veg_type == 20 )  THEN
    2269           IF ( pave_depth == 9999999.9_wp )  THEN
    2270              pave_depth = zs(nzb_soil) 
    2271           ELSE
    2272              pave_depth = MAX( zs(nzb_soil), pave_depth )
    2273           ENDIF
    2274        ENDIF
    2275 
    2276 !
    2277 !--    Map vegetation and soil types to 2D array to allow for heterogeneous
    2278 !--    surfaces via user interface see below
    2279 !--    First for horizontal surfaces
    2280        surf_lsm_h%veg_type_2d  = veg_type
    2281        surf_lsm_h%soil_type_2d = soil_type
     2559           ENDDO
     2560       ENDIF
    22822561
    22832562!
     
    22872566       surf_lsm_h%c_veg             = vegetation_coverage
    22882567       surf_lsm_h%g_d               = canopy_resistance_coefficient
    2289        surf_lsm_h%lambda_surface_s  = lambda_surface_stable
    2290        surf_lsm_h%lambda_surface_u  = lambda_surface_unstable
    22912568       surf_lsm_h%f_sw_in           = f_shortwave_incoming
    2292        surf_lsm_h%z0                = z0_eb
    2293        surf_lsm_h%z0h               = z0h_eb
    2294        surf_lsm_h%z0q               = z0q_eb
    22952569
    22962570!--    Vertical surfaces
    22972571       DO  l = 0, 3
    2298           surf_lsm_v(l)%veg_type_2d  = veg_type
    2299           surf_lsm_v(l)%soil_type_2d = soil_type
    23002572
    23012573!
     
    23052577          surf_lsm_v(l)%c_veg             = vegetation_coverage
    23062578          surf_lsm_v(l)%g_d               = canopy_resistance_coefficient
    2307           surf_lsm_v(l)%lambda_surface_s  = lambda_surface_stable
    2308           surf_lsm_v(l)%lambda_surface_u  = lambda_surface_unstable
    23092579          surf_lsm_v(l)%f_sw_in           = f_shortwave_incoming
    2310           surf_lsm_v(l)%z0                = z0_eb
    2311           surf_lsm_v(l)%z0h               = z0h_eb
    2312           surf_lsm_v(l)%z0q               = z0q_eb
    23132580       ENDDO
    2314 
     2581 
    23152582!
    23162583!--    Possibly do user-defined actions (e.g. define heterogeneous land surface)
    23172584       CALL user_init_land_surface
    23182585
    2319 !
    2320 !--    Set flag parameter if vegetation type was set to a water surface. Also
    2321 !--    set temperature to a constant value in all "soil" layers.
    2322 !--    For now, do not set water surface for vertical surfaces
    2323        DO  m = 1, surf_lsm_h%ns
    2324           i   = surf_lsm_h%i(m)           
    2325           j   = surf_lsm_h%j(m)
    2326 
    2327           IF ( surf_lsm_h%veg_type_2d(m) == 14  .OR.                           &
    2328                surf_lsm_h%veg_type_2d(m) == 15 )  THEN
    2329              surf_lsm_h%water_surface(m) = .TRUE.
    2330              t_soil_h%var_2d(:,m) = t_surface_h%var_1d(m)
    2331           ELSEIF ( surf_lsm_h%veg_type_2d(m) == 20 )  THEN
    2332              surf_lsm_h%pave_surface(m) = .TRUE.
    2333              m_soil_h%var_2d(:,m)            = 0.0_wp
    2334           ENDIF
    2335 
    2336        ENDDO
    23372586
    23382587!
     
    23432592       t_soil_h_p    = t_soil_h
    23442593       m_soil_h_p    = m_soil_h
    2345        m_liq_eb_h_p  = m_liq_eb_h
     2594       m_liq_h_p  = m_liq_h
    23462595       t_surface_h_p = t_surface_h
    23472596
    23482597       t_soil_v_p    = t_soil_v
    23492598       m_soil_v_p    = m_soil_v
    2350        m_liq_eb_v_p  = m_liq_eb_v
     2599       m_liq_v_p  = m_liq_v
    23512600       t_surface_v_p = t_surface_v
    23522601
     
    23772626!
    23782627!--    Allocate global 1D arrays
    2379        ALLOCATE ( ddz_soil(nzb_soil:nzt_soil+1) )
    2380        ALLOCATE ( ddz_soil_stag(nzb_soil:nzt_soil) )
    2381        ALLOCATE ( dz_soil(nzb_soil:nzt_soil+1) )
    2382        ALLOCATE ( dz_soil_stag(nzb_soil:nzt_soil) )
     2628       ALLOCATE ( ddz_soil(nzb_soil:nzt_soil) )
     2629       ALLOCATE ( ddz_soil_layer(nzb_soil:nzt_soil+1) )
     2630       ALLOCATE ( dz_soil(nzb_soil:nzt_soil) )
     2631       ALLOCATE ( dz_soil_layer(nzb_soil:nzt_soil+1) )
    23832632       ALLOCATE ( root_extr(nzb_soil:nzt_soil) )
    23842633 
     
    23932642!
    23942643!--    Horizontal surfaces
    2395        ALLOCATE ( m_liq_eb_h%var_1d(1:surf_lsm_h%ns)                     )
    2396        ALLOCATE ( m_liq_eb_h_p%var_1d(1:surf_lsm_h%ns)                   )
     2644       ALLOCATE ( m_liq_h%var_1d(1:surf_lsm_h%ns)                     )
     2645       ALLOCATE ( m_liq_h_p%var_1d(1:surf_lsm_h%ns)                   )
    23972646       ALLOCATE ( t_surface_h%var_1d(1:surf_lsm_h%ns)                    )
    23982647       ALLOCATE ( t_surface_h_p%var_1d(1:surf_lsm_h%ns)                  )
     
    24012650       ALLOCATE ( t_soil_h%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns)   )
    24022651       ALLOCATE ( t_soil_h_p%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) )
     2652
    24032653!
    24042654!--    Vertical surfaces
    24052655       DO  l = 0, 3
    2406           ALLOCATE ( m_liq_eb_v(l)%var_1d(1:surf_lsm_v(l)%ns)                     )
    2407           ALLOCATE ( m_liq_eb_v(l)_p%var_1d(1:surf_lsm_v(l)%ns)                   )
     2656          ALLOCATE ( m_liq_v(l)%var_1d(1:surf_lsm_v(l)%ns)                     )
     2657          ALLOCATE ( m_liq_v(l)_p%var_1d(1:surf_lsm_v(l)%ns)                   )
    24082658          ALLOCATE ( t_surface_v(l)%var_1d(1:surf_lsm_v(l)%ns)                    )
    24092659          ALLOCATE ( t_surface_v(l)_p%var_1d(1:surf_lsm_v(l)%ns)                  )
     
    24162666!
    24172667!--    Horizontal surfaces
    2418        ALLOCATE ( m_liq_eb_h_1%var_1d(1:surf_lsm_h%ns)                   )
    2419        ALLOCATE ( m_liq_eb_h_2%var_1d(1:surf_lsm_h%ns)                   )
     2668       ALLOCATE ( m_liq_h_1%var_1d(1:surf_lsm_h%ns)                   )
     2669       ALLOCATE ( m_liq_h_2%var_1d(1:surf_lsm_h%ns)                   )
    24202670       ALLOCATE ( t_surface_h_1%var_1d(1:surf_lsm_h%ns)                  )
    24212671       ALLOCATE ( t_surface_h_2%var_1d(1:surf_lsm_h%ns)                  )
     
    24272677!--    Vertical surfaces
    24282678       DO  l = 0, 3
    2429           ALLOCATE ( m_liq_eb_v_1(l)%var_1d(1:surf_lsm_v(l)%ns)                   )
    2430           ALLOCATE ( m_liq_eb_v_2(l)%var_1d(1:surf_lsm_v(l)%ns)                   )
     2679          ALLOCATE ( m_liq_v_1(l)%var_1d(1:surf_lsm_v(l)%ns)                   )
     2680          ALLOCATE ( m_liq_v_2(l)%var_1d(1:surf_lsm_v(l)%ns)                   )
    24312681          ALLOCATE ( t_surface_v_1(l)%var_1d(1:surf_lsm_v(l)%ns)                  )
    24322682          ALLOCATE ( t_surface_v_2(l)%var_1d(1:surf_lsm_v(l)%ns)                  )
     
    24412691!--    Allocate intermediate timestep arrays
    24422692!--    Horizontal surfaces
    2443        ALLOCATE ( tm_liq_eb_h_m%var_1d(1:surf_lsm_h%ns)                  )
     2693       ALLOCATE ( tm_liq_h_m%var_1d(1:surf_lsm_h%ns)                  )
    24442694       ALLOCATE ( tt_surface_h_m%var_1d(1:surf_lsm_h%ns)                 )
    24452695       ALLOCATE ( tm_soil_h_m%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns)  )
     
    24482698!--    Horizontal surfaces
    24492699       DO  l = 0, 3
    2450           ALLOCATE ( tm_liq_eb_v_m(l)%var_1d(1:surf_lsm_v(l)%ns)                  )
     2700          ALLOCATE ( tm_liq_v_m(l)%var_1d(1:surf_lsm_v(l)%ns)                  )
    24512701          ALLOCATE ( tt_surface_v_m(l)%var_1d(1:surf_lsm_v(l)%ns)                 )
    24522702          ALLOCATE ( tm_soil_v_m(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)  )
     
    24632713!--    Allocate 2D vegetation model arrays
    24642714!--    Horizontal surfaces
    2465        ALLOCATE ( surf_lsm_h%alpha_vg(1:surf_lsm_h%ns)         )
    24662715       ALLOCATE ( surf_lsm_h%building_surface(1:surf_lsm_h%ns) )
    24672716       ALLOCATE ( surf_lsm_h%c_liq(1:surf_lsm_h%ns)            )
     2717       ALLOCATE ( surf_lsm_h%c_surface(1:surf_lsm_h%ns)        )
    24682718       ALLOCATE ( surf_lsm_h%c_veg(1:surf_lsm_h%ns)            )
    24692719       ALLOCATE ( surf_lsm_h%f_sw_in(1:surf_lsm_h%ns)          )
    2470        ALLOCATE ( surf_lsm_h%ghf_eb(1:surf_lsm_h%ns)           )
    2471        ALLOCATE ( surf_lsm_h%gamma_w_sat(1:surf_lsm_h%ns)      )
     2720       ALLOCATE ( surf_lsm_h%ghf(1:surf_lsm_h%ns)              )
    24722721       ALLOCATE ( surf_lsm_h%g_d(1:surf_lsm_h%ns)              )
    24732722       ALLOCATE ( surf_lsm_h%lai(1:surf_lsm_h%ns)              )
    2474        ALLOCATE ( surf_lsm_h%l_vg(1:surf_lsm_h%ns)             )
     2723       ALLOCATE ( surf_lsm_h%lambda_h_def(1:surf_lsm_h%ns)     )
    24752724       ALLOCATE ( surf_lsm_h%lambda_surface_u(1:surf_lsm_h%ns) )
    24762725       ALLOCATE ( surf_lsm_h%lambda_surface_s(1:surf_lsm_h%ns) )
    2477        ALLOCATE ( surf_lsm_h%m_fc(1:surf_lsm_h%ns)             )
    2478        ALLOCATE ( surf_lsm_h%m_res(1:surf_lsm_h%ns)            )
    2479        ALLOCATE ( surf_lsm_h%m_sat(1:surf_lsm_h%ns)            )
    2480        ALLOCATE ( surf_lsm_h%m_wilt(1:surf_lsm_h%ns)           )
    2481        ALLOCATE ( surf_lsm_h%n_vg(1:surf_lsm_h%ns)             )
    2482        ALLOCATE ( surf_lsm_h%pave_surface(1:surf_lsm_h%ns)     )
    2483        ALLOCATE ( surf_lsm_h%qsws_eb(1:surf_lsm_h%ns)          )
    2484        ALLOCATE ( surf_lsm_h%qsws_soil_eb(1:surf_lsm_h%ns)     )
    2485        ALLOCATE ( surf_lsm_h%qsws_liq_eb(1:surf_lsm_h%ns)      )
    2486        ALLOCATE ( surf_lsm_h%qsws_veg_eb(1:surf_lsm_h%ns)      )
     2726       ALLOCATE ( surf_lsm_h%vegetation_surface(1:surf_lsm_h%ns)  )
     2727       ALLOCATE ( surf_lsm_h%pavement_surface(1:surf_lsm_h%ns) )
     2728       ALLOCATE ( surf_lsm_h%qsws_soil(1:surf_lsm_h%ns)        )
     2729       ALLOCATE ( surf_lsm_h%qsws_liq(1:surf_lsm_h%ns)         )
     2730       ALLOCATE ( surf_lsm_h%qsws_veg(1:surf_lsm_h%ns)         )
    24872731       ALLOCATE ( surf_lsm_h%rad_net_l(1:surf_lsm_h%ns)        )
     2732       ALLOCATE ( surf_lsm_h%rho_c_def(1:surf_lsm_h%ns)        ) 
    24882733       ALLOCATE ( surf_lsm_h%r_a(1:surf_lsm_h%ns)              )
    24892734       ALLOCATE ( surf_lsm_h%r_canopy(1:surf_lsm_h%ns)         )
     
    24922737       ALLOCATE ( surf_lsm_h%r_s(1:surf_lsm_h%ns)              )
    24932738       ALLOCATE ( surf_lsm_h%r_canopy_min(1:surf_lsm_h%ns)     )
    2494        ALLOCATE ( surf_lsm_h%shf_eb(1:surf_lsm_h%ns)           )
    2495        ALLOCATE ( surf_lsm_h%soil_type_2d(1:surf_lsm_h%ns)     )
    2496        ALLOCATE ( surf_lsm_h%veg_type_2d(1:surf_lsm_h%ns)      )
    24972739       ALLOCATE ( surf_lsm_h%water_surface(1:surf_lsm_h%ns)    )
    24982740
    2499        surf_lsm_h%water_surface = .FALSE.
    2500        surf_lsm_h%pave_surface  = .FALSE.
     2741       surf_lsm_h%water_surface     = .FALSE.
     2742       surf_lsm_h%pavement_surface  = .FALSE.
     2743       surf_lsm_h%vegetation_surface   = .FALSE.
    25012744!
    25022745!--    Vertical surfaces
    25032746       DO  l = 0, 3
    2504           ALLOCATE ( surf_lsm_v(l)%alpha_vg(1:surf_lsm_v(l)%ns)         )
    25052747          ALLOCATE ( surf_lsm_v(l)%building_surface(1:surf_lsm_v(l)%ns) )
    25062748          ALLOCATE ( surf_lsm_v(l)%c_liq(1:surf_lsm_v(l)%ns)            )
     2749          ALLOCATE ( surf_lsm_v(l)%c_surface(1:surf_lsm_v(l)%ns)        )
    25072750          ALLOCATE ( surf_lsm_v(l)%c_veg(1:surf_lsm_v(l)%ns)            )
    25082751          ALLOCATE ( surf_lsm_v(l)%f_sw_in(1:surf_lsm_v(l)%ns)          )
    2509           ALLOCATE ( surf_lsm_v(l)%ghf_eb(1:surf_lsm_v(l)%ns)           )
    2510           ALLOCATE ( surf_lsm_v(l)%gamma_w_sat(1:surf_lsm_v(l)%ns)      )
     2752          ALLOCATE ( surf_lsm_v(l)%ghf(1:surf_lsm_v(l)%ns)               )
    25112753          ALLOCATE ( surf_lsm_v(l)%g_d(1:surf_lsm_v(l)%ns)              )
    25122754          ALLOCATE ( surf_lsm_v(l)%lai(1:surf_lsm_v(l)%ns)              )
    2513           ALLOCATE ( surf_lsm_v(l)%l_vg(1:surf_lsm_v(l)%ns)             )
     2755          ALLOCATE ( surf_lsm_v(l)%lambda_h_def(1:surf_lsm_v(l)%ns)     )
    25142756          ALLOCATE ( surf_lsm_v(l)%lambda_surface_u(1:surf_lsm_v(l)%ns) )
    25152757          ALLOCATE ( surf_lsm_v(l)%lambda_surface_s(1:surf_lsm_v(l)%ns) )
    2516           ALLOCATE ( surf_lsm_v(l)%m_fc(1:surf_lsm_v(l)%ns)             )
    2517           ALLOCATE ( surf_lsm_v(l)%m_res(1:surf_lsm_v(l)%ns)            )
    2518           ALLOCATE ( surf_lsm_v(l)%m_sat(1:surf_lsm_v(l)%ns)            )
    2519           ALLOCATE ( surf_lsm_v(l)%m_wilt(1:surf_lsm_v(l)%ns)           )
    2520           ALLOCATE ( surf_lsm_v(l)%n_vg(1:surf_lsm_v(l)%ns)             )
    2521           ALLOCATE ( surf_lsm_v(l)%pave_surface(1:surf_lsm_v(l)%ns)     )
    2522           ALLOCATE ( surf_lsm_v(l)%qsws_eb(1:surf_lsm_v(l)%ns)          )
    2523           ALLOCATE ( surf_lsm_v(l)%qsws_soil_eb(1:surf_lsm_v(l)%ns)     )
    2524           ALLOCATE ( surf_lsm_v(l)%qsws_liq_eb(1:surf_lsm_v(l)%ns)      )
    2525           ALLOCATE ( surf_lsm_v(l)%qsws_veg_eb(1:surf_lsm_v(l)%ns)      )
     2758          ALLOCATE ( surf_lsm_v(l)%vegetation_surface(1:surf_lsm_v(l)%ns)  )
     2759          ALLOCATE ( surf_lsm_v(l)%pavement_surface(1:surf_lsm_v(l)%ns) )
     2760          ALLOCATE ( surf_lsm_v(l)%qsws_soil(1:surf_lsm_v(l)%ns)        )
     2761          ALLOCATE ( surf_lsm_v(l)%qsws_liq(1:surf_lsm_v(l)%ns)         )
     2762          ALLOCATE ( surf_lsm_v(l)%qsws_veg(1:surf_lsm_v(l)%ns)         )
    25262763          ALLOCATE ( surf_lsm_v(l)%rad_net_l(1:surf_lsm_v(l)%ns)        )
     2764          ALLOCATE ( surf_lsm_v(l)%rho_c_def(1:surf_lsm_h%ns)           ) 
    25272765          ALLOCATE ( surf_lsm_v(l)%r_a(1:surf_lsm_v(l)%ns)              )
    25282766          ALLOCATE ( surf_lsm_v(l)%r_canopy(1:surf_lsm_v(l)%ns)         )
     
    25312769          ALLOCATE ( surf_lsm_v(l)%r_s(1:surf_lsm_v(l)%ns)              )
    25322770          ALLOCATE ( surf_lsm_v(l)%r_canopy_min(1:surf_lsm_v(l)%ns)     )
    2533           ALLOCATE ( surf_lsm_v(l)%shf_eb(1:surf_lsm_v(l)%ns)           )
    2534           ALLOCATE ( surf_lsm_v(l)%soil_type_2d(1:surf_lsm_v(l)%ns)     )
    2535           ALLOCATE ( surf_lsm_v(l)%veg_type_2d(1:surf_lsm_v(l)%ns)      )
    25362771          ALLOCATE ( surf_lsm_v(l)%water_surface(1:surf_lsm_v(l)%ns)    )
    25372772
    2538           surf_lsm_v(l)%water_surface = .FALSE.
    2539           surf_lsm_v(l)%pave_surface  = .FALSE.
     2773          surf_lsm_v(l)%water_surface     = .FALSE.
     2774          surf_lsm_v(l)%pavement_surface  = .FALSE.
     2775          surf_lsm_v(l)%vegetation_surface   = .FALSE.
     2776         
    25402777       ENDDO
    25412778   
     
    25472784       t_surface_h => t_surface_h_1; t_surface_h_p => t_surface_h_2
    25482785       m_soil_h    => m_soil_h_1;    m_soil_h_p    => m_soil_h_2
    2549        m_liq_eb_h  => m_liq_eb_h_1;  m_liq_eb_h_p  => m_liq_eb_h_2
     2786       m_liq_h  => m_liq_h_1;  m_liq_h_p  => m_liq_h_2
    25502787!
    25512788!--    Vertical surfaces
     
    25532790       t_surface_v => t_surface_v_1; t_surface_v_p => t_surface_v_2
    25542791       m_soil_v    => m_soil_v_1;    m_soil_v_p    => m_soil_v_2
    2555        m_liq_eb_v  => m_liq_eb_v_1;  m_liq_eb_v_p  => m_liq_eb_v_2
     2792       m_liq_v  => m_liq_v_1;  m_liq_v_p  => m_liq_v_2
    25562793#endif
    25572794
     
    25742811       NAMELIST /lsm_par/         alpha_vangenuchten, c_surface,               &
    25752812                                  canopy_resistance_coefficient,               &
     2813                                  constant_roughness,                          &
    25762814                                  conserve_water_content,                      &
    25772815                                  f_shortwave_incoming, field_capacity,        &
     
    25812819                                  l_vangenuchten, min_canopy_resistance,       &
    25822820                                  min_soil_resistance, n_vangenuchten,         &
    2583                                   pave_depth, pave_heat_capacity,              &
    2584                                   pave_heat_conductivity,                      &
     2821                                  pavement_depth, pavement_heat_capacity,      &
     2822                                  pavement_heat_conduct, pavement_type,        &
    25852823                                  residual_moisture, root_fraction,            &
    25862824                                  saturation_moisture, skip_time_do_lsm,       &
    25872825                                  soil_moisture, soil_temperature, soil_type,  &
    2588                                   vegetation_coverage, veg_type, wilting_point,&
    2589                                   zs, z0_eb, z0h_eb, z0q_eb
     2826                                  surface_type,                                &
     2827                                  vegetation_coverage, vegetation_type,        &
     2828                                  water_temperature, water_type,               &
     2829                                  wilting_point, zs, z0_vegetation,            &
     2830                                  z0h_vegetation, z0q_vegetation, z0_water,    &
     2831                                  z0h_water, z0q_water, z0_pavement,           &
     2832                                  z0h_pavement, z0q_pavement
    25902833       
    25912834       line = ' '
     
    26682911       DO  m = 1, surf%ns
    26692912
    2670           IF ( surf%pave_surface(m) )  THEN
    2671              surf%rho_c_total(nzb_soil,m) = pave_heat_capacity
    2672              lambda_temp(nzb_soil)          = pave_heat_conductivity
    2673           ENDIF
    2674 
    26752913          IF (  .NOT.  surf%water_surface(m) )  THEN
    26762914             DO  k = nzb_soil, nzt_soil
    26772915
    2678                 IF ( surf%pave_surface(m)  .AND.  zs(k) <= pave_depth )  THEN
     2916                IF ( surf%pavement_surface(m)  .AND.  zs(k) <= pavement_depth )  THEN
    26792917                   
    2680                    surf%rho_c_total(k,m) = pave_heat_capacity
    2681                    lambda_temp(k)        = pave_heat_conductivity 
     2918                   surf%rho_c_total(k,m) = surf%rho_c_def(m)
     2919                   lambda_temp(k)   = surf%lambda_h_def(m)
    26822920
    26832921                ELSE           
     
    26862924!--                into account water content
    26872925                   surf%rho_c_total(k,m) = (rho_c_soil *                 &
    2688                                                ( 1.0_wp - surf%m_sat(m) )&
     2926                                               ( 1.0_wp - surf%m_sat(k,m) )&
    26892927                                               + rho_c_water * surf_m_soil%var_2d(k,m) )
    26902928
     
    26922930!--                Calculate soil heat conductivity at the center of the soil
    26932931!--                layers
    2694                    lambda_h_sat = lambda_h_sm**(1.0_wp - surf%m_sat(m)) *&
     2932                   lambda_h_sat = lambda_h_sm**(1.0_wp - surf%m_sat(k,m)) *      &
    26952933                                  lambda_h_water ** surf_m_soil%var_2d(k,m)
    26962934
    2697                    ke = 1.0_wp + LOG10( MAX( 0.1_wp, surf_m_soil%var_2d(k,m) /             &
    2698                                                      surf%m_sat(m) ) )
     2935                   ke = 1.0_wp + LOG10( MAX( 0.1_wp, surf_m_soil%var_2d(k,m) / &
     2936                                                     surf%m_sat(k,m) ) )
    26992937
    27002938                   lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) +       &
     
    27042942
    27052943!
    2706 !--          Calculate soil heat conductivity (lambda_h) at the _stag level
     2944!--          Calculate soil heat conductivity (lambda_h) at the _layer level
    27072945!--          using linear interpolation. For pavement surface, the
    27082946!--          true pavement depth is considered
    27092947             DO  k = nzb_soil, nzt_soil-1
    2710                 IF ( surf%pave_surface(m)  .AND.  zs(k)   < pave_depth   &
    2711                                                  .AND.  zs(k+1) > pave_depth )  THEN
    2712                    surf%lambda_h(k,m) = ( pave_depth - zs(k) ) / dz_soil(k+1)&
    2713                                      * lambda_temp(k)                          &
    2714                                      + ( 1.0_wp - ( pave_depth - zs(k) )       &
    2715                                      / dz_soil(k+1) ) * lambda_temp(k+1)
     2948                IF ( surf%pavement_surface(m)  .AND.  zs(k) < pavement_depth   &
     2949                     .AND.  zs(k+1) > pavement_depth )                         &
     2950                THEN
     2951                   surf%lambda_h(k,m) = ( pavement_depth - zs(k) ) * ddz_soil(k+1)   &
     2952                                      * lambda_temp(k)                               &
     2953                                      + ( zs(k+1) - pavement_depth ) * ddz_soil(k+1) &
     2954                                      * lambda_temp(k+1)
    27162955                ELSE
    2717                    surf%lambda_h(k,m) = ( lambda_temp(k+1) + lambda_temp(k) )&
    2718                                      * 0.5_wp
     2956                   surf%lambda_h(k,m) = ( lambda_temp(k+1) + lambda_temp(k) )  &
     2957                                        * 0.5_wp
    27192958                ENDIF
     2959
    27202960             ENDDO
    27212961             surf%lambda_h(nzt_soil,m) = lambda_temp(nzt_soil)
     
    27252965             tend(:) = 0.0_wp
    27262966
    2727              tend(nzb_soil) = ( 1.0_wp / surf%rho_c_total(nzb_soil,m) ) *&
     2967             tend(nzb_soil) = ( 1.0_wp / surf%rho_c_total(nzb_soil,m) ) *            &
    27282968                    ( surf%lambda_h(nzb_soil,m) * ( surf_t_soil%var_2d(nzb_soil+1,m) &
    2729                       - surf_t_soil%var_2d(nzb_soil,m) ) * ddz_soil(nzb_soil+1)            &
    2730                       + surf%ghf_eb(m) ) * ddz_soil_stag(nzb_soil)
     2969                      - surf_t_soil%var_2d(nzb_soil,m) ) * ddz_soil(nzb_soil)        &
     2970                      + surf%ghf(m) ) * ddz_soil_layer(nzb_soil)
    27312971
    27322972             DO  k = nzb_soil+1, nzt_soil
    2733                 tend(k) = ( 1.0_wp / surf%rho_c_total(k,m) )             &
    2734                           * (   surf%lambda_h(k,m)                       &
    2735                               * ( surf_t_soil%var_2d(k+1,m) - surf_t_soil%var_2d(k,m) )                &
    2736                               * ddz_soil(k+1)                                  &
    2737                               - surf%lambda_h(k-1,m)                     &
    2738                               * ( surf_t_soil%var_2d(k,m) - surf_t_soil%var_2d(k-1,m) )                &
    2739                               * ddz_soil(k)                                    &
    2740                             ) * ddz_soil_stag(k)
     2973                tend(k) = ( 1.0_wp / surf%rho_c_total(k,m) )                            &
     2974                          * (   surf%lambda_h(k,m)                                      &
     2975                              * ( surf_t_soil%var_2d(k+1,m) - surf_t_soil%var_2d(k,m) ) &
     2976                              * ddz_soil(k)                                             &
     2977                              - surf%lambda_h(k-1,m)                                    &
     2978                              * ( surf_t_soil%var_2d(k,m) - surf_t_soil%var_2d(k-1,m) ) &
     2979                              * ddz_soil(k-1)                                           &
     2980                            ) * ddz_soil_layer(k)
    27412981
    27422982             ENDDO
     
    27693009!
    27703010!--             Calculate soil diffusivity at the center of the soil layers
    2771                 lambda_temp(k) = (- b_ch * surf%gamma_w_sat(m) * psi_sat     &
    2772                                   / surf%m_sat(m) ) * ( MAX( surf_m_soil%var_2d(k,m),    &
    2773                                   surf%m_wilt(m) ) / surf%m_sat(m) )**(&
     3011                lambda_temp(k) = (- b_ch * surf%gamma_w_sat(k,m) * psi_sat     &
     3012                                  / surf%m_sat(k,m) ) * ( MAX( surf_m_soil%var_2d(k,m),    &
     3013                                  surf%m_wilt(k,m) ) / surf%m_sat(k,m) )**(&
    27743014                                                            b_ch + 2.0_wp )
    27753015
    27763016!
    27773017!--             Parametrization of Van Genuchten
    2778                 IF ( soil_type /= 7 )  THEN
    2779 !
    2780 !--                Calculate the hydraulic conductivity after Van Genuchten
    2781 !--                (1980)
    2782                    h_vg = ( ( ( surf%m_res(m) - surf%m_sat(m) ) /  &
    2783                               ( surf%m_res(m) -                          &
    2784                                 MAX( surf_m_soil%var_2d(k,m), surf%m_wilt(m) )       &
    2785                               )                                                &
    2786                             )**(                                               &
    2787                           surf%n_vg(m) / ( surf%n_vg(m) - 1.0_wp ) &
    2788                                ) - 1.0_wp                                      &
    2789                           )**( 1.0_wp / surf%n_vg(m) ) / surf%alpha_vg(m)
    2790 
    2791                    gamma_temp(k) = surf%gamma_w_sat(m) * ( ( ( 1.0_wp +  &
    2792                           ( surf%alpha_vg(m) * h_vg )**surf%n_vg(m)&
     3018!--             Calculate the hydraulic conductivity after Van Genuchten (1980)
     3019                h_vg = ( ( ( surf%m_res(k,m) - surf%m_sat(k,m) ) /  &
     3020                           ( surf%m_res(k,m) -                          &
     3021                             MAX( surf_m_soil%var_2d(k,m), surf%m_wilt(k,m) )       &
     3022                           )                                                &
     3023                         )**(                                               &
     3024                       surf%n_vg(k,m) / ( surf%n_vg(k,m) - 1.0_wp ) &
     3025                            ) - 1.0_wp                                      &
     3026                       )**( 1.0_wp / surf%n_vg(k,m) ) / surf%alpha_vg(k,m)
     3027
     3028                gamma_temp(k) = surf%gamma_w_sat(k,m) * ( ( ( 1.0_wp +  &
     3029                       ( surf%alpha_vg(k,m) * h_vg )**surf%n_vg(k,m)&
    27933030                                                               )**(            &
    2794                               1.0_wp - 1.0_wp / surf%n_vg(m)) - (        &
    2795                           surf%alpha_vg(m) * h_vg )**( surf%n_vg(m)&
    2796                               - 1.0_wp) )**2 )                                 &
    2797                               / ( ( 1.0_wp + ( surf%alpha_vg(m) * h_vg   &
    2798                               )**surf%n_vg(m) )**( ( 1.0_wp  - 1.0_wp    &
    2799                               / surf%n_vg(m) ) *                         &
    2800                               ( surf%l_vg(m) + 2.0_wp) ) )
    2801 !
    2802 !--             Parametrization of Clapp & Hornberger
    2803                 ELSE
    2804                    gamma_temp(k) = surf%gamma_w_sat(m) * ( surf_m_soil%var_2d(k,m)   &
    2805                                    / surf%m_sat(m) )**(2.0_wp * b_ch + 3.0_wp)
    2806                 ENDIF
     3031                           1.0_wp - 1.0_wp / surf%n_vg(k,m)) - (        &
     3032                       surf%alpha_vg(k,m) * h_vg )**( surf%n_vg(k,m)&
     3033                           - 1.0_wp) )**2 )                                 &
     3034                           / ( ( 1.0_wp + ( surf%alpha_vg(k,m) * h_vg   &
     3035                           )**surf%n_vg(k,m) )**( ( 1.0_wp  - 1.0_wp    &
     3036                           / surf%n_vg(k,m) ) *                         &
     3037                           ( surf%l_vg(k,m) + 2.0_wp) ) )
    28073038
    28083039             ENDDO
     
    28123043!--          when humidity is enabled in the atmosphere and the surface type
    28133044!--          is not pavement (implies dry soil below).
    2814              IF ( humidity  .AND.  .NOT.  surf%pave_surface(m) )  THEN
    2815 !
    2816 !--             Calculate soil diffusivity (lambda_w) at the _stag level
     3045             IF ( humidity  .AND.  .NOT.  surf%pavement_surface(m) )  THEN
     3046!
     3047!--             Calculate soil diffusivity (lambda_w) at the _layer level
    28173048!--             using linear interpolation. To do: replace this with
    28183049!--             ECMWF-IFS Eq. 8.81
     
    28373068                ENDIF     
    28383069
    2839 !--             The root extraction (= root_extr * qsws_veg_eb / (rho_l     
     3070!--             The root extraction (= root_extr * qsws_veg / (rho_l     
    28403071!--             * l_v)) ensures the mass conservation for water. The         
    28413072!--             transpiration of plants equals the cumulative withdrawals by
     
    28523083                m_total = 0.0_wp
    28533084                DO  k = nzb_soil, nzt_soil
    2854                     IF ( surf_m_soil%var_2d(k,m) > surf%m_wilt(m) )  THEN
     3085                    IF ( surf_m_soil%var_2d(k,m) > surf%m_wilt(k,m) )  THEN
    28553086                       m_total = m_total + surf%root_fr(k,m) * surf_m_soil%var_2d(k,m)
    28563087                    ENDIF
     
    28583089                 IF ( m_total > 0.0_wp )  THEN
    28593090                   DO  k = nzb_soil, nzt_soil
    2860                       IF ( surf_m_soil%var_2d(k,m) > surf%m_wilt(m) )  THEN
     3091                      IF ( surf_m_soil%var_2d(k,m) > surf%m_wilt(k,m) )  THEN
    28613092                         root_extr(k) = surf%root_fr(k,m) * surf_m_soil%var_2d(k,m)  &
    28623093                                                            / m_total
     
    28723103                tend(nzb_soil) = ( surf%lambda_w(nzb_soil,m) *   (       &
    28733104                         surf_m_soil%var_2d(nzb_soil+1,m) - surf_m_soil%var_2d(nzb_soil,m) )           &
    2874                          * ddz_soil(nzb_soil+1) - surf%gamma_w(nzb_soil,m) - &
     3105                         * ddz_soil(nzb_soil) - surf%gamma_w(nzb_soil,m) - &
    28753106                                                     (                         &
    2876                             root_extr(nzb_soil) * surf%qsws_veg_eb(m)    &
    2877                             + surf%qsws_soil_eb(m) ) * drho_l_lv )       &
    2878                             * ddz_soil_stag(nzb_soil)
     3107                            root_extr(nzb_soil) * surf%qsws_veg(m)    &
     3108                            + surf%qsws_soil(m) ) * drho_l_lv )       &
     3109                            * ddz_soil_layer(nzb_soil)
    28793110
    28803111                DO  k = nzb_soil+1, nzt_soil-1
    28813112                   tend(k) = ( surf%lambda_w(k,m) * ( surf_m_soil%var_2d(k+1,m)      &
    2882                              - surf_m_soil%var_2d(k,m) ) * ddz_soil(k+1)                   &
     3113                             - surf_m_soil%var_2d(k,m) ) * ddz_soil(k)                   &
    28833114                             - surf%gamma_w(k,m)                         &
    28843115                             - surf%lambda_w(k-1,m) * ( surf_m_soil%var_2d(k,m) -    &
    2885                              surf_m_soil%var_2d(k-1,m)) * ddz_soil(k)                      &
     3116                             surf_m_soil%var_2d(k-1,m)) * ddz_soil(k-1)                      &
    28863117                             + surf%gamma_w(k-1,m) - (root_extr(k)       &
    2887                              * surf%qsws_veg_eb(m) * drho_l_lv)          &
    2888                              ) * ddz_soil_stag(k)
     3118                             * surf%qsws_veg(m) * drho_l_lv)          &
     3119                             ) * ddz_soil_layer(k)
    28893120                ENDDO
    28903121                tend(nzt_soil) = ( - surf%gamma_w(nzt_soil,m)            &
     
    28923123                                     * ( surf_m_soil%var_2d(nzt_soil,m)                    &
    28933124                                     - surf_m_soil%var_2d(nzt_soil-1,m))                   &
    2894                                      * ddz_soil(nzt_soil)                      &
     3125                                     * ddz_soil(nzt_soil-1)                      &
    28953126                                     + surf%gamma_w(nzt_soil-1,m) - (    &
    28963127                                       root_extr(nzt_soil)                     &
    2897                                      * surf%qsws_veg_eb(m) * drho_l_lv ) &
    2898                                   ) * ddz_soil_stag(nzt_soil)             
     3128                                     * surf%qsws_veg(m) * drho_l_lv ) &
     3129                                  ) * ddz_soil_layer(nzt_soil)             
    28993130
    29003131                surf_m_soil_p%var_2d(nzb_soil:nzt_soil,m) = surf_m_soil%var_2d(nzb_soil:nzt_soil,m)    &
     
    29503181       IF ( humidity )  THEN
    29513182          m_soil_h    = m_soil_h_p
    2952           m_liq_eb_h  = m_liq_eb_h_p
     3183          m_liq_h  = m_liq_h_p
    29533184       ENDIF
    29543185!
     
    29583189       IF ( humidity )  THEN
    29593190          m_soil_v    = m_soil_v_p
    2960           m_liq_eb_v  = m_liq_eb_v_p
     3191          m_liq_v  = m_liq_v_p
    29613192       ENDIF
    29623193
     
    29723203             IF ( humidity )  THEN
    29733204                m_soil_h    => m_soil_h_1;   m_soil_h_p    => m_soil_h_2
    2974                 m_liq_eb_h  => m_liq_eb_h_1; m_liq_eb_h_p  => m_liq_eb_h_2
     3205                m_liq_h  => m_liq_h_1; m_liq_h_p  => m_liq_h_2
    29753206             ENDIF
    29763207!
     
    29803211             IF ( humidity )  THEN
    29813212                m_soil_v    => m_soil_v_1;   m_soil_v_p    => m_soil_v_2
    2982                 m_liq_eb_v  => m_liq_eb_v_1; m_liq_eb_v_p  => m_liq_eb_v_2
     3213                m_liq_v  => m_liq_v_1; m_liq_v_p  => m_liq_v_2
    29833214             ENDIF
    29843215
     
    29923223             IF ( humidity )  THEN
    29933224                m_soil_h    => m_soil_h_2;   m_soil_h_p    => m_soil_h_1
    2994                 m_liq_eb_h  => m_liq_eb_h_2; m_liq_eb_h_p  => m_liq_eb_h_1
     3225                m_liq_h  => m_liq_h_2; m_liq_h_p  => m_liq_h_1
    29953226             ENDIF
    29963227!
     
    30003231             IF ( humidity )  THEN
    30013232                m_soil_v    => m_soil_v_2;   m_soil_v_p    => m_soil_v_1
    3002                 m_liq_eb_v  => m_liq_eb_v_2; m_liq_eb_v_p  => m_liq_eb_v_1
     3233                m_liq_v  => m_liq_v_2; m_liq_v_p  => m_liq_v_1
    30033234             ENDIF
    30043235
     
    30583289                c_veg_av = 0.0_wp
    30593290
    3060              CASE ( 'ghf_eb*' )
    3061                 IF ( .NOT. ALLOCATED( ghf_eb_av ) )  THEN
    3062                    ALLOCATE( ghf_eb_av(nysg:nyng,nxlg:nxrg) )
     3291             CASE ( 'ghf*' )
     3292                IF ( .NOT. ALLOCATED( ghf_av ) )  THEN
     3293                   ALLOCATE( ghf_av(nysg:nyng,nxlg:nxrg) )
    30633294                ENDIF
    3064                 ghf_eb_av = 0.0_wp
     3295                ghf_av = 0.0_wp
    30653296
    30663297             CASE ( 'lai*' )
     
    30703301                lai_av = 0.0_wp
    30713302
    3072              CASE ( 'm_liq_eb*' )
    3073                 IF ( .NOT. ALLOCATED( m_liq_eb_av ) )  THEN
    3074                    ALLOCATE( m_liq_eb_av(nysg:nyng,nxlg:nxrg) )
     3303             CASE ( 'm_liq*' )
     3304                IF ( .NOT. ALLOCATED( m_liq_av ) )  THEN
     3305                   ALLOCATE( m_liq_av(nysg:nyng,nxlg:nxrg) )
    30753306                ENDIF
    3076                 m_liq_eb_av = 0.0_wp
     3307                m_liq_av = 0.0_wp
    30773308
    30783309             CASE ( 'm_soil' )
     
    30823313                m_soil_av = 0.0_wp
    30833314
    3084              CASE ( 'qsws_eb*' )
    3085                 IF ( .NOT. ALLOCATED( qsws_eb_av ) )  THEN
    3086                    ALLOCATE( qsws_eb_av(nysg:nyng,nxlg:nxrg) )
     3315             CASE ( 'qsws_liq*' )
     3316                IF ( .NOT. ALLOCATED( qsws_liq_av ) )  THEN
     3317                   ALLOCATE( qsws_liq_av(nysg:nyng,nxlg:nxrg) )
    30873318                ENDIF
    3088                 qsws_eb_av = 0.0_wp
    3089 
    3090              CASE ( 'qsws_liq_eb*' )
    3091                 IF ( .NOT. ALLOCATED( qsws_liq_eb_av ) )  THEN
    3092                    ALLOCATE( qsws_liq_eb_av(nysg:nyng,nxlg:nxrg) )
     3319                qsws_liq_av = 0.0_wp
     3320
     3321             CASE ( 'qsws_soil*' )
     3322                IF ( .NOT. ALLOCATED( qsws_soil_av ) )  THEN
     3323                   ALLOCATE( qsws_soil_av(nysg:nyng,nxlg:nxrg) )
    30933324                ENDIF
    3094                 qsws_liq_eb_av = 0.0_wp
    3095 
    3096              CASE ( 'qsws_soil_eb*' )
    3097                 IF ( .NOT. ALLOCATED( qsws_soil_eb_av ) )  THEN
    3098                    ALLOCATE( qsws_soil_eb_av(nysg:nyng,nxlg:nxrg) )
     3325                qsws_soil_av = 0.0_wp
     3326
     3327             CASE ( 'qsws_veg*' )
     3328                IF ( .NOT. ALLOCATED( qsws_veg_av ) )  THEN
     3329                   ALLOCATE( qsws_veg_av(nysg:nyng,nxlg:nxrg) )
    30993330                ENDIF
    3100                 qsws_soil_eb_av = 0.0_wp
    3101 
    3102              CASE ( 'qsws_veg_eb*' )
    3103                 IF ( .NOT. ALLOCATED( qsws_veg_eb_av ) )  THEN
    3104                    ALLOCATE( qsws_veg_eb_av(nysg:nyng,nxlg:nxrg) )
    3105                 ENDIF
    3106                 qsws_veg_eb_av = 0.0_wp
     3331                qsws_veg_av = 0.0_wp
    31073332
    31083333             CASE ( 'r_a*' )
     
    31183343                r_s_av = 0.0_wp
    31193344
    3120              CASE ( 'shf_eb*' )
    3121                 IF ( .NOT. ALLOCATED( shf_eb_av ) )  THEN
    3122                    ALLOCATE( shf_eb_av(nysg:nyng,nxlg:nxrg) )
    3123                 ENDIF
    3124                 shf_eb_av = 0.0_wp
    3125 
    31263345             CASE ( 't_soil' )
    31273346                IF ( .NOT. ALLOCATED( t_soil_av ) )  THEN
     
    31603379             ENDDO
    31613380
    3162           CASE ( 'ghf_eb*' )
     3381          CASE ( 'ghf*' )
    31633382             DO  m = 1, surf_lsm_h%ns
    31643383                i   = surf_lsm_h%i(m)           
    31653384                j   = surf_lsm_h%j(m)
    3166                 ghf_eb_av(j,i) = ghf_eb_av(j,i) + surf_lsm_h%ghf_eb(m)
     3385                ghf_av(j,i) = ghf_av(j,i) + surf_lsm_h%ghf(m)
    31673386             ENDDO
    31683387
     
    31743393             ENDDO
    31753394
    3176           CASE ( 'm_liq_eb*' )
     3395          CASE ( 'm_liq*' )
    31773396             DO  m = 1, surf_lsm_h%ns
    31783397                i   = surf_lsm_h%i(m)           
    31793398                j   = surf_lsm_h%j(m)
    3180                 m_liq_eb_av(j,i) = m_liq_eb_av(j,i) + m_liq_eb_h%var_1d(m)
     3399                m_liq_av(j,i) = m_liq_av(j,i) + m_liq_h%var_1d(m)
    31813400             ENDDO
    31823401
     
    31903409             ENDDO
    31913410
    3192           CASE ( 'qsws_eb*' )
     3411          CASE ( 'qsws_liq*' )
    31933412             DO  m = 1, surf_lsm_h%ns
    31943413                i   = surf_lsm_h%i(m)           
    31953414                j   = surf_lsm_h%j(m)
    3196                 qsws_eb_av(j,i) = qsws_eb_av(j,i) + surf_lsm_h%qsws_eb(m)
    3197              ENDDO
    3198 
    3199           CASE ( 'qsws_liq_eb*' )
     3415                qsws_liq_av(j,i) = qsws_liq_av(j,i) +                    &
     3416                                      surf_lsm_h%qsws_liq(m)
     3417             ENDDO
     3418
     3419          CASE ( 'qsws_soil*' )
    32003420             DO  m = 1, surf_lsm_h%ns
    32013421                i   = surf_lsm_h%i(m)           
    32023422                j   = surf_lsm_h%j(m)
    3203                 qsws_liq_eb_av(j,i) = qsws_liq_eb_av(j,i) +                    &
    3204                                       surf_lsm_h%qsws_liq_eb(m)
    3205              ENDDO
    3206 
    3207           CASE ( 'qsws_soil_eb*' )
     3423                qsws_soil_av(j,i) = qsws_soil_av(j,i) +                  &
     3424                                       surf_lsm_h%qsws_soil(m)
     3425             ENDDO
     3426
     3427          CASE ( 'qsws_veg*' )
    32083428             DO  m = 1, surf_lsm_h%ns
    32093429                i   = surf_lsm_h%i(m)           
    32103430                j   = surf_lsm_h%j(m)
    3211                 qsws_soil_eb_av(j,i) = qsws_soil_eb_av(j,i) +                  &
    3212                                        surf_lsm_h%qsws_soil_eb(m)
    3213              ENDDO
    3214 
    3215           CASE ( 'qsws_veg_eb*' )
    3216              DO  m = 1, surf_lsm_h%ns
    3217                 i   = surf_lsm_h%i(m)           
    3218                 j   = surf_lsm_h%j(m)
    3219                 qsws_veg_eb_av(j,i) = qsws_veg_eb_av(j,i) +                    &
    3220                                       surf_lsm_h%qsws_veg_eb(m)
     3431                qsws_veg_av(j,i) = qsws_veg_av(j,i) +                    &
     3432                                      surf_lsm_h%qsws_veg(m)
    32213433             ENDDO
    32223434
     
    32333445                j   = surf_lsm_h%j(m)
    32343446                r_s_av(j,i) = r_s_av(j,i) + surf_lsm_h%r_s(m)
    3235              ENDDO
    3236 
    3237           CASE ( 'shf_eb*' )
    3238              DO  m = 1, surf_lsm_h%ns
    3239                 i   = surf_lsm_h%i(m)           
    3240                 j   = surf_lsm_h%j(m)
    3241                 shf_eb_av(j,i) = shf_eb_av(j,i) + surf_lsm_h%shf_eb(m)
    32423447             ENDDO
    32433448
     
    32813486             ENDDO
    32823487
    3283           CASE ( 'ghf_eb*' )
     3488          CASE ( 'ghf*' )
    32843489             DO  i = nxl, nxr
    32853490                DO  j = nys, nyn
    3286                    ghf_eb_av(j,i) = ghf_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
     3491                   ghf_av(j,i) = ghf_av(j,i) / REAL( average_count_3d, KIND=wp )
    32873492                ENDDO
    32883493             ENDDO
     
    32953500             ENDDO
    32963501
    3297           CASE ( 'm_liq_eb*' )
     3502          CASE ( 'm_liq*' )
    32983503             DO  i = nxl, nxr
    32993504                DO  j = nys, nyn
    3300                    m_liq_eb_av(j,i) = m_liq_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
     3505                   m_liq_av(j,i) = m_liq_av(j,i) / REAL( average_count_3d, KIND=wp )
    33013506                ENDDO
    33023507             ENDDO
     
    33113516             ENDDO
    33123517
    3313           CASE ( 'qsws_eb*' )
     3518          CASE ( 'qsws_liq*' )
    33143519             DO  i = nxl, nxr
    33153520                DO  j = nys, nyn
    3316                    qsws_eb_av(j,i) = qsws_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
     3521                   qsws_liq_av(j,i) = qsws_liq_av(j,i) / REAL( average_count_3d, KIND=wp )
    33173522                ENDDO
    33183523             ENDDO
    33193524
    3320           CASE ( 'qsws_liq_eb*' )
     3525          CASE ( 'qsws_soil*' )
    33213526             DO  i = nxl, nxr
    33223527                DO  j = nys, nyn
    3323                    qsws_liq_eb_av(j,i) = qsws_liq_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
     3528                   qsws_soil_av(j,i) = qsws_soil_av(j,i) / REAL( average_count_3d, KIND=wp )
    33243529                ENDDO
    33253530             ENDDO
    33263531
    3327           CASE ( 'qsws_soil_eb*' )
     3532          CASE ( 'qsws_veg*' )
    33283533             DO  i = nxl, nxr
    33293534                DO  j = nys, nyn
    3330                    qsws_soil_eb_av(j,i) = qsws_soil_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
    3331                 ENDDO
    3332              ENDDO
    3333 
    3334           CASE ( 'qsws_veg_eb*' )
    3335              DO  i = nxl, nxr
    3336                 DO  j = nys, nyn
    3337                    qsws_veg_eb_av(j,i) = qsws_veg_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
     3535                   qsws_veg_av(j,i) = qsws_veg_av(j,i) / REAL( average_count_3d, KIND=wp )
    33383536                ENDDO
    33393537             ENDDO
     
    35033701          grid = 'zu1'
    35043702
    3505        CASE ( 'ghf_eb*_xy' )        ! 2d-array
     3703       CASE ( 'ghf*_xy' )        ! 2d-array
    35063704          IF ( av == 0 )  THEN
    35073705             DO  m = 1, surf_lsm_h%ns
    35083706                i                   = surf_lsm_h%i(m)           
    35093707                j                   = surf_lsm_h%j(m)
    3510                 local_pf(i,j,nzb+1) = surf_lsm_h%ghf_eb(m)
     3708                local_pf(i,j,nzb+1) = surf_lsm_h%ghf(m)
    35113709             ENDDO
    35123710          ELSE
    35133711             DO  i = nxlg, nxrg
    35143712                DO  j = nysg, nyng
    3515                    local_pf(i,j,nzb+1) = ghf_eb_av(j,i)
     3713                   local_pf(i,j,nzb+1) = ghf_av(j,i)
    35163714                ENDDO
    35173715             ENDDO
     
    35393737          grid = 'zu1'
    35403738
    3541        CASE ( 'm_liq_eb*_xy' )        ! 2d-array
     3739       CASE ( 'm_liq*_xy' )        ! 2d-array
    35423740          IF ( av == 0 )  THEN
    35433741             DO  m = 1, surf_lsm_h%ns
    35443742                i                   = surf_lsm_h%i(m)           
    35453743                j                   = surf_lsm_h%j(m)
    3546                 local_pf(i,j,nzb+1) = m_liq_eb_h%var_1d(m)
     3744                local_pf(i,j,nzb+1) = m_liq_h%var_1d(m)
    35473745             ENDDO
    35483746          ELSE
    35493747             DO  i = nxlg, nxrg
    35503748                DO  j = nysg, nyng
    3551                    local_pf(i,j,nzb+1) = m_liq_eb_av(j,i)
     3749                   local_pf(i,j,nzb+1) = m_liq_av(j,i)
    35523750                ENDDO
    35533751             ENDDO
     
    35813779          IF ( mode == 'xy' ) grid = 'zs'
    35823780
    3583        CASE ( 'qsws_eb*_xy' )        ! 2d-array
     3781       CASE ( 'qsws_liq*_xy' )        ! 2d-array
    35843782          IF ( av == 0 ) THEN
    35853783             DO  m = 1, surf_lsm_h%ns
    35863784                i                   = surf_lsm_h%i(m)           
    35873785                j                   = surf_lsm_h%j(m)
    3588                 local_pf(i,j,nzb+1) =  surf_lsm_h%qsws_eb(m)
     3786                local_pf(i,j,nzb+1) = surf_lsm_h%qsws_liq(m)
    35893787             ENDDO
    35903788          ELSE
    35913789             DO  i = nxlg, nxrg
    35923790                DO  j = nysg, nyng
    3593                    local_pf(i,j,nzb+1) =  qsws_eb_av(j,i)
     3791                   local_pf(i,j,nzb+1) =  qsws_liq_av(j,i)
    35943792                ENDDO
    35953793             ENDDO
     
    35993797          grid = 'zu1'
    36003798
    3601        CASE ( 'qsws_liq_eb*_xy' )        ! 2d-array
     3799       CASE ( 'qsws_soil*_xy' )        ! 2d-array
    36023800          IF ( av == 0 ) THEN
    36033801             DO  m = 1, surf_lsm_h%ns
    36043802                i                   = surf_lsm_h%i(m)           
    36053803                j                   = surf_lsm_h%j(m)
    3606                 local_pf(i,j,nzb+1) = surf_lsm_h%qsws_liq_eb(m)
     3804                local_pf(i,j,nzb+1) =  surf_lsm_h%qsws_soil(m)
    36073805             ENDDO
    36083806          ELSE
    36093807             DO  i = nxlg, nxrg
    36103808                DO  j = nysg, nyng
    3611                    local_pf(i,j,nzb+1) =  qsws_liq_eb_av(j,i)
     3809                   local_pf(i,j,nzb+1) =  qsws_soil_av(j,i)
    36123810                ENDDO
    36133811             ENDDO
     
    36173815          grid = 'zu1'
    36183816
    3619        CASE ( 'qsws_soil_eb*_xy' )        ! 2d-array
     3817       CASE ( 'qsws_veg*_xy' )        ! 2d-array
    36203818          IF ( av == 0 ) THEN
    36213819             DO  m = 1, surf_lsm_h%ns
    36223820                i                   = surf_lsm_h%i(m)           
    36233821                j                   = surf_lsm_h%j(m)
    3624                 local_pf(i,j,nzb+1) =  surf_lsm_h%qsws_soil_eb(m)
     3822                local_pf(i,j,nzb+1) =  surf_lsm_h%qsws_veg(m)
    36253823             ENDDO
    36263824          ELSE
    36273825             DO  i = nxlg, nxrg
    36283826                DO  j = nysg, nyng
    3629                    local_pf(i,j,nzb+1) =  qsws_soil_eb_av(j,i)
    3630                 ENDDO
    3631              ENDDO
    3632           ENDIF
    3633 
    3634           two_d = .TRUE.
    3635           grid = 'zu1'
    3636 
    3637        CASE ( 'qsws_veg_eb*_xy' )        ! 2d-array
    3638           IF ( av == 0 ) THEN
    3639              DO  m = 1, surf_lsm_h%ns
    3640                 i                   = surf_lsm_h%i(m)           
    3641                 j                   = surf_lsm_h%j(m)
    3642                 local_pf(i,j,nzb+1) =  surf_lsm_h%qsws_veg_eb(m)
    3643              ENDDO
    3644           ELSE
    3645              DO  i = nxlg, nxrg
    3646                 DO  j = nysg, nyng
    3647                    local_pf(i,j,nzb+1) =  qsws_veg_eb_av(j,i)
     3827                   local_pf(i,j,nzb+1) =  qsws_veg_av(j,i)
    36483828                ENDDO
    36493829             ENDDO
     
    36833863                DO  j = nysg, nyng
    36843864                   local_pf(i,j,nzb+1) = r_s_av(j,i)
    3685                 ENDDO
    3686              ENDDO
    3687           ENDIF
    3688 
    3689           two_d = .TRUE.
    3690           grid = 'zu1'
    3691 
    3692        CASE ( 'shf_eb*_xy' )        ! 2d-array
    3693           IF ( av == 0 ) THEN
    3694              DO  m = 1, surf_lsm_h%ns
    3695                 i                   = surf_lsm_h%i(m)           
    3696                 j                   = surf_lsm_h%j(m)
    3697                 local_pf(i,j,nzb+1) =  surf_lsm_h%shf_eb(m)
    3698              ENDDO
    3699           ELSE
    3700              DO  i = nxlg, nxrg
    3701                 DO  j = nysg, nyng
    3702                    local_pf(i,j,nzb+1) =  shf_eb_av(j,i)
    37033865                ENDDO
    37043866             ENDDO
     
    38524014          WRITE ( 14 )  'c_veg_av            ';  WRITE ( 14 ) c_veg_av
    38534015       ENDIF
    3854        IF ( ALLOCATED( ghf_eb_av ) )  THEN
    3855           WRITE ( 14 )  'ghf_eb_av           ';  WRITE ( 14 )  ghf_eb_av
     4016       IF ( ALLOCATED( ghf_av ) )  THEN
     4017          WRITE ( 14 )  'ghf_av           ';  WRITE ( 14 )  ghf_av
    38564018       ENDIF
    38574019       IF ( ALLOCATED( lai_av ) )  THEN
    38584020          WRITE ( 14 )  'lai_av              ';  WRITE ( 14 )  lai_av
    38594021       ENDIF
    3860        WRITE ( 14 )     'm_liq_eb            ';  WRITE ( 14 )  m_liq_eb_h%var_1d
    3861        IF ( ALLOCATED( m_liq_eb_av ) )  THEN
    3862           WRITE ( 14 )  'm_liq_eb_av         ';  WRITE ( 14 )  m_liq_eb_av
     4022       WRITE ( 14 )     'm_liq            ';  WRITE ( 14 )  m_liq_h%var_1d
     4023       IF ( ALLOCATED( m_liq_av ) )  THEN
     4024          WRITE ( 14 )  'm_liq_av         ';  WRITE ( 14 )  m_liq_av
    38634025       ENDIF
    38644026       WRITE ( 14 )     'm_soil              ';  WRITE ( 14 )  m_soil_h%var_2d
     
    38664028          WRITE ( 14 )  'm_soil_av           ';  WRITE ( 14 )  m_soil_av
    38674029       ENDIF
    3868        IF ( ALLOCATED( qsws_eb_av ) )  THEN
    3869           WRITE ( 14 )  'qsws_eb_av          ';  WRITE ( 14 )  qsws_eb_av
    3870        ENDIF   
    3871        IF ( ALLOCATED( qsws_liq_eb_av ) )  THEN
    3872           WRITE ( 14 )  'qsws_liq_eb_av      ';  WRITE ( 14 )  qsws_liq_eb_av
     4030       IF ( ALLOCATED( qsws_liq_av ) )  THEN
     4031          WRITE ( 14 )  'qsws_liq_av      ';  WRITE ( 14 )  qsws_liq_av
    38734032       ENDIF 
    3874        IF ( ALLOCATED( qsws_soil_eb_av ) )  THEN
    3875           WRITE ( 14 )  'qsws_soil_eb_av     ';  WRITE ( 14 )  qsws_soil_eb_av
     4033       IF ( ALLOCATED( qsws_soil_av ) )  THEN
     4034          WRITE ( 14 )  'qsws_soil_av     ';  WRITE ( 14 )  qsws_soil_av
    38764035       ENDIF
    3877        IF ( ALLOCATED( qsws_veg_eb_av ) )  THEN
    3878           WRITE ( 14 )  'qsws_veg_eb_av      ';  WRITE ( 14 )  qsws_veg_eb_av
     4036       IF ( ALLOCATED( qsws_veg_av ) )  THEN
     4037          WRITE ( 14 )  'qsws_veg_av      ';  WRITE ( 14 )  qsws_veg_av
    38794038       ENDIF
    3880        IF ( ALLOCATED( shf_eb_av ) )  THEN
    3881           WRITE ( 14 )  'shf_eb_av           ';  WRITE ( 14 )  shf_eb_av
    3882        ENDIF
    38834039       WRITE ( 14 )     't_soil              ';  WRITE ( 14 )  t_soil_h%var_2d
    38844040       IF ( ALLOCATED( t_soil_av ) )  THEN
     
    40024158                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    40034159
    4004                 CASE ( 'ghf_eb_av' )
    4005                    IF ( .NOT. ALLOCATED( ghf_eb_av ) )  THEN
    4006                       ALLOCATE( ghf_eb_av(nysg:nyng,nxlg:nxrg) )
     4160                CASE ( 'ghf_av' )
     4161                   IF ( .NOT. ALLOCATED( ghf_av ) )  THEN
     4162                      ALLOCATE( ghf_av(nysg:nyng,nxlg:nxrg) )
    40074163                   ENDIF
    40084164                   IF ( k == 1 )  READ ( 13 )  tmp_2d
    4009                    ghf_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =        &
     4165                   ghf_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =        &
    40104166                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    40114167
    4012                 CASE ( 'm_liq_eb' )
     4168                CASE ( 'm_liq' )
    40134169                   IF ( k == 1 )  READ ( 13 )  tmp_walltype_1d !tmp_2d
    4014                    m_liq_eb_h%var_1d(1:surf_lsm_h%ns)  =                       &
     4170                   m_liq_h%var_1d(1:surf_lsm_h%ns)  =                       &
    40154171                                 tmp_walltype_1d(1:surf_lsm_h%ns)
    40164172
     
    40234179                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    40244180
    4025                 CASE ( 'm_liq_eb_av' )
    4026                    IF ( .NOT. ALLOCATED( m_liq_eb_av ) )  THEN
    4027                       ALLOCATE( m_liq_eb_av(nysg:nyng,nxlg:nxrg) )
     4181                CASE ( 'm_liq_av' )
     4182                   IF ( .NOT. ALLOCATED( m_liq_av ) )  THEN
     4183                      ALLOCATE( m_liq_av(nysg:nyng,nxlg:nxrg) )
    40284184                   ENDIF
    40294185                   IF ( k == 1 )  READ ( 13 )  tmp_2d
    4030                    m_liq_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =      &
     4186                   m_liq_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =      &
    40314187                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    40324188
     
    40454201                                    -nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    40464202
    4047                 CASE ( 'qsws_eb_av' )
    4048                    IF ( .NOT. ALLOCATED( qsws_eb_av ) )  THEN
    4049                       ALLOCATE( qsws_eb_av(nysg:nyng,nxlg:nxrg) )
     4203                CASE ( 'qsws_liq_av' )
     4204                   IF ( .NOT. ALLOCATED( qsws_liq_av ) )  THEN
     4205                      ALLOCATE( qsws_liq_av(nysg:nyng,nxlg:nxrg) )
    40504206                   ENDIF 
    40514207                   IF ( k == 1 )  READ ( 13 )  tmp_2d
    4052                    qsws_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =      &
     4208                   qsws_liq_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =  &
    40534209                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    4054 
    4055                 CASE ( 'qsws_liq_eb_av' )
    4056                    IF ( .NOT. ALLOCATED( qsws_liq_eb_av ) )  THEN
    4057                       ALLOCATE( qsws_liq_eb_av(nysg:nyng,nxlg:nxrg) )
     4210                CASE ( 'qsws_soil_av' )
     4211                   IF ( .NOT. ALLOCATED( qsws_soil_av ) )  THEN
     4212                      ALLOCATE( qsws_soil_av(nysg:nyng,nxlg:nxrg) )
    40584213                   ENDIF 
    40594214                   IF ( k == 1 )  READ ( 13 )  tmp_2d
    4060                    qsws_liq_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
     4215                   qsws_soil_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
    40614216                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    4062                 CASE ( 'qsws_soil_eb_av' )
    4063                    IF ( .NOT. ALLOCATED( qsws_soil_eb_av ) )  THEN
    4064                       ALLOCATE( qsws_soil_eb_av(nysg:nyng,nxlg:nxrg) )
     4217
     4218                CASE ( 'qsws_veg_av' )
     4219                   IF ( .NOT. ALLOCATED( qsws_veg_av ) )  THEN
     4220                      ALLOCATE( qsws_veg_av(nysg:nyng,nxlg:nxrg) )
    40654221                   ENDIF 
    40664222                   IF ( k == 1 )  READ ( 13 )  tmp_2d
    4067                    qsws_soil_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
     4223                   qsws_veg_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
    40684224                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    4069 
    4070                 CASE ( 'qsws_veg_eb_av' )
    4071                    IF ( .NOT. ALLOCATED( qsws_veg_eb_av ) )  THEN
    4072                       ALLOCATE( qsws_veg_eb_av(nysg:nyng,nxlg:nxrg) )
    4073                    ENDIF 
    4074                    IF ( k == 1 )  READ ( 13 )  tmp_2d
    4075                    qsws_veg_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =  &
    4076                                           tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    4077 
    4078                 CASE ( 'shf_eb_av' )
    4079                    IF ( .NOT. ALLOCATED( shf_eb_av ) )  THEN
    4080                       ALLOCATE( shf_eb_av(nysg:nyng,nxlg:nxrg) )
    4081                    ENDIF
    4082                    IF ( k == 1 )  READ ( 13 )  tmp_2d
    4083                    shf_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =       &
    4084                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    40854225
    40864226                CASE ( 't_soil' )
Note: See TracChangeset for help on using the changeset viewer.