Ignore:
Timestamp:
May 30, 2017 5:47:52 PM (4 years ago)
Author:
suehring
Message:

Adjustments according new topography and surface-modelling concept implemented

File:
1 edited

Legend:

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

    r2150 r2232  
    2020! Current revisions:
    2121! -----------------
    22 !
     22! Adjustments to new topography and surface concept
     23!   - now, also vertical walls are possible
     24!   - for vertical walls, parametrization of r_a (aerodynamic resisistance) is
     25!     implemented.
     26!
     27! Add check for soil moisture, it must not exceed its saturation value.
    2328!
    2429! Former revisions:
     
    159164!> DALES and UCLA-LES models.
    160165!>
     166!> @todo Extensive verification energy-balance solver for vertical surfaces,
     167!>       e.g. parametrization of r_a
     168!> @todo Revise single land-surface processes for vertical surfaces, e.g.
     169!>       treatment of humidity, etc.
    161170!> @todo Consider partial absorption of the net shortwave radiation by the
    162171!>       skin layer.
     
    174183 
    175184    USE arrays_3d,                                                             &
    176         ONLY:  hyp, ol, pt, pt_p, prr, q, q_p, ql, qsws, shf, ts, us, vpt, z0, &
    177                z0h, z0q
     185        ONLY:  hyp, pt, pt_p, prr, q, q_p, ql, vpt, u, v, w
    178186
    179187    USE cloud_parameters,                                                      &
     
    183191        ONLY:  cloud_physics, dt_3d, humidity, intermediate_timestep_count,    &
    184192               initializing_actions, intermediate_timestep_count_max,          &
    185                max_masks, precipitation, pt_surface,                           &
     193               land_surface, max_masks, precipitation, pt_surface,             &
    186194               rho_surface, roughness_length, surface_pressure,                &
    187195               timestep_scheme, tsc, z0h_factor, time_since_reference_point
    188196
    189197    USE indices,                                                               &
    190         ONLY:  nbgp, nxlg, nxrg, nyng, nysg, nzb, nzb_s_inner
     198        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb
    191199
    192200    USE kinds
     
    196204    USE radiation_model_mod,                                                   &
    197205        ONLY:  force_radiation_call, rad_net, rad_sw_in, rad_lw_out,           &
    198                rad_lw_out_change_0, unscheduled_radiation_calls
     206               rad_lw_out_change_0, radiation_scheme, unscheduled_radiation_calls
    199207       
    200208    USE statistics,                                                            &
    201209        ONLY:  hom, statistic_regions
    202210
     211    USE surface_mod,                                                        &
     212        ONLY :  surf_lsm_h, surf_lsm_v, surf_type
     213
    203214    IMPLICIT NONE
    204215
     216    TYPE surf_type_lsm
     217       REAL(wp), DIMENSION(:),   ALLOCATABLE ::  var_1D !< 1D prognostic variable
     218       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  var_2D !< 2D prognostic variable
     219    END TYPE surf_type_lsm
    205220!
    206221!-- LSM model constants
     
    225240                    soil_type = 3    !< NAMELIST soil_type_2d
    226241
    227     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  soil_type_2d, &  !< soil type, 0: user-defined, 1-7: generic (see list)
    228                                                   veg_type_2d      !< vegetation type, 0: user-defined, 1-19: generic (see list)
    229 
    230     LOGICAL, DIMENSION(:,:), ALLOCATABLE :: water_surface,     & !< flag parameter for water surfaces (classes 14+15)
    231                                             pave_surface,      & !< flag parameter for pavements (asphalt etc.) (class 20)
    232                                             building_surface     !< flag parameter indicating that the surface element is covered by buildings (no LSM actions, not implemented yet)
    233 
    234242    LOGICAL :: conserve_water_content = .TRUE.,  & !< open or closed bottom surface for the soil model
    235                force_radiation_call_l = .FALSE., & !< flag parameter for unscheduled radiation model calls
    236                land_surface = .FALSE.              !< flag parameter indicating wheather the lsm is used
     243               force_radiation_call_l = .FALSE., & !< flag to force calling of radiation routine
     244               aero_resist_kray = .TRUE.           !< flag to control parametrization of aerodynamic resistance at vertical surface elements
    237245
    238246!   value 9999999.9_wp -> generic available or user-defined value must be set
     
    288296              dz_soil                                                            !< soil grid spacing (edge-edge)
    289297
     298
    290299#if defined( __nopointer )
    291     REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_surface,   & !< surface temperature (K)
    292                                                      t_surface_p, & !< progn. surface temperature (K)
    293                                                      m_liq_eb,    & !< liquid water reservoir (m)
    294                                                      m_liq_eb_av, & !< liquid water reservoir (m)
    295                                                      m_liq_eb_p     !< progn. liquid water reservoir (m)
     300    TYPE(surf_type_lsm), TARGET  ::  t_soil_h,    & !< Soil temperature (K), horizontal surface elements
     301                                     t_soil_h_p,  & !< Prog. soil temperature (K), horizontal surface elements
     302                                     m_soil_h,    & !< Soil moisture (m3/m3), horizontal surface elements
     303                                     m_soil_h_p     !< Prog. soil moisture (m3/m3), horizontal surface elements
     304
     305    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET  ::  &
     306                                     t_soil_v,    & !< Soil temperature (K), vertical surface elements
     307                                     t_soil_v_p,  & !< Prog. soil temperature (K), vertical surface elements
     308                                     m_soil_v,    & !< Soil moisture (m3/m3), vertical surface elements
     309                                     m_soil_v_p     !< Prog. soil moisture (m3/m3), vertical surface elements
    296310#else
    297     REAL(wp), DIMENSION(:,:), POINTER :: t_surface,      &
    298                                          t_surface_p,    &
    299                                          m_liq_eb,       &
    300                                          m_liq_eb_p
    301 
    302     REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: t_surface_1, t_surface_2, &
    303                                                      m_liq_eb_av,              &
    304                                                      m_liq_eb_1, m_liq_eb_2
     311    TYPE(surf_type_lsm), POINTER ::  t_soil_h,    & !< Soil temperature (K), horizontal surface elements
     312                                     t_soil_h_p,  & !< Prog. soil temperature (K), horizontal surface elements
     313                                     m_soil_h,    & !< Soil moisture (m3/m3), horizontal surface elements
     314                                     m_soil_h_p     !< Prog. soil moisture (m3/m3), horizontal surface elements
     315
     316    TYPE(surf_type_lsm), TARGET  ::  t_soil_h_1,  & !<
     317                                     t_soil_h_2,  & !<
     318                                     m_soil_h_1,  & !<
     319                                     m_soil_h_2     !<
     320
     321    TYPE(surf_type_lsm), DIMENSION(:), POINTER :: &
     322                                     t_soil_v,    & !< Soil temperature (K), vertical surface elements
     323                                     t_soil_v_p,  & !< Prog. soil temperature (K), vertical surface elements
     324                                     m_soil_v,    & !< Soil moisture (m3/m3), vertical surface elements
     325                                     m_soil_v_p     !< Prog. soil moisture (m3/m3), vertical surface elements   
     326
     327    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::&
     328                                     t_soil_v_1,  & !<
     329                                     t_soil_v_2,  & !<
     330                                     m_soil_v_1,  & !<
     331                                     m_soil_v_2     !<
     332#endif   
     333
     334#if defined( __nopointer )
     335    TYPE(surf_type_lsm), TARGET   ::  t_surface_h,    & !< surface temperature (K), horizontal surface elements
     336                                      t_surface_h_p,  & !< progn. surface temperature (K), horizontal surface elements
     337                                      m_liq_eb_h,     & !< liquid water reservoir (m), horizontal surface elements
     338                                      m_liq_eb_h_p      !< progn. liquid water reservoir (m), horizontal surface elements
     339
     340    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET   ::  &
     341                                      t_surface_v,    & !< surface temperature (K), vertical surface elements
     342                                      t_surface_v_p,  & !< progn. surface temperature (K), vertical surface elements
     343                                      m_liq_eb_v,     & !< liquid water reservoir (m), vertical surface elements
     344                                      m_liq_eb_v_p      !< progn. liquid water reservoir (m), vertical surface elements
     345#else
     346    TYPE(surf_type_lsm), POINTER  ::  t_surface_h,    & !< surface temperature (K), horizontal surface elements
     347                                      t_surface_h_p,  & !< progn. surface temperature (K), horizontal surface elements
     348                                      m_liq_eb_h,     & !< liquid water reservoir (m), horizontal surface elements
     349                                      m_liq_eb_h_p      !< progn. liquid water reservoir (m), horizontal surface elements
     350
     351    TYPE(surf_type_lsm), TARGET   ::  t_surface_h_1,  & !<
     352                                      t_surface_h_2,  & !<
     353                                      m_liq_eb_h_1,   & !<
     354                                      m_liq_eb_h_2      !<
     355
     356    TYPE(surf_type_lsm), DIMENSION(:), POINTER  ::    &
     357                                      t_surface_v,    & !< surface temperature (K), vertical surface elements
     358                                      t_surface_v_p,  & !< progn. surface temperature (K), vertical surface elements
     359                                      m_liq_eb_v,     & !< liquid water reservoir (m), vertical surface elements
     360                                      m_liq_eb_v_p      !< progn. liquid water reservoir (m), vertical surface elements
     361
     362    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET   ::  &
     363                                      t_surface_v_1,  & !<
     364                                      t_surface_v_2,  & !<
     365                                      m_liq_eb_v_1,   & !<
     366                                      m_liq_eb_v_2      !<
    305367#endif
    306368
    307 !
    308 !-- Temporal tendencies for time stepping           
    309     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: tt_surface_m,  & !< surface temperature tendency (K)
    310                                              tm_liq_eb_m      !< liquid water reservoir tendency (m)
     369#if defined( __nopointer )
     370    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: m_liq_eb_av
     371#else
     372    REAL(wp), DIMENSION(:,:), ALLOCATABLE, TARGET :: m_liq_eb_av
     373#endif
     374
     375#if defined( __nopointer )
     376    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  t_soil_av, & !< Average of t_soil
     377                                                        m_soil_av    !< Average of m_soil
     378#else
     379    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::  t_soil_av, & !< Average of t_soil
     380                                                        m_soil_av    !< Average of m_soil
     381#endif
     382
     383    TYPE(surf_type_lsm), TARGET ::  tm_liq_eb_h_m   !< liquid water reservoir tendency (m), horizontal surface elements
     384    TYPE(surf_type_lsm), TARGET ::  tt_surface_h_m  !< surface temperature tendency (K), horizontal surface elements
     385    TYPE(surf_type_lsm), TARGET ::  tt_soil_h_m     !< t_soil storage array, horizontal surface elements
     386    TYPE(surf_type_lsm), TARGET ::  tm_soil_h_m     !< m_soil storage array, horizontal surface elements
     387
     388    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tm_liq_eb_v_m   !< liquid water reservoir tendency (m), vertical surface elements
     389    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tt_surface_v_m  !< surface temperature tendency (K), vertical surface elements
     390    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tt_soil_v_m     !< t_soil storage array, vertical surface elements
     391    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET ::  tm_soil_v_m     !< m_soil storage array, vertical surface elements
    311392
    312393!
    313394!-- Energy balance variables               
    314395    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
    315               alpha_vg,         & !< coef. of Van Genuchten
    316               c_liq,            & !< liquid water coverage (of vegetated area)
    317396              c_liq_av,         & !< average of c_liq
    318397              c_soil_av,        & !< average of c_soil
    319               c_veg,            & !< vegetation coverage
    320398              c_veg_av,         & !< average of c_veg
    321               f_sw_in,          & !< fraction of absorbed shortwave radiation by the surface layer (not implemented yet)
    322               ghf_eb,           & !< ground heat flux
    323399              ghf_eb_av,        & !< average of ghf_eb
    324               gamma_w_sat,      & !< hydraulic conductivity at saturation
    325               g_d,              & !< coefficient for dependence of r_canopy on water vapour pressure deficit
    326               lai,              & !< leaf area index
    327400              lai_av,           & !< average of lai
    328               lambda_surface_s, & !< coupling between surface and soil (depends on vegetation type)
    329               lambda_surface_u, & !< coupling between surface and soil (depends on vegetation type)
    330               l_vg,             & !< coef. of Van Genuchten
    331               m_fc,             & !< soil moisture at field capacity (m3/m3)
    332               m_res,            & !< residual soil moisture
    333               m_sat,            & !< saturation soil moisture (m3/m3)
    334               m_wilt,           & !< soil moisture at permanent wilting point (m3/m3)
    335               n_vg,             & !< coef. Van Genuchten 
    336               qsws_eb,          & !< surface flux of latent heat (total)
    337401              qsws_eb_av,       & !< average of qsws_eb
    338               qsws_liq_eb,      & !< surface flux of latent heat (liquid water portion)
    339402              qsws_liq_eb_av,   & !< average of qsws_liq_eb
    340               qsws_soil_eb,     & !< surface flux of latent heat (soil portion)
    341403              qsws_soil_eb_av,  & !< average of qsws_soil_eb
    342               qsws_veg_eb,      & !< surface flux of latent heat (vegetation portion)
    343404              qsws_veg_eb_av,   & !< average of qsws_veg_eb
    344               rad_net_l,        & !< local copy of rad_net (net radiation at surface)
    345               r_a,              & !< aerodynamic resistance
    346405              r_a_av,           & !< average of r_a
    347               r_canopy,         & !< canopy resistance
    348               r_soil,           & !< soil resistance
    349               r_soil_min,       & !< minimum soil resistance
    350               r_s,              & !< total surface resistance (combination of r_soil and r_canopy)
    351406              r_s_av,           & !< average of r_s
    352               r_canopy_min,     & !< minimum canopy (stomatal) resistance
    353               shf_eb,           & !< surface flux of sensible heat
    354407              shf_eb_av           !< average of shf_eb
    355 
    356 
    357     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
    358               lambda_h, &   !< heat conductivity of soil (W/m/K)                           
    359               lambda_w, &   !< hydraulic diffusivity of soil (?)
    360               gamma_w,  &   !< hydraulic conductivity of soil (W/m/K)
    361               rho_c_total   !< volumetric heat capacity of the actual soil matrix (?)
    362 
    363 #if defined( __nopointer )
    364     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
    365               t_soil,    & !< Soil temperature (K)
    366               t_soil_av, & !< Average of t_soil
    367               t_soil_p,  & !< Prog. soil temperature (K)
    368               m_soil,    & !< Soil moisture (m3/m3)
    369               m_soil_av, & !< Average of m_soil
    370               m_soil_p     !< Prog. soil moisture (m3/m3)
    371 #else
    372     REAL(wp), DIMENSION(:,:,:), POINTER ::                                     &
    373               t_soil, t_soil_p, &
    374               m_soil, m_soil_p   
    375 
    376     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                         &
    377               t_soil_av, t_soil_1, t_soil_2,                                   &
    378               m_soil_av, m_soil_1, m_soil_2
    379 #endif
    380 
    381 
    382     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
    383               tt_soil_m, & !< t_soil storage array
    384               tm_soil_m, & !< m_soil storage array
    385               root_fr      !< root fraction (sum=1)
    386408
    387409
     
    571593!
    572594!-- Public parameters, constants and initial values
    573     PUBLIC land_surface, skip_time_do_lsm
     595    PUBLIC aero_resist_kray, skip_time_do_lsm
    574596
    575597!
     
    578600
    579601!
    580 !-- Public 2D output variables
    581     PUBLIC ghf_eb, qsws_eb, qsws_liq_eb, qsws_soil_eb,qsws_veg_eb, r_a, r_s,   &
    582            shf_eb
    583 
    584 !
    585602!-- Public prognostic variables
    586     PUBLIC m_soil, t_soil
     603    PUBLIC m_soil_h, t_soil_h
    587604
    588605
     
    871888    USE control_parameters,                                                    &
    872889        ONLY:  bc_pt_b, bc_q_b, constant_flux_layer, message_string,           &
    873                most_method, topography
     890               most_method
    874891                 
    875892    USE radiation_model_mod,                                                   &
     
    879896    IMPLICIT NONE
    880897
    881  
     898    INTEGER(iwp) ::  k        !< running index, z-dimension
    882899!
    883900!-- Dirichlet boundary conditions are required as the surface fluxes are
     
    897914    ENDIF
    898915
    899     IF ( topography /= 'flat' )  THEN
    900        message_string = 'lsm cannot be used ' //                            &
    901                         'in combination with  topography /= "flat"'
    902        CALL message( 'check_parameters', 'PA0415', 1, 2, 0, 6, 0 )
    903     ENDIF
    904 
    905916    IF ( ( veg_type == 14  .OR.  veg_type == 15 ) .AND.                     &
    906917           most_method == 'lookup' )  THEN
     
    10441055
    10451056    ENDIF
     1057!
     1058!-- Check for proper setting of soil moisture, must not be larger than its
     1059!-- saturation value.
     1060    DO  k = nzb_soil, nzt_soil
     1061       IF ( soil_moisture(k) > m_soil_pars(0,soil_type) )  THEN
     1062          message_string = 'soil_moisture must not exceed its saturation' // &
     1063                           ' value'
     1064          CALL message( 'check_parameters', 'PA0403', 1, 2, 0, 6, 0 )
     1065       ENDIF
     1066    ENDDO
    10461067
    10471068    IF (  .NOT.  radiation )  THEN
     
    10591080!> Solver for the energy balance at the surface.
    10601081!------------------------------------------------------------------------------!
    1061  SUBROUTINE lsm_energy_balance
     1082 SUBROUTINE lsm_energy_balance( horizontal, l )
    10621083
    10631084
     
    10651086
    10661087    INTEGER(iwp) ::  i         !< running index
     1088    INTEGER(iwp) ::  i_off     !< offset to determine index of surface element, seen from atmospheric grid point, for x
    10671089    INTEGER(iwp) ::  j         !< running index
    1068     INTEGER(iwp) ::  k, ks     !< running index
     1090    INTEGER(iwp) ::  j_off     !< offset to determine index of surface element, seen from atmospheric grid point, for y
     1091    INTEGER(iwp) ::  k         !< running index
     1092    INTEGER(iwp) ::  k_off     !< offset to determine index of surface element, seen from atmospheric grid point, for z
     1093    INTEGER(iwp) ::  k_rad     !< index to access radiation array
     1094    INTEGER(iwp) ::  ks        !< running index
     1095    INTEGER(iwp) ::  l         !< surface-facing index
     1096    INTEGER(iwp) ::  m         !< running index concerning wall elements
     1097
     1098    LOGICAL      ::  horizontal !< Flag indicating horizontal or vertical surfaces
    10691099
    10701100    REAL(wp) :: c_surface_tmp,& !< temporary variable for storing the volumetric heat capacity of the surface
     
    10901120                qv1            !< specific humidity at first grid level
    10911121
     1122    TYPE(surf_type_lsm), POINTER ::  surf_t_surface
     1123    TYPE(surf_type_lsm), POINTER ::  surf_t_surface_p
     1124    TYPE(surf_type_lsm), POINTER ::  surf_tt_surface_m
     1125    TYPE(surf_type_lsm), POINTER ::  surf_m_liq_eb
     1126    TYPE(surf_type_lsm), POINTER ::  surf_m_liq_eb_p
     1127    TYPE(surf_type_lsm), POINTER ::  surf_tm_liq_eb_m
     1128
     1129    TYPE(surf_type_lsm), POINTER ::  surf_m_soil
     1130    TYPE(surf_type_lsm), POINTER ::  surf_t_soil
     1131
     1132    TYPE(surf_type), POINTER  ::  surf  !< surface-date type variable
     1133
     1134    IF ( horizontal )  THEN
     1135       surf              => surf_lsm_h
     1136
     1137       surf_t_surface    => t_surface_h
     1138       surf_t_surface_p  => t_surface_h_p
     1139       surf_tt_surface_m => tt_surface_h_m
     1140       surf_m_liq_eb     => m_liq_eb_h
     1141       surf_m_liq_eb_p   => m_liq_eb_h_p
     1142       surf_tm_liq_eb_m  => tm_liq_eb_h_m
     1143       surf_m_soil       => m_soil_h
     1144       surf_t_soil       => t_soil_h
     1145
     1146       k_off     = -1
     1147       j_off     = 0
     1148       i_off     = 0
     1149    ELSE
     1150       surf              => surf_lsm_v(l)
     1151
     1152       surf_t_surface    => t_surface_v(l)
     1153       surf_t_surface_p  => t_surface_v_p(l)
     1154       surf_tt_surface_m => tt_surface_v_m(l)
     1155       surf_m_liq_eb     => m_liq_eb_v(l)
     1156       surf_m_liq_eb_p   => m_liq_eb_v_p(l)
     1157       surf_tm_liq_eb_m  => tm_liq_eb_v_m(l)
     1158       surf_m_soil       => m_soil_v(l)
     1159       surf_t_soil       => t_soil_v(l)
     1160
     1161       k_off = 0
     1162       IF ( l == 0 )  THEN
     1163          j_off = -1
     1164          i_off = 0
     1165       ELSEIF ( l == 1 )  THEN
     1166          j_off = 1
     1167          i_off = 0
     1168       ELSEIF ( l == 2 )  THEN
     1169          j_off = 0
     1170          i_off = -1
     1171       ELSEIF ( l == 3 )  THEN
     1172          j_off = 0
     1173          i_off = 1
     1174       ENDIF
     1175    ENDIF
     1176
    10921177!
    10931178!-- Calculate the exner function for the current time step
    10941179    exn = ( surface_pressure / 1000.0_wp )**0.286_wp
    10951180
    1096     DO  i = nxlg, nxrg
    1097        DO  j = nysg, nyng
    1098           k = nzb_s_inner(j,i)
    1099 
    1100 !
    1101 !--       Set lambda_surface according to stratification between skin layer and soil
    1102           IF (  .NOT.  pave_surface(j,i) )  THEN
    1103 
    1104              c_surface_tmp = c_surface
    1105 
    1106              IF ( t_surface(j,i) >= t_soil(nzb_soil,j,i))  THEN
    1107                 lambda_surface = lambda_surface_s(j,i)
    1108              ELSE
    1109                 lambda_surface = lambda_surface_u(j,i)
    1110              ENDIF
     1181    DO  m = 1, surf%ns
     1182
     1183       i   = surf%i(m)           
     1184       j   = surf%j(m)
     1185       k   = surf%k(m)
     1186!
     1187!--    Determine height index for radiation. Note, in clear-sky case radiation
     1188!--    arrays have rank 0 in first dimensions, so index must be zero. In case
     1189!--    of RRTMG radiation arrays have non-zero rank in first dimension, so that
     1190!--    radiation can be obtained at respective height level
     1191       k_rad = MERGE( 0, k + k_off, radiation_scheme /= 'rrtmg' )
     1192!
     1193!--    Set lambda_surface according to stratification between skin layer and soil
     1194       IF (  .NOT.  surf%pave_surface(m) )  THEN
     1195
     1196          c_surface_tmp = c_surface
     1197
     1198          IF ( surf_t_surface%var_1d(m) >= surf_t_soil%var_2d(nzb_soil,m))  THEN
     1199             lambda_surface = surf%lambda_surface_s(m)
    11111200          ELSE
    1112 
    1113              c_surface_tmp = pave_heat_capacity * dz_soil(nzb_soil) * 0.5_wp
    1114              lambda_surface = pave_heat_conductivity * ddz_soil(nzb_soil)
    1115 
    1116           ENDIF
    1117 
    1118 !
    1119 !--       First step: calculate aerodyamic resistance. As pt, us, ts
    1120 !--       are not available for the prognostic time step, data from the last
    1121 !--       time step is used here. Note that this formulation is the
    1122 !--       equivalent to the ECMWF formulation using drag coefficients
    1123           IF ( cloud_physics )  THEN
    1124              pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
    1125              qv1 = q(k+1,j,i) - ql(k+1,j,i)
    1126           ELSE
    1127              pt1 = pt(k+1,j,i)
    1128              qv1 = q(k+1,j,i)
    1129           ENDIF
    1130 
    1131           r_a(j,i) = (pt1 - pt(k,j,i)) / (ts(j,i) * us(j,i) + 1.0E-20_wp)
    1132 
    1133 !
    1134 !--       Make sure that the resistance does not drop to zero for neutral
    1135 !--       stratification
    1136           IF ( ABS(r_a(j,i)) < 1.0_wp )  r_a(j,i) = 1.0_wp
    1137 
    1138 !
    1139 !--       Second step: calculate canopy resistance r_canopy
    1140 !--       f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
     1201             lambda_surface = surf%lambda_surface_u(m)
     1202          ENDIF
     1203       ELSE
     1204
     1205          c_surface_tmp = pave_heat_capacity * dz_soil(nzb_soil) * 0.5_wp
     1206          lambda_surface = pave_heat_conductivity * ddz_soil(nzb_soil)
     1207
     1208       ENDIF
     1209
     1210!
     1211!--    First step: calculate aerodyamic resistance. As pt, us, ts
     1212!--    are not available for the prognostic time step, data from the last
     1213!--    time step is used here. Note that this formulation is the
     1214!--    equivalent to the ECMWF formulation using drag coefficients
     1215       IF ( cloud_physics )  THEN
     1216          pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
     1217          qv1 = q(k,j,i) - ql(k,j,i)
     1218       ELSE
     1219          pt1 = pt(k,j,i)
     1220          qv1 = q(k,j,i)
     1221       ENDIF
     1222!
     1223!--    Calculate aerodynamical resistance. For horizontal and vertical
     1224!--    surfaces MOST is applied. Moreover, for vertical surfaces, resistance
     1225!--    can be obtain via parameterization of Mason (2000) /
     1226!--    Krayenhoff and Voogt (2006).
     1227!--    To do: detailed investigation which approach is better!
     1228       IF ( horizontal  .OR.  .NOT. aero_resist_kray )  THEN
     1229          surf%r_a(m) = ( pt1 - surf_lsm_h%pt_surface(m) ) /                   &
     1230                        ( surf%ts(m) * surf%us(m) + 1.0E-20_wp )
     1231       ELSE
     1232          surf%r_a(m) = 1.0_wp / ( 11.8_wp + 4.2_wp *                          &
     1233                        SQRT( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 +      &
     1234                              ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 +      &
     1235                              ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2 )      &
     1236                                 )
     1237       ENDIF
     1238!
     1239!--    Make sure that the resistance does not drop to zero for neutral
     1240!--    stratification
     1241       IF ( ABS( surf%r_a(m) ) < 1.0_wp )  surf%r_a(m) = 1.0_wp
     1242!
     1243!--    Second step: calculate canopy resistance r_canopy
     1244!--    f1-f3 here are defined as 1/f1-f3 as in ECMWF documentation
    11411245 
    1142 !--       f1: correction for incoming shortwave radiation (stomata close at
    1143 !--       night)
    1144           f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k,j,i) + 0.05_wp ) /        &
    1145                            (0.81_wp * (0.004_wp * rad_sw_in(k,j,i)             &
    1146                             + 1.0_wp)) )
    1147 
    1148 
    1149 
    1150 !
    1151 !--       f2: correction for soil moisture availability to plants (the
    1152 !--       integrated soil moisture must thus be considered here)
    1153 !--       f2 = 0 for very dry soils
    1154           m_total = 0.0_wp
    1155           DO  ks = nzb_soil, nzt_soil
    1156               m_total = m_total + root_fr(ks,j,i)                              &
    1157                         * MAX(m_soil(ks,j,i),m_wilt(j,i))
    1158           ENDDO
    1159 
    1160           IF ( m_total > m_wilt(j,i)  .AND.  m_total < m_fc(j,i) )  THEN
    1161              f2 = ( m_total - m_wilt(j,i) ) / (m_fc(j,i) - m_wilt(j,i) )
    1162           ELSEIF ( m_total >= m_fc(j,i) )  THEN
    1163              f2 = 1.0_wp
    1164           ELSE
    1165              f2 = 1.0E-20_wp
    1166           ENDIF
    1167 
    1168 !
    1169 !--       Calculate water vapour pressure at saturation
    1170           e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface(j,i)        &
    1171                         - 273.16_wp ) / ( t_surface(j,i) - 35.86_wp ) )
    1172 
    1173 !
    1174 !--       f3: correction for vapour pressure deficit
    1175           IF ( g_d(j,i) /= 0.0_wp )  THEN
    1176 !
    1177 !--          Calculate vapour pressure
    1178              e  = qv1 * surface_pressure / 0.622_wp
    1179              f3 = EXP ( -g_d(j,i) * (e_s - e) )
    1180           ELSE
    1181              f3 = 1.0_wp
    1182           ENDIF
    1183 
    1184 !
    1185 !--       Calculate canopy resistance. In case that c_veg is 0 (bare soils),
    1186 !--       this calculation is obsolete, as r_canopy is not used below.
    1187 !--       To do: check for very dry soil -> r_canopy goes to infinity
    1188           r_canopy(j,i) = r_canopy_min(j,i) / (lai(j,i) * f1 * f2 * f3         &
    1189                                           + 1.0E-20_wp)
    1190 
    1191 !
    1192 !--       Third step: calculate bare soil resistance r_soil. The Clapp &
    1193 !--       Hornberger parametrization does not consider c_veg.
    1194           IF ( soil_type_2d(j,i) /= 7 )  THEN
    1195              m_min = c_veg(j,i) * m_wilt(j,i) + (1.0_wp - c_veg(j,i)) *        &
    1196                      m_res(j,i)
    1197           ELSE
    1198              m_min = m_wilt(j,i)
    1199           ENDIF
    1200 
    1201           f2 = ( m_soil(nzb_soil,j,i) - m_min ) / ( m_fc(j,i) - m_min )
    1202           f2 = MAX(f2,1.0E-20_wp)
    1203           f2 = MIN(f2,1.0_wp)
    1204 
    1205           r_soil(j,i) = r_soil_min(j,i) / f2
    1206 
    1207 !
    1208 !--       Calculate the maximum possible liquid water amount on plants and
    1209 !--       bare surface. For vegetated surfaces, a maximum depth of 0.2 mm is
    1210 !--       assumed, while paved surfaces might hold up 1 mm of water. The
    1211 !--       liquid water fraction for paved surfaces is calculated after
    1212 !--       Noilhan & Planton (1989), while the ECMWF formulation is used for
    1213 !--       vegetated surfaces and bare soils.
    1214           IF ( pave_surface(j,i) )  THEN
    1215              m_liq_eb_max = m_max_depth * 5.0_wp
    1216              c_liq(j,i) = MIN( 1.0_wp, (m_liq_eb(j,i) / m_liq_eb_max)**0.67 )
    1217           ELSE
    1218              m_liq_eb_max = m_max_depth * ( c_veg(j,i) * lai(j,i)              &
    1219                             + (1.0_wp - c_veg(j,i)) )
    1220              c_liq(j,i) = MIN( 1.0_wp, m_liq_eb(j,i) / m_liq_eb_max )
    1221           ENDIF
    1222 
    1223 !
    1224 !--       Calculate saturation specific humidity
    1225           q_s = 0.622_wp * e_s / surface_pressure
    1226 
    1227 !
    1228 !--       In case of dewfall, set evapotranspiration to zero
    1229 !--       All super-saturated water is then removed from the air
    1230           IF ( humidity  .AND.  q_s <= qv1 )  THEN
    1231              r_canopy(j,i) = 0.0_wp
    1232              r_soil(j,i)   = 0.0_wp
    1233           ENDIF
    1234 
    1235 !
    1236 !--       Calculate coefficients for the total evapotranspiration
    1237 !--       In case of water surface, set vegetation and soil fluxes to zero.
    1238 !--       For pavements, only evaporation of liquid water is possible.
    1239           IF ( water_surface(j,i) )  THEN
    1240              f_qsws_veg  = 0.0_wp
    1241              f_qsws_soil = 0.0_wp
    1242              f_qsws_liq  = rho_lv / r_a(j,i)
    1243           ELSEIF ( pave_surface (j,i) )  THEN
    1244              f_qsws_veg  = 0.0_wp
    1245              f_qsws_soil = 0.0_wp
    1246              f_qsws_liq  = rho_lv * c_liq(j,i) / r_a(j,i)
    1247           ELSE
    1248              f_qsws_veg  = rho_lv * c_veg(j,i) * (1.0_wp - c_liq(j,i))/        &
    1249                            (r_a(j,i) + r_canopy(j,i))
    1250              f_qsws_soil = rho_lv * (1.0_wp - c_veg(j,i)) / (r_a(j,i) +        &
    1251                                                              r_soil(j,i))
    1252              f_qsws_liq  = rho_lv * c_veg(j,i) * c_liq(j,i) / r_a(j,i)
    1253           ENDIF
     1246!--    f1: correction for incoming shortwave radiation (stomata close at
     1247!--    night)
     1248       f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k_rad,j+j_off,i+i_off) + 0.05_wp ) /&
     1249                        (0.81_wp * (0.004_wp * rad_sw_in(k_rad,j+j_off,i+i_off)&
     1250                         + 1.0_wp)) )
     1251!
     1252!--    f2: correction for soil moisture availability to plants (the
     1253!--    integrated soil moisture must thus be considered here)
     1254!--    f2 = 0 for very dry soils
     1255       m_total = 0.0_wp
     1256       DO  ks = nzb_soil, nzt_soil
     1257           m_total = m_total + surf%root_fr(ks,m)                              &
     1258                     * MAX( surf_m_soil%var_2d(ks,m), surf%m_wilt(m) )
     1259       ENDDO
     1260
     1261       IF ( m_total > surf%m_wilt(m)  .AND.                                    &
     1262            m_total < surf%m_fc(m) )  THEN
     1263          f2 = ( m_total - surf%m_wilt(m) ) /                                  &
     1264               ( surf%m_fc(m) - surf%m_wilt(m) )
     1265       ELSEIF ( m_total >= surf%m_fc(m) )  THEN
     1266          f2 = 1.0_wp
     1267       ELSE
     1268          f2 = 1.0E-20_wp
     1269       ENDIF
     1270!
     1271!--    Calculate water vapour pressure at saturation
     1272       e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( surf_t_surface%var_1d(m) &
     1273                     - 273.16_wp ) / ( surf_t_surface%var_1d(m) - 35.86_wp ) )
     1274!
     1275!--    f3: correction for vapour pressure deficit
     1276       IF ( surf%g_d(m) /= 0.0_wp )  THEN
     1277!
     1278!--       Calculate vapour pressure
     1279          e  = qv1 * surface_pressure / 0.622_wp
     1280          f3 = EXP ( - surf%g_d(m) * (e_s - e) )
     1281       ELSE
     1282          f3 = 1.0_wp
     1283       ENDIF
     1284!
     1285!--    Calculate canopy resistance. In case that c_veg is 0 (bare soils),
     1286!--    this calculation is obsolete, as r_canopy is not used below.
     1287!--    To do: check for very dry soil -> r_canopy goes to infinity
     1288       surf%r_canopy(m) = surf%r_canopy_min(m) /                               &
     1289                              ( surf%lai(m) * f1 * f2 * f3 + 1.0E-20_wp )
     1290!
     1291!--    Third step: calculate bare soil resistance r_soil. The Clapp &
     1292!--    Hornberger parametrization does not consider c_veg.
     1293       IF ( surf%soil_type_2d(m) /= 7 )  THEN
     1294          m_min = surf%c_veg(m) * surf%m_wilt(m) +                             &
     1295                         ( 1.0_wp - surf%c_veg(m) ) * surf%m_res(m)
     1296       ELSE
     1297          m_min = surf%m_wilt(m)
     1298       ENDIF
     1299
     1300       f2 = ( surf_m_soil%var_2d(nzb_soil,m) - m_min ) / ( surf%m_fc(m) - m_min )
     1301       f2 = MAX( f2, 1.0E-20_wp )
     1302       f2 = MIN( f2, 1.0_wp     )
     1303
     1304       surf%r_soil(m) = surf%r_soil_min(m) / f2
     1305
     1306!
     1307!--    Calculate the maximum possible liquid water amount on plants and
     1308!--    bare surface. For vegetated surfaces, a maximum depth of 0.2 mm is
     1309!--    assumed, while paved surfaces might hold up 1 mm of water. The
     1310!--    liquid water fraction for paved surfaces is calculated after
     1311!--    Noilhan & Planton (1989), while the ECMWF formulation is used for
     1312!--    vegetated surfaces and bare soils.
     1313       IF ( surf%pave_surface(m) )  THEN
     1314          m_liq_eb_max = m_max_depth * 5.0_wp
     1315          surf%c_liq(m) = MIN( 1.0_wp, ( surf_m_liq_eb%var_1d(m) / m_liq_eb_max)**0.67 )
     1316       ELSE
     1317          m_liq_eb_max = m_max_depth * ( surf%c_veg(m) * surf%lai(m)&
     1318                            + ( 1.0_wp - surf%c_veg(m) ) )
     1319          surf%c_liq(m) = MIN( 1.0_wp, surf_m_liq_eb%var_1d(m) / m_liq_eb_max )
     1320       ENDIF
     1321!
     1322!--    Calculate saturation specific humidity
     1323       q_s = 0.622_wp * e_s / surface_pressure
     1324!
     1325!--    In case of dewfall, set evapotranspiration to zero
     1326!--    All super-saturated water is then removed from the air
     1327       IF ( humidity  .AND.  q_s <= qv1 )  THEN
     1328          surf%r_canopy(m) = 0.0_wp
     1329          surf%r_soil(m)   = 0.0_wp
     1330       ENDIF
     1331
     1332!
     1333!--    Calculate coefficients for the total evapotranspiration
     1334!--    In case of water surface, set vegetation and soil fluxes to zero.
     1335!--    For pavements, only evaporation of liquid water is possible.
     1336       IF ( surf%water_surface(m) )  THEN
     1337          f_qsws_veg  = 0.0_wp
     1338          f_qsws_soil = 0.0_wp
     1339          f_qsws_liq  = rho_lv / surf%r_a(m)
     1340       ELSEIF ( surf%pave_surface (m) )  THEN
     1341          f_qsws_veg  = 0.0_wp
     1342          f_qsws_soil = 0.0_wp
     1343          f_qsws_liq  = rho_lv * surf%c_liq(m) / surf%r_a(m)
     1344       ELSE
     1345          f_qsws_veg  = rho_lv * surf%c_veg(m) *                               &
     1346                            ( 1.0_wp        - surf%c_liq(m)    ) /             &
     1347                            ( surf%r_a(m) + surf%r_canopy(m) )
     1348          f_qsws_soil = rho_lv * (1.0_wp    - surf%c_veg(m)    ) /             &
     1349                            ( surf%r_a(m) + surf%r_soil(m)   )
     1350          f_qsws_liq  = rho_lv * surf%c_veg(m) * surf%c_liq(m)   /             &
     1351                              surf%r_a(m)
     1352       ENDIF
    12541353!
    12551354!--       If soil moisture is below wilting point, plants do no longer
     
    12591358!           ENDIF
    12601359
    1261           f_shf  = rho_cp / r_a(j,i)
    1262           f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
    1263 
    1264 !
    1265 !--       Calculate derivative of q_s for Taylor series expansion
    1266           e_s_dt = e_s * ( 17.269_wp / (t_surface(j,i) - 35.86_wp) -           &
    1267                            17.269_wp*(t_surface(j,i) - 273.16_wp)              &
    1268                            / (t_surface(j,i) - 35.86_wp)**2 )
    1269 
    1270           dq_s_dt = 0.622_wp * e_s_dt / surface_pressure
    1271 
    1272 !
    1273 !--       Add LW up so that it can be removed in prognostic equation
    1274           rad_net_l(j,i) = rad_net(j,i) + rad_lw_out(nzb,j,i)
    1275 
    1276 !
    1277 !--       Calculate new skin temperature
    1278           IF ( humidity )  THEN
    1279 
    1280 !
    1281 !--          Numerator of the prognostic equation
    1282              coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)                &
    1283                       * t_surface(j,i) - rad_lw_out(nzb,j,i)                   &
    1284                       + f_shf * pt1 + f_qsws * ( qv1 - q_s                     &
    1285                       + dq_s_dt * t_surface(j,i) ) + lambda_surface            &
    1286                       * t_soil(nzb_soil,j,i)
    1287 
    1288 !
    1289 !--          Denominator of the prognostic equation
    1290              coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt              &
    1291                       + lambda_surface + f_shf / exn
    1292           ELSE
    1293 
    1294 !
    1295 !--          Numerator of the prognostic equation
    1296              coef_1 = rad_net_l(j,i) + rad_lw_out_change_0(j,i)                &
    1297                       * t_surface(j,i) - rad_lw_out(nzb,j,i)                   &
    1298                       + f_shf * pt1  + lambda_surface                          &
    1299                       * t_soil(nzb_soil,j,i)
    1300 
    1301 !
    1302 !--          Denominator of the prognostic equation
    1303              coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn
    1304 
    1305           ENDIF
    1306 
    1307           tend = 0.0_wp
    1308 
    1309 !
    1310 !--       Implicit solution when the surface layer has no heat capacity,
    1311 !--       otherwise use RK3 scheme.
    1312           t_surface_p(j,i) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp *       &
    1313                              t_surface(j,i) ) / ( c_surface_tmp + coef_2       &
    1314                                 * dt_3d * tsc(2) )
    1315 
    1316 !
    1317 !--       Add RK3 term
    1318           IF ( c_surface_tmp /= 0.0_wp )  THEN
    1319 
    1320              t_surface_p(j,i) = t_surface_p(j,i) + dt_3d * tsc(3)              &
    1321                                 * tt_surface_m(j,i)
    1322 
    1323 !
    1324 !--          Calculate true tendency
    1325              tend = (t_surface_p(j,i) - t_surface(j,i) - dt_3d * tsc(3)        &
    1326                     * tt_surface_m(j,i)) / (dt_3d  * tsc(2))
    1327 !
    1328 !--          Calculate t_surface tendencies for the next Runge-Kutta step
    1329              IF ( timestep_scheme(1:5) == 'runge' )  THEN
    1330                 IF ( intermediate_timestep_count == 1 )  THEN
    1331                    tt_surface_m(j,i) = tend
    1332                 ELSEIF ( intermediate_timestep_count <                         &
    1333                          intermediate_timestep_count_max )  THEN
    1334                    tt_surface_m(j,i) = -9.5625_wp * tend + 5.3125_wp           &
    1335                                        * tt_surface_m(j,i)
    1336                 ENDIF
     1360       f_shf  = rho_cp / surf%r_a(m)
     1361       f_qsws = f_qsws_veg + f_qsws_soil + f_qsws_liq
     1362!
     1363!--    Calculate derivative of q_s for Taylor series expansion
     1364       e_s_dt = e_s * ( 17.269_wp / ( surf_t_surface%var_1d(m) - 35.86_wp) -   &
     1365                        17.269_wp*( surf_t_surface%var_1d(m) - 273.16_wp)      &
     1366                       / ( surf_t_surface%var_1d(m) - 35.86_wp)**2 )
     1367
     1368       dq_s_dt = 0.622_wp * e_s_dt / surface_pressure
     1369!
     1370!--    Add LW up so that it can be removed in prognostic equation
     1371       surf%rad_net_l(m) = rad_net(j,i) + rad_lw_out(nzb,j,i)
     1372!
     1373!--    Calculate new skin temperature
     1374       IF ( humidity )  THEN
     1375!
     1376!--       Numerator of the prognostic equation
     1377          coef_1 = surf%rad_net_l(m) + rad_lw_out_change_0(j,i)                &
     1378                   * surf_t_surface%var_1d(m) - rad_lw_out(nzb,j,i)            &
     1379                   + f_shf * pt1 + f_qsws * ( qv1 - q_s                        &
     1380                   + dq_s_dt * surf_t_surface%var_1d(m) ) + lambda_surface     &
     1381                   * surf_t_soil%var_2d(nzb_soil,m)
     1382!
     1383!--       Denominator of the prognostic equation
     1384          coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt                 &
     1385                   + lambda_surface + f_shf / exn
     1386       ELSE
     1387!
     1388!--       Numerator of the prognostic equation
     1389          coef_1 = surf%rad_net_l(m) + rad_lw_out_change_0(j,i)                &
     1390                   * surf_t_surface%var_1d(m) - rad_lw_out(nzb,j,i)            &
     1391                   + f_shf * pt1  + lambda_surface                             &
     1392                   * surf_t_soil%var_2d(nzb_soil,m)
     1393!
     1394!--       Denominator of the prognostic equation
     1395          coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn
     1396
     1397       ENDIF
     1398
     1399       tend = 0.0_wp
     1400
     1401!
     1402!--    Implicit solution when the surface layer has no heat capacity,
     1403!--    otherwise use RK3 scheme.
     1404       surf_t_surface_p%var_1d(m) = ( coef_1 * dt_3d * tsc(2) + c_surface_tmp *&
     1405                          surf_t_surface%var_1d(m) ) / ( c_surface_tmp + coef_2&
     1406                                             * dt_3d * tsc(2) )
     1407
     1408!
     1409!--    Add RK3 term
     1410       IF ( c_surface_tmp /= 0.0_wp )  THEN
     1411
     1412          surf_t_surface_p%var_1d(m) = surf_t_surface_p%var_1d(m) + dt_3d *    &
     1413                                       tsc(3) * surf_tt_surface_m%var_1d(m)
     1414
     1415!
     1416!--       Calculate true tendency
     1417          tend = ( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) -     &
     1418                   dt_3d * tsc(3) * surf_tt_surface_m%var_1d(m)) / (dt_3d  * tsc(2))
     1419!
     1420!--       Calculate t_surface tendencies for the next Runge-Kutta step
     1421          IF ( timestep_scheme(1:5) == 'runge' )  THEN
     1422             IF ( intermediate_timestep_count == 1 )  THEN
     1423                surf_tt_surface_m%var_1d(m) = tend
     1424             ELSEIF ( intermediate_timestep_count <                            &
     1425                      intermediate_timestep_count_max )  THEN
     1426                surf_tt_surface_m%var_1d(m) = -9.5625_wp * tend +              &
     1427                                               5.3125_wp * surf_tt_surface_m%var_1d(m)
    13371428             ENDIF
    13381429          ENDIF
    1339 
    1340 !
    1341 !--       In case of fast changes in the skin temperature, it is possible to
    1342 !--       update the radiative fluxes independently from the prescribed
    1343 !--       radiation call frequency. This effectively prevents oscillations,
    1344 !--       especially when setting skip_time_do_radiation /= 0. The threshold
    1345 !--       value of 0.2 used here is just a first guess. This method should be
    1346 !--       revised in the future as tests have shown that the threshold is
    1347 !--       often reached, when no oscillations would occur (causes immense
    1348 !--       computing time for the radiation code).
    1349           IF ( ABS( t_surface_p(j,i) - t_surface(j,i) ) > 0.2_wp  .AND.        &
    1350                unscheduled_radiation_calls )  THEN
    1351              force_radiation_call_l = .TRUE.
    1352           ENDIF
    1353 
    1354           pt(k,j,i) = t_surface_p(j,i) / exn
    1355 
    1356 !
    1357 !--       Calculate fluxes
    1358           rad_net_l(j,i)   = rad_net_l(j,i) + rad_lw_out_change_0(j,i)         &
    1359                              * t_surface(j,i) - rad_lw_out(nzb,j,i)            &
    1360                              - rad_lw_out_change_0(j,i) * t_surface_p(j,i)
    1361 
    1362           rad_net(j,i) = rad_net_l(j,i)
    1363           rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i) + rad_lw_out_change_0(j,i) &
    1364                                 * ( t_surface_p(j,i) - t_surface(j,i) )
    1365 
    1366           ghf_eb(j,i)    = lambda_surface * (t_surface_p(j,i)                  &
    1367                            - t_soil(nzb_soil,j,i))
    1368 
    1369           shf_eb(j,i)    = - f_shf * ( pt1 - pt(k,j,i) )
    1370 
    1371           shf(j,i) = shf_eb(j,i) / rho_cp
    1372 
    1373           IF ( humidity )  THEN
    1374              qsws_eb(j,i)  = - f_qsws    * ( qv1 - q_s + dq_s_dt               &
    1375                              * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )
    1376 
    1377              qsws(j,i) = qsws_eb(j,i) / rho_lv
    1378 
    1379              qsws_veg_eb(j,i)  = - f_qsws_veg  * ( qv1 - q_s                   &
    1380                                  + dq_s_dt * t_surface(j,i) - dq_s_dt          &
    1381                                  * t_surface_p(j,i) )
    1382 
    1383              qsws_soil_eb(j,i) = - f_qsws_soil * ( qv1 - q_s                   &
    1384                                  + dq_s_dt * t_surface(j,i) - dq_s_dt          &
    1385                                  * t_surface_p(j,i) )
    1386 
    1387              qsws_liq_eb(j,i)  = - f_qsws_liq * ( qv1 - q_s                   &
    1388                                  + dq_s_dt * t_surface(j,i) - dq_s_dt          &
    1389                                  * t_surface_p(j,i) )
    1390           ENDIF
    1391 
    1392 !
    1393 !--       Calculate the true surface resistance
    1394           IF ( qsws_eb(j,i) == 0.0_wp )  THEN
    1395              r_s(j,i) = 1.0E10_wp
    1396           ELSE
    1397              r_s(j,i) = - rho_lv * ( qv1 - q_s + dq_s_dt                       &
    1398                         * t_surface(j,i) - dq_s_dt * t_surface_p(j,i) )        &
    1399                         / qsws_eb(j,i) - r_a(j,i)
    1400           ENDIF
    1401 
    1402 !
    1403 !--       Calculate change in liquid water reservoir due to dew fall or
    1404 !--       evaporation of liquid water
    1405           IF ( humidity )  THEN
    1406 !
    1407 !--          If precipitation is activated, add rain water to qsws_liq_eb
    1408 !--          and qsws_soil_eb according the the vegetation coverage.
    1409 !--          precipitation_rate is given in mm.
    1410              IF ( precipitation )  THEN
    1411 
    1412 !
    1413 !--             Add precipitation to liquid water reservoir, if possible.
    1414 !--             Otherwise, add the water to soil. In case of
    1415 !--             pavements, the exceeding water amount is implicitely removed
    1416 !--             as runoff as qsws_soil_eb is then not used in the soil model
    1417                 IF ( m_liq_eb(j,i) /= m_liq_eb_max )  THEN
    1418                    qsws_liq_eb(j,i) = qsws_liq_eb(j,i)                         &
    1419                                     + c_veg(j,i) * prr(k,j,i) * hyrho(k)       &
    1420                                     * 0.001_wp * rho_l * l_v
    1421                 ELSE
    1422                    qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                       &
    1423                                      + c_veg(j,i) * prr(k,j,i) * hyrho(k)      &
    1424                                      * 0.001_wp * rho_l * l_v
    1425                 ENDIF
    1426 
    1427 !--             Add precipitation to bare soil according to the bare soil
    1428 !--             coverage.
    1429                 qsws_soil_eb(j,i) = qsws_soil_eb(j,i) + (1.0_wp                &
    1430                                     - c_veg(j,i)) * prr(k,j,i) * hyrho(k)      &
    1431                                     * 0.001_wp * rho_l * l_v
     1430       ENDIF
     1431
     1432!
     1433!--    In case of fast changes in the skin temperature, it is possible to
     1434!--    update the radiative fluxes independently from the prescribed
     1435!--    radiation call frequency. This effectively prevents oscillations,
     1436!--    especially when setting skip_time_do_radiation /= 0. The threshold
     1437!--    value of 0.2 used here is just a first guess. This method should be
     1438!--    revised in the future as tests have shown that the threshold is
     1439!--    often reached, when no oscillations would occur (causes immense
     1440!--    computing time for the radiation code).
     1441       IF ( ABS( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) ) > 0.2_wp  .AND. &
     1442            unscheduled_radiation_calls )  THEN
     1443          force_radiation_call_l = .TRUE.
     1444       ENDIF
     1445
     1446       pt(k+k_off,j+j_off,i+i_off) = surf_t_surface_p%var_1d(m) / exn  !is actually no air temperature
     1447       surf%pt_surface(m)          = surf_t_surface_p%var_1d(m) / exn
     1448
     1449!
     1450!--    Calculate fluxes
     1451       surf%rad_net_l(m) = surf%rad_net_l(m) +                                 &
     1452                            rad_lw_out_change_0(j,i)                           &
     1453                          * surf_t_surface%var_1d(m) - rad_lw_out(nzb,j,i)     &
     1454                          - rad_lw_out_change_0(j,i) * surf_t_surface_p%var_1d(m)
     1455
     1456       rad_net(j,i) = surf%rad_net_l(m)
     1457       rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i) + rad_lw_out_change_0(j,i) *  &
     1458                     ( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) )
     1459
     1460       surf%ghf_eb(m) = lambda_surface * ( surf_t_surface_p%var_1d(m)          &
     1461                                             - surf_t_soil%var_2d(nzb_soil,m) )
     1462
     1463       surf%shf_eb(m) = - f_shf * ( pt1 - surf%pt_surface(m) )
     1464
     1465       surf%shf(m) = surf%shf_eb(m) / rho_cp
     1466
     1467       IF ( humidity )  THEN
     1468          surf%qsws_eb(m)  = - f_qsws * ( qv1 - q_s + dq_s_dt                  &
     1469                          * surf_t_surface%var_1d(m) - dq_s_dt *               &
     1470                            surf_t_surface_p%var_1d(m) )
     1471
     1472          surf%qsws(m) = surf%qsws_eb(m) / rho_lv
     1473
     1474          surf%qsws_veg_eb(m)  = - f_qsws_veg * ( qv1 - q_s                   &
     1475                              + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt   &
     1476                              * surf_t_surface_p%var_1d(m) )
     1477
     1478          surf%qsws_soil_eb(m) = - f_qsws_soil * ( qv1 - q_s                   &
     1479                              + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt   &
     1480                              * surf_t_surface_p%var_1d(m) )
     1481
     1482          surf%qsws_liq_eb(m)  = - f_qsws_liq  * ( qv1 - q_s                   &
     1483                              + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt   &
     1484                              * surf_t_surface_p%var_1d(m) )
     1485       ENDIF
     1486
     1487!
     1488!--    Calculate the true surface resistance
     1489       IF ( surf%qsws_eb(m) == 0.0_wp )  THEN
     1490          surf%r_s(m) = 1.0E10_wp
     1491       ELSE
     1492          surf%r_s(m) = - rho_lv * ( qv1 - q_s + dq_s_dt                       &
     1493                          * surf_t_surface%var_1d(m) - dq_s_dt *               &
     1494                            surf_t_surface_p%var_1d(m) ) /                     &
     1495                            surf%qsws_eb(m) - surf%r_a(m)
     1496       ENDIF
     1497
     1498!
     1499!--    Calculate change in liquid water reservoir due to dew fall or
     1500!--    evaporation of liquid water
     1501       IF ( humidity )  THEN
     1502!
     1503!--       If precipitation is activated, add rain water to qsws_liq_eb
     1504!--       and qsws_soil_eb according the the vegetation coverage.
     1505!--       precipitation_rate is given in mm.
     1506          IF ( precipitation )  THEN
     1507
     1508!
     1509!--          Add precipitation to liquid water reservoir, if possible.
     1510!--          Otherwise, add the water to soil. In case of
     1511!--          pavements, the exceeding water amount is implicitely removed
     1512!--          as runoff as qsws_soil_eb is then not used in the soil model
     1513             IF ( surf_m_liq_eb%var_1d(m) /= m_liq_eb_max )  THEN
     1514                surf%qsws_liq_eb(m) = surf%qsws_liq_eb(m)                      &
     1515                                 + surf%c_veg(m) * prr(k+k_off,j+j_off,i+i_off)&
     1516                                 * hyrho(k+k_off)                              &
     1517                                 * 0.001_wp * rho_l * l_v
     1518             ELSE
     1519                surf%qsws_soil_eb(m) = surf%qsws_soil_eb(m)                    &
     1520                                 + surf%c_veg(m) * prr(k+k_off,j+j_off,i+i_off)&
     1521                                 * hyrho(k+k_off)                              &
     1522                                 * 0.001_wp * rho_l * l_v
    14321523             ENDIF
    14331524
    1434 !
    1435 !--          If the air is saturated, check the reservoir water level
    1436              IF ( qsws_eb(j,i) < 0.0_wp )  THEN
    1437 
    1438 !
    1439 !--             Check if reservoir is full (avoid values > m_liq_eb_max)
    1440 !--             In that case, qsws_liq_eb goes to qsws_soil_eb. In this
    1441 !--             case qsws_veg_eb is zero anyway (because c_liq = 1),       
    1442 !--             so that tend is zero and no further check is needed
    1443                 IF ( m_liq_eb(j,i) == m_liq_eb_max )  THEN
    1444                    qsws_soil_eb(j,i) = qsws_soil_eb(j,i)                       &
    1445                                         + qsws_liq_eb(j,i)
    1446 
    1447                    qsws_liq_eb(j,i)  = 0.0_wp
    1448                 ENDIF
    1449 
    1450 !
    1451 !--             In case qsws_veg_eb becomes negative (unphysical behavior),
    1452 !--             let the water enter the liquid water reservoir as dew on the
    1453 !--             plant
    1454                 IF ( qsws_veg_eb(j,i) < 0.0_wp )  THEN
    1455                    qsws_liq_eb(j,i) = qsws_liq_eb(j,i) + qsws_veg_eb(j,i)
    1456                    qsws_veg_eb(j,i) = 0.0_wp
    1457                 ENDIF
    1458              ENDIF                   
    1459  
    1460              tend = - qsws_liq_eb(j,i) * drho_l_lv
    1461              m_liq_eb_p(j,i) = m_liq_eb(j,i) + dt_3d * ( tsc(2) * tend         &
    1462                                                 + tsc(3) * tm_liq_eb_m(j,i) )
    1463 
    1464 !
    1465 !--          Check if reservoir is overfull -> reduce to maximum
    1466 !--          (conservation of water is violated here)
    1467              m_liq_eb_p(j,i) = MIN(m_liq_eb_p(j,i),m_liq_eb_max)
    1468 
    1469 !
    1470 !--          Check if reservoir is empty (avoid values < 0.0)
    1471 !--          (conservation of water is violated here)
    1472              m_liq_eb_p(j,i) = MAX(m_liq_eb_p(j,i),0.0_wp)
    1473 
    1474 
    1475 !
    1476 !--          Calculate m_liq_eb tendencies for the next Runge-Kutta step
    1477              IF ( timestep_scheme(1:5) == 'runge' )  THEN
    1478                 IF ( intermediate_timestep_count == 1 )  THEN
    1479                    tm_liq_eb_m(j,i) = tend
    1480                 ELSEIF ( intermediate_timestep_count <                         &
    1481                          intermediate_timestep_count_max )  THEN
    1482                    tm_liq_eb_m(j,i) = -9.5625_wp * tend + 5.3125_wp            &
    1483                                     * tm_liq_eb_m(j,i)
    1484                 ENDIF
     1525!--          Add precipitation to bare soil according to the bare soil
     1526!--          coverage.
     1527             surf%qsws_soil_eb(m) = surf%qsws_soil_eb(m) + ( 1.0_wp            &
     1528                               - surf%c_veg(m) ) * prr(k+k_off,j+j_off,i+i_off)&
     1529                               * hyrho(k+k_off)                                &
     1530                               * 0.001_wp * rho_l * l_v
     1531          ENDIF
     1532
     1533!
     1534!--       If the air is saturated, check the reservoir water level
     1535          IF ( surf%qsws_eb(m) < 0.0_wp )  THEN
     1536!
     1537!--          Check if reservoir is full (avoid values > m_liq_eb_max)
     1538!--          In that case, qsws_liq_eb goes to qsws_soil_eb. In this
     1539!--          case qsws_veg_eb is zero anyway (because c_liq = 1),       
     1540!--          so that tend is zero and no further check is needed
     1541             IF ( surf_m_liq_eb%var_1d(m) == m_liq_eb_max )  THEN
     1542                surf%qsws_soil_eb(m) = surf%qsws_soil_eb(m) + surf%qsws_liq_eb(m)
     1543
     1544                surf%qsws_liq_eb(m)  = 0.0_wp
    14851545             ENDIF
    14861546
    1487           ENDIF
    1488 
    1489        ENDDO
     1547!
     1548!--          In case qsws_veg_eb becomes negative (unphysical behavior),
     1549!--          let the water enter the liquid water reservoir as dew on the
     1550!--          plant
     1551             IF ( surf%qsws_veg_eb(m) < 0.0_wp )  THEN
     1552                surf%qsws_liq_eb(m) = surf%qsws_liq_eb(m) + surf%qsws_veg_eb(m)
     1553                surf%qsws_veg_eb(m) = 0.0_wp
     1554             ENDIF
     1555          ENDIF                   
     1556 
     1557          tend = - surf%qsws_liq_eb(m) * drho_l_lv
     1558          surf_m_liq_eb_p%var_1d(m) = surf_m_liq_eb%var_1d(m) + dt_3d *        &
     1559                                        ( tsc(2) * tend +                      &
     1560                                          tsc(3) * surf_tm_liq_eb_m%var_1d(m) )
     1561!
     1562!--       Check if reservoir is overfull -> reduce to maximum
     1563!--       (conservation of water is violated here)
     1564          surf_m_liq_eb_p%var_1d(m) = MIN( surf_m_liq_eb_p%var_1d(m),m_liq_eb_max )
     1565
     1566!
     1567!--       Check if reservoir is empty (avoid values < 0.0)
     1568!--       (conservation of water is violated here)
     1569          surf_m_liq_eb_p%var_1d(m) = MAX( surf_m_liq_eb_p%var_1d(m), 0.0_wp )
     1570!
     1571!--       Calculate m_liq_eb tendencies for the next Runge-Kutta step
     1572          IF ( timestep_scheme(1:5) == 'runge' )  THEN
     1573             IF ( intermediate_timestep_count == 1 )  THEN
     1574                surf_tm_liq_eb_m%var_1d(m) = tend
     1575             ELSEIF ( intermediate_timestep_count <                            &
     1576                      intermediate_timestep_count_max )  THEN
     1577                surf_tm_liq_eb_m%var_1d(m) = -9.5625_wp * tend +               &
     1578                                              5.3125_wp * surf_tm_liq_eb_m%var_1d(m)
     1579             ENDIF
     1580          ENDIF
     1581
     1582       ENDIF
     1583
    14901584    ENDDO
    14911585
     
    14971591#if defined( __parallel )
    14981592       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1499        CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call,    &
     1593       CALL MPI_ALLREDUCE( force_radiation_call_l, force_radiation_call,       &
    15001594                           1, MPI_LOGICAL, MPI_LOR, comm2d, ierr )
    15011595#else
     
    15131607!
    15141608!-- Calculate new roughness lengths (for water surfaces only)
    1515     CALL calc_z0_water_surface
     1609    IF ( horizontal )  CALL calc_z0_water_surface
     1610
     1611    CONTAINS
     1612!------------------------------------------------------------------------------!
     1613! Description:
     1614! ------------
     1615!> Calculation of specific humidity of the skin layer (surface). It is assumend
     1616!> that the skin is always saturated.
     1617!------------------------------------------------------------------------------!
     1618    SUBROUTINE calc_q_surface
     1619
     1620       IMPLICIT NONE
     1621
     1622       REAL(wp) :: resistance    !< aerodynamic and soil resistance term
     1623
     1624       DO  m = 1, surf%ns
     1625
     1626          i   = surf%i(m)           
     1627          j   = surf%j(m)
     1628          k   = surf%k(m)
     1629
     1630!
     1631!--       Calculate water vapour pressure at saturation
     1632          e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp *                         &
     1633                                 ( surf_t_surface_p%var_1d(m) - 273.16_wp ) /  &
     1634                                 ( surf_t_surface_p%var_1d(m) - 35.86_wp  )    &
     1635                                         )
     1636
     1637!
     1638!--       Calculate specific humidity at saturation
     1639          q_s = 0.622_wp * e_s / surface_pressure
     1640
     1641          resistance = surf%r_a(m) / ( surf%r_a(m) + surf%r_s(m) )
     1642
     1643!
     1644!--       Calculate specific humidity at surface
     1645          IF ( cloud_physics )  THEN
     1646             q(k+k_off,j+j_off,i+i_off) = resistance * q_s +                   &
     1647                                        ( 1.0_wp - resistance ) *              &
     1648                                        ( q(k,j,i) - ql(k,j,i) )
     1649          ELSE
     1650             q(k+k_off,j+j_off,i+i_off) = resistance * q_s +                   &
     1651                                        ( 1.0_wp - resistance ) *              &
     1652                                          q(k,j,i)
     1653          ENDIF
     1654
     1655!
     1656!--       Update virtual potential temperature
     1657          vpt(k+k_off,j+j_off,i+i_off) = pt(k+k_off,j+j_off,i+i_off) *         &
     1658                     ( 1.0_wp + 0.61_wp * q(k+k_off,j+j_off,i+i_off) )
     1659
     1660       ENDDO
     1661
     1662    END SUBROUTINE calc_q_surface
     1663
    15161664
    15171665
     
    16071755       IMPLICIT NONE
    16081756
    1609        INTEGER(iwp) ::  i !< running index
    1610        INTEGER(iwp) ::  j !< running index
    1611        INTEGER(iwp) ::  k !< running index
     1757       INTEGER(iwp) ::  i     !< running index
     1758       INTEGER(iwp) ::  i_off !< index offset of surface element, seen from atmospheric grid point
     1759       INTEGER(iwp) ::  j     !< running index
     1760       INTEGER(iwp) ::  j_off !< index offset of surface element, seen from atmospheric grid point
     1761       INTEGER(iwp) ::  k     !< running index
     1762       INTEGER(iwp) ::  l     !< running index surface facing
     1763       INTEGER(iwp) ::  m     !< running index
    16121764
    16131765       REAL(wp) :: pt1   !< potential temperature at first grid level
    1614 
    16151766
    16161767!
    16171768!--    Calculate Exner function
    16181769       exn = ( surface_pressure / 1000.0_wp )**0.286_wp
    1619 
    1620 
    16211770!
    16221771!--    If no cloud physics is used, rho_surface has not been calculated before
     
    16341783!
    16351784!--    Set inital values for prognostic quantities
    1636        tt_surface_m = 0.0_wp
    1637        tt_soil_m    = 0.0_wp
    1638        tm_soil_m    = 0.0_wp
    1639        tm_liq_eb_m  = 0.0_wp
    1640        c_liq        = 0.0_wp
    1641 
    1642        ghf_eb = 0.0_wp
    1643        shf_eb = rho_cp * shf
     1785!--    Horizontal surfaces
     1786       tt_surface_h_m%var_1d = 0.0_wp
     1787       tt_soil_h_m%var_2d    = 0.0_wp
     1788       tm_soil_h_m%var_2d    = 0.0_wp
     1789       tm_liq_eb_h_m%var_1d  = 0.0_wp
     1790       surf_lsm_h%c_liq = 0.0_wp
     1791
     1792       surf_lsm_h%ghf_eb = 0.0_wp
     1793       surf_lsm_h%shf_eb = rho_cp * surf_lsm_h%shf
    16441794
    16451795       IF ( humidity )  THEN
    1646           qsws_eb = rho_lv * qsws
     1796          surf_lsm_h%qsws_eb = rho_lv * surf_lsm_h%qsws
    16471797       ELSE
    1648           qsws_eb = 0.0_wp
    1649        ENDIF
    1650 
    1651        qsws_liq_eb  = 0.0_wp
    1652        qsws_soil_eb = 0.0_wp
    1653        qsws_veg_eb  = 0.0_wp
    1654 
    1655        r_a        = 50.0_wp
    1656        r_s        = 50.0_wp
    1657        r_canopy   = 0.0_wp
    1658        r_soil     = 0.0_wp
     1798          surf_lsm_h%qsws_eb = 0.0_wp
     1799       ENDIF
     1800
     1801       surf_lsm_h%qsws_liq_eb  = 0.0_wp
     1802       surf_lsm_h%qsws_soil_eb = 0.0_wp
     1803       surf_lsm_h%qsws_veg_eb  = 0.0_wp
     1804
     1805       surf_lsm_h%r_a        = 50.0_wp
     1806       surf_lsm_h%r_s        = 50.0_wp
     1807       surf_lsm_h%r_canopy   = 0.0_wp
     1808       surf_lsm_h%r_soil     = 0.0_wp
     1809!
     1810!--    Do the same for vertical surfaces
     1811       DO  l = 0, 3
     1812          tt_surface_v_m(l)%var_1d = 0.0_wp
     1813          tt_soil_v_m(l)%var_2d    = 0.0_wp
     1814          tm_soil_v_m(l)%var_2d    = 0.0_wp
     1815          tm_liq_eb_v_m(l)%var_1d  = 0.0_wp
     1816          surf_lsm_v(l)%c_liq = 0.0_wp
     1817
     1818          surf_lsm_v(l)%ghf_eb = 0.0_wp
     1819          surf_lsm_v(l)%shf_eb = rho_cp * surf_lsm_v(l)%shf
     1820
     1821          IF ( humidity )  THEN
     1822             surf_lsm_v(l)%qsws_eb = rho_lv * surf_lsm_v(l)%qsws
     1823          ELSE
     1824             surf_lsm_v(l)%qsws_eb = 0.0_wp
     1825          ENDIF
     1826
     1827          surf_lsm_v(l)%qsws_liq_eb  = 0.0_wp
     1828          surf_lsm_v(l)%qsws_soil_eb = 0.0_wp
     1829          surf_lsm_v(l)%qsws_veg_eb  = 0.0_wp
     1830
     1831          surf_lsm_v(l)%r_a        = 50.0_wp
     1832          surf_lsm_v(l)%r_s        = 50.0_wp
     1833          surf_lsm_v(l)%r_canopy   = 0.0_wp
     1834          surf_lsm_v(l)%r_soil     = 0.0_wp
     1835       ENDDO
    16591836
    16601837!
    16611838!--    Allocate 3D soil model arrays
    1662        ALLOCATE ( root_fr(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
    1663        ALLOCATE ( lambda_h(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
    1664        ALLOCATE ( rho_c_total(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
    1665 
    1666        lambda_h = 0.0_wp
     1839!--    First, for horizontal surfaces
     1840       ALLOCATE ( surf_lsm_h%root_fr(nzb_soil:nzt_soil,1:surf_lsm_h%ns)     )
     1841       ALLOCATE ( surf_lsm_h%lambda_h(nzb_soil:nzt_soil,1:surf_lsm_h%ns)    )
     1842       ALLOCATE ( surf_lsm_h%rho_c_total(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
     1843
     1844       surf_lsm_h%lambda_h = 0.0_wp
    16671845!
    16681846!--    If required, allocate humidity-related variables for the soil model
    16691847       IF ( humidity )  THEN
    1670           ALLOCATE ( lambda_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
    1671           ALLOCATE ( gamma_w(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )   
    1672 
    1673           lambda_w = 0.0_wp
    1674        ENDIF
     1848          ALLOCATE ( surf_lsm_h%lambda_w(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
     1849          ALLOCATE ( surf_lsm_h%gamma_w(nzb_soil:nzt_soil,1:surf_lsm_h%ns)  ) 
     1850
     1851          surf_lsm_h%lambda_w = 0.0_wp
     1852       ENDIF
     1853!
     1854!--    For vertical surfaces
     1855       DO  l = 0, 3
     1856          ALLOCATE ( surf_lsm_v(l)%root_fr(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)     )
     1857          ALLOCATE ( surf_lsm_v(l)%lambda_h(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)    )
     1858          ALLOCATE ( surf_lsm_v(l)%rho_c_total(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
     1859
     1860          surf_lsm_v(l)%lambda_h = 0.0_wp
     1861!
     1862!--       If required, allocate humidity-related variables for the soil model
     1863          IF ( humidity )  THEN
     1864             ALLOCATE ( surf_lsm_v(l)%lambda_w(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
     1865             ALLOCATE ( surf_lsm_v(l)%gamma_w(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)  ) 
     1866
     1867             surf_lsm_v(l)%lambda_w = 0.0_wp
     1868          ENDIF     
     1869       ENDDO
    16751870
    16761871!
     
    17351930!
    17361931!--    Map values to the respective 2D arrays
    1737        alpha_vg      = alpha_vangenuchten
    1738        l_vg          = l_vangenuchten
    1739        n_vg          = n_vangenuchten
    1740        gamma_w_sat   = hydraulic_conductivity
    1741        m_sat         = saturation_moisture
    1742        m_fc          = field_capacity
    1743        m_wilt        = wilting_point
    1744        m_res         = residual_moisture
    1745        r_soil_min    = min_soil_resistance
     1932!--    Horizontal surfaces
     1933       surf_lsm_h%alpha_vg      = alpha_vangenuchten
     1934       surf_lsm_h%l_vg          = l_vangenuchten
     1935       surf_lsm_h%n_vg          = n_vangenuchten
     1936       surf_lsm_h%gamma_w_sat   = hydraulic_conductivity
     1937       surf_lsm_h%m_sat         = saturation_moisture
     1938       surf_lsm_h%m_fc          = field_capacity
     1939       surf_lsm_h%m_wilt        = wilting_point
     1940       surf_lsm_h%m_res         = residual_moisture
     1941       surf_lsm_h%r_soil_min    = min_soil_resistance
     1942!
     1943!--    Vertical surfaces
     1944       DO  l = 0, 3
     1945          surf_lsm_v(l)%alpha_vg      = alpha_vangenuchten
     1946          surf_lsm_v(l)%l_vg          = l_vangenuchten
     1947          surf_lsm_v(l)%n_vg          = n_vangenuchten
     1948          surf_lsm_v(l)%gamma_w_sat   = hydraulic_conductivity
     1949          surf_lsm_v(l)%m_sat         = saturation_moisture
     1950          surf_lsm_v(l)%m_fc          = field_capacity
     1951          surf_lsm_v(l)%m_wilt        = wilting_point
     1952          surf_lsm_v(l)%m_res         = residual_moisture
     1953          surf_lsm_v(l)%r_soil_min    = min_soil_resistance
     1954       ENDDO
     1955
     1956
    17461957
    17471958!
    17481959!--    Initial run actions
    17491960       IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
    1750 
    1751           t_soil    = 0.0_wp
    1752           m_liq_eb  = 0.0_wp
    1753           m_soil    = 0.0_wp
     1961!
     1962!--       First, for horizontal surfaces
     1963          t_soil_h%var_2d    = 0.0_wp
     1964          m_soil_h%var_2d    = 0.0_wp
     1965          m_liq_eb_h%var_1d  = 0.0_wp
    17541966
    17551967!
     
    17581970!--       wilting point) -> problems with devision by zero)
    17591971          DO  k = nzb_soil, nzt_soil
    1760              t_soil(k,:,:)    = soil_temperature(k)
    1761              m_soil(k,:,:)    = MAX(soil_moisture(k),m_wilt(:,:))
     1972             t_soil_h%var_2d(k,:)      = soil_temperature(k)
     1973             m_soil_h%var_2d(k,:)      = MAX(soil_moisture(k),surf_lsm_h%m_wilt(:))
    17621974             soil_moisture(k) = MAX(soil_moisture(k),wilting_point)
    17631975          ENDDO
    1764           t_soil(nzt_soil+1,:,:) = soil_temperature(nzt_soil+1)
     1976          t_soil_h%var_2d(nzt_soil+1,:) = soil_temperature(nzt_soil+1)
    17651977
    17661978!
    17671979!--       Calculate surface temperature
    1768           t_surface   = pt_surface * exn
     1980          t_surface_h%var_1d(:)    = pt_surface * exn
     1981          surf_lsm_h%pt_surface(:) = pt_surface
    17691982
    17701983!
    17711984!--       Set artifical values for ts and us so that r_a has its initial value
    1772 !--       for the first time step
    1773           DO  i = nxlg, nxrg
    1774              DO  j = nysg, nyng
    1775                 k = nzb_s_inner(j,i)
     1985!--       for the first time step. Only for interior core domain, not for ghost points
     1986          DO  m = 1, surf_lsm_h%ns
     1987             i   = surf_lsm_h%i(m)           
     1988             j   = surf_lsm_h%j(m)
     1989             k   = surf_lsm_h%k(m)
     1990
     1991             IF ( cloud_physics )  THEN
     1992                pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
     1993             ELSE
     1994                pt1 = pt(k,j,i)
     1995             ENDIF
     1996
     1997!
     1998!--          Assure that r_a cannot be zero at model start
     1999             IF ( pt1 == surf_lsm_h%pt_surface(m) )  pt1 = pt1 + 1.0E-10_wp
     2000
     2001             surf_lsm_h%us(m)   = 0.1_wp
     2002             surf_lsm_h%ts(m)   = ( pt1 - surf_lsm_h%pt_surface(m) ) / surf_lsm_h%r_a(m)
     2003             surf_lsm_h%shf(m)  = - surf_lsm_h%us(m) * surf_lsm_h%ts(m)
     2004          ENDDO
     2005!
     2006!--       Vertical surfaces
     2007          DO  l = 0, 3
     2008
     2009             t_soil_v(l)%var_2d    = 0.0_wp
     2010             m_soil_v(l)%var_2d    = 0.0_wp
     2011             m_liq_eb_v(l)%var_1d  = 0.0_wp
     2012
     2013!
     2014!--          Map user settings of T and q for each soil layer
     2015!--          (make sure that the soil moisture does not drop below the permanent
     2016!--          wilting point) -> problems with devision by zero)
     2017             DO  k = nzb_soil, nzt_soil
     2018                t_soil_v(l)%var_2d(k,:)      = soil_temperature(k)
     2019                m_soil_v(l)%var_2d(k,:)      = MAX(soil_moisture(k),surf_lsm_v(l)%m_wilt(:))
     2020                soil_moisture(k) = MAX(soil_moisture(k),wilting_point)
     2021             ENDDO
     2022             t_soil_v(l)%var_2d(nzt_soil+1,:) = soil_temperature(nzt_soil+1)
     2023
     2024!
     2025!--          Calculate surface temperature
     2026             t_surface_v(l)%var_1d(:)   = pt_surface * exn
     2027             surf_lsm_h%pt_surface(:)   = pt_surface
     2028
     2029!
     2030!--          Set artifical values for ts and us so that r_a has its initial value
     2031!--          for the first time step. Only for interior core domain, not for ghost points
     2032             DO  m = 1, surf_lsm_v(l)%ns
     2033                i   = surf_lsm_v(l)%i(m)           
     2034                j   = surf_lsm_v(l)%j(m)
     2035                k   = surf_lsm_v(l)%k(m)
    17762036
    17772037                IF ( cloud_physics )  THEN
    1778                    pt1 = pt(k+1,j,i) + l_d_cp * pt_d_t(k+1) * ql(k+1,j,i)
     2038                   pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
    17792039                ELSE
    1780                    pt1 = pt(k+1,j,i)
     2040                   pt1 = pt(k,j,i)
    17812041                ENDIF
    17822042
    17832043!
    17842044!--             Assure that r_a cannot be zero at model start
    1785                 IF ( pt1 == pt(k,j,i) )  pt1 = pt1 + 1.0E-10_wp
    1786 
    1787                 us(j,i)  = 0.1_wp
    1788                 ts(j,i)  = (pt1 - pt(k,j,i)) / r_a(j,i)
    1789                 shf(j,i) = - us(j,i) * ts(j,i)
     2045                IF ( pt1 == surf_lsm_h%pt_surface(m) )  pt1 = pt1 + 1.0E-10_wp
     2046
     2047                surf_lsm_v(l)%us(m)   = 0.1_wp
     2048                surf_lsm_v(l)%ts(m)   = ( pt1 - surf_lsm_h%pt_surface(m) ) / surf_lsm_v(l)%r_a(m)
     2049                surf_lsm_v(l)%shf(m)  = - surf_lsm_v(l)%us(m) * surf_lsm_v(l)%ts(m)
    17902050             ENDDO
    17912051          ENDDO
     
    17942054!--    Actions for restart runs
    17952055       ELSE
    1796 
    1797           DO  i = nxlg, nxrg
    1798              DO  j = nysg, nyng
    1799                 k = nzb_s_inner(j,i)               
    1800                 t_surface(j,i) = pt(k,j,i) * exn
    1801              ENDDO
     2056!
     2057!--       Horizontal surfaces
     2058          DO  m = 1, surf_lsm_h%ns
     2059             i   = surf_lsm_h%i(m)           
     2060             j   = surf_lsm_h%j(m)
     2061             k   = surf_lsm_h%k(m)         
     2062             t_surface_h%var_1d(m) = pt(k-1,j,i) * exn
    18022063          ENDDO
    1803 
    1804        ENDIF
    1805 
    1806        DO  k = nzb_soil, nzt_soil
    1807           root_fr(k,:,:) = root_fraction(k)
     2064!
     2065!--       Vertical surfaces
     2066          DO  l = 0, 3
     2067!
     2068!--          Set index offset of surface element, seen from atmospheric grid point
     2069             IF ( l == 0 )  THEN
     2070                j_off = -1
     2071                i_off = 0
     2072             ELSEIF ( l == 1 )  THEN
     2073                j_off = 1
     2074                i_off = 0
     2075             ELSEIF ( l == 2 )  THEN
     2076                j_off = 0
     2077                i_off = -1
     2078             ELSEIF ( l == 3 )  THEN
     2079                j_off = 0
     2080                i_off = 1
     2081             ENDIF
     2082             DO  m = 1, surf_lsm_v(l)%ns
     2083                i   = surf_lsm_v(l)%i(m)           
     2084                j   = surf_lsm_v(l)%j(m)
     2085                k   = surf_lsm_v(l)%k(m)         
     2086                t_surface_v(l)%var_1d(m) = pt(k,j+j_off,i+i_off) * exn
     2087             ENDDO
     2088          ENDDO
     2089
     2090       ENDIF
     2091!
     2092!--    Initialize root fraction
     2093!--    Horizontal surfaces
     2094       DO  m = 1, surf_lsm_h%ns
     2095          i   = surf_lsm_h%i(m)           
     2096          j   = surf_lsm_h%j(m)
     2097
     2098          DO  k = nzb_soil, nzt_soil
     2099             surf_lsm_h%root_fr(k,m) = root_fraction(k)
     2100          ENDDO
     2101       ENDDO
     2102!
     2103!--    Vertical surfaces
     2104       DO  l = 0, 3
     2105          DO  m = 1, surf_lsm_v(l)%ns
     2106             i   = surf_lsm_v(l)%i(m)           
     2107             j   = surf_lsm_v(l)%j(m)
     2108
     2109             DO  k = nzb_soil, nzt_soil
     2110                surf_lsm_v(l)%root_fr(k,m) = root_fraction(k)
     2111             ENDDO
     2112          ENDDO
    18082113       ENDDO
    18092114
     
    18432148
    18442149          IF ( ANY( root_fraction == 9999999.9_wp ) )  THEN
    1845              DO  k = nzb_soil, nzt_soil
    1846                 root_fr(k,:,:) = root_distribution(k,veg_type)
    1847                 root_fraction(k) = root_distribution(k,veg_type)
     2150             DO  m = 1, surf_lsm_h%ns
     2151                i   = surf_lsm_h%i(m)           
     2152                j   = surf_lsm_h%j(m)
     2153
     2154                DO  k = nzb_soil, nzt_soil
     2155                   surf_lsm_h%root_fr(k,m) = root_distribution(k,veg_type)
     2156                   root_fraction(k)    = root_distribution(k,veg_type)
     2157                ENDDO
     2158             ENDDO
     2159             DO  l = 0, 3
     2160                DO  m = 1, surf_lsm_v(l)%ns
     2161                   i   = surf_lsm_v(l)%i(m)           
     2162                   j   = surf_lsm_v(l)%j(m)
     2163
     2164                   DO  k = nzb_soil, nzt_soil
     2165                      surf_lsm_v(l)%root_fr(k,m) = root_distribution(k,veg_type)
     2166                      root_fraction(k)    = root_distribution(k,veg_type)
     2167                   ENDDO
     2168                ENDDO
    18482169             ENDDO
    18492170          ENDIF
     
    18772198!--    Map vegetation and soil types to 2D array to allow for heterogeneous
    18782199!--    surfaces via user interface see below
    1879        veg_type_2d = veg_type
    1880        soil_type_2d = soil_type
     2200!--    First for horizontal surfaces
     2201       surf_lsm_h%veg_type_2d  = veg_type
     2202       surf_lsm_h%soil_type_2d = soil_type
    18812203
    18822204!
    18832205!--    Map vegetation parameters to the respective 2D arrays
    1884        r_canopy_min         = min_canopy_resistance
    1885        lai                  = leaf_area_index
    1886        c_veg                = vegetation_coverage
    1887        g_d                  = canopy_resistance_coefficient
    1888        lambda_surface_s     = lambda_surface_stable
    1889        lambda_surface_u     = lambda_surface_unstable
    1890        f_sw_in              = f_shortwave_incoming
    1891        z0                   = z0_eb
    1892        z0h                  = z0h_eb
    1893        z0q                  = z0q_eb
     2206       surf_lsm_h%r_canopy_min      = min_canopy_resistance
     2207       surf_lsm_h%lai               = leaf_area_index
     2208       surf_lsm_h%c_veg             = vegetation_coverage
     2209       surf_lsm_h%g_d               = canopy_resistance_coefficient
     2210       surf_lsm_h%lambda_surface_s  = lambda_surface_stable
     2211       surf_lsm_h%lambda_surface_u  = lambda_surface_unstable
     2212       surf_lsm_h%f_sw_in           = f_shortwave_incoming
     2213       surf_lsm_h%z0                = z0_eb
     2214       surf_lsm_h%z0h               = z0h_eb
     2215       surf_lsm_h%z0q               = z0q_eb
     2216
     2217!--    Vertical surfaces
     2218       DO  l = 0, 3
     2219          surf_lsm_v(l)%veg_type_2d  = veg_type
     2220          surf_lsm_v(l)%soil_type_2d = soil_type
     2221
     2222!
     2223!--       Map vegetation parameters to the respective 2D arrays
     2224          surf_lsm_v(l)%r_canopy_min      = min_canopy_resistance
     2225          surf_lsm_v(l)%lai               = leaf_area_index
     2226          surf_lsm_v(l)%c_veg             = vegetation_coverage
     2227          surf_lsm_v(l)%g_d               = canopy_resistance_coefficient
     2228          surf_lsm_v(l)%lambda_surface_s  = lambda_surface_stable
     2229          surf_lsm_v(l)%lambda_surface_u  = lambda_surface_unstable
     2230          surf_lsm_v(l)%f_sw_in           = f_shortwave_incoming
     2231          surf_lsm_v(l)%z0                = z0_eb
     2232          surf_lsm_v(l)%z0h               = z0h_eb
     2233          surf_lsm_v(l)%z0q               = z0q_eb
     2234       ENDDO
    18942235
    18952236!
     
    19002241!--    Set flag parameter if vegetation type was set to a water surface. Also
    19012242!--    set temperature to a constant value in all "soil" layers.
    1902        DO  i = nxlg, nxrg
    1903           DO  j = nysg, nyng
    1904              IF ( veg_type_2d(j,i) == 14  .OR.  veg_type_2d(j,i) == 15 )  THEN
    1905                 water_surface(j,i) = .TRUE.
    1906                 t_soil(:,j,i) = t_surface(j,i)
    1907              ELSEIF ( veg_type_2d(j,i) == 20 )  THEN
    1908                 pave_surface(j,i) = .TRUE.
    1909                 m_soil(:,j,i) = 0.0_wp
    1910              ENDIF
    1911 
    1912           ENDDO
     2243!--    For now, do not set water surface for vertical surfaces
     2244       DO  m = 1, surf_lsm_h%ns
     2245          i   = surf_lsm_h%i(m)           
     2246          j   = surf_lsm_h%j(m)
     2247
     2248          IF ( surf_lsm_h%veg_type_2d(m) == 14  .OR.                           &
     2249               surf_lsm_h%veg_type_2d(m) == 15 )  THEN
     2250             surf_lsm_h%water_surface(m) = .TRUE.
     2251             t_soil_h%var_2d(:,m) = t_surface_h%var_1d(m)
     2252          ELSEIF ( surf_lsm_h%veg_type_2d(m) == 20 )  THEN
     2253             surf_lsm_h%pave_surface(m) = .TRUE.
     2254             m_soil_h%var_2d(:,m)            = 0.0_wp
     2255          ENDIF
     2256
    19132257       ENDDO
    19142258
    19152259!
    1916 !--    Calculate new roughness lengths (for water surfaces only)
     2260!--    Calculate new roughness lengths (for water surfaces only, i.e. only
     2261!-     horizontal surfaces)
    19172262       CALL calc_z0_water_surface
    19182263
    1919        t_soil_p    = t_soil
    1920        m_soil_p    = m_soil
    1921        m_liq_eb_p  = m_liq_eb
    1922        t_surface_p = t_surface
     2264       t_soil_h_p    = t_soil_h
     2265       m_soil_h_p    = m_soil_h
     2266       m_liq_eb_h_p  = m_liq_eb_h
     2267       t_surface_h_p = t_surface_h
     2268
     2269       t_soil_v_p    = t_soil_v
     2270       m_soil_v_p    = m_soil_v
     2271       m_liq_eb_v_p  = m_liq_eb_v
     2272       t_surface_v_p = t_surface_v
    19232273
    19242274
     
    19262276!--    Store initial profiles of t_soil and m_soil (assuming they are
    19272277!--    horizontally homogeneous on this PE)
    1928        hom(nzb_soil:nzt_soil,1,90,:)  = SPREAD( t_soil(nzb_soil:nzt_soil,      &
    1929                                                 nysg,nxlg), 2,                 &
    1930                                                 statistic_regions+1 )
    1931        hom(nzb_soil:nzt_soil,1,92,:)  = SPREAD( m_soil(nzb_soil:nzt_soil,      &
    1932                                                 nysg,nxlg), 2,                 &
    1933                                                 statistic_regions+1 )
     2278       hom(nzb_soil:nzt_soil,1,90,:)  = SPREAD( t_soil_h%var_2d(nzb_soil:nzt_soil,1),  &
     2279                                                2, statistic_regions+1 )
     2280       hom(nzb_soil:nzt_soil,1,92,:)  = SPREAD( m_soil_h%var_2d(nzb_soil:nzt_soil,1),  &
     2281                                                2, statistic_regions+1 )
    19342282
    19352283    END SUBROUTINE lsm_init
     
    19462294       IMPLICIT NONE
    19472295
    1948 !
    1949 !--    Allocate surface and soil temperature / humidity
     2296       INTEGER(iwp) ::  l !< index indicating facing of surface array
     2297
     2298!
     2299!--    Allocate surface and soil temperature / humidity. Please note,
     2300!--    these arrays are allocated according to surface-data structure,
     2301!--    even if they do not belong to the data type due to the
     2302!--    pointer arithmetric (TARGET attribute is not allowed in a data-type).
    19502303#if defined( __nopointer )
    1951        ALLOCATE ( m_liq_eb(nysg:nyng,nxlg:nxrg) )
    1952        ALLOCATE ( m_liq_eb_p(nysg:nyng,nxlg:nxrg) )
    1953        ALLOCATE ( m_soil(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
    1954        ALLOCATE ( m_soil_p(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
    1955        ALLOCATE ( t_surface(nysg:nyng,nxlg:nxrg) )
    1956        ALLOCATE ( t_surface_p(nysg:nyng,nxlg:nxrg) )
    1957        ALLOCATE ( t_soil(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
    1958        ALLOCATE ( t_soil_p(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
     2304!
     2305!--    Horizontal surfaces
     2306       ALLOCATE ( m_liq_eb_h%var_1d(1:surf_lsm_h%ns)                     )
     2307       ALLOCATE ( m_liq_eb_h_p%var_1d(1:surf_lsm_h%ns)                   )
     2308       ALLOCATE ( t_surface_h%var_1d(1:surf_lsm_h%ns)                    )
     2309       ALLOCATE ( t_surface_h_p%var_1d(1:surf_lsm_h%ns)                  )
     2310       ALLOCATE ( m_soil_h%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns)     )
     2311       ALLOCATE ( m_soil_h_p%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns)   )
     2312       ALLOCATE ( t_soil_h%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns)   )
     2313       ALLOCATE ( t_soil_h_p%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) )
     2314!
     2315!--    Vertical surfaces
     2316       DO  l = 0, 3
     2317          ALLOCATE ( m_liq_eb_v(l)%var_1d(1:surf_lsm_v(l)%ns)                     )
     2318          ALLOCATE ( m_liq_eb_v(l)_p%var_1d(1:surf_lsm_v(l)%ns)                   )
     2319          ALLOCATE ( t_surface_v(l)%var_1d(1:surf_lsm_v(l)%ns)                    )
     2320          ALLOCATE ( t_surface_v(l)_p%var_1d(1:surf_lsm_v(l)%ns)                  )
     2321          ALLOCATE ( m_soil_v(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)     )
     2322          ALLOCATE ( m_soil_v(l)_p%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)   )
     2323          ALLOCATE ( t_soil_v(l)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(l)%ns)   )
     2324          ALLOCATE ( t_soil_v(l)_p%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(l)%ns) )
     2325       ENDDO
    19592326#else
    1960        ALLOCATE ( m_liq_eb_1(nysg:nyng,nxlg:nxrg) )
    1961        ALLOCATE ( m_liq_eb_2(nysg:nyng,nxlg:nxrg) )
    1962        ALLOCATE ( m_soil_1(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
    1963        ALLOCATE ( m_soil_2(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
    1964        ALLOCATE ( t_surface_1(nysg:nyng,nxlg:nxrg) )
    1965        ALLOCATE ( t_surface_2(nysg:nyng,nxlg:nxrg) )
    1966        ALLOCATE ( t_soil_1(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
    1967        ALLOCATE ( t_soil_2(nzb_soil:nzt_soil+1,nysg:nyng,nxlg:nxrg) )
     2327!
     2328!--    Horizontal surfaces
     2329       ALLOCATE ( m_liq_eb_h_1%var_1d(1:surf_lsm_h%ns)                   )
     2330       ALLOCATE ( m_liq_eb_h_2%var_1d(1:surf_lsm_h%ns)                   )
     2331       ALLOCATE ( t_surface_h_1%var_1d(1:surf_lsm_h%ns)                  )
     2332       ALLOCATE ( t_surface_h_2%var_1d(1:surf_lsm_h%ns)                  )
     2333       ALLOCATE ( m_soil_h_1%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns)   )
     2334       ALLOCATE ( m_soil_h_2%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns)   )
     2335       ALLOCATE ( t_soil_h_1%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) )
     2336       ALLOCATE ( t_soil_h_2%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) )
     2337!
     2338!--    Vertical surfaces
     2339       DO  l = 0, 3
     2340          ALLOCATE ( m_liq_eb_v_1(l)%var_1d(1:surf_lsm_v(l)%ns)                   )
     2341          ALLOCATE ( m_liq_eb_v_2(l)%var_1d(1:surf_lsm_v(l)%ns)                   )
     2342          ALLOCATE ( t_surface_v_1(l)%var_1d(1:surf_lsm_v(l)%ns)                  )
     2343          ALLOCATE ( t_surface_v_2(l)%var_1d(1:surf_lsm_v(l)%ns)                  )
     2344          ALLOCATE ( m_soil_v_1(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)   )
     2345          ALLOCATE ( m_soil_v_2(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)   )
     2346          ALLOCATE ( t_soil_v_1(l)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(l)%ns) )
     2347          ALLOCATE ( t_soil_v_2(l)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(l)%ns) )
     2348       ENDDO
    19682349#endif
    19692350
    19702351!
    19712352!--    Allocate intermediate timestep arrays
    1972        ALLOCATE ( tm_liq_eb_m(nysg:nyng,nxlg:nxrg) )
    1973        ALLOCATE ( tm_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
    1974        ALLOCATE ( tt_surface_m(nysg:nyng,nxlg:nxrg) )
    1975        ALLOCATE ( tt_soil_m(nzb_soil:nzt_soil,nysg:nyng,nxlg:nxrg) )
     2353!--    Horizontal surfaces
     2354       ALLOCATE ( tm_liq_eb_h_m%var_1d(1:surf_lsm_h%ns)                  )
     2355       ALLOCATE ( tt_surface_h_m%var_1d(1:surf_lsm_h%ns)                 )
     2356       ALLOCATE ( tm_soil_h_m%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns)  )
     2357       ALLOCATE ( tt_soil_h_m%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns)  )
     2358!
     2359!--    Horizontal surfaces
     2360       DO  l = 0, 3
     2361          ALLOCATE ( tm_liq_eb_v_m(l)%var_1d(1:surf_lsm_v(l)%ns)                  )
     2362          ALLOCATE ( tt_surface_v_m(l)%var_1d(1:surf_lsm_v(l)%ns)                 )
     2363          ALLOCATE ( tm_soil_v_m(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)  )
     2364          ALLOCATE ( tt_soil_v_m(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)  )
     2365       ENDDO
     2366!
     2367!--    Allocate skin-surface temperature
     2368       ALLOCATE ( surf_lsm_h%pt_surface(1:surf_lsm_h%ns) )
     2369       DO  l = 0, 3
     2370          ALLOCATE ( surf_lsm_v(l)%pt_surface(1:surf_lsm_v(l)%ns) )
     2371       ENDDO
    19762372
    19772373!
    19782374!--    Allocate 2D vegetation model arrays
    1979        ALLOCATE ( alpha_vg(nysg:nyng,nxlg:nxrg) )
    1980        ALLOCATE ( building_surface(nysg:nyng,nxlg:nxrg) )
    1981        ALLOCATE ( c_liq(nysg:nyng,nxlg:nxrg) )
    1982        ALLOCATE ( c_veg(nysg:nyng,nxlg:nxrg) )
    1983        ALLOCATE ( f_sw_in(nysg:nyng,nxlg:nxrg) )
    1984        ALLOCATE ( ghf_eb(nysg:nyng,nxlg:nxrg) )
    1985        ALLOCATE ( gamma_w_sat(nysg:nyng,nxlg:nxrg) )
    1986        ALLOCATE ( g_d(nysg:nyng,nxlg:nxrg) )
    1987        ALLOCATE ( lai(nysg:nyng,nxlg:nxrg) )
    1988        ALLOCATE ( l_vg(nysg:nyng,nxlg:nxrg) )
    1989        ALLOCATE ( lambda_surface_u(nysg:nyng,nxlg:nxrg) )
    1990        ALLOCATE ( lambda_surface_s(nysg:nyng,nxlg:nxrg) )
    1991        ALLOCATE ( m_fc(nysg:nyng,nxlg:nxrg) )
    1992        ALLOCATE ( m_res(nysg:nyng,nxlg:nxrg) )
    1993        ALLOCATE ( m_sat(nysg:nyng,nxlg:nxrg) )
    1994        ALLOCATE ( m_wilt(nysg:nyng,nxlg:nxrg) )
    1995        ALLOCATE ( n_vg(nysg:nyng,nxlg:nxrg) )
    1996        ALLOCATE ( pave_surface(nysg:nyng,nxlg:nxrg) )
    1997        ALLOCATE ( qsws_eb(nysg:nyng,nxlg:nxrg) )
    1998        ALLOCATE ( qsws_soil_eb(nysg:nyng,nxlg:nxrg) )
    1999        ALLOCATE ( qsws_liq_eb(nysg:nyng,nxlg:nxrg) )
    2000        ALLOCATE ( qsws_veg_eb(nysg:nyng,nxlg:nxrg) )
    2001        ALLOCATE ( rad_net_l(nysg:nyng,nxlg:nxrg) )
    2002        ALLOCATE ( r_a(nysg:nyng,nxlg:nxrg) )
    2003        ALLOCATE ( r_canopy(nysg:nyng,nxlg:nxrg) )
    2004        ALLOCATE ( r_soil(nysg:nyng,nxlg:nxrg) )
    2005        ALLOCATE ( r_soil_min(nysg:nyng,nxlg:nxrg) )
    2006        ALLOCATE ( r_s(nysg:nyng,nxlg:nxrg) )
    2007        ALLOCATE ( r_canopy_min(nysg:nyng,nxlg:nxrg) )
    2008        ALLOCATE ( shf_eb(nysg:nyng,nxlg:nxrg) )
    2009        ALLOCATE ( soil_type_2d(nysg:nyng,nxlg:nxrg) )
    2010        ALLOCATE ( veg_type_2d(nysg:nyng,nxlg:nxrg) )
    2011        ALLOCATE ( water_surface(nysg:nyng,nxlg:nxrg) )
    2012 
    2013        water_surface = .FALSE.
    2014        pave_surface  = .FALSE.
    2015 
     2375!--    Horizontal surfaces
     2376       ALLOCATE ( surf_lsm_h%alpha_vg(1:surf_lsm_h%ns)         )
     2377       ALLOCATE ( surf_lsm_h%building_surface(1:surf_lsm_h%ns) )
     2378       ALLOCATE ( surf_lsm_h%c_liq(1:surf_lsm_h%ns)            )
     2379       ALLOCATE ( surf_lsm_h%c_veg(1:surf_lsm_h%ns)            )
     2380       ALLOCATE ( surf_lsm_h%f_sw_in(1:surf_lsm_h%ns)          )
     2381       ALLOCATE ( surf_lsm_h%ghf_eb(1:surf_lsm_h%ns)           )
     2382       ALLOCATE ( surf_lsm_h%gamma_w_sat(1:surf_lsm_h%ns)      )
     2383       ALLOCATE ( surf_lsm_h%g_d(1:surf_lsm_h%ns)              )
     2384       ALLOCATE ( surf_lsm_h%lai(1:surf_lsm_h%ns)              )
     2385       ALLOCATE ( surf_lsm_h%l_vg(1:surf_lsm_h%ns)             )
     2386       ALLOCATE ( surf_lsm_h%lambda_surface_u(1:surf_lsm_h%ns) )
     2387       ALLOCATE ( surf_lsm_h%lambda_surface_s(1:surf_lsm_h%ns) )
     2388       ALLOCATE ( surf_lsm_h%m_fc(1:surf_lsm_h%ns)             )
     2389       ALLOCATE ( surf_lsm_h%m_res(1:surf_lsm_h%ns)            )
     2390       ALLOCATE ( surf_lsm_h%m_sat(1:surf_lsm_h%ns)            )
     2391       ALLOCATE ( surf_lsm_h%m_wilt(1:surf_lsm_h%ns)           )
     2392       ALLOCATE ( surf_lsm_h%n_vg(1:surf_lsm_h%ns)             )
     2393       ALLOCATE ( surf_lsm_h%pave_surface(1:surf_lsm_h%ns)     )
     2394       ALLOCATE ( surf_lsm_h%qsws_eb(1:surf_lsm_h%ns)          )
     2395       ALLOCATE ( surf_lsm_h%qsws_soil_eb(1:surf_lsm_h%ns)     )
     2396       ALLOCATE ( surf_lsm_h%qsws_liq_eb(1:surf_lsm_h%ns)      )
     2397       ALLOCATE ( surf_lsm_h%qsws_veg_eb(1:surf_lsm_h%ns)      )
     2398       ALLOCATE ( surf_lsm_h%rad_net_l(1:surf_lsm_h%ns)        )
     2399       ALLOCATE ( surf_lsm_h%r_a(1:surf_lsm_h%ns)              )
     2400       ALLOCATE ( surf_lsm_h%r_canopy(1:surf_lsm_h%ns)         )
     2401       ALLOCATE ( surf_lsm_h%r_soil(1:surf_lsm_h%ns)           )
     2402       ALLOCATE ( surf_lsm_h%r_soil_min(1:surf_lsm_h%ns)       )
     2403       ALLOCATE ( surf_lsm_h%r_s(1:surf_lsm_h%ns)              )
     2404       ALLOCATE ( surf_lsm_h%r_canopy_min(1:surf_lsm_h%ns)     )
     2405       ALLOCATE ( surf_lsm_h%shf_eb(1:surf_lsm_h%ns)           )
     2406       ALLOCATE ( surf_lsm_h%soil_type_2d(1:surf_lsm_h%ns)     )
     2407       ALLOCATE ( surf_lsm_h%veg_type_2d(1:surf_lsm_h%ns)      )
     2408       ALLOCATE ( surf_lsm_h%water_surface(1:surf_lsm_h%ns)    )
     2409
     2410       surf_lsm_h%water_surface = .FALSE.
     2411       surf_lsm_h%pave_surface  = .FALSE.
     2412!
     2413!--    Vertical surfaces
     2414       DO  l = 0, 3
     2415          ALLOCATE ( surf_lsm_v(l)%alpha_vg(1:surf_lsm_v(l)%ns)         )
     2416          ALLOCATE ( surf_lsm_v(l)%building_surface(1:surf_lsm_v(l)%ns) )
     2417          ALLOCATE ( surf_lsm_v(l)%c_liq(1:surf_lsm_v(l)%ns)            )
     2418          ALLOCATE ( surf_lsm_v(l)%c_veg(1:surf_lsm_v(l)%ns)            )
     2419          ALLOCATE ( surf_lsm_v(l)%f_sw_in(1:surf_lsm_v(l)%ns)          )
     2420          ALLOCATE ( surf_lsm_v(l)%ghf_eb(1:surf_lsm_v(l)%ns)           )
     2421          ALLOCATE ( surf_lsm_v(l)%gamma_w_sat(1:surf_lsm_v(l)%ns)      )
     2422          ALLOCATE ( surf_lsm_v(l)%g_d(1:surf_lsm_v(l)%ns)              )
     2423          ALLOCATE ( surf_lsm_v(l)%lai(1:surf_lsm_v(l)%ns)              )
     2424          ALLOCATE ( surf_lsm_v(l)%l_vg(1:surf_lsm_v(l)%ns)             )
     2425          ALLOCATE ( surf_lsm_v(l)%lambda_surface_u(1:surf_lsm_v(l)%ns) )
     2426          ALLOCATE ( surf_lsm_v(l)%lambda_surface_s(1:surf_lsm_v(l)%ns) )
     2427          ALLOCATE ( surf_lsm_v(l)%m_fc(1:surf_lsm_v(l)%ns)             )
     2428          ALLOCATE ( surf_lsm_v(l)%m_res(1:surf_lsm_v(l)%ns)            )
     2429          ALLOCATE ( surf_lsm_v(l)%m_sat(1:surf_lsm_v(l)%ns)            )
     2430          ALLOCATE ( surf_lsm_v(l)%m_wilt(1:surf_lsm_v(l)%ns)           )
     2431          ALLOCATE ( surf_lsm_v(l)%n_vg(1:surf_lsm_v(l)%ns)             )
     2432          ALLOCATE ( surf_lsm_v(l)%pave_surface(1:surf_lsm_v(l)%ns)     )
     2433          ALLOCATE ( surf_lsm_v(l)%qsws_eb(1:surf_lsm_v(l)%ns)          )
     2434          ALLOCATE ( surf_lsm_v(l)%qsws_soil_eb(1:surf_lsm_v(l)%ns)     )
     2435          ALLOCATE ( surf_lsm_v(l)%qsws_liq_eb(1:surf_lsm_v(l)%ns)      )
     2436          ALLOCATE ( surf_lsm_v(l)%qsws_veg_eb(1:surf_lsm_v(l)%ns)      )
     2437          ALLOCATE ( surf_lsm_v(l)%rad_net_l(1:surf_lsm_v(l)%ns)        )
     2438          ALLOCATE ( surf_lsm_v(l)%r_a(1:surf_lsm_v(l)%ns)              )
     2439          ALLOCATE ( surf_lsm_v(l)%r_canopy(1:surf_lsm_v(l)%ns)         )
     2440          ALLOCATE ( surf_lsm_v(l)%r_soil(1:surf_lsm_v(l)%ns)           )
     2441          ALLOCATE ( surf_lsm_v(l)%r_soil_min(1:surf_lsm_v(l)%ns)       )
     2442          ALLOCATE ( surf_lsm_v(l)%r_s(1:surf_lsm_v(l)%ns)              )
     2443          ALLOCATE ( surf_lsm_v(l)%r_canopy_min(1:surf_lsm_v(l)%ns)     )
     2444          ALLOCATE ( surf_lsm_v(l)%shf_eb(1:surf_lsm_v(l)%ns)           )
     2445          ALLOCATE ( surf_lsm_v(l)%soil_type_2d(1:surf_lsm_v(l)%ns)     )
     2446          ALLOCATE ( surf_lsm_v(l)%veg_type_2d(1:surf_lsm_v(l)%ns)      )
     2447          ALLOCATE ( surf_lsm_v(l)%water_surface(1:surf_lsm_v(l)%ns)    )
     2448
     2449          surf_lsm_v(l)%water_surface = .FALSE.
     2450          surf_lsm_v(l)%pave_surface  = .FALSE.
     2451       ENDDO
     2452   
    20162453#if ! defined( __nopointer )
    20172454!
    20182455!--    Initial assignment of the pointers
    2019        t_soil    => t_soil_1;    t_soil_p    => t_soil_2
    2020        t_surface => t_surface_1; t_surface_p => t_surface_2
    2021        m_soil    => m_soil_1;    m_soil_p    => m_soil_2
    2022        m_liq_eb  => m_liq_eb_1;  m_liq_eb_p  => m_liq_eb_2
     2456!--    Horizontal surfaces
     2457       t_soil_h    => t_soil_h_1;    t_soil_h_p    => t_soil_h_2
     2458       t_surface_h => t_surface_h_1; t_surface_h_p => t_surface_h_2
     2459       m_soil_h    => m_soil_h_1;    m_soil_h_p    => m_soil_h_2
     2460       m_liq_eb_h  => m_liq_eb_h_1;  m_liq_eb_h_p  => m_liq_eb_h_2
     2461!
     2462!--    Vertical surfaces
     2463       t_soil_v    => t_soil_v_1;    t_soil_v_p    => t_soil_v_2
     2464       t_surface_v => t_surface_v_1; t_surface_v_p => t_surface_v_2
     2465       m_soil_v    => m_soil_v_1;    m_soil_v_p    => m_soil_v_2
     2466       m_liq_eb_v  => m_liq_eb_v_1;  m_liq_eb_v_p  => m_liq_eb_v_2
    20232467#endif
    20242468
     
    20432487                                  conserve_water_content,                      &
    20442488                                  f_shortwave_incoming, field_capacity,        &
    2045                                   hydraulic_conductivity,                      &
     2489                                  aero_resist_kray, hydraulic_conductivity,                      &
    20462490                                  lambda_surface_stable,                       &
    20472491                                  lambda_surface_unstable, leaf_area_index,    &
     
    20872531!> temperature and water content.
    20882532!------------------------------------------------------------------------------!
    2089     SUBROUTINE lsm_soil_model
     2533    SUBROUTINE lsm_soil_model( horizontal, l )
    20902534
    20912535
    20922536       IMPLICIT NONE
    20932537
    2094        INTEGER(iwp) ::  i   !< running index
    2095        INTEGER(iwp) ::  j   !< running index
    2096        INTEGER(iwp) ::  k   !< running index
    2097 
    2098        REAL(wp)     :: h_vg !< Van Genuchten coef. h
     2538       INTEGER(iwp) ::  k       !< running index
     2539       INTEGER(iwp) ::  l       !< surface-data type index indication facing
     2540       INTEGER(iwp) ::  m       !< running index
     2541
     2542       LOGICAL      ::  horizontal !< flag indication horizontal wall, required to set pointer accordingly
     2543
     2544       REAL(wp)     ::  h_vg !< Van Genuchten coef. h
    20992545
    21002546       REAL(wp), DIMENSION(nzb_soil:nzt_soil) :: gamma_temp,  & !< temp. gamma
     
    21022548                                                 tend           !< tendency
    21032549
    2104        DO  i = nxlg, nxrg
    2105           DO  j = nysg, nyng
    2106 
    2107              IF ( pave_surface(j,i) )  THEN
    2108                 rho_c_total(nzb_soil,j,i) = pave_heat_capacity
    2109                 lambda_temp(nzb_soil)     = pave_heat_conductivity
     2550       TYPE(surf_type_lsm), POINTER ::  surf_m_soil
     2551       TYPE(surf_type_lsm), POINTER ::  surf_m_soil_p
     2552       TYPE(surf_type_lsm), POINTER ::  surf_t_soil
     2553       TYPE(surf_type_lsm), POINTER ::  surf_t_soil_p
     2554       TYPE(surf_type_lsm), POINTER ::  surf_tm_soil_m
     2555       TYPE(surf_type_lsm), POINTER ::  surf_tt_soil_m
     2556
     2557       TYPE(surf_type), POINTER  ::  surf  !< surface-date type variable
     2558
     2559       IF ( horizontal )  THEN
     2560          surf           => surf_lsm_h
     2561
     2562          surf_m_soil    => m_soil_h
     2563          surf_m_soil_p  => m_soil_h_p
     2564          surf_t_soil    => t_soil_h
     2565          surf_t_soil_p  => t_soil_h_p
     2566          surf_tm_soil_m => tm_soil_h_m
     2567          surf_tt_soil_m => tt_soil_h_m
     2568       ELSE
     2569          surf           => surf_lsm_v(l)
     2570
     2571          surf_m_soil    => m_soil_v(l)
     2572          surf_m_soil_p  => m_soil_v_p(l)
     2573          surf_t_soil    => t_soil_v(l)
     2574          surf_t_soil_p  => t_soil_v_p(l)
     2575          surf_tm_soil_m => tm_soil_v_m(l)
     2576          surf_tt_soil_m => tt_soil_v_m(l)
     2577       ENDIF
     2578
     2579       DO  m = 1, surf%ns
     2580
     2581          IF ( surf%pave_surface(m) )  THEN
     2582             surf%rho_c_total(nzb_soil,m) = pave_heat_capacity
     2583             lambda_temp(nzb_soil)          = pave_heat_conductivity
     2584          ENDIF
     2585
     2586          IF (  .NOT.  surf%water_surface(m) )  THEN
     2587             DO  k = nzb_soil, nzt_soil
     2588
     2589                IF ( surf%pave_surface(m)  .AND.  zs(k) <= pave_depth )  THEN
     2590                   
     2591                   surf%rho_c_total(k,m) = pave_heat_capacity
     2592                   lambda_temp(k)        = pave_heat_conductivity   
     2593
     2594                ELSE           
     2595!
     2596!--                Calculate volumetric heat capacity of the soil, taking
     2597!--                into account water content
     2598                   surf%rho_c_total(k,m) = (rho_c_soil *                 &
     2599                                               ( 1.0_wp - surf%m_sat(m) )&
     2600                                               + rho_c_water * surf_m_soil%var_2d(k,m) )
     2601
     2602!
     2603!--                Calculate soil heat conductivity at the center of the soil
     2604!--                layers
     2605                   lambda_h_sat = lambda_h_sm**(1.0_wp - surf%m_sat(m)) *&
     2606                                  lambda_h_water ** surf_m_soil%var_2d(k,m)
     2607
     2608                   ke = 1.0_wp + LOG10( MAX( 0.1_wp, surf_m_soil%var_2d(k,m) /             &
     2609                                                     surf%m_sat(m) ) )
     2610
     2611                   lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) +       &
     2612                                    lambda_h_dry
     2613                ENDIF
     2614             ENDDO
     2615
     2616!
     2617!--          Calculate soil heat conductivity (lambda_h) at the _stag level
     2618!--          using linear interpolation. For pavement surface, the
     2619!--          true pavement depth is considered
     2620             DO  k = nzb_soil, nzt_soil-1
     2621                IF ( surf%pave_surface(m)  .AND.  zs(k)   < pave_depth   &
     2622                                                 .AND.  zs(k+1) > pave_depth )  THEN
     2623                   surf%lambda_h(k,m) = ( pave_depth - zs(k) ) / dz_soil(k+1)&
     2624                                     * lambda_temp(k)                          &
     2625                                     + ( 1.0_wp - ( pave_depth - zs(k) )       &
     2626                                     / dz_soil(k+1) ) * lambda_temp(k+1)
     2627                ELSE
     2628                   surf%lambda_h(k,m) = ( lambda_temp(k+1) + lambda_temp(k) )&
     2629                                     * 0.5_wp
     2630                ENDIF
     2631             ENDDO
     2632             surf%lambda_h(nzt_soil,m) = lambda_temp(nzt_soil)
     2633
     2634!
     2635!--          Prognostic equation for soil temperature t_soil
     2636             tend(:) = 0.0_wp
     2637
     2638             tend(nzb_soil) = ( 1.0_wp / surf%rho_c_total(nzb_soil,m) ) *&
     2639                    ( surf%lambda_h(nzb_soil,m) * ( surf_t_soil%var_2d(nzb_soil+1,m) &
     2640                      - surf_t_soil%var_2d(nzb_soil,m) ) * ddz_soil(nzb_soil+1)            &
     2641                      + surf%ghf_eb(m) ) * ddz_soil_stag(nzb_soil)
     2642
     2643             DO  k = nzb_soil+1, nzt_soil
     2644                tend(k) = ( 1.0_wp / surf%rho_c_total(k,m) )             &
     2645                          * (   surf%lambda_h(k,m)                       &
     2646                              * ( surf_t_soil%var_2d(k+1,m) - surf_t_soil%var_2d(k,m) )                &
     2647                              * ddz_soil(k+1)                                  &
     2648                              - surf%lambda_h(k-1,m)                     &
     2649                              * ( surf_t_soil%var_2d(k,m) - surf_t_soil%var_2d(k-1,m) )                &
     2650                              * ddz_soil(k)                                    &
     2651                            ) * ddz_soil_stag(k)
     2652
     2653             ENDDO
     2654
     2655             surf_t_soil_p%var_2d(nzb_soil:nzt_soil,m) = surf_t_soil%var_2d(nzb_soil:nzt_soil,m)       &
     2656                                               + dt_3d * ( tsc(2)              &
     2657                                               * tend(nzb_soil:nzt_soil)       &
     2658                                               + tsc(3)                        &
     2659                                               * surf_tt_soil_m%var_2d(nzb_soil:nzt_soil,m) )
     2660
     2661!
     2662!--          Calculate t_soil tendencies for the next Runge-Kutta step
     2663             IF ( timestep_scheme(1:5) == 'runge' )  THEN
     2664                IF ( intermediate_timestep_count == 1 )  THEN
     2665                   DO  k = nzb_soil, nzt_soil
     2666                      surf_tt_soil_m%var_2d(k,m) = tend(k)
     2667                   ENDDO
     2668                ELSEIF ( intermediate_timestep_count <                         &
     2669                         intermediate_timestep_count_max )  THEN
     2670                   DO  k = nzb_soil, nzt_soil
     2671                      surf_tt_soil_m%var_2d(k,m) = -9.5625_wp * tend(k) + 5.3125_wp        &
     2672                                         * surf_tt_soil_m%var_2d(k,m)
     2673                   ENDDO
     2674                ENDIF
    21102675             ENDIF
    21112676
    2112              IF (  .NOT.  water_surface(j,i) )  THEN
     2677
     2678             DO  k = nzb_soil, nzt_soil
     2679
     2680!
     2681!--             Calculate soil diffusivity at the center of the soil layers
     2682                lambda_temp(k) = (- b_ch * surf%gamma_w_sat(m) * psi_sat     &
     2683                                  / surf%m_sat(m) ) * ( MAX( surf_m_soil%var_2d(k,m),    &
     2684                                  surf%m_wilt(m) ) / surf%m_sat(m) )**(&
     2685                                                            b_ch + 2.0_wp )
     2686
     2687!
     2688!--             Parametrization of Van Genuchten
     2689                IF ( soil_type /= 7 )  THEN
     2690!
     2691!--                Calculate the hydraulic conductivity after Van Genuchten
     2692!--                (1980)
     2693                   h_vg = ( ( ( surf%m_res(m) - surf%m_sat(m) ) /  &
     2694                              ( surf%m_res(m) -                          &
     2695                                MAX( surf_m_soil%var_2d(k,m), surf%m_wilt(m) )       &
     2696                              )                                                &
     2697                            )**(                                               &
     2698                          surf%n_vg(m) / ( surf%n_vg(m) - 1.0_wp ) &
     2699                               ) - 1.0_wp                                      &
     2700                          )**( 1.0_wp / surf%n_vg(m) ) / surf%alpha_vg(m)
     2701
     2702                   gamma_temp(k) = surf%gamma_w_sat(m) * ( ( ( 1.0_wp +  &
     2703                          ( surf%alpha_vg(m) * h_vg )**surf%n_vg(m)&
     2704                                                               )**(            &
     2705                              1.0_wp - 1.0_wp / surf%n_vg(m)) - (        &
     2706                          surf%alpha_vg(m) * h_vg )**( surf%n_vg(m)&
     2707                              - 1.0_wp) )**2 )                                 &
     2708                              / ( ( 1.0_wp + ( surf%alpha_vg(m) * h_vg   &
     2709                              )**surf%n_vg(m) )**( ( 1.0_wp  - 1.0_wp    &
     2710                              / surf%n_vg(m) ) *                         &
     2711                              ( surf%l_vg(m) + 2.0_wp) ) )
     2712!
     2713!--             Parametrization of Clapp & Hornberger
     2714                ELSE
     2715                   gamma_temp(k) = surf%gamma_w_sat(m) * ( surf_m_soil%var_2d(k,m)   &
     2716                                   / surf%m_sat(m) )**(2.0_wp * b_ch + 3.0_wp)
     2717                ENDIF
     2718
     2719             ENDDO
     2720
     2721!
     2722!--          Prognostic equation for soil moisture content. Only performed,
     2723!--          when humidity is enabled in the atmosphere and the surface type
     2724!--          is not pavement (implies dry soil below).
     2725             IF ( humidity  .AND.  .NOT.  surf%pave_surface(m) )  THEN
     2726!
     2727!--             Calculate soil diffusivity (lambda_w) at the _stag level
     2728!--             using linear interpolation. To do: replace this with
     2729!--             ECMWF-IFS Eq. 8.81
     2730                DO  k = nzb_soil, nzt_soil-1
     2731                   
     2732                   surf%lambda_w(k,m) = ( lambda_temp(k+1) +             &
     2733                                                lambda_temp(k) ) * 0.5_wp
     2734                   surf%gamma_w(k,m)  = ( gamma_temp(k+1)  +             &
     2735                                                gamma_temp(k) )  * 0.5_wp
     2736
     2737                ENDDO
     2738
     2739!
     2740!
     2741!--             In case of a closed bottom (= water content is conserved),
     2742!--             set hydraulic conductivity to zero to that no water will be
     2743!--             lost in the bottom layer.
     2744                IF ( conserve_water_content )  THEN
     2745                   surf%gamma_w(nzt_soil,m) = 0.0_wp
     2746                ELSE
     2747                   surf%gamma_w(nzt_soil,m) = gamma_temp(nzt_soil)
     2748                ENDIF     
     2749
     2750!--             The root extraction (= root_extr * qsws_veg_eb / (rho_l     
     2751!--             * l_v)) ensures the mass conservation for water. The         
     2752!--             transpiration of plants equals the cumulative withdrawals by
     2753!--             the roots in the soil. The scheme takes into account the
     2754!--             availability of water in the soil layers as well as the root
     2755!--             fraction in the respective layer. Layer with moisture below
     2756!--             wilting point will not contribute, which reflects the
     2757!--             preference of plants to take water from moister layers.
     2758!
     2759!--             Calculate the root extraction (ECMWF 7.69, the sum of
     2760!--             root_extr = 1). The energy balance solver guarantees a
     2761!--             positive transpiration, so that there is no need for an
     2762!--             additional check.
     2763                m_total = 0.0_wp
    21132764                DO  k = nzb_soil, nzt_soil
    2114 
    2115 
    2116                    IF ( pave_surface(j,i)  .AND.  zs(k) <= pave_depth )  THEN
    2117                    
    2118                       rho_c_total(k,j,i) = pave_heat_capacity
    2119                       lambda_temp(k)     = pave_heat_conductivity   
    2120 
    2121                    ELSE           
    2122 !
    2123 !--                   Calculate volumetric heat capacity of the soil, taking
    2124 !--                   into account water content
    2125                       rho_c_total(k,j,i) = (rho_c_soil * (1.0_wp - m_sat(j,i)) &
    2126                                            + rho_c_water * m_soil(k,j,i))
    2127 
    2128 !
    2129 !--                   Calculate soil heat conductivity at the center of the soil
    2130 !--                   layers
    2131                       lambda_h_sat = lambda_h_sm ** (1.0_wp - m_sat(j,i)) *    &
    2132                                      lambda_h_water ** m_soil(k,j,i)
    2133 
    2134                       ke = 1.0_wp + LOG10(MAX(0.1_wp,m_soil(k,j,i)             &
    2135                                                      / m_sat(j,i)))
    2136 
    2137                       lambda_temp(k) = ke * (lambda_h_sat - lambda_h_dry) +    &
    2138                                        lambda_h_dry
    2139                    ENDIF
    2140 
     2765                    IF ( surf_m_soil%var_2d(k,m) > surf%m_wilt(m) )  THEN
     2766                       m_total = m_total + surf%root_fr(k,m) * surf_m_soil%var_2d(k,m)
     2767                    ENDIF
     2768                ENDDO 
     2769                 IF ( m_total > 0.0_wp )  THEN
     2770                   DO  k = nzb_soil, nzt_soil
     2771                      IF ( surf_m_soil%var_2d(k,m) > surf%m_wilt(m) )  THEN
     2772                         root_extr(k) = surf%root_fr(k,m) * surf_m_soil%var_2d(k,m)  &
     2773                                                            / m_total
     2774                      ELSE
     2775                         root_extr(k) = 0.0_wp
     2776                      ENDIF
     2777                   ENDDO
     2778                ENDIF
     2779!
     2780!--             Prognostic equation for soil water content m_soil_h.
     2781                tend(:) = 0.0_wp
     2782
     2783                tend(nzb_soil) = ( surf%lambda_w(nzb_soil,m) *   (       &
     2784                         surf_m_soil%var_2d(nzb_soil+1,m) - surf_m_soil%var_2d(nzb_soil,m) )           &
     2785                         * ddz_soil(nzb_soil+1) - surf%gamma_w(nzb_soil,m) - &
     2786                                                     (                         &
     2787                            root_extr(nzb_soil) * surf%qsws_veg_eb(m)    &
     2788                            + surf%qsws_soil_eb(m) ) * drho_l_lv )       &
     2789                            * ddz_soil_stag(nzb_soil)
     2790
     2791                DO  k = nzb_soil+1, nzt_soil-1
     2792                   tend(k) = ( surf%lambda_w(k,m) * ( surf_m_soil%var_2d(k+1,m)      &
     2793                             - surf_m_soil%var_2d(k,m) ) * ddz_soil(k+1)                   &
     2794                             - surf%gamma_w(k,m)                         &
     2795                             - surf%lambda_w(k-1,m) * ( surf_m_soil%var_2d(k,m) -    &
     2796                             surf_m_soil%var_2d(k-1,m)) * ddz_soil(k)                      &
     2797                             + surf%gamma_w(k-1,m) - (root_extr(k)       &
     2798                             * surf%qsws_veg_eb(m) * drho_l_lv)          &
     2799                             ) * ddz_soil_stag(k)
    21412800                ENDDO
    2142 
    2143 !
    2144 !--             Calculate soil heat conductivity (lambda_h) at the _stag level
    2145 !--             using linear interpolation. For pavement surface, the
    2146 !--             true pavement depth is considered
    2147                 DO  k = nzb_soil, nzt_soil-1
    2148                    IF ( pave_surface(j,i)  .AND.  zs(k)   < pave_depth         &
    2149                                            .AND.  zs(k+1) > pave_depth )  THEN
    2150                       lambda_h(k,j,i) = ( pave_depth - zs(k) ) / dz_soil(k+1)  &
    2151                                         * lambda_temp(k)                       &
    2152                                         + ( 1.0_wp - ( pave_depth - zs(k) )    &
    2153                                         / dz_soil(k+1) ) * lambda_temp(k+1)
    2154                    ELSE
    2155                       lambda_h(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) )  &
    2156                                         * 0.5_wp
    2157                    ENDIF
     2801                tend(nzt_soil) = ( - surf%gamma_w(nzt_soil,m)            &
     2802                                     - surf%lambda_w(nzt_soil-1,m)       &
     2803                                     * ( surf_m_soil%var_2d(nzt_soil,m)                    &
     2804                                     - surf_m_soil%var_2d(nzt_soil-1,m))                   &
     2805                                     * ddz_soil(nzt_soil)                      &
     2806                                     + surf%gamma_w(nzt_soil-1,m) - (    &
     2807                                       root_extr(nzt_soil)                     &
     2808                                     * surf%qsws_veg_eb(m) * drho_l_lv ) &
     2809                                  ) * ddz_soil_stag(nzt_soil)             
     2810
     2811                surf_m_soil_p%var_2d(nzb_soil:nzt_soil,m) = surf_m_soil%var_2d(nzb_soil:nzt_soil,m)    &
     2812                                                + dt_3d * ( tsc(2) * tend(:)   &
     2813                                                + tsc(3) * surf_tm_soil_m%var_2d(:,m) )   
     2814   
     2815!
     2816!--             Account for dry soils (find a better solution here!)
     2817                DO  k = nzb_soil, nzt_soil
     2818                   IF ( surf_m_soil_p%var_2d(k,m) < 0.0_wp )  surf_m_soil_p%var_2d(k,m) = 0.0_wp
    21582819                ENDDO
    2159                 lambda_h(nzt_soil,j,i) = lambda_temp(nzt_soil)
    2160 
    2161 
    2162 
    2163 
    2164 !
    2165 !--             Prognostic equation for soil temperature t_soil
    2166                 tend(:) = 0.0_wp
    2167 
    2168                 tend(nzb_soil) = (1.0_wp/rho_c_total(nzb_soil,j,i)) *          &
    2169                           ( lambda_h(nzb_soil,j,i) * ( t_soil(nzb_soil+1,j,i)  &
    2170                             - t_soil(nzb_soil,j,i) ) * ddz_soil(nzb_soil+1)    &
    2171                             + ghf_eb(j,i) ) * ddz_soil_stag(nzb_soil)
    2172 
    2173                 DO  k = nzb_soil+1, nzt_soil
    2174                    tend(k) = (1.0_wp/rho_c_total(k,j,i))                       &
    2175                              * (   lambda_h(k,j,i)                             &
    2176                                  * ( t_soil(k+1,j,i) - t_soil(k,j,i) )         &
    2177                                  * ddz_soil(k+1)                               &
    2178                                  - lambda_h(k-1,j,i)                           &
    2179                                  * ( t_soil(k,j,i) - t_soil(k-1,j,i) )         &
    2180                                  * ddz_soil(k)                                 &
    2181                                ) * ddz_soil_stag(k)
    2182 
    2183                 ENDDO
    2184 
    2185                 t_soil_p(nzb_soil:nzt_soil,j,i) = t_soil(nzb_soil:nzt_soil,j,i)&
    2186                                                   + dt_3d * ( tsc(2)           &
    2187                                                   * tend(nzb_soil:nzt_soil)    &
    2188                                                   + tsc(3)                     &
    2189                                                   * tt_soil_m(:,j,i) )   
    2190 
    2191 !
    2192 !--             Calculate t_soil tendencies for the next Runge-Kutta step
     2820
     2821!
     2822!--             Calculate m_soil tendencies for the next Runge-Kutta step
    21932823                IF ( timestep_scheme(1:5) == 'runge' )  THEN
    21942824                   IF ( intermediate_timestep_count == 1 )  THEN
    21952825                      DO  k = nzb_soil, nzt_soil
    2196                          tt_soil_m(k,j,i) = tend(k)
     2826                         surf_tm_soil_m%var_2d(k,m) = tend(k)
    21972827                      ENDDO
    21982828                   ELSEIF ( intermediate_timestep_count <                      &
    21992829                            intermediate_timestep_count_max )  THEN
    22002830                      DO  k = nzb_soil, nzt_soil
    2201                          tt_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp   &
    2202                                          * tt_soil_m(k,j,i)
     2831                         surf_tm_soil_m%var_2d(k,m) = -9.5625_wp * tend(k) + 5.3125_wp     &
     2832                                  * surf_tm_soil_m%var_2d(k,m)
    22032833                      ENDDO
    22042834                   ENDIF
    22052835                ENDIF
    2206 
    2207 
    2208                 DO  k = nzb_soil, nzt_soil
    2209 
    2210 !
    2211 !--                Calculate soil diffusivity at the center of the soil layers
    2212                    lambda_temp(k) = (- b_ch * gamma_w_sat(j,i) * psi_sat       &
    2213                                      / m_sat(j,i) ) * ( MAX( m_soil(k,j,i),    &
    2214                                      m_wilt(j,i) ) / m_sat(j,i) )**(           &
    2215                                      b_ch + 2.0_wp )
    2216 
    2217 !
    2218 !--                Parametrization of Van Genuchten
    2219                    IF ( soil_type /= 7 )  THEN
    2220 !
    2221 !--                   Calculate the hydraulic conductivity after Van Genuchten
    2222 !--                   (1980)
    2223                       h_vg = ( ( (m_res(j,i) - m_sat(j,i)) / ( m_res(j,i) -    &
    2224                                  MAX( m_soil(k,j,i), m_wilt(j,i) ) ) )**(      &
    2225                                  n_vg(j,i) / (n_vg(j,i) - 1.0_wp ) ) - 1.0_wp  &
    2226                              )**( 1.0_wp / n_vg(j,i) ) / alpha_vg(j,i)
    2227 
    2228 
    2229                       gamma_temp(k) = gamma_w_sat(j,i) * ( ( (1.0_wp +         &
    2230                                       ( alpha_vg(j,i) * h_vg )**n_vg(j,i))**(  &
    2231                                       1.0_wp - 1.0_wp / n_vg(j,i) ) - (        &
    2232                                       alpha_vg(j,i) * h_vg )**( n_vg(j,i)      &
    2233                                       - 1.0_wp) )**2 )                         &
    2234                                       / ( ( 1.0_wp + ( alpha_vg(j,i) * h_vg    &
    2235                                       )**n_vg(j,i) )**( ( 1.0_wp  - 1.0_wp     &
    2236                                       / n_vg(j,i) ) *( l_vg(j,i) + 2.0_wp) ) )
    2237 
    2238 !
    2239 !--                Parametrization of Clapp & Hornberger
    2240                    ELSE
    2241                       gamma_temp(k) = gamma_w_sat(j,i) * ( m_soil(k,j,i)       &
    2242                                       / m_sat(j,i) )**(2.0_wp * b_ch + 3.0_wp)
    2243                    ENDIF
    2244 
    2245                 ENDDO
    2246 
    2247 !
    2248 !--             Prognostic equation for soil moisture content. Only performed,
    2249 !--             when humidity is enabled in the atmosphere and the surface type
    2250 !--             is not pavement (implies dry soil below).
    2251                 IF ( humidity  .AND.  .NOT.  pave_surface(j,i) )  THEN
    2252 !
    2253 !--                Calculate soil diffusivity (lambda_w) at the _stag level
    2254 !--                using linear interpolation. To do: replace this with
    2255 !--                ECMWF-IFS Eq. 8.81
    2256                    DO  k = nzb_soil, nzt_soil-1
    2257                      
    2258                       lambda_w(k,j,i) = ( lambda_temp(k+1) + lambda_temp(k) )  &
    2259                                         * 0.5_wp
    2260                       gamma_w(k,j,i)  = ( gamma_temp(k+1) + gamma_temp(k) )    &
    2261                                         * 0.5_wp
    2262 
    2263                    ENDDO
    2264 
    2265 !
    2266 !
    2267 !--                In case of a closed bottom (= water content is conserved),
    2268 !--                set hydraulic conductivity to zero to that no water will be
    2269 !--                lost in the bottom layer.
    2270                    IF ( conserve_water_content )  THEN
    2271                       gamma_w(nzt_soil,j,i) = 0.0_wp
    2272                    ELSE
    2273                       gamma_w(nzt_soil,j,i) = gamma_temp(nzt_soil)
    2274                    ENDIF     
    2275 
    2276 !--                The root extraction (= root_extr * qsws_veg_eb / (rho_l     
    2277 !--                * l_v)) ensures the mass conservation for water. The         
    2278 !--                transpiration of plants equals the cumulative withdrawals by
    2279 !--                the roots in the soil. The scheme takes into account the
    2280 !--                availability of water in the soil layers as well as the root
    2281 !--                fraction in the respective layer. Layer with moisture below
    2282 !--                wilting point will not contribute, which reflects the
    2283 !--                preference of plants to take water from moister layers.
    2284 
    2285 !
    2286 !--                Calculate the root extraction (ECMWF 7.69, the sum of
    2287 !--                root_extr = 1). The energy balance solver guarantees a
    2288 !--                positive transpiration, so that there is no need for an
    2289 !--                additional check.
    2290                    m_total = 0.0_wp
    2291                    DO  k = nzb_soil, nzt_soil
    2292                        IF ( m_soil(k,j,i) > m_wilt(j,i) )  THEN
    2293                           m_total = m_total + root_fr(k,j,i) * m_soil(k,j,i)
    2294                        ENDIF
    2295                    ENDDO 
    2296 
    2297                    IF ( m_total > 0.0_wp )  THEN
    2298                       DO  k = nzb_soil, nzt_soil
    2299                          IF ( m_soil(k,j,i) > m_wilt(j,i) )  THEN
    2300                             root_extr(k) = root_fr(k,j,i) * m_soil(k,j,i)      &
    2301                                                             / m_total
    2302                          ELSE
    2303                             root_extr(k) = 0.0_wp
    2304                          ENDIF
    2305                       ENDDO
    2306                    ENDIF
    2307 
    2308 !
    2309 !--                Prognostic equation for soil water content m_soil.
    2310                    tend(:) = 0.0_wp
    2311 
    2312                    tend(nzb_soil) = ( lambda_w(nzb_soil,j,i) * (               &
    2313                             m_soil(nzb_soil+1,j,i) - m_soil(nzb_soil,j,i) )    &
    2314                             * ddz_soil(nzb_soil+1) - gamma_w(nzb_soil,j,i) - ( &
    2315                                root_extr(nzb_soil) * qsws_veg_eb(j,i)          &
    2316                                + qsws_soil_eb(j,i) ) * drho_l_lv )             &
    2317                                * ddz_soil_stag(nzb_soil)
    2318 
    2319                    DO  k = nzb_soil+1, nzt_soil-1
    2320                       tend(k) = ( lambda_w(k,j,i) * ( m_soil(k+1,j,i)          &
    2321                                 - m_soil(k,j,i) ) * ddz_soil(k+1)              &
    2322                                 - gamma_w(k,j,i)                               &
    2323                                 - lambda_w(k-1,j,i) * (m_soil(k,j,i) -         &
    2324                                 m_soil(k-1,j,i)) * ddz_soil(k)                 &
    2325                                 + gamma_w(k-1,j,i) - (root_extr(k)             &
    2326                                 * qsws_veg_eb(j,i) * drho_l_lv)                &
    2327                                 ) * ddz_soil_stag(k)
    2328 
    2329                    ENDDO
    2330                    tend(nzt_soil) = ( - gamma_w(nzt_soil,j,i)                  &
    2331                                            - lambda_w(nzt_soil-1,j,i)          &
    2332                                            * (m_soil(nzt_soil,j,i)             &
    2333                                            - m_soil(nzt_soil-1,j,i))           &
    2334                                            * ddz_soil(nzt_soil)                &
    2335                                            + gamma_w(nzt_soil-1,j,i) - (       &
    2336                                              root_extr(nzt_soil)               &
    2337                                            * qsws_veg_eb(j,i) * drho_l_lv  )   &
    2338                                      ) * ddz_soil_stag(nzt_soil)             
    2339 
    2340                    m_soil_p(nzb_soil:nzt_soil,j,i) = m_soil(nzb_soil:nzt_soil,j,i)&
    2341                                                    + dt_3d * ( tsc(2) * tend(:)   &
    2342                                                    + tsc(3) * tm_soil_m(:,j,i) )   
    2343    
    2344 !
    2345 !--                Account for dry soils (find a better solution here!)
    2346                    DO  k = nzb_soil, nzt_soil
    2347                       IF ( m_soil_p(k,j,i) < 0.0_wp )  m_soil_p(k,j,i) = 0.0_wp
    2348                    ENDDO
    2349 
    2350 !
    2351 !--                Calculate m_soil tendencies for the next Runge-Kutta step
    2352                    IF ( timestep_scheme(1:5) == 'runge' )  THEN
    2353                       IF ( intermediate_timestep_count == 1 )  THEN
    2354                          DO  k = nzb_soil, nzt_soil
    2355                             tm_soil_m(k,j,i) = tend(k)
    2356                          ENDDO
    2357                       ELSEIF ( intermediate_timestep_count <                   &
    2358                                intermediate_timestep_count_max )  THEN
    2359                          DO  k = nzb_soil, nzt_soil
    2360                             tm_soil_m(k,j,i) = -9.5625_wp * tend(k) + 5.3125_wp&
    2361                                      * tm_soil_m(k,j,i)
    2362                          ENDDO
    2363                       ENDIF
    2364                    ENDIF
    2365 
    2366                 ENDIF
    2367 
    23682836             ENDIF
    23692837
    2370           ENDDO
     2838          ENDIF
     2839
    23712840       ENDDO
    23722841
     
    23862855
    23872856#if defined( __nopointer )
    2388 
    2389        t_surface    = t_surface_p
    2390        t_soil       = t_soil_p
     2857!
     2858!--    Horizontal surfaces
     2859       t_surface_h  = t_surface_h_p
     2860       t_soil_h     = t_soil_h_p
    23912861       IF ( humidity )  THEN
    2392           m_soil    = m_soil_p
    2393           m_liq_eb  = m_liq_eb_p
     2862          m_soil_h    = m_soil_h_p
     2863          m_liq_eb_h  = m_liq_eb_h_p
     2864       ENDIF
     2865!
     2866!--    Vertical surfaces
     2867       t_surface_v  = t_surface_v_p
     2868       t_soil_v     = t_soil_v_p
     2869       IF ( humidity )  THEN
     2870          m_soil_v    = m_soil_v_p
     2871          m_liq_eb_v  = m_liq_eb_v_p
    23942872       ENDIF
    23952873
     
    23992877
    24002878          CASE ( 0 )
    2401 
    2402              t_surface  => t_surface_1; t_surface_p  => t_surface_2
    2403              t_soil     => t_soil_1;    t_soil_p     => t_soil_2
     2879!
     2880!--          Horizontal surfaces
     2881             t_surface_h  => t_surface_h_1; t_surface_h_p  => t_surface_h_2
     2882             t_soil_h     => t_soil_h_1;    t_soil_h_p     => t_soil_h_2
    24042883             IF ( humidity )  THEN
    2405                 m_soil    => m_soil_1;   m_soil_p    => m_soil_2
    2406                 m_liq_eb  => m_liq_eb_1; m_liq_eb_p  => m_liq_eb_2
     2884                m_soil_h    => m_soil_h_1;   m_soil_h_p    => m_soil_h_2
     2885                m_liq_eb_h  => m_liq_eb_h_1; m_liq_eb_h_p  => m_liq_eb_h_2
    24072886             ENDIF
     2887!
     2888!--          Vertical surfaces
     2889             t_surface_v  => t_surface_v_1; t_surface_v_p  => t_surface_v_2
     2890             t_soil_v     => t_soil_v_1;    t_soil_v_p     => t_soil_v_2
     2891             IF ( humidity )  THEN
     2892                m_soil_v    => m_soil_v_1;   m_soil_v_p    => m_soil_v_2
     2893                m_liq_eb_v  => m_liq_eb_v_1; m_liq_eb_v_p  => m_liq_eb_v_2
     2894             ENDIF
     2895
    24082896
    24092897
    24102898          CASE ( 1 )
    2411 
    2412              t_surface  => t_surface_2; t_surface_p  => t_surface_1
    2413              t_soil     => t_soil_2;    t_soil_p     => t_soil_1
     2899!
     2900!--          Horizontal surfaces
     2901             t_surface_h  => t_surface_h_2; t_surface_h_p  => t_surface_h_1
     2902             t_soil_h     => t_soil_h_2;    t_soil_h_p     => t_soil_h_1
    24142903             IF ( humidity )  THEN
    2415                 m_soil    => m_soil_2;   m_soil_p    => m_soil_1
    2416                 m_liq_eb  => m_liq_eb_2; m_liq_eb_p  => m_liq_eb_1
     2904                m_soil_h    => m_soil_h_2;   m_soil_h_p    => m_soil_h_1
     2905                m_liq_eb_h  => m_liq_eb_h_2; m_liq_eb_h_p  => m_liq_eb_h_1
     2906             ENDIF
     2907!
     2908!--          Vertical surfaces
     2909             t_surface_v  => t_surface_v_2; t_surface_v_p  => t_surface_v_1
     2910             t_soil_v     => t_soil_v_2;    t_soil_v_p     => t_soil_v_1
     2911             IF ( humidity )  THEN
     2912                m_soil_v    => m_soil_v_2;   m_soil_v_p    => m_soil_v_1
     2913                m_liq_eb_v  => m_liq_eb_v_2; m_liq_eb_v_p  => m_liq_eb_v_1
    24172914             ENDIF
    24182915
     
    24452942    CHARACTER (LEN=*) :: variable !<
    24462943
    2447     INTEGER(iwp) ::  i !<
    2448     INTEGER(iwp) ::  j !<
    2449     INTEGER(iwp) ::  k !<
     2944    INTEGER(iwp) ::  i       !<
     2945    INTEGER(iwp) ::  j       !<
     2946    INTEGER(iwp) ::  k       !<
     2947    INTEGER(iwp) ::  m       !< running index
    24502948
    24512949    IF ( mode == 'allocate' )  THEN
     
    25533051
    25543052          CASE ( 'c_liq*' )
    2555              DO  i = nxlg, nxrg
    2556                 DO  j = nysg, nyng
    2557                    c_liq_av(j,i) = c_liq_av(j,i) + c_liq(j,i)
     3053             DO  m = 1, surf_lsm_h%ns
     3054                i   = surf_lsm_h%i(m)           
     3055                j   = surf_lsm_h%j(m)
     3056                c_liq_av(j,i) = c_liq_av(j,i) + surf_lsm_h%c_liq(m)
     3057             ENDDO
     3058
     3059          CASE ( 'c_soil*' )
     3060             DO  m = 1, surf_lsm_h%ns
     3061                i   = surf_lsm_h%i(m)           
     3062                j   = surf_lsm_h%j(m)
     3063                c_soil_av(j,i) = c_soil_av(j,i) + (1.0 - surf_lsm_h%c_veg(m))
     3064             ENDDO
     3065
     3066          CASE ( 'c_veg*' )
     3067             DO  m = 1, surf_lsm_h%ns
     3068                i   = surf_lsm_h%i(m)           
     3069                j   = surf_lsm_h%j(m)
     3070                c_veg_av(j,i) = c_veg_av(j,i) + surf_lsm_h%c_veg(m)
     3071             ENDDO
     3072
     3073          CASE ( 'ghf_eb*' )
     3074             DO  m = 1, surf_lsm_h%ns
     3075                i   = surf_lsm_h%i(m)           
     3076                j   = surf_lsm_h%j(m)
     3077                ghf_eb_av(j,i) = ghf_eb_av(j,i) + surf_lsm_h%ghf_eb(m)
     3078             ENDDO
     3079
     3080          CASE ( 'lai*' )
     3081             DO  m = 1, surf_lsm_h%ns
     3082                i   = surf_lsm_h%i(m)           
     3083                j   = surf_lsm_h%j(m)
     3084                lai_av(j,i) = lai_av(j,i) + surf_lsm_h%lai(m)
     3085             ENDDO
     3086
     3087          CASE ( 'm_liq_eb*' )
     3088             DO  m = 1, surf_lsm_h%ns
     3089                i   = surf_lsm_h%i(m)           
     3090                j   = surf_lsm_h%j(m)
     3091                m_liq_eb_av(j,i) = m_liq_eb_av(j,i) + m_liq_eb_h%var_1d(m)
     3092             ENDDO
     3093
     3094          CASE ( 'm_soil' )
     3095             DO  m = 1, surf_lsm_h%ns
     3096                i   = surf_lsm_h%i(m)           
     3097                j   = surf_lsm_h%j(m)
     3098                DO  k = nzb_soil, nzt_soil
     3099                   m_soil_av(k,j,i) = m_soil_av(k,j,i) + m_soil_h%var_2d(k,m)
    25583100                ENDDO
    25593101             ENDDO
    25603102
    2561           CASE ( 'c_soil*' )
    2562              DO  i = nxlg, nxrg
    2563                 DO  j = nysg, nyng
    2564                    c_soil_av(j,i) = c_soil_av(j,i) + (1.0 - c_veg(j,i))
    2565                 ENDDO
    2566              ENDDO
    2567 
    2568           CASE ( 'c_veg*' )
    2569              DO  i = nxlg, nxrg
    2570                 DO  j = nysg, nyng
    2571                    c_veg_av(j,i) = c_veg_av(j,i) + c_veg(j,i)
    2572                 ENDDO
    2573              ENDDO
    2574 
    2575           CASE ( 'ghf_eb*' )
    2576              DO  i = nxlg, nxrg
    2577                 DO  j = nysg, nyng
    2578                    ghf_eb_av(j,i) = ghf_eb_av(j,i) + ghf_eb(j,i)
    2579                 ENDDO
    2580              ENDDO
    2581 
    2582           CASE ( 'lai*' )
    2583              DO  i = nxlg, nxrg
    2584                 DO  j = nysg, nyng
    2585                    lai_av(j,i) = lai_av(j,i) + lai(j,i)
    2586                 ENDDO
    2587              ENDDO
    2588 
    2589           CASE ( 'm_liq_eb*' )
    2590              DO  i = nxlg, nxrg
    2591                 DO  j = nysg, nyng
    2592                    m_liq_eb_av(j,i) = m_liq_eb_av(j,i) + m_liq_eb(j,i)
    2593                 ENDDO
    2594              ENDDO
    2595 
    2596           CASE ( 'm_soil' )
    2597              DO  i = nxlg, nxrg
    2598                 DO  j = nysg, nyng
    2599                    DO  k = nzb_soil, nzt_soil
    2600                       m_soil_av(k,j,i) = m_soil_av(k,j,i) + m_soil(k,j,i)
    2601                    ENDDO
    2602                 ENDDO
    2603              ENDDO
    2604 
    26053103          CASE ( 'qsws_eb*' )
    2606              DO  i = nxlg, nxrg
    2607                 DO  j = nysg, nyng
    2608                    qsws_eb_av(j,i) = qsws_eb_av(j,i) + qsws_eb(j,i)
    2609                 ENDDO
     3104             DO  m = 1, surf_lsm_h%ns
     3105                i   = surf_lsm_h%i(m)           
     3106                j   = surf_lsm_h%j(m)
     3107                qsws_eb_av(j,i) = qsws_eb_av(j,i) + surf_lsm_h%qsws_eb(m)
    26103108             ENDDO
    26113109
    26123110          CASE ( 'qsws_liq_eb*' )
    2613              DO  i = nxlg, nxrg
    2614                 DO  j = nysg, nyng
    2615                    qsws_liq_eb_av(j,i) = qsws_liq_eb_av(j,i) + qsws_liq_eb(j,i)
    2616                 ENDDO
     3111             DO  m = 1, surf_lsm_h%ns
     3112                i   = surf_lsm_h%i(m)           
     3113                j   = surf_lsm_h%j(m)
     3114                qsws_liq_eb_av(j,i) = qsws_liq_eb_av(j,i) +                    &
     3115                                      surf_lsm_h%qsws_liq_eb(m)
    26173116             ENDDO
    26183117
    26193118          CASE ( 'qsws_soil_eb*' )
    2620              DO  i = nxlg, nxrg
    2621                 DO  j = nysg, nyng
    2622                    qsws_soil_eb_av(j,i) = qsws_soil_eb_av(j,i) + qsws_soil_eb(j,i)
    2623                 ENDDO
     3119             DO  m = 1, surf_lsm_h%ns
     3120                i   = surf_lsm_h%i(m)           
     3121                j   = surf_lsm_h%j(m)
     3122                qsws_soil_eb_av(j,i) = qsws_soil_eb_av(j,i) +                  &
     3123                                       surf_lsm_h%qsws_soil_eb(m)
    26243124             ENDDO
    26253125
    26263126          CASE ( 'qsws_veg_eb*' )
    2627              DO  i = nxlg, nxrg
    2628                 DO  j = nysg, nyng
    2629                    qsws_veg_eb_av(j,i) = qsws_veg_eb_av(j,i) + qsws_veg_eb(j,i)
    2630                 ENDDO
     3127             DO  m = 1, surf_lsm_h%ns
     3128                i   = surf_lsm_h%i(m)           
     3129                j   = surf_lsm_h%j(m)
     3130                qsws_veg_eb_av(j,i) = qsws_veg_eb_av(j,i) +                    &
     3131                                      surf_lsm_h%qsws_veg_eb(m)
    26313132             ENDDO
    26323133
    26333134          CASE ( 'r_a*' )
    2634              DO  i = nxlg, nxrg
    2635                 DO  j = nysg, nyng
    2636                    r_a_av(j,i) = r_a_av(j,i) + r_a(j,i)
    2637                 ENDDO
     3135             DO  m = 1, surf_lsm_h%ns
     3136                i   = surf_lsm_h%i(m)           
     3137                j   = surf_lsm_h%j(m)
     3138                r_a_av(j,i) = r_a_av(j,i) + surf_lsm_h%r_a(m)
    26383139             ENDDO
    26393140
    26403141          CASE ( 'r_s*' )
    2641              DO  i = nxlg, nxrg
    2642                 DO  j = nysg, nyng
    2643                    r_s_av(j,i) = r_s_av(j,i) + r_s(j,i)
    2644                 ENDDO
     3142             DO  m = 1, surf_lsm_h%ns
     3143                i   = surf_lsm_h%i(m)           
     3144                j   = surf_lsm_h%j(m)
     3145                r_s_av(j,i) = r_s_av(j,i) + surf_lsm_h%r_s(m)
    26453146             ENDDO
    26463147
    26473148          CASE ( 'shf_eb*' )
    2648              DO  i = nxlg, nxrg
    2649                 DO  j = nysg, nyng
    2650                    shf_eb_av(j,i) = shf_eb_av(j,i) + shf_eb(j,i)
    2651                 ENDDO
     3149             DO  m = 1, surf_lsm_h%ns
     3150                i   = surf_lsm_h%i(m)           
     3151                j   = surf_lsm_h%j(m)
     3152                shf_eb_av(j,i) = shf_eb_av(j,i) + surf_lsm_h%shf_eb(m)
    26523153             ENDDO
    26533154
    26543155          CASE ( 't_soil' )
    2655              DO  i = nxlg, nxrg
    2656                 DO  j = nysg, nyng
    2657                    DO  k = nzb_soil, nzt_soil
    2658                       t_soil_av(k,j,i) = t_soil_av(k,j,i) + t_soil(k,j,i)
    2659                    ENDDO
     3156             DO  m = 1, surf_lsm_h%ns
     3157                i   = surf_lsm_h%i(m)           
     3158                j   = surf_lsm_h%j(m)
     3159                DO  k = nzb_soil, nzt_soil
     3160                   t_soil_av(k,j,i) = t_soil_av(k,j,i) + t_soil_h%var_2d(k,m)
    26603161                ENDDO
    26613162             ENDDO
     
    26713172
    26723173          CASE ( 'c_liq*' )
    2673              DO  i = nxlg, nxrg
    2674                 DO  j = nysg, nyng
     3174             DO  i = nxl, nxr
     3175                DO  j = nys, nyn
    26753176                   c_liq_av(j,i) = c_liq_av(j,i) / REAL( average_count_3d, KIND=wp )
    26763177                ENDDO
     
    26783179
    26793180          CASE ( 'c_soil*' )
    2680              DO  i = nxlg, nxrg
    2681                 DO  j = nysg, nyng
     3181             DO  i = nxl, nxr
     3182                DO  j = nys, nyn
    26823183                   c_soil_av(j,i) = c_soil_av(j,i) / REAL( average_count_3d, KIND=wp )
    26833184                ENDDO
     
    26853186
    26863187          CASE ( 'c_veg*' )
    2687              DO  i = nxlg, nxrg
    2688                 DO  j = nysg, nyng
     3188             DO  i = nxl, nxr
     3189                DO  j = nys, nyn
    26893190                   c_veg_av(j,i) = c_veg_av(j,i) / REAL( average_count_3d, KIND=wp )
    26903191                ENDDO
     
    26923193
    26933194          CASE ( 'ghf_eb*' )
    2694              DO  i = nxlg, nxrg
    2695                 DO  j = nysg, nyng
     3195             DO  i = nxl, nxr
     3196                DO  j = nys, nyn
    26963197                   ghf_eb_av(j,i) = ghf_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
    26973198                ENDDO
     
    26993200
    27003201         CASE ( 'lai*' )
    2701              DO  i = nxlg, nxrg
    2702                 DO  j = nysg, nyng
     3202             DO  i = nxl, nxr
     3203                DO  j = nys, nyn
    27033204                   lai_av(j,i) = lai_av(j,i) / REAL( average_count_3d, KIND=wp )
    27043205                ENDDO
     
    27063207
    27073208          CASE ( 'm_liq_eb*' )
    2708              DO  i = nxlg, nxrg
    2709                 DO  j = nysg, nyng
     3209             DO  i = nxl, nxr
     3210                DO  j = nys, nyn
    27103211                   m_liq_eb_av(j,i) = m_liq_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
    27113212                ENDDO
     
    27133214
    27143215          CASE ( 'm_soil' )
    2715              DO  i = nxlg, nxrg
    2716                 DO  j = nysg, nyng
     3216             DO  i = nxl, nxr
     3217                DO  j = nys, nyn
    27173218                   DO  k = nzb_soil, nzt_soil
    27183219                      m_soil_av(k,j,i) = m_soil_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     
    27223223
    27233224          CASE ( 'qsws_eb*' )
    2724              DO  i = nxlg, nxrg
    2725                 DO  j = nysg, nyng
     3225             DO  i = nxl, nxr
     3226                DO  j = nys, nyn
    27263227                   qsws_eb_av(j,i) = qsws_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
    27273228                ENDDO
     
    27293230
    27303231          CASE ( 'qsws_liq_eb*' )
    2731              DO  i = nxlg, nxrg
    2732                 DO  j = nysg, nyng
     3232             DO  i = nxl, nxr
     3233                DO  j = nys, nyn
    27333234                   qsws_liq_eb_av(j,i) = qsws_liq_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
    27343235                ENDDO
     
    27363237
    27373238          CASE ( 'qsws_soil_eb*' )
    2738              DO  i = nxlg, nxrg
    2739                 DO  j = nysg, nyng
     3239             DO  i = nxl, nxr
     3240                DO  j = nys, nyn
    27403241                   qsws_soil_eb_av(j,i) = qsws_soil_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
    27413242                ENDDO
     
    27433244
    27443245          CASE ( 'qsws_veg_eb*' )
    2745              DO  i = nxlg, nxrg
    2746                 DO  j = nysg, nyng
     3246             DO  i = nxl, nxr
     3247                DO  j = nys, nyn
    27473248                   qsws_veg_eb_av(j,i) = qsws_veg_eb_av(j,i) / REAL( average_count_3d, KIND=wp )
    27483249                ENDDO
     
    27503251
    27513252          CASE ( 'r_a*' )
    2752              DO  i = nxlg, nxrg
    2753                 DO  j = nysg, nyng
     3253             DO  i = nxl, nxr
     3254                DO  j = nys, nyn
    27543255                   r_a_av(j,i) = r_a_av(j,i) / REAL( average_count_3d, KIND=wp )
    27553256                ENDDO
     
    27573258
    27583259          CASE ( 'r_s*' )
    2759              DO  i = nxlg, nxrg
    2760                 DO  j = nysg, nyng
     3260             DO  i = nxl, nxr
     3261                DO  j = nys, nyn
    27613262                   r_s_av(j,i) = r_s_av(j,i) / REAL( average_count_3d, KIND=wp )
    27623263                ENDDO
     
    27643265
    27653266          CASE ( 't_soil' )
    2766              DO  i = nxlg, nxrg
    2767                 DO  j = nysg, nyng
     3267             DO  i = nxl, nxr
     3268                DO  j = nys, nyn
    27683269                   DO  k = nzb_soil, nzt_soil
    27693270                      t_soil_av(k,j,i) = t_soil_av(k,j,i) / REAL( average_count_3d, KIND=wp )
     
    27713272                ENDDO
    27723273             ENDDO
     3274!
     3275!--
    27733276
    27743277       END SELECT
     
    28373340    CHARACTER (LEN=*) ::  variable !<
    28383341
    2839     INTEGER(iwp) ::  av !<
    2840     INTEGER(iwp) ::  i  !<
    2841     INTEGER(iwp) ::  j  !<
    2842     INTEGER(iwp) ::  k  !<
     3342    INTEGER(iwp) ::  av      !<
     3343    INTEGER(iwp) ::  i       !< running index
     3344    INTEGER(iwp) ::  j       !< running index
     3345    INTEGER(iwp) ::  k       !< running index
     3346    INTEGER(iwp) ::  m       !< running index
    28433347    INTEGER(iwp) ::  nzb_do  !<
    28443348    INTEGER(iwp) ::  nzt_do  !<
     
    28493353    REAL(wp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb:nzt+1) ::  local_pf !<
    28503354
     3355
    28513356    found = .TRUE.
    28523357
    28533358    SELECT CASE ( TRIM( variable ) )
    2854 
    2855 
     3359!
     3360!--    Before data is transfered to local_pf, transfer is it 2D dummy variable and exchange ghost points therein.
     3361!--    However, at this point this is only required for instantaneous arrays, time-averaged quantities are already exchanged.
    28563362       CASE ( 'c_liq*_xy' )        ! 2d-array
    28573363          IF ( av == 0 )  THEN
    2858              DO  i = nxlg, nxrg
    2859                 DO  j = nysg, nyng
    2860                    local_pf(i,j,nzb+1) = c_liq(j,i) * c_veg(j,i)
    2861                 ENDDO
     3364             DO  m = 1, surf_lsm_h%ns
     3365                i                   = surf_lsm_h%i(m)           
     3366                j                   = surf_lsm_h%j(m)
     3367                local_pf(i,j,nzb+1) = surf_lsm_h%c_liq(m) * surf_lsm_h%c_veg(m)
    28623368             ENDDO
    28633369          ELSE
     
    28743380       CASE ( 'c_soil*_xy' )        ! 2d-array
    28753381          IF ( av == 0 )  THEN
    2876              DO  i = nxlg, nxrg
    2877                 DO  j = nysg, nyng
    2878                    local_pf(i,j,nzb+1) = 1.0_wp - c_veg(j,i)
    2879                 ENDDO
     3382             DO  m = 1, surf_lsm_h%ns
     3383                i                   = surf_lsm_h%i(m)           
     3384                j                   = surf_lsm_h%j(m)
     3385                local_pf(i,j,nzb+1) = 1.0_wp - surf_lsm_h%c_veg(m)
    28803386             ENDDO
    28813387          ELSE
     
    28923398       CASE ( 'c_veg*_xy' )        ! 2d-array
    28933399          IF ( av == 0 )  THEN
    2894              DO  i = nxlg, nxrg
    2895                 DO  j = nysg, nyng
    2896                    local_pf(i,j,nzb+1) = c_veg(j,i)
    2897                 ENDDO
     3400             DO  m = 1, surf_lsm_h%ns
     3401                i                   = surf_lsm_h%i(m)           
     3402                j                   = surf_lsm_h%j(m)
     3403                local_pf(i,j,nzb+1) = surf_lsm_h%c_veg(m)
    28983404             ENDDO
    28993405          ELSE
     
    29103416       CASE ( 'ghf_eb*_xy' )        ! 2d-array
    29113417          IF ( av == 0 )  THEN
    2912              DO  i = nxlg, nxrg
    2913                 DO  j = nysg, nyng
    2914                    local_pf(i,j,nzb+1) = ghf_eb(j,i)
    2915                 ENDDO
     3418             DO  m = 1, surf_lsm_h%ns
     3419                i                   = surf_lsm_h%i(m)           
     3420                j                   = surf_lsm_h%j(m)
     3421                local_pf(i,j,nzb+1) = surf_lsm_h%ghf_eb(m)
    29163422             ENDDO
    29173423          ELSE
     
    29283434       CASE ( 'lai*_xy' )        ! 2d-array
    29293435          IF ( av == 0 )  THEN
    2930              DO  i = nxlg, nxrg
    2931                 DO  j = nysg, nyng
    2932                    local_pf(i,j,nzb+1) = lai(j,i)
    2933                 ENDDO
     3436             DO  m = 1, surf_lsm_h%ns
     3437                i                   = surf_lsm_h%i(m)           
     3438                j                   = surf_lsm_h%j(m)
     3439                local_pf(i,j,nzb+1) = surf_lsm_h%lai(m)
    29343440             ENDDO
    29353441          ELSE
     
    29463452       CASE ( 'm_liq_eb*_xy' )        ! 2d-array
    29473453          IF ( av == 0 )  THEN
    2948              DO  i = nxlg, nxrg
    2949                 DO  j = nysg, nyng
    2950                    local_pf(i,j,nzb+1) = m_liq_eb(j,i)
    2951                 ENDDO
     3454             DO  m = 1, surf_lsm_h%ns
     3455                i                   = surf_lsm_h%i(m)           
     3456                j                   = surf_lsm_h%j(m)
     3457                local_pf(i,j,nzb+1) = m_liq_eb_h%var_1d(m)
    29523458             ENDDO
    29533459          ELSE
     
    29643470       CASE ( 'm_soil_xy', 'm_soil_xz', 'm_soil_yz' )
    29653471          IF ( av == 0 )  THEN
    2966              DO  i = nxlg, nxrg
    2967                 DO  j = nysg, nyng
    2968                    DO k = nzb_soil, nzt_soil
    2969                       local_pf(i,j,k) = m_soil(k,j,i)
    2970                    ENDDO
     3472             DO  m = 1, surf_lsm_h%ns
     3473                i   = surf_lsm_h%i(m)           
     3474                j   = surf_lsm_h%j(m)
     3475                DO k = nzb_soil, nzt_soil
     3476                   local_pf(i,j,k) = m_soil_h%var_2d(k,m)
    29713477                ENDDO
    29723478             ENDDO
     
    29883494       CASE ( 'qsws_eb*_xy' )        ! 2d-array
    29893495          IF ( av == 0 ) THEN
    2990              DO  i = nxlg, nxrg
    2991                 DO  j = nysg, nyng
    2992                    local_pf(i,j,nzb+1) =  qsws_eb(j,i)
    2993                 ENDDO
     3496             DO  m = 1, surf_lsm_h%ns
     3497                i                   = surf_lsm_h%i(m)           
     3498                j                   = surf_lsm_h%j(m)
     3499                local_pf(i,j,nzb+1) =  surf_lsm_h%qsws_eb(m)
    29943500             ENDDO
    29953501          ELSE
     
    30063512       CASE ( 'qsws_liq_eb*_xy' )        ! 2d-array
    30073513          IF ( av == 0 ) THEN
    3008              DO  i = nxlg, nxrg
    3009                 DO  j = nysg, nyng
    3010                    local_pf(i,j,nzb+1) =  qsws_liq_eb(j,i)
    3011                 ENDDO
     3514             DO  m = 1, surf_lsm_h%ns
     3515                i                   = surf_lsm_h%i(m)           
     3516                j                   = surf_lsm_h%j(m)
     3517                local_pf(i,j,nzb+1) = surf_lsm_h%qsws_liq_eb(m)
    30123518             ENDDO
    30133519          ELSE
     
    30243530       CASE ( 'qsws_soil_eb*_xy' )        ! 2d-array
    30253531          IF ( av == 0 ) THEN
    3026              DO  i = nxlg, nxrg
    3027                 DO  j = nysg, nyng
    3028                    local_pf(i,j,nzb+1) =  qsws_soil_eb(j,i)
    3029                 ENDDO
     3532             DO  m = 1, surf_lsm_h%ns
     3533                i                   = surf_lsm_h%i(m)           
     3534                j                   = surf_lsm_h%j(m)
     3535                local_pf(i,j,nzb+1) =  surf_lsm_h%qsws_soil_eb(m)
    30303536             ENDDO
    30313537          ELSE
     
    30423548       CASE ( 'qsws_veg_eb*_xy' )        ! 2d-array
    30433549          IF ( av == 0 ) THEN
    3044              DO  i = nxlg, nxrg
    3045                 DO  j = nysg, nyng
    3046                    local_pf(i,j,nzb+1) =  qsws_veg_eb(j,i)
    3047                 ENDDO
     3550             DO  m = 1, surf_lsm_h%ns
     3551                i                   = surf_lsm_h%i(m)           
     3552                j                   = surf_lsm_h%j(m)
     3553                local_pf(i,j,nzb+1) =  surf_lsm_h%qsws_veg_eb(m)
    30483554             ENDDO
    30493555          ELSE
     
    30613567       CASE ( 'r_a*_xy' )        ! 2d-array
    30623568          IF ( av == 0 )  THEN
    3063              DO  i = nxlg, nxrg
    3064                 DO  j = nysg, nyng
    3065                    local_pf(i,j,nzb+1) = r_a(j,i)
    3066                 ENDDO
     3569             DO  m = 1, surf_lsm_h%ns
     3570                i                   = surf_lsm_h%i(m)           
     3571                j                   = surf_lsm_h%j(m)
     3572                local_pf(i,j,nzb+1) = surf_lsm_h%r_a(m)
    30673573             ENDDO
    30683574          ELSE
     
    30793585       CASE ( 'r_s*_xy' )        ! 2d-array
    30803586          IF ( av == 0 )  THEN
    3081              DO  i = nxlg, nxrg
    3082                 DO  j = nysg, nyng
    3083                    local_pf(i,j,nzb+1) = r_s(j,i)
    3084                 ENDDO
     3587             DO  m = 1, surf_lsm_h%ns
     3588                i                   = surf_lsm_h%i(m)           
     3589                j                   = surf_lsm_h%j(m)
     3590                local_pf(i,j,nzb+1) = surf_lsm_h%r_s(m)
    30853591             ENDDO
    30863592          ELSE
     
    30973603       CASE ( 'shf_eb*_xy' )        ! 2d-array
    30983604          IF ( av == 0 ) THEN
    3099              DO  i = nxlg, nxrg
    3100                 DO  j = nysg, nyng
    3101                    local_pf(i,j,nzb+1) =  shf_eb(j,i)
    3102                 ENDDO
     3605             DO  m = 1, surf_lsm_h%ns
     3606                i                   = surf_lsm_h%i(m)           
     3607                j                   = surf_lsm_h%j(m)
     3608                local_pf(i,j,nzb+1) =  surf_lsm_h%shf_eb(m)
    31033609             ENDDO
    31043610          ELSE
     
    31153621       CASE ( 't_soil_xy', 't_soil_xz', 't_soil_yz' )
    31163622          IF ( av == 0 )  THEN
    3117              DO  i = nxlg, nxrg
    3118                 DO  j = nysg, nyng
    3119                    DO k = nzb_soil, nzt_soil
    3120                       local_pf(i,j,k) = t_soil(k,j,i)
    3121                    ENDDO
     3623             DO  m = 1, surf_lsm_h%ns
     3624                i   = surf_lsm_h%i(m)           
     3625                j   = surf_lsm_h%j(m)
     3626                DO k = nzb_soil, nzt_soil
     3627                   local_pf(i,j,k) = t_soil_h%var_2d(k,m)
    31223628                ENDDO
    31233629             ENDDO
     
    31683674    INTEGER(iwp) ::  j     !<
    31693675    INTEGER(iwp) ::  k     !<
     3676    INTEGER(iwp) ::  m     !< running index
    31703677
    31713678    LOGICAL      ::  found !<
     
    31783685
    31793686    SELECT CASE ( TRIM( variable ) )
    3180 
     3687!
     3688!--   Requires 3D exchange
    31813689
    31823690      CASE ( 'm_soil' )
    31833691
    31843692         IF ( av == 0 )  THEN
    3185             DO  i = nxlg, nxrg
    3186                DO  j = nysg, nyng
    3187                   DO  k = nzb_soil, nzt_soil
    3188                      local_pf(i,j,k) = m_soil(k,j,i)
    3189                   ENDDO
    3190                ENDDO
     3693            DO  m = 1, surf_lsm_h%ns
     3694                i   = surf_lsm_h%i(m)           
     3695                j   = surf_lsm_h%j(m)
     3696                DO  k = nzb_soil, nzt_soil
     3697                   local_pf(i,j,k) = m_soil_h%var_2d(k,m)
     3698                ENDDO
    31913699            ENDDO
    31923700         ELSE
     
    32033711
    32043712         IF ( av == 0 )  THEN
    3205             DO  i = nxlg, nxrg
    3206                DO  j = nysg, nyng
    3207                   DO  k = nzb_soil, nzt_soil
    3208                      local_pf(i,j,k) = t_soil(k,j,i)
    3209                   ENDDO
     3713            DO  m = 1, surf_lsm_h%ns
     3714               i   = surf_lsm_h%i(m)           
     3715               j   = surf_lsm_h%j(m)
     3716               DO  k = nzb_soil, nzt_soil
     3717                  local_pf(i,j,k) = t_soil_h%var_2d(k,m)
    32103718               ENDDO
    32113719            ENDDO
     
    32613769          WRITE ( 14 )  'lai_av              ';  WRITE ( 14 )  lai_av
    32623770       ENDIF
    3263        WRITE ( 14 )  'm_liq_eb            ';  WRITE ( 14 )  m_liq_eb
     3771       WRITE ( 14 )     'm_liq_eb            ';  WRITE ( 14 )  m_liq_eb_h
    32643772       IF ( ALLOCATED( m_liq_eb_av ) )  THEN
    32653773          WRITE ( 14 )  'm_liq_eb_av         ';  WRITE ( 14 )  m_liq_eb_av
    32663774       ENDIF
    3267        WRITE ( 14 )  'm_soil              ';  WRITE ( 14 )  m_soil
     3775       WRITE ( 14 )     'm_soil              ';  WRITE ( 14 )  m_soil_h
    32683776       IF ( ALLOCATED( m_soil_av ) )  THEN
    32693777          WRITE ( 14 )  'm_soil_av           ';  WRITE ( 14 )  m_soil_av
     
    32843792          WRITE ( 14 )  'shf_eb_av           ';  WRITE ( 14 )  shf_eb_av
    32853793       ENDIF
    3286        WRITE ( 14 )  't_soil              ';  WRITE ( 14 )  t_soil
     3794       WRITE ( 14 )     't_soil              ';  WRITE ( 14 )  t_soil_h
    32873795       IF ( ALLOCATED( t_soil_av ) )  THEN
    32883796          WRITE ( 14 )  't_soil_av           ';  WRITE ( 14 )  t_soil_av
     
    33493857          tmp_3d2   !<
    33503858
     3859    REAL(wp),                                                                  &
     3860       DIMENSION(1:surf_lsm_h%ns) ::                                  &
     3861          tmp_walltype_1d   !<
     3862
     3863    REAL(wp),                                                                  &
     3864       DIMENSION(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) ::              &
     3865          tmp_walltype_2d   !<
     3866
     3867    REAL(wp),                                                                  &
     3868       DIMENSION(nzb_soil:nzt_soil,1:surf_lsm_h%ns) ::                &
     3869          tmp_walltype_2d2  !<
     3870
    33513871
    33523872   IF ( initializing_actions == 'read_restart_data' )  THEN
     
    33663886            nync = nynfa(i,k) + offset_ya(i,k)
    33673887
    3368 
    33693888            SELECT CASE ( TRIM( field_char ) )
     3889
    33703890
    33713891                CASE ( 'c_liq_av' )
     
    33743894                   ENDIF
    33753895                   IF ( k == 1 )  READ ( 13 )  tmp_2d
    3376                    c_liq_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3896                   c_liq_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =         &
    33773897                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    33783898
     
    33823902                   ENDIF
    33833903                   IF ( k == 1 )  READ ( 13 )  tmp_2d
    3384                    c_soil_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3904                   c_soil_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =        &
    33853905                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    33863906
     
    33903910                   ENDIF
    33913911                   IF ( k == 1 )  READ ( 13 )  tmp_2d
    3392                    c_veg_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3912                   c_veg_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =         &
    33933913                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    33943914
     
    33983918                   ENDIF
    33993919                   IF ( k == 1 )  READ ( 13 )  tmp_2d
    3400                    ghf_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3920                   ghf_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =        &
    34013921                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    34023922
    34033923                CASE ( 'm_liq_eb' )
    3404                    IF ( k == 1 )  READ ( 13 )  tmp_2d
    3405                    m_liq_eb(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =        &
    3406                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     3924                   IF ( k == 1 )  READ ( 13 )  tmp_walltype_1d !tmp_2d
     3925                   m_liq_eb_h%var_1d(1:surf_lsm_h%ns)  =                       &
     3926                                 tmp_walltype_1d(1:surf_lsm_h%ns)
    34073927
    34083928                CASE ( 'lai_av' )
     
    34113931                   ENDIF
    34123932                   IF ( k == 1 )  READ ( 13 )  tmp_2d
    3413                    lai_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     3933                   lai_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =           &
    34143934                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    34153935
     
    34233943
    34243944                CASE ( 'm_soil' )
    3425                    IF ( k == 1 )  READ ( 13 )  tmp_3d2(:,:,:)
    3426                    m_soil(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =         &
    3427                           tmp_3d2(nzb_soil:nzt_soil,nysf-nbgp:nynf             &
    3428                           +nbgp,nxlf-nbgp:nxrf+nbgp)
     3945                   IF ( k == 1 )  READ ( 13 )  tmp_walltype_2d2(:,:)
     3946                   m_soil_h%var_2d(:,1:surf_lsm_h%ns) =                        &
     3947                          tmp_walltype_2d2(:,1:surf_lsm_h%ns)   
    34293948
    34303949                CASE ( 'm_soil_av' )
     
    34423961                   ENDIF 
    34433962                   IF ( k == 1 )  READ ( 13 )  tmp_2d
    3444                    qsws_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
     3963                   qsws_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =      &
    34453964                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    34463965
     
    34503969                   ENDIF 
    34513970                   IF ( k == 1 )  READ ( 13 )  tmp_2d
    3452                    qsws_liq_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
     3971                   qsws_liq_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =  &
    34533972                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    34543973                CASE ( 'qsws_soil_eb_av' )
     
    34733992                   ENDIF
    34743993                   IF ( k == 1 )  READ ( 13 )  tmp_2d
    3475                    shf_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  = &
     3994                   shf_eb_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =       &
    34763995                         tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    34773996
    34783997                CASE ( 't_soil' )
    3479                    IF ( k == 1 )  READ ( 13 )  tmp_3d
    3480                    t_soil(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =         &
    3481                                    tmp_3d(:,nysf-nbgp:nynf+nbgp,               &
    3482                                                 nxlf-nbgp:nxrf+nbgp)
     3998                   IF ( k == 1 )  READ ( 13 )  tmp_walltype_2d(:,:)
     3999                   t_soil_h%var_2d(:,1:surf_lsm_h%ns) =                        &
     4000                         tmp_walltype_2d(:,1:surf_lsm_h%ns)     
    34834001
    34844002                CASE ( 't_soil_av' )
     
    34884006                   IF ( k == 1 )  READ ( 13 )  tmp_3d2(:,:,:)
    34894007                   t_soil_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =      &
    3490                                     tmp_3d(:,nysf-nbgp:nynf+nbgp,             &
     4008                                    tmp_3d2(:,nysf-nbgp:nynf+nbgp,             &
    34914009                                    nxlf-nbgp:nxrf+nbgp)
    34924010
     
    35254043       IMPLICIT NONE
    35264044
    3527        INTEGER :: i  !< running index
    3528        INTEGER :: j  !< running index
     4045       INTEGER(iwp) ::  i       !< running index
     4046       INTEGER(iwp) ::  j       !< running index
     4047       INTEGER(iwp) ::  m       !< running index
    35294048
    35304049       REAL(wp), PARAMETER :: alpha_ch  = 0.018_wp !< Charnock constant (0.01-0.11). Use 0.01 for FLake and 0.018 for ECMWF
     
    35334052!       REAL(wp) :: re_0 !< near-surface roughness Reynolds number
    35344053
    3535 
    3536        DO  i = nxlg, nxrg   
    3537           DO  j = nysg, nyng
    3538              IF ( water_surface(j,i) )  THEN
    3539 
    3540 !
    3541 !--             Disabled: FLake parameterization. Ideally, the Charnock
    3542 !--             coefficient should depend on the water depth and the fetch
    3543 !--             length
    3544 !                re_0 = z0(j,i) * us(j,i) / molecular_viscosity
     4054       DO  m = 1, surf_lsm_h%ns
     4055
     4056          i   = surf_lsm_h%i(m)           
     4057          j   = surf_lsm_h%j(m)
     4058         
     4059          IF ( surf_lsm_h%water_surface(m) )  THEN
     4060
     4061!
     4062!--          Disabled: FLake parameterization. Ideally, the Charnock
     4063!--          coefficient should depend on the water depth and the fetch
     4064!--          length
     4065!             re_0 = z0(j,i) * us(j,i) / molecular_viscosity
    35454066!       
    3546 !                z0(j,i) = MAX( 0.1_wp * molecular_viscosity / us(j,i),            &
    3547 !                              alpha_ch * us(j,i) / g )
    3548 !
    3549 !                z0h(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 3.2_wp ) )
    3550 !                z0q(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 4.2_wp ) )
    3551 
    3552 !
    3553 !--              Set minimum roughness length for u* > 0.2
    3554 !                IF ( us(j,i) > 0.2_wp )  THEN
    3555 !                   z0h(j,i) = MAX( 1.0E-5_wp, z0h(j,i) )
    3556 !                   z0q(j,i) = MAX( 1.0E-5_wp, z0q(j,i) )
    3557 !                ENDIF
    3558 
    3559 !
    3560 !--             ECMWF IFS model parameterization after Beljaars (1994). At low
    3561 !--             wind speed, the sea surface becomes aerodynamically smooth and
    3562 !--             the roughness scales with the viscosity. At high wind speed, the
    3563 !--             Charnock relation is used.
    3564                 z0(j,i) =   ( 0.11_wp * molecular_viscosity / us(j,i) )        &
    3565                           + ( alpha_ch * us(j,i)**2 / g )
    3566 
    3567                 z0h(j,i) = 0.40_wp * molecular_viscosity / us(j,i)
    3568                 z0q(j,i) = 0.62_wp * molecular_viscosity / us(j,i)
    3569 
    3570              ENDIF
    3571           ENDDO
     4067!             z0(j,i) = MAX( 0.1_wp * molecular_viscosity / us(j,i),            &
     4068!                           alpha_ch * us(j,i) / g )
     4069!
     4070!             z0h(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 3.2_wp ) )
     4071!             z0q(j,i) = z0(j,i) * EXP( - kappa / pr_number * ( 4.0_wp * SQRT( re_0 ) - 4.2_wp ) )
     4072
     4073!
     4074!--           Set minimum roughness length for u* > 0.2
     4075!             IF ( us(j,i) > 0.2_wp )  THEN
     4076!                z0h(j,i) = MAX( 1.0E-5_wp, z0h(j,i) )
     4077!                z0q(j,i) = MAX( 1.0E-5_wp, z0q(j,i) )
     4078!             ENDIF
     4079
     4080!
     4081!--          ECMWF IFS model parameterization after Beljaars (1994). At low
     4082!--          wind speed, the sea surface becomes aerodynamically smooth and
     4083!--          the roughness scales with the viscosity. At high wind speed, the
     4084!--          Charnock relation is used.
     4085             surf_lsm_h%z0(m)  = ( 0.11_wp * molecular_viscosity /             &
     4086                                 surf_lsm_h%us(m) )  &
     4087                               + ( alpha_ch * surf_lsm_h%us(m)**2 / g )
     4088
     4089             surf_lsm_h%z0h(m) = 0.40_wp * molecular_viscosity /               &
     4090                                 surf_lsm_h%us(m)
     4091             surf_lsm_h%z0q(m) = 0.62_wp * molecular_viscosity /               &
     4092                                 surf_lsm_h%us(m)
     4093
     4094          ENDIF
    35724095       ENDDO
    35734096
     
    35754098
    35764099
    3577 !------------------------------------------------------------------------------!
    3578 ! Description:
    3579 ! ------------
    3580 !> Calculation of specific humidity of the skin layer (surface). It is assumend
    3581 !> that the skin is always saturated.
    3582 !------------------------------------------------------------------------------!
    3583     SUBROUTINE calc_q_surface
    3584 
    3585        IMPLICIT NONE
    3586 
    3587        INTEGER :: i              !< running index
    3588        INTEGER :: j              !< running index
    3589        INTEGER :: k              !< running index
    3590 
    3591        REAL(wp) :: resistance    !< aerodynamic and soil resistance term
    3592 
    3593        DO  i = nxlg, nxrg   
    3594           DO  j = nysg, nyng
    3595              k = nzb_s_inner(j,i)
    3596 
    3597 !
    3598 !--          Calculate water vapour pressure at saturation
    3599              e_s = 0.01_wp * 610.78_wp * EXP( 17.269_wp * ( t_surface_p(j,i)   &
    3600                    - 273.16_wp ) / ( t_surface_p(j,i) - 35.86_wp ) )
    3601 
    3602 !
    3603 !--          Calculate specific humidity at saturation
    3604              q_s = 0.622_wp * e_s / surface_pressure
    3605 
    3606              resistance = r_a(j,i) / (r_a(j,i) + r_s(j,i))
    3607 
    3608 !
    3609 !--          Calculate specific humidity at surface
    3610              IF ( cloud_physics )  THEN
    3611                 q(k,j,i) = resistance * q_s + (1.0_wp - resistance)            &
    3612                              * ( q(k+1,j,i) - ql(k+1,j,i) )
    3613              ELSE
    3614                 q(k,j,i) = resistance * q_s + (1.0_wp - resistance)            &
    3615                              * q(k+1,j,i)
    3616              ENDIF
    3617 
    3618 !
    3619 !--          Update virtual potential temperature
    3620              vpt(k,j,i) = pt(k,j,i) * ( 1.0_wp + 0.61_wp * q(k,j,i) )
    3621 
    3622           ENDDO
    3623        ENDDO
    3624 
    3625     END SUBROUTINE calc_q_surface
    3626 
    36274100
    36284101 END MODULE land_surface_model_mod
Note: See TracChangeset for help on using the changeset viewer.