Ignore:
Timestamp:
Dec 14, 2017 5:12:51 PM (4 years ago)
Author:
kanani
Message:

Merge of branch palm4u into trunk

Location:
palm/trunk
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk

  • palm/trunk/SOURCE

  • palm/trunk/SOURCE/land_surface_model_mod.f90

    r2608 r2696  
    11!> @file land_surface_model_mod.f90
    22!------------------------------------------------------------------------------!
    3 ! This file is part of PALM.
     3! This file is part of the PALM model system.
    44!
    55! PALM is free software: you can redistribute it and/or modify it under the
     
    2525! -----------------
    2626! $Id$
     27! Bugfix: missing USE statement for calc_mean_profile
     28! do not write surface temperatures onto pt array as this might cause
     29! problems with nesting (MS)
     30! Revised calculation of pt1 and qv1 (now done in surface_layer_fluxes). Bugfix
     31! in calculation of surface density (cannot be done via an surface non-air
     32! temperature) (BM)
     33! Bugfix: g_d was NaN for non-vegetaed surface types (BM)
     34! Bugfix initialization of c_veg and lai
     35! Revise data output to enable _FillValues
     36! Bugfix in calcultion of r_a and rad_net_l (MS)
     37! Bugfix: rad_net is not updated in case of radiation_interaction and must thu
     38! be calculated again from the radiative fluxes
     39! Temporary fix for cases where no soil model is used on some PEs (BM)
     40! Revised input and initialization of soil and surface paramters
     41! pavement_depth is variable for each surface element
     42! radiation quantities belong to surface type now
     43! surface fractions initialized
     44! Rename lsm_last_actions into lsm_write_restart_data (MS)
     45!
     46! 2608 2017-11-13 14:04:26Z schwenkel
    2747! Calculation of magnus equation in external module (diagnostic_quantities_mod).
    2848! Adjust calculation of vapor pressure and saturation mixing ratio that it is
     
    265285!> DALES and UCLA-LES models.
    266286!>
    267 !> @todo Restart data for vertical natural land-surfaces
    268287!> @todo Extensive verification energy-balance solver for vertical surfaces,
    269288!>       e.g. parametrization of r_a
     
    282301!> @note No time step criterion is required as long as the soil layers do not
    283302!>       become too thin.
     303!> @todo Attention, pavement_subpars_1/2 are hardcoded to 8 levels, in case
     304!>       more levels are used this may cause an potential bug
     305!> @todo Routine calc_q_surface required?
    284306!------------------------------------------------------------------------------!
    285307 MODULE land_surface_model_mod
     
    287309    USE arrays_3d,                                                             &
    288310        ONLY:  hyp, pt, prr, q, q_p, ql, vpt, u, v, w
     311
     312    USE calc_mean_profile_mod,                                                 &
     313        ONLY:  calc_mean_profile
    289314
    290315    USE cloud_parameters,                                                      &
     
    301326
    302327    USE indices,                                                               &
    303         ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb
     328        ONLY:  nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb
     329
     330    USE netcdf_data_input_mod,                                                 &
     331        ONLY :  building_type_f, init_3d, input_pids_static,                   &
     332                netcdf_data_input_interpolate,                                 &
     333                pavement_pars_f, pavement_subsurface_pars_f, pavement_type_f,  &
     334                root_area_density_lsm_f, soil_pars_f, soil_type_f,             &
     335                surface_fraction_f, vegetation_pars_f, vegetation_type_f,      &
     336                water_pars_f, water_type_f
    304337
    305338    USE kinds
     
    308341
    309342    USE radiation_model_mod,                                                   &
    310         ONLY:  albedo, albedo_type, emissivity, force_radiation_call, rad_net, &
    311                rad_sw_in, rad_lw_out, rad_lw_out_change_0, radiation_scheme,   &
    312                unscheduled_radiation_calls
     343        ONLY:  albedo, albedo_type, emissivity, force_radiation_call,          &
     344               radiation_scheme, unscheduled_radiation_calls
    313345       
    314346    USE statistics,                                                            &
     
    316348
    317349    USE surface_mod,                                                           &
    318         ONLY :  surf_lsm_h, surf_lsm_v, surf_type
     350        ONLY :  surf_lsm_h, surf_lsm_v, surf_type, surface_restore_elements
    319351
    320352    IMPLICIT NONE
    321353
    322354    TYPE surf_type_lsm
    323        REAL(wp), DIMENSION(:),   ALLOCATABLE ::  var_1D !< 1D prognostic variable
    324        REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  var_2D !< 2D prognostic variable
     355       REAL(wp), DIMENSION(:),   ALLOCATABLE ::  var_1d !< 1D prognostic variable
     356       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  var_2d !< 2D prognostic variable
    325357    END TYPE surf_type_lsm
    326358
     
    355387!
    356388!-- LSM variables
    357     CHARACTER(10) :: surface_type = 'vegetation'  !< general classification. Allowed are:
     389    CHARACTER(10) :: surface_type = 'netcdf'      !< general classification. Allowed are:
    358390                                                  !< 'vegetation', 'pavement', ('building'),
    359391                                                  !< 'water', and 'netcdf'
     
    447479                                     m_soil_v,       & !< Soil moisture (m3/m3), vertical surface elements
    448480                                     m_soil_v_p        !< Prog. soil moisture (m3/m3), vertical surface elements
     481
    449482#else
    450483    TYPE(surf_type_lsm), POINTER ::  t_soil_h,    & !< Soil temperature (K), horizontal surface elements
     
    474507    TYPE(surf_type_lsm), TARGET   ::  t_surface_h,    & !< surface temperature (K), horizontal surface elements
    475508                                      t_surface_h_p,  & !< progn. surface temperature (K), horizontal surface elements
    476                                       m_liq_h,     & !< liquid water reservoir (m), horizontal surface elements
    477                                       m_liq_h_p      !< progn. liquid water reservoir (m), horizontal surface elements
     509                                      m_liq_h,        & !< liquid water reservoir (m), horizontal surface elements
     510                                      m_liq_h_p         !< progn. liquid water reservoir (m), horizontal surface elements
    478511
    479512    TYPE(surf_type_lsm), DIMENSION(0:3), TARGET   ::  &
     
    533566!-- Energy balance variables               
    534567    REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
    535               c_liq_av,     & !< average of c_liq
    536               c_soil_av,    & !< average of c_soil
    537               c_veg_av,     & !< average of c_veg
    538               ghf_av,       & !< average of ghf
    539               lai_av,       & !< average of lai
    540               qsws_liq_av,  & !< average of qsws_liq
    541               qsws_soil_av, & !< average of qsws_soil
    542               qsws_veg_av,  & !< average of qsws_veg
    543               r_a_av,       & !< average of r_a
    544               r_s_av          !< average of r_s
     568              c_liq_av,         & !< average of c_liq
     569              c_soil_av,        & !< average of c_soil
     570              c_veg_av,         & !< average of c_veg
     571              ghf_av,           & !< average of ghf
     572              lai_av,           & !< average of lai
     573              qsws_liq_av,      & !< average of qsws_liq
     574              qsws_soil_av,     & !< average of qsws_soil
     575              qsws_veg_av,      & !< average of qsws_veg
     576              r_a_av,           & !< average of r_a
     577              r_s_av              !< average of r_s
    545578                   
    546579
     
    548581!-- Predefined Land surface classes (vegetation_type)
    549582    CHARACTER(26), DIMENSION(0:18), PARAMETER :: vegetation_type_name = (/ &
    550                                    'user defined              ', & !  0
    551                                    'bare soil                 ', & !  1                           
    552                                    'crops, mixed farming      ', & !  2
    553                                    'short grass               ', & !  3
    554                                    'evergreen needleleaf trees', & !  4
    555                                    'deciduous needleleaf trees', & !  5
    556                                    'evergreen broadleaf trees ', & !  6
    557                                    'deciduous broadleaf trees ', & !  7
    558                                    'tall grass                ', & !  8
    559                                    'desert                    ', & !  9
    560                                    'tundra                    ', & ! 10
    561                                    'irrigated crops           ', & ! 11
    562                                    'semidesert                ', & ! 12
    563                                    'ice caps and glaciers     ', & ! 13
    564                                    'bogs and marshes          ', & ! 14
    565                                    'evergreen shrubs          ', & ! 15
    566                                    'deciduous shrubs          ', & ! 16
    567                                    'mixed forest/woodland     ', & ! 17
    568                                    'interrupted forest        '  & ! 18
     583                                   'user defined              ',           & !  0
     584                                   'bare soil                 ',           & !  1                           
     585                                   'crops, mixed farming      ',           & !  2
     586                                   'short grass               ',           & !  3
     587                                   'evergreen needleleaf trees',           & !  4
     588                                   'deciduous needleleaf trees',           & !  5
     589                                   'evergreen broadleaf trees ',           & !  6
     590                                   'deciduous broadleaf trees ',           & !  7
     591                                   'tall grass                ',           & !  8
     592                                   'desert                    ',           & !  9
     593                                   'tundra                    ',           & ! 10
     594                                   'irrigated crops           ',           & ! 11
     595                                   'semidesert                ',           & ! 12
     596                                   'ice caps and glaciers     ',           & ! 13
     597                                   'bogs and marshes          ',           & ! 14
     598                                   'evergreen shrubs          ',           & ! 15
     599                                   'deciduous shrubs          ',           & ! 16
     600                                   'mixed forest/woodland     ',           & ! 17
     601                                   'interrupted forest        '            & ! 18
    569602                                                                 /)
    570603
     
    572605!-- Soil model classes (soil_type)
    573606    CHARACTER(12), DIMENSION(0:6), PARAMETER :: soil_type_name = (/ &
    574                                    'user defined', & ! 0
    575                                    'coarse      ', & ! 1
    576                                    'medium      ', & ! 2
    577                                    'medium-fine ', & ! 3
    578                                    'fine        ', & ! 4
    579                                    'very fine   ', & ! 5
    580                                    'organic     '  & ! 6
     607                                   'user defined',                  & ! 0
     608                                   'coarse      ',                  & ! 1
     609                                   'medium      ',                  & ! 2
     610                                   'medium-fine ',                  & ! 3
     611                                   'fine        ',                  & ! 4
     612                                   'very fine   ',                  & ! 5
     613                                   'organic     '                   & ! 6
    581614                                                                 /)
    582615
     
    600633                                   'artifical turf (sports)      ', & ! 14
    601634                                   'clay (sports)                '  & ! 15
    602                                                                  /)                                                                 
    603    
    604 
    605 
    606 
    607                                                              
     635                                                                 /)                                                             
     636                                                                 
    608637!
    609638!-- Water classes
    610639    CHARACTER(12), DIMENSION(0:5), PARAMETER :: water_type_name = (/ &
    611                                    'user defined', & ! 0
    612                                    'lake        ', & ! 1
    613                                    'river       ', & ! 2
    614                                    'ocean       ', & ! 3
    615                                    'pond        ', & ! 4
    616                                    'fountain    '  & ! 5
     640                                   'user defined',                   & ! 0
     641                                   'lake        ',                   & ! 1
     642                                   'river       ',                   & ! 2
     643                                   'ocean       ',                   & ! 3
     644                                   'pond        ',                   & ! 4
     645                                   'fountain    '                    & ! 5
    617646                                                                  /)                                                                 
    618                                                                  
    619                                                                  !
     647!
     648!-- IDs for vegetation, pavement and water surfaces
     649    INTEGER(iwp) ::  ind_veg = 0    !< index for vegetation surfaces, used to assess surface-fraction, albedo, etc.     
     650    INTEGER(iwp) ::  ind_pav = 1    !< index for pavement surfaces, used to assess surface-fraction, albedo, etc.       
     651    INTEGER(iwp) ::  ind_wat = 2    !< index for vegetation surfaces, used to assess surface-fraction, albedo, etc.                         
     652                   
     653!
    620654!-- Land surface parameters according to the respective classes (vegetation_type)
    621 
     655    INTEGER(iwp) ::  ind_v_rc_min = 0    !< index for r_canopy_min in vegetation_pars
     656    INTEGER(iwp) ::  ind_v_rc_lai = 1    !< index for LAI in vegetation_pars
     657    INTEGER(iwp) ::  ind_v_c_veg   = 2   !< index for c_veg in vegetation_pars
     658    INTEGER(iwp) ::  ind_v_gd  = 3       !< index for g_d in vegetation_pars
     659    INTEGER(iwp) ::  ind_v_z0 = 4        !< index for z0 in vegetation_pars
     660    INTEGER(iwp) ::  ind_v_z0qh = 5      !< index for z0h / z0q in vegetation_pars
     661    INTEGER(iwp) ::  ind_v_lambda_s = 6  !< index for lambda_s_s in vegetation_pars
     662    INTEGER(iwp) ::  ind_v_lambda_u = 7  !< index for lambda_s_u in vegetation_pars
     663    INTEGER(iwp) ::  ind_v_f_sw_in = 8   !< index for f_sw_in in vegetation_pars
     664    INTEGER(iwp) ::  ind_v_c_surf = 9    !< index for c_surface in vegetation_pars
     665    INTEGER(iwp) ::  ind_v_at = 10       !< index for albedo_type in vegetation_pars
     666    INTEGER(iwp) ::  ind_v_emis = 11     !< index for emissivity in vegetation_pars
     667
     668    INTEGER(iwp) ::  ind_w_temp     = 0    !< index for temperature in water_pars
     669    INTEGER(iwp) ::  ind_w_z0       = 1    !< index for z0 in water_pars
     670    INTEGER(iwp) ::  ind_w_z0h      = 2    !< index for z0h in water_pars
     671    INTEGER(iwp) ::  ind_w_lambda_s = 3    !< index for lambda_s_s in water_pars
     672    INTEGER(iwp) ::  ind_w_lambda_u = 4    !< index for lambda_s_u in water_pars
     673    INTEGER(iwp) ::  ind_w_at       = 5    !< index for albedo type in water_pars
     674    INTEGER(iwp) ::  ind_w_emis     = 6    !< index for emissivity in water_pars
     675
     676    INTEGER(iwp) ::  ind_p_z0       = 0    !< index for z0 in pavement_pars
     677    INTEGER(iwp) ::  ind_p_z0h      = 1    !< index for z0h in pavement_pars
     678    INTEGER(iwp) ::  ind_p_at       = 2    !< index for albedo type in pavement_pars
     679    INTEGER(iwp) ::  ind_p_emis     = 3    !< index for emissivity in pavement_pars
     680    INTEGER(iwp) ::  ind_p_lambda_h = 0    !< index for lambda_h in pavement_subsurface_pars
     681    INTEGER(iwp) ::  ind_p_rho_c    = 1    !< index for rho_c in pavement_pars
    622682!
    623683!-- Land surface parameters
     
    649709!--                                level 1 - level 4 according to zs_ref
    650710    REAL(wp), DIMENSION(0:3,1:18), PARAMETER :: root_distribution = RESHAPE( (/ &
    651                                  1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & !  1
    652                                  0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & !  2
    653                                  0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp, & !  3
    654                                  0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp, & !  4
    655                                  0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp, & !  5
    656                                  0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp, & !  6
    657                                  0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp, & !  7
    658                                  0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp, & !  8
    659                                  1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & !  9
    660                                  0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp, & ! 10
    661                                  0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp, & ! 11
    662                                  0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp, & ! 12
    663                                  0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp, & ! 13
    664                                  0.25_wp, 0.34_wp, 0.27_wp, 0.11_wp, & ! 14
    665                                  0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 15
    666                                  0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp, & ! 16
    667                                  0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp, & ! 17
    668                                  0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp  & ! 18
     711                                 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & !  1
     712                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & !  2
     713                                 0.35_wp, 0.38_wp, 0.23_wp, 0.04_wp,            & !  3
     714                                 0.26_wp, 0.39_wp, 0.29_wp, 0.06_wp,            & !  4
     715                                 0.26_wp, 0.38_wp, 0.29_wp, 0.07_wp,            & !  5
     716                                 0.24_wp, 0.38_wp, 0.31_wp, 0.07_wp,            & !  6
     717                                 0.25_wp, 0.34_wp, 0.27_wp, 0.14_wp,            & !  7
     718                                 0.27_wp, 0.27_wp, 0.27_wp, 0.09_wp,            & !  8
     719                                 1.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & !  9
     720                                 0.47_wp, 0.45_wp, 0.08_wp, 0.00_wp,            & ! 10
     721                                 0.24_wp, 0.41_wp, 0.31_wp, 0.04_wp,            & ! 11
     722                                 0.17_wp, 0.31_wp, 0.33_wp, 0.19_wp,            & ! 12
     723                                 0.00_wp, 0.00_wp, 0.00_wp, 0.00_wp,            & ! 13
     724                                 0.25_wp, 0.34_wp, 0.27_wp, 0.11_wp,            & ! 14
     725                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 15
     726                                 0.23_wp, 0.36_wp, 0.30_wp, 0.11_wp,            & ! 16
     727                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp,            & ! 17
     728                                 0.19_wp, 0.35_wp, 0.36_wp, 0.10_wp             & ! 18
    669729                                 /), (/ 4, 18 /) )
    670730
     
    682742                      1.30_wp,  0.400_wp, 1.20_wp,  0.93E-6_wp, 0.766_wp, 0.663_wp, 0.267_wp, 0.010_wp & ! 6
    683743                                                                     /), (/ 8, 6 /) )
     744
    684745
    685746!
     
    747808!
    748809!-- TO BE FILLED
    749 !-- Water parameters  temperature,     z0,      z0h, albedo_type, emissivity,
    750     REAL(wp), DIMENSION(0:4,1:5), PARAMETER :: water_pars = RESHAPE( (/ &
    751                          283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp, & ! 1
    752                          283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp, & ! 2
    753                          283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp, & ! 3
    754                          283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp, & ! 4
    755                          283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp  & ! 5
    756                                                                 /), (/ 5, 5 /) )                                                                     
     810!-- Water parameters                    temperature,     z0,      z0h, albedo_type, emissivity,
     811    REAL(wp), DIMENSION(0:6,1:5), PARAMETER :: water_pars = RESHAPE( (/ &
     812       283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 1
     813       283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 2
     814       283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 3
     815       283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp, & ! 4
     816       283.0_wp, 0.01_wp, 0.001_wp, 1.0E10_wp, 1.0E10_wp, 1.0_wp, 0.99_wp  & ! 5
     817                                                                     /), (/ 7, 5 /) )                                                                   
    757818                                                                                                                                     
    758819    SAVE
     
    764825!
    765826!-- Public functions
    766     PUBLIC lsm_check_data_output, lsm_check_data_output_pr,                    &
     827    PUBLIC lsm_boundary_condition, lsm_check_data_output,                      &
     828           lsm_check_data_output_pr,                                           &
    767829           lsm_check_parameters, lsm_define_netcdf_grid, lsm_3d_data_averaging,&
    768830           lsm_data_output_2d, lsm_data_output_3d, lsm_energy_balance,         &
    769831           lsm_header, lsm_init, lsm_init_arrays, lsm_parin, lsm_soil_model,   &
    770            lsm_swap_timelevel, lsm_read_restart_data, lsm_last_actions
     832           lsm_swap_timelevel, lsm_read_restart_data, lsm_write_restart_data
    771833! !vegetat
    772834!-- Public parameters, constants and initial values
     
    781843    PUBLIC m_soil_h, t_soil_h
    782844
     845    INTERFACE lsm_boundary_condition
     846       MODULE PROCEDURE lsm_boundary_condition
     847    END INTERFACE lsm_boundary_condition
    783848
    784849    INTERFACE lsm_check_data_output
     
    842907    END INTERFACE lsm_read_restart_data
    843908
    844     INTERFACE lsm_last_actions
    845        MODULE PROCEDURE lsm_last_actions
    846     END INTERFACE lsm_last_actions
     909    INTERFACE lsm_write_restart_data
     910       MODULE PROCEDURE lsm_write_restart_data
     911    END INTERFACE lsm_write_restart_data
    847912
    848913 CONTAINS
     914
     915
     916!------------------------------------------------------------------------------!
     917! Description:
     918! ------------
     919!> Set internal Neumann boundary condition at outer soil grid points
     920!> for temperature and humidity.
     921!------------------------------------------------------------------------------!
     922 SUBROUTINE lsm_boundary_condition
     923 
     924    IMPLICIT NONE
     925
     926    INTEGER(iwp) :: i      !< grid index x-direction
     927    INTEGER(iwp) :: ioff   !< offset index x-direction indicating location of soil grid point
     928    INTEGER(iwp) :: j      !< grid index y-direction
     929    INTEGER(iwp) :: joff   !< offset index x-direction indicating location of soil grid point
     930    INTEGER(iwp) :: k      !< grid index z-direction
     931    INTEGER(iwp) :: koff   !< offset index x-direction indicating location of soil grid point
     932    INTEGER(iwp) :: l      !< running index surface-orientation
     933    INTEGER(iwp) :: m      !< running index surface elements
     934
     935    koff = surf_lsm_h%koff
     936    DO  m = 1, surf_lsm_h%ns
     937       i = surf_lsm_h%i(m)
     938       j = surf_lsm_h%j(m)
     939       k = surf_lsm_h%k(m)
     940       pt(k+koff,j,i) = pt(k,j,i)
     941    ENDDO
     942
     943    DO  l = 0, 3
     944       ioff = surf_lsm_v(l)%ioff
     945       joff = surf_lsm_v(l)%joff
     946       DO  m = 1, surf_lsm_v(l)%ns
     947          i = surf_lsm_v(l)%i(m)
     948          j = surf_lsm_v(l)%j(m)
     949          k = surf_lsm_v(l)%k(m)
     950          pt(k,j+joff,i+ioff) = pt(k,j,i)
     951       ENDDO
     952    ENDDO
     953!
     954!-- In case of humidity, set boundary conditions also for q and vpt.
     955    IF ( humidity )  THEN
     956       koff = surf_lsm_h%koff
     957       DO  m = 1, surf_lsm_h%ns
     958          i = surf_lsm_h%i(m)
     959          j = surf_lsm_h%j(m)
     960          k = surf_lsm_h%k(m)
     961          q(k+koff,j,i)   = q(k,j,i)
     962          vpt(k+koff,j,i) = vpt(k,j,i)
     963       ENDDO
     964
     965       DO  l = 0, 3
     966          ioff = surf_lsm_v(l)%ioff
     967          joff = surf_lsm_v(l)%joff
     968          DO  m = 1, surf_lsm_v(l)%ns
     969             i = surf_lsm_v(l)%i(m)
     970             j = surf_lsm_v(l)%j(m)
     971             k = surf_lsm_v(l)%k(m)
     972             q(k,j+joff,i+ioff)   = q(k,j,i)
     973             vpt(k,j+joff,i+ioff) = vpt(k,j,i)
     974          ENDDO
     975       ENDDO
     976    ENDIF
     977
     978 END SUBROUTINE lsm_boundary_condition
    849979
    850980!------------------------------------------------------------------------------!
     
    9751105
    9761106 END SUBROUTINE lsm_check_data_output
     1107
    9771108
    9781109
     
    10601191        ONLY:  bc_pt_b, bc_q_b, constant_flux_layer, message_string,           &
    10611192               most_method
    1062                  
    1063     USE radiation_model_mod,                                                   &
    1064         ONLY:  radiation
    1065    
     1193                     
    10661194   
    10671195    IMPLICIT NONE
     
    10751203         TRIM( surface_type ) /= 'pavement'    .AND.                           &
    10761204         TRIM( surface_type ) /= 'water'       .AND.                           &
    1077          TRIM( surface_type ) /= 'netcdf' )  THEN
    1078           message_string = 'unknown surface type surface_type = "' //          &
     1205         TRIM( surface_type ) /= 'netcdf' )  THEN 
     1206       message_string = 'unknown surface type surface_type = "' //             &
    10791207                        TRIM( surface_type ) // '"'
    1080                 CALL message( 'check_parameters', 'PA0019', 1, 2, 0, 6, 0 )
     1208       CALL message( 'check_parameters', 'PA0019', 1, 2, 0, 6, 0 )
    10811209    ENDIF
    10821210
     
    11641292          ENDIF
    11651293       ENDIF
    1166      
     1294
    11671295       IF ( vegetation_type == 1 )  THEN
    11681296          IF ( vegetation_coverage /= 9999999.9_wp  .AND.  vegetation_coverage &
     
    12601388   
    12611389    ENDIF
     1390
     1391    IF ( TRIM( surface_type ) == 'netcdf' )  THEN
     1392       IF ( ANY( water_type_f%var /= water_type_f%fill )  .AND.                &
     1393            TRIM( most_method ) == 'lookup' )  THEN   
     1394          WRITE( message_string, * ) 'water-surfaces are not allowed in ' //   &
     1395                                     'combination with most_method = ',        &
     1396                                     TRIM( most_method )
     1397          CALL message( 'check_parameters', 'PA0999', 2, 2, 0, 6, 0 )
     1398       ENDIF
     1399!
     1400!--    MS: Some problme here, after calling message everythings stucks at
     1401!--        MPI_FINALIZE call.
     1402       IF ( ANY( pavement_type_f%var /= pavement_type_f%fill )  .AND.           &
     1403            ANY( dz_soil /= 9999999.9_wp ) )  THEN
     1404          message_string = 'pavement-surfaces are not allowed in ' //           &
     1405                           'combination with a non-default setting of dz_soil'
     1406          CALL message( 'check_parameters', 'PA0999', 2, 2, 0, 6, 0 )
     1407       ENDIF
     1408    ENDIF
    12621409   
    12631410!
    12641411!-- Temporary message as long as NetCDF input is not available
    1265     IF ( TRIM( surface_type ) == 'netcdf' )  THEN
    1266              message_string = 'surface_type = netcdf '//                       &
    1267                               'is not supported at the moment'
    1268              CALL message( 'check_parameters', 'PA0465', 1, 2, 0, 6, 0 )
     1412    IF ( TRIM( surface_type ) == 'netcdf'  .AND.  .NOT.  input_pids_static )   &
     1413    THEN
     1414       message_string = 'surface_type = netcdf requires static input file.'
     1415       CALL message( 'check_parameters', 'PA0465', 1, 2, 0, 6, 0 )
    12691416    ENDIF
    12701417
     
    13621509    DO  k = nzb_soil, nzt_soil
    13631510       IF ( soil_moisture(k) > saturation_moisture )  THEN
    1364           message_string = 'soil_moisture must not exceed its saturation' //   &
     1511          message_string = 'soil_moisture must not exceed its saturation' //    &
    13651512                            ' value'
    13661513          CALL message( 'check_parameters', 'PA0458', 1, 2, 0, 6, 0 )
     
    13781525    ALLOCATE ( dz_soil_center(nzb_soil:nzt_soil) )
    13791526    ALLOCATE ( zs(nzb_soil:nzt_soil+1) )
     1527
    13801528
    13811529    zs(nzb_soil) = 0.5_wp * dz_soil(nzb_soil)
     
    14281576    INTEGER(iwp) ::  k         !< running index
    14291577    INTEGER(iwp) ::  k_off     !< offset to determine index of surface element, seen from atmospheric grid point, for z
    1430     INTEGER(iwp) ::  k_rad     !< index to access radiation array
    14311578    INTEGER(iwp) ::  ks        !< running index
    14321579    INTEGER(iwp) ::  l         !< surface-facing index
     
    14541601                lambda_soil, & !< Thermal conductivity of the uppermost soil layer (W/m2/K)
    14551602                lambda_surface, & !< Current value of lambda_surface (W/m2/K)
    1456                 m_liq_max,   & !< maxmimum value of the liq. water reservoir
    1457                 pt1,         & !< potential temperature at first grid level
    1458                 qv1            !< specific humidity at first grid level
     1603                m_liq_max      !< maxmimum value of the liq. water reservoir
    14591604
    14601605    TYPE(surf_type_lsm), POINTER ::  surf_t_surface
     
    14721617    IF ( horizontal )  THEN
    14731618       surf              => surf_lsm_h
     1619
    14741620       surf_t_surface    => t_surface_h
    14751621       surf_t_surface_p  => t_surface_h_p
     
    14801626       surf_m_soil       => m_soil_h
    14811627       surf_t_soil       => t_soil_h
    1482 
    1483        k_off     = -1
    1484        j_off     = 0
    1485        i_off     = 0
    14861628    ELSE
    14871629       surf              => surf_lsm_v(l)
     1630
    14881631       surf_t_surface    => t_surface_v(l)
    14891632       surf_t_surface_p  => t_surface_v_p(l)
     
    14941637       surf_m_soil       => m_soil_v(l)
    14951638       surf_t_soil       => t_soil_v(l)
    1496 
    1497        k_off = 0
    1498        IF ( l == 0 )  THEN
    1499           j_off = -1
    1500           i_off = 0
    1501        ELSEIF ( l == 1 )  THEN
    1502           j_off = 1
    1503           i_off = 0
    1504        ELSEIF ( l == 2 )  THEN
    1505           j_off = 0
    1506           i_off = -1
    1507        ELSEIF ( l == 3 )  THEN
    1508           j_off = 0
    1509           i_off = 1
    1510        ENDIF
    15111639    ENDIF
     1640
     1641!
     1642!-- Index offset of surface element point with respect to adjoining
     1643!-- atmospheric grid point
     1644    k_off = surf%koff
     1645    j_off = surf%joff
     1646    i_off = surf%ioff
    15121647
    15131648!
     
    15201655       j   = surf%j(m)
    15211656       k   = surf%k(m)
    1522 !
    1523 !--    Determine height index for radiation. Note, in clear-sky case radiation
    1524 !--    arrays have rank 0 in first dimensions, so index must be zero. In case
    1525 !--    of RRTMG radiation arrays have non-zero rank in first dimension, so that
    1526 !--    radiation can be obtained at respective height level
    1527        k_rad = MERGE( 0, k + k_off, radiation_scheme /= 'rrtmg' )
    15281657
    15291658!
     
    15891718!--    time step is used here. Note that this formulation is the
    15901719!--    equivalent to the ECMWF formulation using drag coefficients
    1591        IF ( cloud_physics )  THEN
    1592           pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
    1593           qv1 = q(k,j,i) - ql(k,j,i)
    1594        ELSEIF ( cloud_droplets ) THEN
    1595           pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
    1596           qv1 = q(k,j,i)
    1597        ELSE
    1598           pt1 = pt(k,j,i)
    1599           IF ( humidity )  THEN
    1600              qv1 = q(k,j,i)
    1601           ELSE
    1602              qv1 = 0.0_wp
    1603           ENDIF
    1604        ENDIF
     1720!        IF ( cloud_physics )  THEN
     1721!           pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
     1722!           qv1 = q(k,j,i) - ql(k,j,i)
     1723!        ELSEIF ( cloud_droplets ) THEN
     1724!           pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
     1725!           qv1 = q(k,j,i)
     1726!        ELSE
     1727!           pt1 = pt(k,j,i)
     1728!           IF ( humidity )  THEN
     1729!              qv1 = q(k,j,i)
     1730!           ELSE
     1731!              qv1 = 0.0_wp
     1732!           ENDIF
     1733!        ENDIF
    16051734!
    16061735!--    Calculate aerodynamical resistance. For horizontal and vertical
     
    16081737!--    can be obtain via parameterization of Mason (2000) /
    16091738!--    Krayenhoff and Voogt (2006).
    1610 !--    To do: detailed investigation which approach is better!
     1739!--    To do: detailed investigation which approach gives more reliable results!
     1740!--    Please note, in case of very small friction velocity, e.g. in little
     1741!--    holes, the resistance can become negative. For this reason, limit r_a
     1742!--    to positive values.
    16111743       IF ( horizontal  .OR.  .NOT. aero_resist_kray )  THEN
    1612           surf%r_a(m) = ( pt1 - surf_lsm_h%pt_surface(m) ) /                   &
    1613                         ( surf%ts(m) * surf%us(m) + 1.0E-20_wp )
     1744          surf%r_a(m) = ABS( ( surf%pt1(m) - surf%pt_surface(m) ) /            &
     1745                             ( surf%ts(m) * surf%us(m) + 1.0E-20_wp ) )
    16141746       ELSE
    16151747          surf%r_a(m) = 1.0_wp / ( 11.8_wp + 4.2_wp *                          &
    1616                         SQRT( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 +      &
    1617                               ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 +      &
    1618                               ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2 )      &
     1748                        SQRT( MAX( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 + &
     1749                                   ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 + &
     1750                                   ( ( w(k,j,i) + w(k-1,j,i) ) * 0.5_wp )**2,  &
     1751                              0.01_wp ) )                                      &
    16191752                                 )
    16201753       ENDIF
    16211754!
    16221755!--    Make sure that the resistance does not drop to zero for neutral
    1623 !--    stratification
    1624        IF ( ABS( surf%r_a(m) ) < 1.0_wp )  surf%r_a(m) = 1.0_wp
     1756!--    stratification.
     1757       IF ( surf%r_a(m) < 1.0_wp )  surf%r_a(m) = 1.0_wp
    16251758!
    16261759!--    Second step: calculate canopy resistance r_canopy
     
    16291762!--    f1: correction for incoming shortwave radiation (stomata close at
    16301763!--    night)
    1631        f1 = MIN( 1.0_wp, ( 0.004_wp * rad_sw_in(k_rad,j+j_off,i+i_off)         &
    1632                  + 0.05_wp ) / (0.81_wp * (0.004_wp                            &
    1633                  * rad_sw_in(k_rad,j+j_off,i+i_off) + 1.0_wp)) )
     1764       f1 = MIN( 1.0_wp, ( 0.004_wp * surf%rad_sw_in(m) + 0.05_wp ) /          &
     1765                        (0.81_wp * (0.004_wp * surf%rad_sw_in(m)               &
     1766                         + 1.0_wp)) )
    16341767!
    16351768!--    f2: correction for soil moisture availability to plants (the
     
    16651798!
    16661799!--       Calculate vapour pressure
    1667           e  = qv1 * surface_pressure / ( qv1 + 0.622_wp )
     1800          e  = surf%qv1(m) * surface_pressure / ( surf%qv1(m) + 0.622_wp )
    16681801          f3 = EXP ( - surf%g_d(m) * (e_s - e) )
    16691802       ELSE
     
    16821815
    16831816
    1684        f2 = ( surf_m_soil%var_2d(nzb_soil,m) - m_min )                         &
    1685             / ( surf%m_fc(nzb_soil,m) - m_min )
     1817       f2 = ( surf_m_soil%var_2d(nzb_soil,m) - m_min ) /                       &
     1818            ( surf%m_fc(nzb_soil,m) - m_min )
    16861819       f2 = MAX( f2, 1.0E-20_wp )
    16871820       f2 = MIN( f2, 1.0_wp     )
     
    16981831       IF ( surf%pavement_surface(m) )  THEN
    16991832          m_liq_max = m_max_depth * 5.0_wp
    1700           surf%c_liq(m) = MIN( 1.0_wp, (surf_m_liq%var_1d(m) / m_liq_max)**0.67 )
     1833          surf%c_liq(m) = MIN( 1.0_wp, ( surf_m_liq%var_1d(m) / m_liq_max)**0.67 )
    17011834       ELSE
    17021835          m_liq_max = m_max_depth * ( surf%c_veg(m) * surf%lai(m)              &
     
    17101843!--    In case of dewfall, set evapotranspiration to zero
    17111844!--    All super-saturated water is then removed from the air
    1712        IF ( humidity  .AND.  q_s <= qv1 )  THEN
     1845       IF ( humidity  .AND.  q_s <= surf%qv1(m) )  THEN
    17131846          surf%r_canopy(m) = 0.0_wp
    17141847          surf%r_soil(m)   = 0.0_wp
     
    17471880       dq_s_dt = 0.622_wp * e_s_dt / ( surface_pressure - e_s_dt )
    17481881!
    1749 !--    Add LW up so that it can be removed in prognostic equation
    1750        surf%rad_net_l(m) = rad_net(j,i) + rad_lw_out(nzb,j,i)
     1882!--    Calculate net radiation radiation without longwave outgoing flux because
     1883!--    it has a dependency on surface temperature and thus enters the prognostic
     1884!--    equations directly
     1885       surf%rad_net_l(m) = surf%rad_sw_in(m) - surf%rad_sw_out(m)              &
     1886                           + surf%rad_lw_in(m)
    17511887!
    17521888!--    Calculate new skin temperature
     
    17541890!
    17551891!--       Numerator of the prognostic equation
    1756           coef_1 = surf%rad_net_l(m) + rad_lw_out_change_0(j,i)                &
    1757                    * surf_t_surface%var_1d(m) - rad_lw_out(nzb,j,i)            &
    1758                    + f_shf * pt1 + f_qsws * ( qv1 - q_s                        &
     1892          coef_1 = surf%rad_net_l(m) + surf%rad_lw_out_change_0(m)             &
     1893                   * surf_t_surface%var_1d(m) - surf%rad_lw_out(m)             &
     1894                   + f_shf * surf%pt1(m) + f_qsws * ( surf%qv1(m) - q_s        &
    17591895                   + dq_s_dt * surf_t_surface%var_1d(m) ) + lambda_surface     &
    17601896                   * surf_t_soil%var_2d(nzb_soil,m)
     1897
    17611898!
    17621899!--       Denominator of the prognostic equation
    1763           coef_2 = rad_lw_out_change_0(j,i) + f_qsws * dq_s_dt                 &
     1900          coef_2 = surf%rad_lw_out_change_0(m) + f_qsws * dq_s_dt              &
    17641901                   + lambda_surface + f_shf / exn
    17651902       ELSE
    17661903!
    17671904!--       Numerator of the prognostic equation
    1768           coef_1 = surf%rad_net_l(m) + rad_lw_out_change_0(j,i)                &
    1769                    * surf_t_surface%var_1d(m) - rad_lw_out(nzb,j,i)            &
    1770                    + f_shf * pt1  + lambda_surface                             &
     1905          coef_1 = surf%rad_net_l(m) + surf%rad_lw_out_change_0(m)             &
     1906                   * surf_t_surface%var_1d(m) - surf%rad_lw_out(m)             &
     1907                   + f_shf * surf%pt1(m)  + lambda_surface                     &
    17711908                   * surf_t_soil%var_2d(nzb_soil,m)
    17721909!
    17731910!--       Denominator of the prognostic equation
    1774           coef_2 = rad_lw_out_change_0(j,i) + lambda_surface + f_shf / exn
     1911          coef_2 = surf%rad_lw_out_change_0(m) + lambda_surface + f_shf / exn
    17751912
    17761913       ENDIF
     1914
    17771915
    17781916       tend = 0.0_wp
     
    18191957!--    computing time for the radiation code).
    18201958       IF ( ABS( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) )       &
    1821             > 0.2_wp  .AND. unscheduled_radiation_calls )  THEN
     1959            > 0.2_wp  .AND. &
     1960            unscheduled_radiation_calls )  THEN
    18221961          force_radiation_call_l = .TRUE.
    18231962       ENDIF
    18241963
    1825        pt(k+k_off,j+j_off,i+i_off) = surf_t_surface_p%var_1d(m) / exn  !is actually no air temperature
     1964
     1965!        IF ( surf_t_surface_p%var_1d(m) < 270.0_wp .OR. surf_t_surface_p%var_1d(m) > 350.0_wp )  THEN
     1966!           PRINT*, "myid: ", myid
     1967!           PRINT*, "m: ", m
     1968!           PRINT*, "i,j,k: ", i, j, k
     1969!           PRINT*, "pt_p: ", surf_t_surface_p%var_1d(m)
     1970!           PRINT*, "f_shf: ", f_shf
     1971!           PRINT*, "f_qsws: ", f_qsws
     1972!         ENDIF
     1973
     1974
     1975!        pt(k+k_off,j+j_off,i+i_off) = surf_t_surface_p%var_1d(m) / exn  !is actually no air temperature
    18261976       surf%pt_surface(m)          = surf_t_surface_p%var_1d(m) / exn
    18271977
     
    18291979!--    Calculate fluxes
    18301980       surf%rad_net_l(m) = surf%rad_net_l(m) +                                 &
    1831                             rad_lw_out_change_0(j,i)                           &
    1832                           * surf_t_surface%var_1d(m) - rad_lw_out(nzb,j,i)     &
    1833                           - rad_lw_out_change_0(j,i) * surf_t_surface_p%var_1d(m)
    1834 
    1835        rad_net(j,i) = surf%rad_net_l(m)
    1836        rad_lw_out(nzb,j,i) = rad_lw_out(nzb,j,i) + rad_lw_out_change_0(j,i) * &
     1981                            surf%rad_lw_out_change_0(m)                        &
     1982                          * surf_t_surface%var_1d(m) - surf%rad_lw_out(m)      &
     1983                          - surf%rad_lw_out_change_0(m) * surf_t_surface_p%var_1d(m)
     1984
     1985       surf%rad_net(m) = surf%rad_net_l(m)
     1986       surf%rad_lw_out(m) = surf%rad_lw_out(m) + surf%rad_lw_out_change_0(m) * &
    18371987                     ( surf_t_surface_p%var_1d(m) - surf_t_surface%var_1d(m) )
    18381988
     
    18401990                                             - surf_t_soil%var_2d(nzb_soil,m) )
    18411991
    1842        surf%shf(m) = - f_shf * ( pt1 - surf%pt_surface(m) ) / cp
    1843 
     1992       surf%shf(m) = - f_shf * ( surf%pt1(m) - surf%pt_surface(m) ) / cp
    18441993
    18451994       IF ( humidity )  THEN
    1846           surf%qsws(m)  = - f_qsws * ( qv1 - q_s + dq_s_dt                     &
     1995          surf%qsws(m)  = - f_qsws * ( surf%qv1(m) - q_s + dq_s_dt             &
    18471996                          * surf_t_surface%var_1d(m) - dq_s_dt *               &
    18481997                            surf_t_surface_p%var_1d(m) )
    18491998
    1850           surf%qsws_veg(m)  = - f_qsws_veg  * ( qv1 - q_s                      &
     1999          surf%qsws_veg(m)  = - f_qsws_veg  * ( surf%qv1(m) - q_s              &
    18512000                              + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt   &
    18522001                              * surf_t_surface_p%var_1d(m) )
    18532002
    1854           surf%qsws_soil(m) = - f_qsws_soil * ( qv1 - q_s                      &
     2003          surf%qsws_soil(m) = - f_qsws_soil * ( surf%qv1(m) - q_s              &
    18552004                              + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt   &
    18562005                              * surf_t_surface_p%var_1d(m) )
    18572006
    1858           surf%qsws_liq(m)  = - f_qsws_liq  * ( qv1 - q_s                      &
     2007          surf%qsws_liq(m)  = - f_qsws_liq  * ( surf%qv1(m) - q_s              &
    18592008                              + dq_s_dt * surf_t_surface%var_1d(m) - dq_s_dt   &
    18602009                              * surf_t_surface_p%var_1d(m) )
     
    18662015          surf%r_s(m) = 1.0E10_wp
    18672016       ELSE
    1868           surf%r_s(m) = - rho_lv * ( qv1 - q_s + dq_s_dt                       &
     2017          surf%r_s(m) = - rho_lv * ( surf%qv1(m) - q_s + dq_s_dt               &
    18692018                          * surf_t_surface%var_1d(m) - dq_s_dt *               &
    18702019                            surf_t_surface_p%var_1d(m) ) /                     &
     
    19542103                      intermediate_timestep_count_max )  THEN
    19552104                surf_tm_liq_m%var_1d(m) = -9.5625_wp * tend +                  &
    1956                                               5.3125_wp * surf_tm_liq_m%var_1d(m)
     2105                                           5.3125_wp * surf_tm_liq_m%var_1d(m)
    19572106             ENDIF
    19582107          ENDIF
     
    20152164          q_s = 0.622_wp * e_s / ( surface_pressure - e_s )
    20162165
    2017           resistance = surf%r_a(m) / ( surf%r_a(m) + surf%r_s(m) )
     2166          resistance = surf%r_a(m) / ( surf%r_a(m) + surf%r_s(m) + 1E-5_wp )
    20182167
    20192168!
     
    20282177                                          q(k,j,i)
    20292178          ENDIF
    2030 
    20312179!
    20322180!--       Update virtual potential temperature
     
    20982246       ENDIF
    20992247
    2100        WRITE( io, 4 ) TRIM( vegetation_type_name(vegetation_type) ),           &
    2101                         TRIM (soil_type_name(soil_type) )
    2102        WRITE( io, 5 ) TRIM( soil_depth_chr ), TRIM( t_soil_chr ),              &
     2248       IF ( vegetation_type_f%from_file )  THEN
     2249          WRITE( io, 5 )
     2250       ELSE
     2251          WRITE( io, 4 ) TRIM( vegetation_type_name(vegetation_type) ),        &
     2252                         TRIM (soil_type_name(soil_type) )
     2253       ENDIF
     2254       WRITE( io, 6 ) TRIM( soil_depth_chr ), TRIM( t_soil_chr ),              &
    21032255                        TRIM( m_soil_chr ), TRIM( roots_chr ),                 &
    21042256                        TRIM( vertical_index_chr )
     
    211122634   FORMAT ('    --> Land surface type  : ',A,/                                &
    21122264            '    --> Soil porosity type : ',A)
    2113 5   FORMAT (/'    Initial soil temperature and moisture profile:'//            &
     22655   FORMAT ('    --> Land surface type  : read from file' /                    &
     2266            '    --> Soil porosity type : read from file' )
     22676   FORMAT (/'    Initial soil temperature and moisture profile:'//            &
    21142268            '       Height:        ',A,'  m'/                                  &
    21152269            '       Temperature:   ',A,'  K'/                                  &
     
    21172271            '       Root fraction: ',A,'  '/                                   &
    21182272            '       Grid point:    ',A)
     2273
    21192274
    21202275    END SUBROUTINE lsm_header
     
    21332288       IMPLICIT NONE
    21342289
    2135        INTEGER(iwp) ::  i     !< running index
    2136        INTEGER(iwp) ::  i_off !< index offset of surface element, seen from atmospheric grid point
    2137        INTEGER(iwp) ::  j     !< running index
    2138        INTEGER(iwp) ::  j_off !< index offset of surface element, seen from atmospheric grid point
    2139        INTEGER(iwp) ::  k     !< running index
    2140        INTEGER(iwp) ::  kn    !< running index
    2141        INTEGER(iwp) ::  ko    !< running index
    2142        INTEGER(iwp) ::  kroot !< running index
    2143        INTEGER(iwp) ::  kzs   !< running index
    2144        INTEGER(iwp) ::  l     !< running index surface facing
    2145        INTEGER(iwp) ::  m     !< running index
     2290       INTEGER(iwp) ::  i                       !< running index
     2291       INTEGER(iwp) ::  i_off                   !< index offset of surface element, seen from atmospheric grid point
     2292       INTEGER(iwp) ::  j                       !< running index
     2293       INTEGER(iwp) ::  j_off                   !< index offset of surface element, seen from atmospheric grid point
     2294       INTEGER(iwp) ::  k                       !< running index
     2295       INTEGER(iwp) ::  kn                      !< running index
     2296       INTEGER(iwp) ::  ko                      !< running index
     2297       INTEGER(iwp) ::  kroot                   !< running index
     2298       INTEGER(iwp) ::  kzs                     !< running index
     2299       INTEGER(iwp) ::  l                       !< running index surface facing
     2300       INTEGER(iwp) ::  m                       !< running index
     2301       INTEGER(iwp) ::  st                      !< soil-type index
    21462302       INTEGER(iwp) ::  n_soil_layers_total     !< temperature variable, stores the total number of soil layers + 4
    2147 
    2148 
    2149        REAL(wp) :: pt1   !< potential temperature at first grid level
     2303       INTEGER(iwp) ::  n_surf                  !< number of surface types of given surface element
    21502304
    21512305       REAL(wp), DIMENSION(:), ALLOCATABLE ::  bound, bound_root_fr  !< temporary arrays for storing index bounds
     
    21572311!--    If no cloud physics is used, rho_surface has not been calculated before
    21582312       IF (  .NOT.  cloud_physics )  THEN
    2159           rho_surface = surface_pressure * 100.0_wp / ( r_d * pt_surface * exn )
     2313          CALL calc_mean_profile( pt, 4 )
     2314          rho_surface = surface_pressure * 100.0_wp / ( r_d * hom(nzb+1,1,4,0) * exn )
    21602315       ENDIF
    21612316
     
    22062361          surf_lsm_v(l)%r_soil     = 0.0_wp
    22072362       ENDDO
    2208    
     2363
     2364!
     2365!--    Set initial values for prognostic soil quantities
     2366       IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
     2367          t_soil_h%var_2d = 0.0_wp
     2368          m_soil_h%var_2d = 0.0_wp
     2369          m_liq_h%var_1d  = 0.0_wp
     2370
     2371          DO  l = 0, 3
     2372             t_soil_v(l)%var_2d = 0.0_wp
     2373             m_soil_v(l)%var_2d = 0.0_wp
     2374             m_liq_v(l)%var_1d  = 0.0_wp
     2375          ENDDO
     2376       ENDIF
    22092377!
    22102378!--    Allocate 3D soil model arrays
     
    22612429          ENDIF     
    22622430       ENDDO
    2263 
     2431!
     2432!--    Allocate albedo type and emissivity for vegetation, water and pavement
     2433!--    fraction.
     2434!--    Set default values at each surface element.
     2435       ALLOCATE ( surf_lsm_h%albedo_type(0:2,1:surf_lsm_h%ns) )
     2436       ALLOCATE ( surf_lsm_h%emissivity(0:2,1:surf_lsm_h%ns) )
     2437       surf_lsm_h%albedo_type = 0     
     2438       surf_lsm_h%emissivity  = emissivity
     2439       DO  l = 0, 3
     2440          ALLOCATE ( surf_lsm_v(l)%albedo_type(0:2,1:surf_lsm_v(l)%ns) )
     2441          ALLOCATE ( surf_lsm_v(l)%emissivity(0:2,1:surf_lsm_v(l)%ns)  )
     2442          surf_lsm_v(l)%albedo_type = 0
     2443          surf_lsm_v(l)%emissivity  = emissivity
     2444       ENDDO
     2445!
     2446!--    Allocate arrays for relative surface fraction.
     2447!--    0 - vegetation fraction, 2 - water fraction, 1 - pavement fraction
     2448       ALLOCATE( surf_lsm_h%frac(0:2,1:surf_lsm_h%ns) )
     2449       surf_lsm_h%frac = 0.0_wp
     2450       DO  l = 0, 3
     2451          ALLOCATE( surf_lsm_v(l)%frac(0:2,1:surf_lsm_v(l)%ns) )
     2452          surf_lsm_v(l)%frac = 0.0_wp
     2453       ENDDO
     2454!
     2455!--    For vertical walls only - allocate special flag indicating if any building is on
     2456!--    top of any natural surfaces. Used for initialization only.
     2457       DO  l = 0, 3
     2458          ALLOCATE( surf_lsm_v(l)%building_covered(1:surf_lsm_v(l)%ns) )
     2459       ENDDO
    22642460!
    22652461!--    Set flag parameter for the prescribed surface type depending on user
    2266 !--    input.
    2267        DO  m = 1, surf_lsm_h%ns
    2268 
    2269           SELECT CASE ( TRIM( surface_type ) )
     2462!--    input. Set surface fraction to 1 for the respective type.
     2463       SELECT CASE ( TRIM( surface_type ) )
    22702464         
    2271              CASE ( 'vegetation' )
     2465          CASE ( 'vegetation' )
    22722466         
    2273                 surf_lsm_h%vegetation_surface(m)    = .TRUE.
     2467             surf_lsm_h%vegetation_surface = .TRUE.
     2468             surf_lsm_h%frac(ind_veg,:) = 1.0_wp
     2469             DO  l = 0, 3
     2470                surf_lsm_v(l)%vegetation_surface = .TRUE.
     2471                surf_lsm_v(l)%frac(ind_veg,:) = 1.0_wp
     2472             ENDDO
    22742473   
    2275              CASE ( 'water' )
     2474          CASE ( 'water' )
    22762475             
    2277                 surf_lsm_h%water_surface(m)      = .TRUE.
    2278 
    2279              CASE ( 'pavement' )
     2476             surf_lsm_h%water_surface = .TRUE.
     2477             surf_lsm_h%frac(ind_wat,:) = 1.0_wp
     2478!
     2479!--          Note, vertical water surface does not really make sense.
     2480             DO  l = 0, 3 
     2481                surf_lsm_v(l)%water_surface   = .TRUE.
     2482                surf_lsm_v(l)%frac(ind_wat,:) = 1.0_wp
     2483             ENDDO
     2484
     2485          CASE ( 'pavement' )
    22802486             
    2281                 surf_lsm_h%pavement_surface(m)   = .TRUE.
    2282    
    2283           END SELECT
    2284  
    2285        ENDDO
    2286        
    2287        DO  l = 0, 3
    2288           SELECT CASE ( TRIM( surface_type ) )
    2289          
    2290              CASE ( 'vegetation' )
    2291          
    2292                 surf_lsm_v(l)%vegetation_surface = .TRUE.
    2293    
    2294              CASE ( 'water' )
    2295 
    2296                 surf_lsm_v(l)%water_surface   = .TRUE.
    2297 
    2298              CASE ( 'pavement' )
    2299              
    2300                 surf_lsm_h%pavement_surface   = .TRUE.
    2301 
    2302           END SELECT   
    2303 
    2304        ENDDO
    2305 !
    2306 !--    Initialize standard soil types. It is possible to overwrite each
    2307 !--    parameter by setting the respecticy NAMELIST variable to a
    2308 !--    value /= 9999999.9.
     2487             surf_lsm_h%pavement_surface = .TRUE.
     2488                surf_lsm_h%frac(ind_pav,:) = 1.0_wp
     2489             DO  l = 0, 3
     2490                surf_lsm_v(l)%pavement_surface   = .TRUE.
     2491                surf_lsm_v(l)%frac(ind_pav,:) = 1.0_wp
     2492             ENDDO
     2493
     2494          CASE ( 'netcdf' )
     2495
     2496             DO  m = 1, surf_lsm_h%ns
     2497                i = surf_lsm_h%i(m)
     2498                j = surf_lsm_h%j(m)
     2499                IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill )        &
     2500                   surf_lsm_h%vegetation_surface(m) = .TRUE.
     2501                IF ( pavement_type_f%var(j,i)   /= pavement_type_f%fill )          &
     2502                   surf_lsm_h%pavement_surface(m) = .TRUE.
     2503                IF ( water_type_f%var(j,i)      /= water_type_f%fill )             &
     2504                   surf_lsm_h%water_surface(m) = .TRUE.
     2505             ENDDO
     2506             DO  l = 0, 3
     2507                DO  m = 1, surf_lsm_v(l)%ns
     2508!
     2509!--                Only for vertical surfaces. Check if natural walls at reference
     2510!--                grid point are covered by any building. This case, problems
     2511!--                with initialization will aris if index offsets are used.
     2512!--                In order to deal with this, set special flag.
     2513                   surf_lsm_v(l)%building_covered(m) = .FALSE.
     2514                   IF ( building_type_f%from_file )  THEN
     2515                      i = surf_lsm_v(l)%i(m) + surf_lsm_v(l)%ioff
     2516                      j = surf_lsm_v(l)%j(m) + surf_lsm_v(l)%joff
     2517                      IF ( building_type_f%var(j,i) /= 0 )                     &
     2518                         surf_lsm_v(l)%building_covered(m) = .TRUE.
     2519                   ENDIF
     2520!
     2521!--                Normally proceed with setting surface types.
     2522                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,     &
     2523                                            surf_lsm_v(l)%building_covered(m) )
     2524                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,     &
     2525                                            surf_lsm_v(l)%building_covered(m) )
     2526                   IF ( vegetation_type_f%var(j,i) /= vegetation_type_f%fill ) &
     2527                      surf_lsm_v(l)%vegetation_surface(m) = .TRUE.
     2528                   IF ( pavement_type_f%var(j,i)   /= pavement_type_f%fill )   &
     2529                      surf_lsm_v(l)%pavement_surface(m) = .TRUE.
     2530                   IF ( water_type_f%var(j,i)      /= water_type_f%fill )      &
     2531                      surf_lsm_v(l)%water_surface(m) = .TRUE.
     2532                ENDDO
     2533             ENDDO
     2534
     2535       END SELECT
     2536!
     2537!--    In case of netcdf input file, further initialize surface fractions.
     2538!--    At the moment only 1 surface is given at a location, so that the fraction
     2539!--    is either 0 or 1. This will be revised later.
     2540       IF ( input_pids_static )  THEN
     2541          DO  m = 1, surf_lsm_h%ns
     2542             i = surf_lsm_h%i(m)
     2543             j = surf_lsm_h%j(m)
     2544!
     2545!--          0 - vegetation fraction, 1 - pavement fraction, 2 - water fraction             
     2546             surf_lsm_h%frac(ind_veg,m) = surface_fraction_f%frac(ind_veg,j,i)         
     2547             surf_lsm_h%frac(ind_pav,m) = surface_fraction_f%frac(ind_pav,j,i)       
     2548             surf_lsm_h%frac(ind_wat,m) = surface_fraction_f%frac(ind_wat,j,i)
     2549
     2550          ENDDO
     2551          DO  l = 0, 3
     2552             DO  m = 1, surf_lsm_v(l)%ns
     2553                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
     2554                                                surf_lsm_v(l)%building_covered(m) )
     2555                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
     2556                                                surf_lsm_v(l)%building_covered(m) )
     2557!
     2558!--             0 - vegetation fraction, 1 - pavement fraction, 2 - water fraction       
     2559                surf_lsm_v(l)%frac(ind_veg,m) = surface_fraction_f%frac(ind_veg,j,i)         
     2560                surf_lsm_v(l)%frac(ind_pav,m) = surface_fraction_f%frac(ind_pav,j,i)       
     2561                surf_lsm_v(l)%frac(ind_wat,m) = surface_fraction_f%frac(ind_wat,j,i)
     2562
     2563             ENDDO
     2564          ENDDO
     2565       ENDIF
     2566!
     2567!--    Level 1, initialization of soil parameters.
     2568!--    It is possible to overwrite each parameter by setting the respecticy
     2569!--    NAMELIST variable to a value /= 9999999.9.
    23092570       IF ( soil_type /= 0 )  THEN 
    23102571 
     
    23392600          IF ( residual_moisture == 9999999.9_wp )  THEN
    23402601             residual_moisture = soil_pars(7,soil_type)       
    2341           ENDIF
    2342 
    2343        ENDIF   
    2344  
    2345 !
    2346 !--    Check whether parameters from the lookup tables are to be used
    2347        IF ( vegetation_type /= 0 )  THEN
    2348 
    2349           IF ( min_canopy_resistance == 9999999.9_wp )  THEN
    2350              min_canopy_resistance = vegetation_pars(0,vegetation_type)
    2351           ENDIF
    2352 
    2353           IF ( leaf_area_index == 9999999.9_wp )  THEN
    2354              leaf_area_index = vegetation_pars(1,vegetation_type)         
    2355           ENDIF
    2356 
    2357           IF ( vegetation_coverage == 9999999.9_wp )  THEN
    2358              vegetation_coverage = vegetation_pars(2,vegetation_type)     
    2359           ENDIF
    2360 
    2361           IF ( canopy_resistance_coefficient == 9999999.9_wp )  THEN
    2362               canopy_resistance_coefficient= vegetation_pars(3,vegetation_type)     
    2363           ENDIF
    2364 
    2365           IF ( z0_vegetation == 9999999.9_wp )  THEN
    2366              z0_vegetation  = vegetation_pars(4,vegetation_type)
    2367           ENDIF
    2368 
    2369           IF ( z0h_vegetation == 9999999.9_wp )  THEN
    2370              z0h_vegetation = vegetation_pars(5,vegetation_type)
    2371           ENDIF   
    2372          
    2373           IF ( lambda_surface_stable == 9999999.9_wp )  THEN
    2374              lambda_surface_stable = vegetation_pars(6,vegetation_type)
    2375           ENDIF
    2376 
    2377           IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
    2378              lambda_surface_unstable = vegetation_pars(7,vegetation_type)           
    2379           ENDIF
    2380 
    2381           IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
    2382              f_shortwave_incoming = vegetation_pars(8,vegetation_type)       
    2383           ENDIF
    2384 
    2385           IF ( c_surface == 9999999.9_wp )  THEN
    2386              c_surface = vegetation_pars(9,vegetation_type)       
    2387           ENDIF
    2388  
    2389           IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
    2390              albedo_type = INT(vegetation_pars(10,vegetation_type))       
    2391           ENDIF
    2392    
    2393           IF ( emissivity == 9999999.9_wp )  THEN
    2394              emissivity = vegetation_pars(11,vegetation_type)     
    2395           ENDIF
    2396          
    2397        ENDIF
    2398 
    2399        IF ( water_type /= 0 )  THEN
    2400 
    2401           IF ( water_temperature == 9999999.9_wp )  THEN
    2402              water_temperature = water_pars(0,water_type)       
    2403           ENDIF
    2404 
    2405           IF ( z0_water == 9999999.9_wp )  THEN
    2406              z0_water = water_pars(1,water_type)
    2407           ENDIF       
    2408 
    2409           IF ( z0h_water == 9999999.9_wp )  THEN
    2410              z0h_water = water_pars(2,water_type)   
    2411           ENDIF 
    2412 
    2413           IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
    2414              albedo_type = INT(water_pars(3,water_type))       
    2415           ENDIF
    2416    
    2417           IF ( emissivity == 9999999.9_wp )  THEN
    2418              emissivity = water_pars(4,water_type)       
    2419           ENDIF
    2420 
    2421        ENDIF
    2422        
    2423        IF ( pavement_type /= 0 )  THEN 
    2424 
    2425 !
    2426 !--       When a pavement_type is used, overwrite a possible setting of
    2427 !--       the pavement depth as it is already defined by the pavement type
    2428           pavement_depth_level = 0
    2429 
    2430           IF ( z0_pavement == 9999999.9_wp )  THEN
    2431              z0_pavement  = pavement_pars(0,pavement_type)
    2432           ENDIF
    2433 
    2434           IF ( z0h_pavement == 9999999.9_wp )  THEN
    2435              z0h_pavement = pavement_pars(1,pavement_type)
    2436           ENDIF
    2437 
    2438           IF ( pavement_heat_conduct == 9999999.9_wp )  THEN
    2439              pavement_heat_conduct = pavement_subsurface_pars_1(0,pavement_type)
    2440           ENDIF
    2441 
    2442           IF ( pavement_heat_capacity == 9999999.9_wp )  THEN
    2443              pavement_heat_capacity = pavement_subsurface_pars_2(0,pavement_type)
    2444           ENDIF   
    2445    
    2446           IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
    2447              albedo_type = INT(pavement_pars(2,pavement_type))       
    2448           ENDIF
    2449    
    2450           IF ( emissivity == 9999999.9_wp )  THEN
    2451              emissivity = pavement_pars(3,pavement_type)       
    2452           ENDIF
    2453 
    2454 !
    2455 !--       If the depth level of the pavement is not set, determine it from
    2456 !--       lookup table.
    2457           IF ( pavement_depth_level == 0 )  THEN
    2458              DO  k = nzb_soil, nzt_soil 
    2459                 IF ( pavement_subsurface_pars_1(k,pavement_type) == 9999999.9_wp &
    2460                    .OR. pavement_subsurface_pars_2(k,pavement_type)              &
    2461                    == 9999999.9_wp)  THEN
    2462                    nzt_pavement = k-1
    2463                    EXIT
    2464                 ENDIF
    2465              ENDDO
    2466           ELSE
    2467              nzt_pavement = pavement_depth_level
    24682602          ENDIF
    24692603
     
    24942628          surf_lsm_v(l)%r_soil_min    = min_soil_resistance
    24952629       ENDDO
    2496 
     2630!
     2631!--    Level 2, initialization of soil parameters via soil_type read from file.
     2632!--    Soil parameters are initialized for each (y,x)-grid point
     2633!--    individually using default paramter settings according to the given
     2634!--    soil type.
     2635       IF ( soil_type_f%from_file )  THEN
     2636!
     2637!--       Level of detail = 1, i.e. a homogeneous soil distribution along the
     2638!--       vertical dimension is assumed.
     2639          IF ( soil_type_f%lod == 1 )  THEN
     2640!
     2641!--          Horizontal surfaces
     2642             DO  m = 1, surf_lsm_h%ns
     2643                i = surf_lsm_h%i(m)
     2644                j = surf_lsm_h%j(m)
     2645             
     2646                st = soil_type_f%var_2d(j,i)
     2647                IF ( st /= soil_type_f%fill )  THEN
     2648                   surf_lsm_h%alpha_vg(:,m)    = soil_pars(0,st)
     2649                   surf_lsm_h%l_vg(:,m)        = soil_pars(1,st)
     2650                   surf_lsm_h%n_vg(:,m)        = soil_pars(2,st)
     2651                   surf_lsm_h%gamma_w_sat(:,m) = soil_pars(3,st)
     2652                   surf_lsm_h%m_sat(:,m)       = soil_pars(4,st)
     2653                   surf_lsm_h%m_fc(:,m)        = soil_pars(5,st)
     2654                   surf_lsm_h%m_wilt(:,m)      = soil_pars(6,st)
     2655                   surf_lsm_h%m_res(:,m)       = soil_pars(7,st)
     2656                ENDIF
     2657             ENDDO
     2658!
     2659!--          Vertical surfaces ( assumes the soil type given at respective (x,y)
     2660             DO  l = 0, 3
     2661                DO  m = 1, surf_lsm_v(l)%ns
     2662                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
     2663                                                   surf_lsm_v(l)%building_covered(m) )
     2664                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
     2665                                                   surf_lsm_v(l)%building_covered(m) )
     2666
     2667                   st = soil_type_f%var_2d(j,i)
     2668                   IF ( st /= soil_type_f%fill )  THEN
     2669                      surf_lsm_v(l)%alpha_vg(:,m)    = soil_pars(0,st)
     2670                      surf_lsm_v(l)%l_vg(:,m)        = soil_pars(1,st)
     2671                      surf_lsm_v(l)%n_vg(:,m)        = soil_pars(2,st)
     2672                      surf_lsm_v(l)%gamma_w_sat(:,m) = soil_pars(3,st)
     2673                      surf_lsm_v(l)%m_sat(:,m)       = soil_pars(4,st)
     2674                      surf_lsm_v(l)%m_fc(:,m)        = soil_pars(5,st)
     2675                      surf_lsm_v(l)%m_wilt(:,m)      = soil_pars(6,st)
     2676                      surf_lsm_v(l)%m_res(:,m)       = soil_pars(7,st)
     2677                   ENDIF
     2678                ENDDO
     2679             ENDDO
     2680!
     2681!--       Level of detail = 2, i.e. soil type and thus the soil parameters
     2682!--       can be heterogeneous along the vertical dimension.
     2683          ELSE
     2684!
     2685!--          Horizontal surfaces
     2686             DO  m = 1, surf_lsm_h%ns
     2687                i = surf_lsm_h%i(m)
     2688                j = surf_lsm_h%j(m)
     2689             
     2690                DO  k = nzb_soil, nzt_soil
     2691                   st = soil_type_f%var_3d(k,j,i)
     2692                   IF ( st /= soil_type_f%fill )  THEN
     2693                      surf_lsm_h%alpha_vg(k,m)    = soil_pars(0,st)
     2694                      surf_lsm_h%l_vg(k,m)        = soil_pars(1,st)
     2695                      surf_lsm_h%n_vg(k,m)        = soil_pars(2,st)
     2696                      surf_lsm_h%gamma_w_sat(k,m) = soil_pars(3,st)
     2697                      surf_lsm_h%m_sat(k,m)       = soil_pars(4,st)
     2698                      surf_lsm_h%m_fc(k,m)        = soil_pars(5,st)
     2699                      surf_lsm_h%m_wilt(k,m)      = soil_pars(6,st)
     2700                      surf_lsm_h%m_res(k,m)       = soil_pars(7,st)
     2701                   ENDIF
     2702                ENDDO
     2703             ENDDO
     2704!
     2705!--          Vertical surfaces ( assumes the soil type given at respective (x,y)
     2706             DO  l = 0, 3
     2707                DO  m = 1, surf_lsm_v(l)%ns
     2708                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
     2709                                                   surf_lsm_v(l)%building_covered(m) )
     2710                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
     2711                                                   surf_lsm_v(l)%building_covered(m) )
     2712
     2713                   DO  k = nzb_soil, nzt_soil
     2714                      st = soil_type_f%var_3d(k,j,i)
     2715                      IF ( st /= soil_type_f%fill )  THEN
     2716                         surf_lsm_v(l)%alpha_vg(k,m)    = soil_pars(0,st)
     2717                         surf_lsm_v(l)%l_vg(k,m)        = soil_pars(1,st)
     2718                         surf_lsm_v(l)%n_vg(k,m)        = soil_pars(2,st)
     2719                         surf_lsm_v(l)%gamma_w_sat(k,m) = soil_pars(3,st)
     2720                         surf_lsm_v(l)%m_sat(k,m)       = soil_pars(4,st)
     2721                         surf_lsm_v(l)%m_fc(k,m)        = soil_pars(5,st)
     2722                         surf_lsm_v(l)%m_wilt(k,m)      = soil_pars(6,st)
     2723                         surf_lsm_v(l)%m_res(k,m)       = soil_pars(7,st)
     2724                      ENDIF
     2725                   ENDDO
     2726                ENDDO
     2727             ENDDO
     2728          ENDIF
     2729       ENDIF
     2730!
     2731!--    Level 3, initialization of single soil parameters at single z,x,y
     2732!--    position via soil_pars read from file.
     2733       IF ( soil_pars_f%from_file )  THEN
     2734!
     2735!--       Level of detail = 1, i.e. a homogeneous vertical distribution of soil
     2736!--       parameters is assumed.
     2737!--       Horizontal surfaces
     2738          IF ( soil_pars_f%lod == 1 )  THEN
     2739!
     2740!--          Horizontal surfaces
     2741             DO  m = 1, surf_lsm_h%ns
     2742                i = surf_lsm_h%i(m)
     2743                j = surf_lsm_h%j(m)
     2744
     2745                IF ( soil_pars_f%pars_xy(0,j,i) /= soil_pars_f%fill )              &
     2746                   surf_lsm_h%alpha_vg(:,m)    = soil_pars_f%pars_xy(0,j,i)
     2747                IF ( soil_pars_f%pars_xy(1,j,i) /= soil_pars_f%fill )              &
     2748                   surf_lsm_h%l_vg(:,m)        = soil_pars_f%pars_xy(1,j,i)
     2749                IF ( soil_pars_f%pars_xy(2,j,i) /= soil_pars_f%fill )              &
     2750                   surf_lsm_h%n_vg(:,m)        = soil_pars_f%pars_xy(2,j,i)
     2751                IF ( soil_pars_f%pars_xy(3,j,i) /= soil_pars_f%fill )              &
     2752                   surf_lsm_h%gamma_w_sat(:,m) = soil_pars_f%pars_xy(3,j,i)
     2753                IF ( soil_pars_f%pars_xy(4,j,i) /= soil_pars_f%fill )              &
     2754                   surf_lsm_h%m_sat(:,m)       = soil_pars_f%pars_xy(4,j,i)
     2755                IF ( soil_pars_f%pars_xy(5,j,i) /= soil_pars_f%fill )              &
     2756                   surf_lsm_h%m_fc(:,m)        = soil_pars_f%pars_xy(5,j,i)
     2757                IF ( soil_pars_f%pars_xy(6,j,i) /= soil_pars_f%fill )              &
     2758                   surf_lsm_h%m_wilt(:,m)      = soil_pars_f%pars_xy(6,j,i)
     2759                IF ( soil_pars_f%pars_xy(7,j,i) /= soil_pars_f%fill )              &
     2760                   surf_lsm_h%m_res(:,m)       = soil_pars_f%pars_xy(7,j,i)
     2761
     2762             ENDDO
     2763!
     2764!--          Vertical surfaces
     2765             DO  l = 0, 3
     2766                DO  m = 1, surf_lsm_v(l)%ns
     2767                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
     2768                                                   surf_lsm_v(l)%building_covered(m) )
     2769                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
     2770                                                   surf_lsm_v(l)%building_covered(m) )
     2771
     2772                   IF ( soil_pars_f%pars_xy(0,j,i) /= soil_pars_f%fill )           &
     2773                      surf_lsm_v(l)%alpha_vg(:,m)    = soil_pars_f%pars_xy(0,j,i)
     2774                   IF ( soil_pars_f%pars_xy(1,j,i) /= soil_pars_f%fill )           &
     2775                      surf_lsm_v(l)%l_vg(:,m)        = soil_pars_f%pars_xy(1,j,i)
     2776                   IF ( soil_pars_f%pars_xy(2,j,i) /= soil_pars_f%fill )           &
     2777                      surf_lsm_v(l)%n_vg(:,m)        = soil_pars_f%pars_xy(2,j,i)
     2778                   IF ( soil_pars_f%pars_xy(3,j,i) /= soil_pars_f%fill )           &
     2779                      surf_lsm_v(l)%gamma_w_sat(:,m) = soil_pars_f%pars_xy(3,j,i)
     2780                   IF ( soil_pars_f%pars_xy(4,j,i) /= soil_pars_f%fill )           &
     2781                      surf_lsm_v(l)%m_sat(:,m)       = soil_pars_f%pars_xy(4,j,i)
     2782                   IF ( soil_pars_f%pars_xy(5,j,i) /= soil_pars_f%fill )           &
     2783                      surf_lsm_v(l)%m_fc(:,m)        = soil_pars_f%pars_xy(5,j,i)
     2784                   IF ( soil_pars_f%pars_xy(6,j,i) /= soil_pars_f%fill )           &
     2785                      surf_lsm_v(l)%m_wilt(:,m)      = soil_pars_f%pars_xy(6,j,i)
     2786                   IF ( soil_pars_f%pars_xy(7,j,i) /= soil_pars_f%fill )           &
     2787                      surf_lsm_v(l)%m_res(:,m)       = soil_pars_f%pars_xy(7,j,i)
     2788
     2789                ENDDO
     2790             ENDDO
     2791!
     2792!--       Level of detail = 2, i.e. soil parameters can be set at each soil
     2793!--       layer individually.
     2794          ELSE
     2795!
     2796!--          Horizontal surfaces
     2797             DO  m = 1, surf_lsm_h%ns
     2798                i = surf_lsm_h%i(m)
     2799                j = surf_lsm_h%j(m)
     2800
     2801                DO  k = nzb_soil, nzt_soil
     2802                   IF ( soil_pars_f%pars_xyz(0,k,j,i) /= soil_pars_f%fill )        &
     2803                      surf_lsm_h%alpha_vg(k,m)    = soil_pars_f%pars_xyz(0,k,j,i)
     2804                   IF ( soil_pars_f%pars_xyz(1,k,j,i) /= soil_pars_f%fill )        &
     2805                      surf_lsm_h%l_vg(k,m)        = soil_pars_f%pars_xyz(1,k,j,i)
     2806                   IF ( soil_pars_f%pars_xyz(2,k,j,i) /= soil_pars_f%fill )        &
     2807                      surf_lsm_h%n_vg(k,m)        = soil_pars_f%pars_xyz(2,k,j,i)
     2808                   IF ( soil_pars_f%pars_xyz(3,k,j,i) /= soil_pars_f%fill )        &
     2809                      surf_lsm_h%gamma_w_sat(k,m) = soil_pars_f%pars_xyz(3,k,j,i)
     2810                   IF ( soil_pars_f%pars_xyz(4,k,j,i) /= soil_pars_f%fill )        &
     2811                      surf_lsm_h%m_sat(k,m)       = soil_pars_f%pars_xyz(4,k,j,i)
     2812                   IF ( soil_pars_f%pars_xyz(5,k,j,i) /= soil_pars_f%fill )        &
     2813                      surf_lsm_h%m_fc(k,m)        = soil_pars_f%pars_xyz(5,k,j,i)
     2814                   IF ( soil_pars_f%pars_xyz(6,k,j,i) /= soil_pars_f%fill )        &
     2815                      surf_lsm_h%m_wilt(k,m)      = soil_pars_f%pars_xyz(6,k,j,i)
     2816                   IF ( soil_pars_f%pars_xyz(7,k,j,i) /= soil_pars_f%fill )        &
     2817                      surf_lsm_h%m_res(k,m)       = soil_pars_f%pars_xyz(7,k,j,i)
     2818                ENDDO
     2819
     2820             ENDDO
     2821!
     2822!--          Vertical surfaces
     2823             DO  l = 0, 3
     2824                DO  m = 1, surf_lsm_v(l)%ns
     2825                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
     2826                                                   surf_lsm_v(l)%building_covered(m) )
     2827                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
     2828                                                   surf_lsm_v(l)%building_covered(m) )
     2829
     2830                   DO  k = nzb_soil, nzt_soil
     2831                      IF ( soil_pars_f%pars_xyz(0,k,j,i) /= soil_pars_f%fill )        &
     2832                         surf_lsm_v(l)%alpha_vg(k,m)    = soil_pars_f%pars_xyz(0,k,j,i)
     2833                      IF ( soil_pars_f%pars_xyz(1,k,j,i) /= soil_pars_f%fill )        &
     2834                         surf_lsm_v(l)%l_vg(k,m)        = soil_pars_f%pars_xyz(1,k,j,i)
     2835                      IF ( soil_pars_f%pars_xyz(2,k,j,i) /= soil_pars_f%fill )        &
     2836                         surf_lsm_v(l)%n_vg(k,m)        = soil_pars_f%pars_xyz(2,k,j,i)
     2837                      IF ( soil_pars_f%pars_xyz(3,k,j,i) /= soil_pars_f%fill )        &
     2838                         surf_lsm_v(l)%gamma_w_sat(k,m) = soil_pars_f%pars_xyz(3,k,j,i)
     2839                      IF ( soil_pars_f%pars_xyz(4,k,j,i) /= soil_pars_f%fill )        &
     2840                         surf_lsm_v(l)%m_sat(k,m)       = soil_pars_f%pars_xyz(4,k,j,i)
     2841                      IF ( soil_pars_f%pars_xyz(5,k,j,i) /= soil_pars_f%fill )        &
     2842                         surf_lsm_v(l)%m_fc(k,m)        = soil_pars_f%pars_xyz(5,k,j,i)
     2843                      IF ( soil_pars_f%pars_xyz(6,k,j,i) /= soil_pars_f%fill )        &
     2844                         surf_lsm_v(l)%m_wilt(k,m)      = soil_pars_f%pars_xyz(6,k,j,i)
     2845                      IF ( soil_pars_f%pars_xyz(7,k,j,i) /= soil_pars_f%fill )        &
     2846                         surf_lsm_v(l)%m_res(k,m)       = soil_pars_f%pars_xyz(7,k,j,i)
     2847                   ENDDO
     2848
     2849                ENDDO
     2850             ENDDO
     2851
     2852          ENDIF
     2853       ENDIF
     2854
     2855!
     2856!--    Level 1, initialization of vegetation parameters. A horizontally
     2857!--    homogeneous distribution is assumed here.
     2858       IF ( vegetation_type /= 0 )  THEN
     2859
     2860          IF ( min_canopy_resistance == 9999999.9_wp )  THEN
     2861             min_canopy_resistance = vegetation_pars(ind_v_rc_min,vegetation_type)
     2862          ENDIF
     2863
     2864          IF ( leaf_area_index == 9999999.9_wp )  THEN
     2865             leaf_area_index = vegetation_pars(ind_v_rc_lai,vegetation_type)         
     2866          ENDIF
     2867
     2868          IF ( vegetation_coverage == 9999999.9_wp )  THEN
     2869             vegetation_coverage = vegetation_pars(ind_v_c_veg,vegetation_type)     
     2870          ENDIF
     2871
     2872          IF ( canopy_resistance_coefficient == 9999999.9_wp )  THEN
     2873              canopy_resistance_coefficient= vegetation_pars(ind_v_gd,vegetation_type)     
     2874          ENDIF
     2875
     2876          IF ( z0_vegetation == 9999999.9_wp )  THEN
     2877             z0_vegetation  = vegetation_pars(ind_v_z0,vegetation_type)
     2878          ENDIF
     2879
     2880          IF ( z0h_vegetation == 9999999.9_wp )  THEN
     2881             z0h_vegetation = vegetation_pars(ind_v_z0qh,vegetation_type)
     2882          ENDIF   
     2883         
     2884          IF ( lambda_surface_stable == 9999999.9_wp )  THEN
     2885             lambda_surface_stable = vegetation_pars(ind_v_lambda_s,vegetation_type)
     2886          ENDIF
     2887
     2888          IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
     2889             lambda_surface_unstable = vegetation_pars(ind_v_lambda_u,vegetation_type)           
     2890          ENDIF
     2891
     2892          IF ( f_shortwave_incoming == 9999999.9_wp )  THEN
     2893             f_shortwave_incoming = vegetation_pars(ind_v_f_sw_in,vegetation_type)       
     2894          ENDIF
     2895
     2896          IF ( c_surface == 9999999.9_wp )  THEN
     2897             c_surface = vegetation_pars(ind_v_c_surf,vegetation_type)       
     2898          ENDIF
     2899
     2900          IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
     2901             albedo_type = INT(vegetation_pars(ind_v_at,vegetation_type))       
     2902          ENDIF
     2903   
     2904          IF ( emissivity == 9999999.9_wp )  THEN
     2905             emissivity = vegetation_pars(ind_v_emis,vegetation_type)     
     2906          ENDIF
     2907
     2908       ENDIF
     2909!
     2910!--    Map values onto horizontal elemements
     2911       DO  m = 1, surf_lsm_h%ns
     2912          IF ( surf_lsm_h%vegetation_surface(m) )  THEN
     2913             surf_lsm_h%r_canopy_min(m)     = min_canopy_resistance
     2914             surf_lsm_h%lai(m)              = leaf_area_index
     2915             surf_lsm_h%c_veg(m)            = vegetation_coverage
     2916             surf_lsm_h%g_d(m)              = canopy_resistance_coefficient
     2917             surf_lsm_h%z0(m)               = z0_vegetation
     2918             surf_lsm_h%z0h(m)              = z0h_vegetation
     2919             surf_lsm_h%z0q(m)              = z0h_vegetation
     2920             surf_lsm_h%lambda_surface_s(m) = lambda_surface_stable
     2921             surf_lsm_h%lambda_surface_u(m) = lambda_surface_unstable
     2922             surf_lsm_h%f_sw_in(m)          = f_shortwave_incoming
     2923             surf_lsm_h%c_surface(m)        = c_surface
     2924             surf_lsm_h%albedo_type(ind_veg,m) = albedo_type
     2925             surf_lsm_h%emissivity(ind_veg,m)  = emissivity
     2926          ELSE
     2927             surf_lsm_h%lai(m)   = 0.0_wp
     2928             surf_lsm_h%c_veg(m) = 0.0_wp
     2929             surf_lsm_h%g_d(m)   = 0.0_wp
     2930          ENDIF
     2931 
     2932       ENDDO
     2933!
     2934!--    Map values onto vertical elements, even though this does not make
     2935!--    much sense.
     2936       DO  l = 0, 3
     2937          DO  m = 1, surf_lsm_v(l)%ns
     2938             IF ( surf_lsm_v(l)%vegetation_surface(m) )  THEN
     2939                surf_lsm_v(l)%r_canopy_min(m)     = min_canopy_resistance
     2940                surf_lsm_v(l)%lai(m)              = leaf_area_index
     2941                surf_lsm_v(l)%c_veg(m)            = vegetation_coverage
     2942                surf_lsm_v(l)%g_d(m)              = canopy_resistance_coefficient
     2943                surf_lsm_v(l)%z0(m)               = z0_vegetation
     2944                surf_lsm_v(l)%z0h(m)              = z0h_vegetation
     2945                surf_lsm_v(l)%z0q(m)              = z0h_vegetation
     2946                surf_lsm_v(l)%lambda_surface_s(m) = lambda_surface_stable
     2947                surf_lsm_v(l)%lambda_surface_u(m) = lambda_surface_unstable
     2948                surf_lsm_v(l)%f_sw_in(m)          = f_shortwave_incoming
     2949                surf_lsm_v(l)%c_surface(m)        = c_surface
     2950                surf_lsm_v(l)%albedo_type(ind_veg,m) = albedo_type
     2951                surf_lsm_v(l)%emissivity(ind_veg,m)  = emissivity
     2952             ELSE
     2953                surf_lsm_v(l)%lai(m)   = 0.0_wp
     2954                surf_lsm_v(l)%c_veg(m) = 0.0_wp
     2955                surf_lsm_v(l)%g_d(m)   = 0.0_wp
     2956             ENDIF
     2957          ENDDO
     2958       ENDDO
     2959
     2960!
     2961!--    Level 2, initialization of vegation parameters via vegetation_type read
     2962!--    from file. Vegetation parameters are initialized for each (y,x)-grid point
     2963!--    individually using default paramter settings according to the given
     2964!--    vegetation type.
     2965       IF ( vegetation_type_f%from_file )  THEN
     2966!
     2967!--       Horizontal surfaces
     2968          DO  m = 1, surf_lsm_h%ns
     2969             i = surf_lsm_h%i(m)
     2970             j = surf_lsm_h%j(m)
     2971             
     2972             st = vegetation_type_f%var(j,i)
     2973             IF ( st /= vegetation_type_f%fill  .AND.  st /= 0 )  THEN
     2974                surf_lsm_h%r_canopy_min(m)     = vegetation_pars(ind_v_rc_min,st)
     2975                surf_lsm_h%lai(m)              = vegetation_pars(ind_v_rc_lai,st)
     2976                surf_lsm_h%c_veg(m)            = vegetation_pars(ind_v_c_veg,st)
     2977                surf_lsm_h%g_d(m)              = vegetation_pars(ind_v_gd,st)
     2978                surf_lsm_h%z0(m)               = vegetation_pars(ind_v_z0,st)
     2979                surf_lsm_h%z0h(m)              = vegetation_pars(ind_v_z0qh,st)
     2980                surf_lsm_h%z0q(m)              = vegetation_pars(ind_v_z0qh,st)
     2981                surf_lsm_h%lambda_surface_s(m) = vegetation_pars(ind_v_lambda_s,st)
     2982                surf_lsm_h%lambda_surface_u(m) = vegetation_pars(ind_v_lambda_u,st)
     2983                surf_lsm_h%f_sw_in(m)          = vegetation_pars(ind_v_f_sw_in,st)
     2984                surf_lsm_h%c_surface(m)        = vegetation_pars(ind_v_c_surf,st)
     2985                surf_lsm_h%albedo_type(ind_veg,m) = INT( vegetation_pars(ind_v_at,st) )
     2986                surf_lsm_h%emissivity(ind_veg,m)  = vegetation_pars(ind_v_emis,st)
     2987             ENDIF
     2988          ENDDO
     2989!
     2990!--       Vertical surfaces
     2991          DO  l = 0, 3
     2992             DO  m = 1, surf_lsm_v(l)%ns
     2993                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
     2994                                                surf_lsm_v(l)%building_covered(m) )
     2995                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
     2996                                                   surf_lsm_v(l)%building_covered(m) )
     2997             
     2998                st = vegetation_type_f%var(j,i)
     2999                IF ( st /= vegetation_type_f%fill  .AND.  st /= 0 )  THEN
     3000                   surf_lsm_v(l)%r_canopy_min(m)     = vegetation_pars(ind_v_rc_min,st)
     3001                   surf_lsm_v(l)%lai(m)              = vegetation_pars(ind_v_rc_lai,st)
     3002                   surf_lsm_v(l)%c_veg(m)            = vegetation_pars(ind_v_c_veg,st)
     3003                   surf_lsm_v(l)%g_d(m)              = vegetation_pars(ind_v_gd,st)
     3004                   surf_lsm_v(l)%z0(m)               = vegetation_pars(ind_v_z0,st)
     3005                   surf_lsm_v(l)%z0h(m)              = vegetation_pars(ind_v_z0qh,st)
     3006                   surf_lsm_v(l)%z0q(m)              = vegetation_pars(ind_v_z0qh,st)
     3007                   surf_lsm_v(l)%lambda_surface_s(m) = vegetation_pars(ind_v_lambda_s,st)
     3008                   surf_lsm_v(l)%lambda_surface_u(m) = vegetation_pars(ind_v_lambda_u,st)
     3009                   surf_lsm_v(l)%f_sw_in(m)          = vegetation_pars(ind_v_f_sw_in,st)
     3010                   surf_lsm_v(l)%c_surface(m)        = vegetation_pars(ind_v_c_surf,st)
     3011                   surf_lsm_v(l)%albedo_type(ind_veg,m) = INT( vegetation_pars(ind_v_at,st) )
     3012                   surf_lsm_v(l)%emissivity(ind_veg,m)  = vegetation_pars(ind_v_emis,st)
     3013                ENDIF
     3014             ENDDO
     3015          ENDDO
     3016       ENDIF
     3017!
     3018!--    Level 3, initialization of vegation parameters at single (x,y)
     3019!--    position via vegetation_pars read from file.
     3020       IF ( vegetation_pars_f%from_file )  THEN
     3021!
     3022!--       Horizontal surfaces
     3023          DO  m = 1, surf_lsm_h%ns
     3024
     3025             i = surf_lsm_h%i(m)
     3026             j = surf_lsm_h%j(m)
     3027!
     3028!--          If surface element is not a vegetation surface and any value in
     3029!--          vegetation_pars is given, neglect this information and give an
     3030!--          informative message that this value will not be used.   
     3031             IF ( .NOT. surf_lsm_h%vegetation_surface(m)  .AND.                &
     3032                   ANY( vegetation_pars_f%pars_xy(:,j,i) /=                    &
     3033                   vegetation_pars_f%fill ) )  THEN
     3034                WRITE( message_string, * )                                     &
     3035                                 'surface element at grid point (j,i) = (',    &
     3036                                 j, i, ') is not a vegation surface, ',        &
     3037                                 'so that information given in ',              &
     3038                                 'vegetation_pars at this point is neglected.'
     3039                CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 )
     3040             ELSE
     3041
     3042                IF ( vegetation_pars_f%pars_xy(ind_v_rc_min,j,i) /=            &
     3043                     vegetation_pars_f%fill )                                  &
     3044                   surf_lsm_h%r_canopy_min(m)  =                               &
     3045                                   vegetation_pars_f%pars_xy(ind_v_rc_min,j,i)
     3046                IF ( vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i) /=            &
     3047                     vegetation_pars_f%fill )                                  &
     3048                   surf_lsm_h%lai(m)           =                               &
     3049                                   vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i)
     3050                IF ( vegetation_pars_f%pars_xy(ind_v_c_veg,j,i) /=             &
     3051                     vegetation_pars_f%fill )                                  &
     3052                   surf_lsm_h%c_veg(m)         =                               &
     3053                                   vegetation_pars_f%pars_xy(ind_v_c_veg,j,i)
     3054                IF ( vegetation_pars_f%pars_xy(ind_v_gd,j,i) /=                &
     3055                     vegetation_pars_f%fill )                                  &
     3056                   surf_lsm_h%g_d(m)           =                               &
     3057                                   vegetation_pars_f%pars_xy(ind_v_gd,j,i)
     3058                IF ( vegetation_pars_f%pars_xy(ind_v_z0,j,i) /=                &
     3059                     vegetation_pars_f%fill )                                  &
     3060                   surf_lsm_h%z0(m)            =                               &
     3061                                   vegetation_pars_f%pars_xy(ind_v_z0,j,i)
     3062                IF ( vegetation_pars_f%pars_xy(ind_v_z0qh,j,i) /=              &
     3063                     vegetation_pars_f%fill )  THEN
     3064                   surf_lsm_h%z0h(m)           =                               &
     3065                                   vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
     3066                   surf_lsm_h%z0q(m)           =                               &
     3067                                   vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
     3068                ENDIF
     3069                IF ( vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i) /=          &
     3070                     vegetation_pars_f%fill )                                  &
     3071                   surf_lsm_h%lambda_surface_s(m) =                            &
     3072                                   vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i)
     3073                IF ( vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i) /=          &
     3074                     vegetation_pars_f%fill )                                  &
     3075                   surf_lsm_h%lambda_surface_u(m) =                            &
     3076                                   vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i)
     3077                IF ( vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i) /=           &
     3078                     vegetation_pars_f%fill )                                  &
     3079                   surf_lsm_h%f_sw_in(m)          =                            &
     3080                                   vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i)
     3081                IF ( vegetation_pars_f%pars_xy(ind_v_c_surf,j,i) /=            &
     3082                     vegetation_pars_f%fill )                                  &
     3083                   surf_lsm_h%c_surface(m)        =                            &
     3084                                   vegetation_pars_f%pars_xy(ind_v_c_surf,j,i)
     3085                IF ( vegetation_pars_f%pars_xy(ind_v_at,j,i) /=                &
     3086                     vegetation_pars_f%fill )                                  &
     3087                   surf_lsm_h%albedo_type(ind_veg,m) =                         &
     3088                                   INT( vegetation_pars_f%pars_xy(ind_v_at,j,i) )
     3089                IF ( vegetation_pars_f%pars_xy(ind_v_emis,j,i) /=              &
     3090                     vegetation_pars_f%fill )                                  &
     3091                   surf_lsm_h%emissivity(ind_veg,m)  =                         &
     3092                                   vegetation_pars_f%pars_xy(ind_v_emis,j,i)
     3093             ENDIF
     3094          ENDDO
     3095!
     3096!--       Vertical surfaces
     3097          DO  l = 0, 3
     3098             DO  m = 1, surf_lsm_v(l)%ns
     3099                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
     3100                                                surf_lsm_v(l)%building_covered(m) )
     3101                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
     3102                                                surf_lsm_v(l)%building_covered(m) )
     3103!
     3104!--             If surface element is not a vegetation surface and any value in
     3105!--             vegetation_pars is given, neglect this information and give an
     3106!--             informative message that this value will not be used.   
     3107                IF ( .NOT. surf_lsm_v(l)%vegetation_surface(m)  .AND.          &
     3108                      ANY( vegetation_pars_f%pars_xy(:,j,i) /=                 &
     3109                      vegetation_pars_f%fill ) )  THEN
     3110                   WRITE( message_string, * )                                  &
     3111                                 'surface element at grid point (j,i) = (',    &
     3112                                 j, i, ') is not a vegation surface, ',        &
     3113                                 'so that information given in ',              &
     3114                                 'vegetation_pars at this point is neglected.'
     3115                   CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 )
     3116                ELSE
     3117
     3118                   IF ( vegetation_pars_f%pars_xy(ind_v_rc_min,j,i) /=         &
     3119                        vegetation_pars_f%fill )                               &
     3120                      surf_lsm_v(l)%r_canopy_min(m)  =                         &
     3121                                   vegetation_pars_f%pars_xy(ind_v_rc_min,j,i)
     3122                   IF ( vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i) /=         &
     3123                        vegetation_pars_f%fill )                               &
     3124                      surf_lsm_v(l)%lai(m)           =                         &
     3125                                   vegetation_pars_f%pars_xy(ind_v_rc_lai,j,i)
     3126                   IF ( vegetation_pars_f%pars_xy(ind_v_c_veg,j,i) /=          &
     3127                        vegetation_pars_f%fill )                               &
     3128                      surf_lsm_v(l)%c_veg(m)         =                         &
     3129                                   vegetation_pars_f%pars_xy(ind_v_c_veg,j,i)
     3130                   IF ( vegetation_pars_f%pars_xy(ind_v_gd,j,i) /=             &
     3131                        vegetation_pars_f%fill )                               &
     3132                     surf_lsm_v(l)%g_d(m)            =                         &
     3133                                   vegetation_pars_f%pars_xy(ind_v_gd,j,i)
     3134                   IF ( vegetation_pars_f%pars_xy(ind_v_z0,j,i) /=             &
     3135                        vegetation_pars_f%fill )                               &
     3136                      surf_lsm_v(l)%z0(m)            =                         &
     3137                                   vegetation_pars_f%pars_xy(ind_v_z0,j,i)
     3138                   IF ( vegetation_pars_f%pars_xy(ind_v_z0qh,j,i) /=           &
     3139                        vegetation_pars_f%fill )  THEN
     3140                      surf_lsm_v(l)%z0h(m)           =                         &
     3141                                   vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
     3142                      surf_lsm_v(l)%z0q(m)           =                         &
     3143                                   vegetation_pars_f%pars_xy(ind_v_z0qh,j,i)
     3144                   ENDIF
     3145                   IF ( vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i) /=       &
     3146                        vegetation_pars_f%fill )                               &
     3147                      surf_lsm_v(l)%lambda_surface_s(m)  =                     &
     3148                                   vegetation_pars_f%pars_xy(ind_v_lambda_s,j,i)
     3149                   IF ( vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i) /=       &
     3150                        vegetation_pars_f%fill )                               &
     3151                      surf_lsm_v(l)%lambda_surface_u(m)  =                     &
     3152                                   vegetation_pars_f%pars_xy(ind_v_lambda_u,j,i)
     3153                   IF ( vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i) /=        &
     3154                        vegetation_pars_f%fill )                               &
     3155                      surf_lsm_v(l)%f_sw_in(m)           =                     &
     3156                                   vegetation_pars_f%pars_xy(ind_v_f_sw_in,j,i)
     3157                   IF ( vegetation_pars_f%pars_xy(ind_v_c_surf,j,i) /=         &
     3158                        vegetation_pars_f%fill )                               &
     3159                      surf_lsm_v(l)%c_surface(m)         =                     &
     3160                                   vegetation_pars_f%pars_xy(ind_v_c_surf,j,i)
     3161                   IF ( vegetation_pars_f%pars_xy(ind_v_at,j,i) /=             &
     3162                        vegetation_pars_f%fill )                               &
     3163                      surf_lsm_v(l)%albedo_type(ind_veg,m) =                   &
     3164                                   INT( vegetation_pars_f%pars_xy(ind_v_at,j,i) )
     3165                   IF ( vegetation_pars_f%pars_xy(ind_v_emis,j,i) /=           &
     3166                        vegetation_pars_f%fill )                               &
     3167                      surf_lsm_v(l)%emissivity(ind_veg,m)  =                   &
     3168                                   vegetation_pars_f%pars_xy(ind_v_emis,j,i)
     3169                ENDIF
     3170
     3171             ENDDO
     3172          ENDDO
     3173       ENDIF
     3174
     3175!
     3176!--    Level 1, initialization of water parameters. A horizontally
     3177!--    homogeneous distribution is assumed here.
     3178       IF ( water_type /= 0 )  THEN
     3179
     3180          IF ( water_temperature == 9999999.9_wp )  THEN
     3181             water_temperature = water_pars(ind_w_temp,water_type)       
     3182          ENDIF
     3183
     3184          IF ( z0_water == 9999999.9_wp )  THEN
     3185             z0_water = water_pars(ind_w_z0,water_type)       
     3186          ENDIF       
     3187
     3188          IF ( z0h_water == 9999999.9_wp )  THEN
     3189             z0h_water = water_pars(ind_w_z0h,water_type)       
     3190          ENDIF 
     3191
     3192          IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
     3193             albedo_type = INT(water_pars(ind_w_at,water_type))       
     3194          ENDIF
     3195   
     3196          IF ( emissivity == 9999999.9_wp )  THEN
     3197             emissivity = water_pars(ind_w_emis,water_type)       
     3198          ENDIF
     3199
     3200       ENDIF
     3201!
     3202!--    Map values onto horizontal elemements
     3203       DO  m = 1, surf_lsm_h%ns
     3204          IF ( surf_lsm_h%water_surface(m) )  THEN
     3205             IF ( TRIM( initializing_actions ) /= 'read_restart_data' )        &
     3206                t_soil_h%var_2d(:,m)           = water_temperature
     3207             surf_lsm_h%z0(m)               = z0_water
     3208             surf_lsm_h%z0h(m)              = z0h_water
     3209             surf_lsm_h%z0q(m)              = z0h_water
     3210             surf_lsm_h%lambda_surface_s(m) = 1.0E10_wp
     3211             surf_lsm_h%lambda_surface_u(m) = 1.0E10_wp               
     3212             surf_lsm_h%c_surface(m)        = 0.0_wp
     3213             surf_lsm_h%albedo_type(ind_wat,m) = albedo_type
     3214             surf_lsm_h%emissivity(ind_wat,m)  = emissivity
     3215          ENDIF
     3216       ENDDO
     3217!
     3218!--    Map values onto vertical elements, even though this does not make
     3219!--    much sense.
     3220       DO  l = 0, 3
     3221          DO  m = 1, surf_lsm_v(l)%ns
     3222             IF ( surf_lsm_v(l)%water_surface(m) )  THEN
     3223                IF ( TRIM( initializing_actions ) /= 'read_restart_data' )     &
     3224                   t_soil_v(l)%var_2d(:,m)           = water_temperature
     3225                surf_lsm_v(l)%z0(m)               = z0_water
     3226                surf_lsm_v(l)%z0h(m)              = z0h_water
     3227                surf_lsm_v(l)%z0q(m)              = z0h_water
     3228                surf_lsm_v(l)%lambda_surface_s(m) = 1.0E10_wp
     3229                surf_lsm_v(l)%lambda_surface_u(m) = 1.0E10_wp               
     3230                surf_lsm_v(l)%c_surface(m)        = 0.0_wp
     3231                surf_lsm_v(l)%albedo_type(ind_wat,m) = albedo_type
     3232                surf_lsm_v(l)%emissivity(ind_wat,m)  = emissivity
     3233             ENDIF
     3234          ENDDO
     3235       ENDDO
     3236!
     3237!
     3238!--    Level 2, initialization of water parameters via water_type read
     3239!--    from file. Water surfaces are initialized for each (y,x)-grid point
     3240!--    individually using default paramter settings according to the given
     3241!--    water type.
     3242!--    Note, parameter 3/4 of water_pars are albedo and emissivity,
     3243!--    whereas paramter 3/4 of water_pars_f are heat conductivities!
     3244       IF ( water_type_f%from_file )  THEN
     3245!
     3246!--       Horizontal surfaces
     3247          DO  m = 1, surf_lsm_h%ns
     3248             i = surf_lsm_h%i(m)
     3249             j = surf_lsm_h%j(m)
     3250             
     3251             st = water_type_f%var(j,i)
     3252             IF ( st /= water_type_f%fill  .AND.  st /= 0 )  THEN
     3253                IF ( TRIM( initializing_actions ) /= 'read_restart_data' )     &
     3254                   t_soil_h%var_2d(:,m) = water_pars(ind_w_temp,st)
     3255                surf_lsm_h%z0(m)     = water_pars(ind_w_z0,st)
     3256                surf_lsm_h%z0h(m)    = water_pars(ind_w_z0h,st)
     3257                surf_lsm_h%z0q(m)    = water_pars(ind_w_z0h,st)
     3258                surf_lsm_h%lambda_surface_s(m) = water_pars(ind_w_lambda_s,st)
     3259                surf_lsm_h%lambda_surface_u(m) = water_pars(ind_w_lambda_u,st)             
     3260                surf_lsm_h%c_surface(m)        = 0.0_wp
     3261                surf_lsm_h%albedo_type(ind_wat,m) = INT( water_pars(ind_w_at,st) )
     3262                surf_lsm_h%emissivity(ind_wat,m)  = water_pars(ind_w_emis,st)
     3263             ENDIF
     3264          ENDDO
     3265!
     3266!--       Vertical surfaces
     3267          DO  l = 0, 3
     3268             DO  m = 1, surf_lsm_v(l)%ns
     3269                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
     3270                                                surf_lsm_v(l)%building_covered(m) )
     3271                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
     3272                                                surf_lsm_v(l)%building_covered(m) )
     3273             
     3274                st = water_type_f%var(j,i)
     3275                IF ( st /= water_type_f%fill  .AND.  st /= 0 )  THEN
     3276                   IF ( TRIM( initializing_actions ) /= 'read_restart_data' )  &
     3277                      t_soil_v(l)%var_2d(:,m) = water_pars(ind_w_temp,st)
     3278                   surf_lsm_v(l)%z0(m)     = water_pars(ind_w_z0,st)
     3279                   surf_lsm_v(l)%z0h(m)    = water_pars(ind_w_z0h,st)
     3280                   surf_lsm_v(l)%z0q(m)    = water_pars(ind_w_z0h,st)
     3281                   surf_lsm_v(l)%lambda_surface_s(m) =                         &
     3282                                                   water_pars(ind_w_lambda_s,st)
     3283                   surf_lsm_v(l)%lambda_surface_u(m) =                         &
     3284                                                   water_pars(ind_w_lambda_u,st)           
     3285                   surf_lsm_v(l)%c_surface(m)     = 0.0_wp
     3286                   surf_lsm_v(l)%albedo_type(ind_wat,m) =                      &
     3287                                                  INT( water_pars(ind_w_at,st) )
     3288                   surf_lsm_v(l)%emissivity(ind_wat,m)  =                      &
     3289                                                  water_pars(ind_w_emis,st)
     3290                ENDIF
     3291             ENDDO
     3292          ENDDO
     3293       ENDIF     
     3294
     3295!
     3296!--    Level 3, initialization of water parameters at single (x,y)
     3297!--    position via water_pars read from file.
     3298       IF ( water_pars_f%from_file )  THEN
     3299!
     3300!--       Horizontal surfaces
     3301          DO  m = 1, surf_lsm_h%ns
     3302             i = surf_lsm_h%i(m)
     3303             j = surf_lsm_h%j(m)
     3304!
     3305!--          If surface element is not a water surface and any value in
     3306!--          water_pars is given, neglect this information and give an
     3307!--          informative message that this value will not be used.   
     3308             IF ( .NOT. surf_lsm_h%water_surface(m)  .AND.                     &
     3309                   ANY( water_pars_f%pars_xy(:,j,i) /= water_pars_f%fill ) )  THEN
     3310                WRITE( message_string, * )                                     &
     3311                              'surface element at grid point (j,i) = (',       &
     3312                              j, i, ') is not a water surface, ',              &
     3313                              'so that information given in ',                 &
     3314                              'water_pars at this point is neglected.'
     3315                CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 )
     3316             ELSE
     3317                IF ( water_pars_f%pars_xy(ind_w_temp,j,i) /=                   &
     3318                     water_pars_f%fill  .AND.                                  &
     3319                     TRIM( initializing_actions ) /= 'read_restart_data' )     &
     3320                      t_soil_h%var_2d(:,m) = water_pars_f%pars_xy(ind_w_temp,j,i)
     3321
     3322                IF ( water_pars_f%pars_xy(ind_w_z0,j,i) /= water_pars_f%fill ) &
     3323                   surf_lsm_h%z0(m)     = water_pars_f%pars_xy(ind_w_z0,j,i)
     3324
     3325                IF ( water_pars_f%pars_xy(ind_w_z0h,j,i) /= water_pars_f%fill )&
     3326                THEN
     3327                   surf_lsm_h%z0h(m)    = water_pars_f%pars_xy(ind_w_z0h,j,i)
     3328                   surf_lsm_h%z0q(m)    = water_pars_f%pars_xy(ind_w_z0h,j,i)
     3329                ENDIF
     3330
     3331                IF ( water_pars_f%pars_xy(ind_w_lambda_s,j,i) /=               &
     3332                     water_pars_f%fill )                                       &
     3333                   surf_lsm_h%lambda_surface_s(m) =                            &
     3334                                        water_pars_f%pars_xy(ind_w_lambda_s,j,i)
     3335
     3336                IF ( water_pars_f%pars_xy(ind_w_lambda_u,j,i) /=               &
     3337                      water_pars_f%fill )                                      &
     3338                   surf_lsm_h%lambda_surface_u(m) =                            &
     3339                                        water_pars_f%pars_xy(ind_w_lambda_u,j,i)     
     3340       
     3341                IF ( water_pars_f%pars_xy(ind_w_at,j,i) /=                     &
     3342                     water_pars_f%fill )                                       &
     3343                   surf_lsm_h%albedo_type(ind_wat,m) =                         &
     3344                                       INT( water_pars_f%pars_xy(ind_w_at,j,i) )
     3345
     3346                IF ( water_pars_f%pars_xy(ind_w_emis,j,i) /=                   &
     3347                     water_pars_f%fill )                                       &
     3348                   surf_lsm_h%emissivity(ind_wat,m) =                          &
     3349                   water_pars_f%pars_xy(ind_w_emis,j,i) 
     3350             ENDIF
     3351          ENDDO
     3352!
     3353!--       Vertical surfaces
     3354          DO  l = 0, 3
     3355             DO  m = 1, surf_lsm_v(l)%ns
     3356                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
     3357                                                surf_lsm_v(l)%building_covered(m) )
     3358                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
     3359                                                surf_lsm_v(l)%building_covered(m) )
     3360!
     3361!--             If surface element is not a water surface and any value in
     3362!--             water_pars is given, neglect this information and give an
     3363!--             informative message that this value will not be used.   
     3364                IF ( .NOT. surf_lsm_v(l)%water_surface(m)  .AND.               &
     3365                      ANY( water_pars_f%pars_xy(:,j,i) /=                      &
     3366                      water_pars_f%fill ) )  THEN
     3367                   WRITE( message_string, * )                                  &
     3368                              'surface element at grid point (j,i) = (',       &
     3369                              j, i, ') is not a water surface, ',              &
     3370                              'so that information given in ',                 &
     3371                              'water_pars at this point is neglected.'
     3372                   CALL message( 'land_surface_model_mod', 'PA0999',           &
     3373                                  0, 0, 0, 6, 0 )
     3374                ELSE
     3375
     3376                   IF ( water_pars_f%pars_xy(ind_w_temp,j,i) /=                &
     3377                     water_pars_f%fill  .AND.                                  &
     3378                     TRIM( initializing_actions ) /= 'read_restart_data' )     &
     3379                      t_soil_v(l)%var_2d(:,m) = water_pars_f%pars_xy(ind_w_temp,j,i)
     3380
     3381                   IF ( water_pars_f%pars_xy(ind_w_z0,j,i) /=                  &
     3382                        water_pars_f%fill )                                    &
     3383                      surf_lsm_v(l)%z0(m)   = water_pars_f%pars_xy(ind_w_z0,j,i)
     3384
     3385                   IF ( water_pars_f%pars_xy(ind_w_z0h,j,i) /=                 &
     3386                       water_pars_f%fill )  THEN
     3387                      surf_lsm_v(l)%z0h(m)  = water_pars_f%pars_xy(ind_w_z0h,j,i)
     3388                      surf_lsm_v(l)%z0q(m)  = water_pars_f%pars_xy(ind_w_z0h,j,i)
     3389                   ENDIF
     3390
     3391                   IF ( water_pars_f%pars_xy(ind_w_lambda_s,j,i) /=            &
     3392                        water_pars_f%fill )                                    &
     3393                      surf_lsm_v(l)%lambda_surface_s(m) =                      &
     3394                                      water_pars_f%pars_xy(ind_w_lambda_s,j,i)
     3395
     3396                   IF ( water_pars_f%pars_xy(ind_w_lambda_u,j,i) /=            &
     3397                        water_pars_f%fill )                                    &
     3398                      surf_lsm_v(l)%lambda_surface_u(m) =                      &
     3399                                      water_pars_f%pars_xy(ind_w_lambda_u,j,i)   
     3400 
     3401                   IF ( water_pars_f%pars_xy(ind_w_at,j,i) /=                  &
     3402                        water_pars_f%fill )                                    &
     3403                      surf_lsm_v(l)%albedo_type(ind_wat,m)    =                &
     3404                                      INT( water_pars_f%pars_xy(ind_w_at,j,i) )
     3405
     3406                   IF ( water_pars_f%pars_xy(ind_w_emis,j,i) /=                &
     3407                        water_pars_f%fill )                                    &
     3408                      surf_lsm_v(l)%emissivity(ind_wat,m)     =                &
     3409                                      water_pars_f%pars_xy(ind_w_emis,j,i) 
     3410                ENDIF
     3411             ENDDO
     3412          ENDDO
     3413
     3414       ENDIF
     3415!
     3416!--    Initialize pavement-type surfaces, level 1
     3417       IF ( pavement_type /= 0 )  THEN 
     3418
     3419!
     3420!--       When a pavement_type is used, overwrite a possible setting of
     3421!--       the pavement depth as it is already defined by the pavement type
     3422          pavement_depth_level = 0
     3423
     3424          IF ( z0_pavement == 9999999.9_wp )  THEN
     3425             z0_pavement  = pavement_pars(ind_p_z0,pavement_type)
     3426          ENDIF
     3427
     3428          IF ( z0h_pavement == 9999999.9_wp )  THEN
     3429             z0h_pavement = pavement_pars(ind_p_z0h,pavement_type)
     3430          ENDIF
     3431
     3432          IF ( pavement_heat_conduct == 9999999.9_wp )  THEN
     3433             pavement_heat_conduct = pavement_subsurface_pars_1(0,pavement_type)
     3434          ENDIF
     3435
     3436          IF ( pavement_heat_capacity == 9999999.9_wp )  THEN
     3437             pavement_heat_capacity = pavement_subsurface_pars_2(0,pavement_type)
     3438          ENDIF   
     3439   
     3440          IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
     3441             albedo_type = INT(pavement_pars(ind_p_at,pavement_type))       
     3442          ENDIF
     3443   
     3444          IF ( emissivity == 9999999.9_wp )  THEN
     3445             emissivity = pavement_pars(ind_p_emis,pavement_type)       
     3446          ENDIF
     3447
     3448!
     3449!--       If the depth level of the pavement is not set, determine it from
     3450!--       lookup table.
     3451          IF ( pavement_depth_level == 0 )  THEN
     3452             DO  k = nzb_soil, nzt_soil 
     3453                IF ( pavement_subsurface_pars_1(k,pavement_type) == 9999999.9_wp &
     3454                .OR. pavement_subsurface_pars_2(k,pavement_type) == 9999999.9_wp)&
     3455                THEN
     3456                   nzt_pavement = k-1
     3457                   EXIT
     3458                ENDIF
     3459             ENDDO
     3460          ELSE
     3461             nzt_pavement = pavement_depth_level
     3462          ENDIF
     3463
     3464       ENDIF
     3465!
     3466!--    Level 1 initialization of pavement type surfaces. Horizontally
     3467!--    homogeneous characteristics are assumed
     3468       surf_lsm_h%nzt_pavement = pavement_depth_level
     3469       DO  m = 1, surf_lsm_h%ns
     3470          IF ( surf_lsm_h%pavement_surface(m) )  THEN
     3471             surf_lsm_h%nzt_pavement(m)        = nzt_pavement
     3472             surf_lsm_h%z0(m)                  = z0_pavement
     3473             surf_lsm_h%z0h(m)                 = z0h_pavement
     3474             surf_lsm_h%z0q(m)                 = z0h_pavement
     3475             surf_lsm_h%lambda_surface_s(m)    = pavement_heat_conduct         &
     3476                                                  * ddz_soil(nzb_soil)         &
     3477                                                  * 2.0_wp   
     3478             surf_lsm_h%lambda_surface_u(m)    = pavement_heat_conduct         &
     3479                                                  * ddz_soil(nzb_soil)         &
     3480                                                  * 2.0_wp           
     3481             surf_lsm_h%c_surface(m)           = pavement_heat_capacity        &
     3482                                                        * dz_soil(nzb_soil)    &
     3483                                                        * 0.25_wp                                   
     3484
     3485             surf_lsm_h%albedo_type(ind_pav,m) = albedo_type
     3486             surf_lsm_h%emissivity(ind_pav,m)  = emissivity     
     3487     
     3488             IF ( pavement_type /= 0 )  THEN
     3489                DO  k = nzb_soil, surf_lsm_h%nzt_pavement(m)
     3490                   surf_lsm_h%lambda_h_def(k,m)    =                           &
     3491                                     pavement_subsurface_pars_1(k,pavement_type)                       
     3492                   surf_lsm_h%rho_c_total_def(k,m) =                           &
     3493                                     pavement_subsurface_pars_2(k,pavement_type)
     3494                ENDDO
     3495             ELSE
     3496                surf_lsm_v(l)%lambda_h_def(:,m)     = pavement_heat_conduct
     3497                surf_lsm_v(l)%rho_c_total_def(:,m)  = pavement_heat_capacity
     3498             ENDIF       
     3499          ENDIF
     3500       ENDDO                               
     3501
     3502       DO  l = 0, 3
     3503          surf_lsm_v(l)%nzt_pavement = pavement_depth_level
     3504          DO  m = 1, surf_lsm_v(l)%ns
     3505             IF ( surf_lsm_v(l)%pavement_surface(m) )  THEN
     3506                surf_lsm_v(l)%nzt_pavement(m)        = nzt_pavement
     3507                surf_lsm_v(l)%z0(m)                  = z0_pavement
     3508                surf_lsm_v(l)%z0h(m)                 = z0h_pavement
     3509                surf_lsm_v(l)%z0q(m)                 = z0h_pavement
     3510                surf_lsm_v(l)%lambda_surface_s(m)    = pavement_heat_conduct   &
     3511                                                  * ddz_soil(nzb_soil)         &
     3512                                                  * 2.0_wp   
     3513                surf_lsm_v(l)%lambda_surface_u(m)    = pavement_heat_conduct   &
     3514                                                  * ddz_soil(nzb_soil)         &
     3515                                                  * 2.0_wp           
     3516                surf_lsm_v(l)%c_surface(m)           = pavement_heat_capacity  &
     3517                                                        * dz_soil(nzb_soil)    &
     3518                                                        * 0.25_wp                                     
     3519
     3520                surf_lsm_v(l)%albedo_type(ind_pav,m) = albedo_type
     3521                surf_lsm_v(l)%emissivity(ind_pav,m)  = emissivity
     3522
     3523                IF ( pavement_type /= 0 )  THEN
     3524                   DO  k = nzb_soil, surf_lsm_v(l)%nzt_pavement(m)
     3525                      surf_lsm_v(l)%lambda_h_def(k,m)    =                     &
     3526                                     pavement_subsurface_pars_1(k,pavement_type)                       
     3527                      surf_lsm_v(l)%rho_c_total_def(k,m) =                     &
     3528                                     pavement_subsurface_pars_2(k,pavement_type)
     3529                   ENDDO
     3530                ELSE
     3531                   surf_lsm_v(l)%lambda_h_def(:,m)     = pavement_heat_conduct
     3532                   surf_lsm_v(l)%rho_c_total_def(:,m)  = pavement_heat_capacity
     3533                ENDIF     
     3534             ENDIF
     3535          ENDDO
     3536       ENDDO
     3537!
     3538!--    Level 2 initialization of pavement type surfaces via pavement_type read
     3539!--    from file. Pavement surfaces are initialized for each (y,x)-grid point
     3540!--    individually.
     3541       IF ( pavement_type_f%from_file )  THEN
     3542!
     3543!--       Horizontal surfaces
     3544          DO  m = 1, surf_lsm_h%ns
     3545             i = surf_lsm_h%i(m)
     3546             j = surf_lsm_h%j(m)
     3547             
     3548             st = pavement_type_f%var(j,i)
     3549             IF ( st /= pavement_type_f%fill  .AND.  st /= 0 )  THEN
     3550!
     3551!--             Determine deepmost index of pavement layer
     3552                DO  k = nzb_soil, nzt_soil 
     3553                   IF ( pavement_subsurface_pars_1(k,st) == 9999999.9_wp       &
     3554                   .OR. pavement_subsurface_pars_2(k,st) == 9999999.9_wp)      &
     3555                   THEN
     3556                      surf_lsm_h%nzt_pavement(m) = k-1
     3557                      EXIT
     3558                   ENDIF
     3559                ENDDO
     3560
     3561                surf_lsm_h%z0(m)                = pavement_pars(ind_p_z0,st)
     3562                surf_lsm_h%z0h(m)               = pavement_pars(ind_p_z0h,st)
     3563                surf_lsm_h%z0q(m)               = pavement_pars(ind_p_z0h,st)
     3564
     3565                surf_lsm_h%lambda_surface_s(m)  =                              &
     3566                                              pavement_subsurface_pars_1(0,st) &
     3567                                                  * ddz_soil(nzb_soil)         &
     3568                                                  * 2.0_wp   
     3569                surf_lsm_h%lambda_surface_u(m)  =                              &
     3570                                              pavement_subsurface_pars_1(0,st) &
     3571                                                  * ddz_soil(nzb_soil)         &
     3572                                                  * 2.0_wp       
     3573                surf_lsm_h%c_surface(m)         =                              &
     3574                                               pavement_subsurface_pars_2(0,st)&
     3575                                                        * dz_soil(nzb_soil)    &
     3576                                                        * 0.25_wp                               
     3577                surf_lsm_h%albedo_type(ind_pav,m) = INT( pavement_pars(ind_p_at,st) )
     3578                surf_lsm_h%emissivity(ind_pav,m)  = pavement_pars(ind_p_emis,st) 
     3579
     3580                DO  k = nzb_soil, surf_lsm_h%nzt_pavement(m)
     3581                   surf_lsm_h%lambda_h_def(k,m)    =                           &
     3582                                     pavement_subsurface_pars_1(k,pavement_type)                       
     3583                   surf_lsm_h%rho_c_total_def(k,m) =                           &
     3584                                     pavement_subsurface_pars_2(k,pavement_type)
     3585                ENDDO   
     3586             ENDIF
     3587          ENDDO
     3588!
     3589!--       Vertical surfaces
     3590          DO  l = 0, 3
     3591             DO  m = 1, surf_lsm_v(l)%ns
     3592                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
     3593                                                surf_lsm_v(l)%building_covered(m) )
     3594                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
     3595                                                surf_lsm_v(l)%building_covered(m) )
     3596             
     3597                st = pavement_type_f%var(j,i)
     3598                IF ( st /= pavement_type_f%fill  .AND.  st /= 0 )  THEN
     3599!
     3600!--                Determine deepmost index of pavement layer
     3601                   DO  k = nzb_soil, nzt_soil 
     3602                      IF ( pavement_subsurface_pars_1(k,st) == 9999999.9_wp    &
     3603                      .OR. pavement_subsurface_pars_2(k,st) == 9999999.9_wp)   &
     3604                      THEN
     3605                         surf_lsm_v(l)%nzt_pavement(m) = k-1
     3606                         EXIT
     3607                      ENDIF
     3608                   ENDDO
     3609
     3610                   surf_lsm_v(l)%z0(m)  = pavement_pars(ind_p_z0,st)
     3611                   surf_lsm_v(l)%z0h(m) = pavement_pars(ind_p_z0h,st)
     3612                   surf_lsm_v(l)%z0q(m) = pavement_pars(ind_p_z0h,st)
     3613
     3614                   surf_lsm_v(l)%lambda_surface_s(m)  =                        &
     3615                                              pavement_subsurface_pars_1(0,st) &
     3616                                                  * ddz_soil(nzb_soil)         &
     3617                                                  * 2.0_wp   
     3618                   surf_lsm_v(l)%lambda_surface_u(m)  =                        &
     3619                                              pavement_subsurface_pars_1(0,st) &
     3620                                                  * ddz_soil(nzb_soil)         &
     3621                                                  * 2.0_wp     
     3622
     3623                   surf_lsm_v(l)%c_surface(m)    =                             &
     3624                                           pavement_subsurface_pars_2(0,st)    &
     3625                                                        * dz_soil(nzb_soil)    &
     3626                                                        * 0.25_wp                                   
     3627                   surf_lsm_v(l)%albedo_type(ind_pav,m) =                      &
     3628                                              INT( pavement_pars(ind_p_at,st) )
     3629                   surf_lsm_v(l)%emissivity(ind_pav,m)  =                      &
     3630                                              pavement_pars(ind_p_emis,st)   
     3631
     3632                   DO  k = nzb_soil, surf_lsm_h%nzt_pavement(m)
     3633                      surf_lsm_v(l)%lambda_h_def(k,m)    =                    &
     3634                                    pavement_subsurface_pars_1(k,pavement_type)                       
     3635                      surf_lsm_v(l)%rho_c_total_def(k,m) =                    &
     3636                                    pavement_subsurface_pars_2(k,pavement_type)
     3637                   ENDDO   
     3638                ENDIF
     3639             ENDDO
     3640          ENDDO
     3641       ENDIF
     3642!
     3643!--    Level 3, initialization of pavement parameters at single (x,y)
     3644!--    position via pavement_pars read from file.
     3645       IF ( pavement_pars_f%from_file )  THEN
     3646!
     3647!--       Horizontal surfaces
     3648          DO  m = 1, surf_lsm_h%ns
     3649             i = surf_lsm_h%i(m)
     3650             j = surf_lsm_h%j(m)
     3651!
     3652!--          If surface element is not a pavement surface and any value in
     3653!--          pavement_pars is given, neglect this information and give an
     3654!--          informative message that this value will not be used.   
     3655             IF ( .NOT. surf_lsm_h%pavement_surface(m)  .AND.                  &
     3656                   ANY( pavement_pars_f%pars_xy(:,j,i) /=                      &
     3657                   pavement_pars_f%fill ) )  THEN
     3658                WRITE( message_string, * )                                     &
     3659                              'surface element at grid point (j,i) = (',       &
     3660                              j, i, ') is not a pavement surface, ',           &
     3661                              'so that information given in ',                 &
     3662                              'pavement_pars at this point is neglected.'
     3663                CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 )
     3664             ELSE
     3665                IF ( pavement_pars_f%pars_xy(ind_p_z0,j,i) /=                  &
     3666                     pavement_pars_f%fill )                                    &
     3667                   surf_lsm_h%z0(m)  = pavement_pars_f%pars_xy(ind_p_z0,j,i)
     3668                IF ( pavement_pars_f%pars_xy(ind_p_z0h,j,i) /=                 &
     3669                     pavement_pars_f%fill )  THEN
     3670                   surf_lsm_h%z0h(m) = pavement_pars_f%pars_xy(ind_p_z0h,j,i)
     3671                   surf_lsm_h%z0q(m) = pavement_pars_f%pars_xy(ind_p_z0h,j,i)
     3672                ENDIF
     3673                IF ( pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i) &
     3674                     /= pavement_subsurface_pars_f%fill )  THEN
     3675                   surf_lsm_h%lambda_surface_s(m)  =                           &
     3676                      pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)&
     3677                                                  * ddz_soil(nzb_soil)         &
     3678                                                  * 2.0_wp
     3679                   surf_lsm_h%lambda_surface_u(m)  =                           &
     3680                      pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)&
     3681                                                  * ddz_soil(nzb_soil)         &
     3682                                                  * 2.0_wp   
     3683                ENDIF
     3684                IF ( pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,0,j,i) /= &
     3685                     pavement_subsurface_pars_f%fill )  THEN
     3686                   surf_lsm_h%c_surface(m)     =                               &
     3687                      pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,0,j,i)   &
     3688                                                  * dz_soil(nzb_soil)          &
     3689                                                  * 0.25_wp                                   
     3690                ENDIF
     3691                IF ( pavement_pars_f%pars_xy(ind_p_at,j,i) /=                  &
     3692                     pavement_pars_f%fill )                                    &
     3693                   surf_lsm_h%albedo_type(ind_pav,m) =                         &
     3694                                              INT( pavement_pars(ind_p_at,st) )
     3695                IF ( pavement_pars_f%pars_xy(ind_p_emis,j,i) /=                &
     3696                     pavement_pars_f%fill )                                    &
     3697                   surf_lsm_h%emissivity(ind_pav,m)  =                         &
     3698                                              pavement_pars(ind_p_emis,st)
     3699             ENDIF
     3700
     3701          ENDDO
     3702!
     3703!--       Vertical surfaces
     3704          DO  l = 0, 3
     3705             DO  m = 1, surf_lsm_v(l)%ns
     3706                i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,         &
     3707                                                surf_lsm_v(l)%building_covered(m) )
     3708                j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,         &
     3709                                                surf_lsm_v(l)%building_covered(m) )
     3710!
     3711!--             If surface element is not a pavement surface and any value in
     3712!--             pavement_pars is given, neglect this information and give an
     3713!--             informative message that this value will not be used.   
     3714                IF ( .NOT. surf_lsm_v(l)%pavement_surface(m)  .AND.            &
     3715                      ANY( pavement_pars_f%pars_xy(:,j,i) /=                   &
     3716                      pavement_pars_f%fill ) )  THEN
     3717                   WRITE( message_string, * )                                  &
     3718                                 'surface element at grid point (j,i) = (',    &
     3719                                 j, i, ') is not a pavement surface, ',        &
     3720                                 'so that information given in ',              &
     3721                                 'pavement_pars at this point is neglected.'
     3722                   CALL message( 'land_surface_model_mod', 'PA0999', 0, 0, 0, 6, 0 )
     3723                ELSE
     3724
     3725                   IF ( pavement_pars_f%pars_xy(ind_p_z0,j,i) /=               &
     3726                        pavement_pars_f%fill )                                 &
     3727                      surf_lsm_v(l)%z0(m) = pavement_pars_f%pars_xy(ind_p_z0,j,i)
     3728                   IF ( pavement_pars_f%pars_xy(ind_p_z0h,j,i) /=              &
     3729                        pavement_pars_f%fill )  THEN
     3730                      surf_lsm_v(l)%z0h(m) = pavement_pars_f%pars_xy(ind_p_z0h,j,i)
     3731                      surf_lsm_v(l)%z0q(m) = pavement_pars_f%pars_xy(ind_p_z0h,j,i)
     3732                   ENDIF
     3733                   IF ( pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)&
     3734                        /= pavement_subsurface_pars_f%fill )  THEN
     3735                      surf_lsm_v(l)%lambda_surface_s(m) =                      &
     3736                      pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)&
     3737                                                  * ddz_soil(nzb_soil)         &
     3738                                                  * 2.0_wp
     3739                      surf_lsm_v(l)%lambda_surface_u(m) =                      &
     3740                      pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,0,j,i)&
     3741                                                  * ddz_soil(nzb_soil)         &
     3742                                                  * 2.0_wp   
     3743                   ENDIF
     3744                   IF ( pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,0,j,i) &
     3745                        /= pavement_subsurface_pars_f%fill )  THEN
     3746                      surf_lsm_v(l)%c_surface(m)    =                          &
     3747                         pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,0,j,i)&
     3748                                                  * dz_soil(nzb_soil)          &
     3749                                                  * 0.25_wp                                 
     3750                   ENDIF
     3751                   IF ( pavement_pars_f%pars_xy(ind_p_at,j,i) /=               &
     3752                        pavement_pars_f%fill )                                 &
     3753                      surf_lsm_v(l)%albedo_type(ind_pav,m) =                   &
     3754                                            INT( pavement_pars(ind_p_at,st) )
     3755
     3756                   IF ( pavement_pars_f%pars_xy(ind_p_emis,j,i) /=             &
     3757                        pavement_pars_f%fill )                                 &
     3758                      surf_lsm_v(l)%emissivity(ind_pav,m)  =                   &
     3759                                            pavement_pars(ind_p_emis,st) 
     3760                ENDIF
     3761             ENDDO
     3762          ENDDO
     3763       ENDIF
     3764!
     3765!--    Moreover, for grid points which are flagged with pavement-type 0 or whre
     3766!--    pavement_subsurface_pars_f is provided, soil heat conductivity and
     3767!--    capacity are initialized with parameters given in       
     3768!--    pavement_subsurface_pars read from file.
     3769       IF ( pavement_subsurface_pars_f%from_file )  THEN
     3770!
     3771!--       Set pavement depth to nzt_soil. Please note, this is just a
     3772!--       workaround at the moment.
     3773          DO  m = 1, surf_lsm_h%ns
     3774             IF ( surf_lsm_h%pavement_surface(m) )  THEN
     3775
     3776                i = surf_lsm_h%i(m)
     3777                j = surf_lsm_h%j(m)
     3778
     3779                surf_lsm_h%nzt_pavement(m) = nzt_soil
     3780
     3781                DO  k = nzb_soil, nzt_soil 
     3782                   surf_lsm_h%lambda_h_def(k,m) =                              &
     3783                       pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,k,j,i)
     3784                   surf_lsm_h%rho_c_total_def(k,m) =                           &
     3785                       pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,k,j,i)
     3786                ENDDO
     3787
     3788             ENDIF
     3789          ENDDO
     3790          DO  l = 0, 3
     3791             DO  m = 1, surf_lsm_v(l)%ns
     3792                IF ( surf_lsm_v(l)%pavement_surface(m) )  THEN
     3793
     3794                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
     3795                                                surf_lsm_v(l)%building_covered(m) )
     3796                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
     3797                                                surf_lsm_v(l)%building_covered(m) )
     3798
     3799                   surf_lsm_v(l)%nzt_pavement(m) = nzt_soil
     3800
     3801                   DO  k = nzb_soil, nzt_soil 
     3802                      surf_lsm_v(l)%lambda_h_def(k,m) =                        &
     3803                       pavement_subsurface_pars_f%pars_xyz(ind_p_lambda_h,k,j,i)
     3804                      surf_lsm_v(l)%rho_c_total_def(k,m) =                     &
     3805                       pavement_subsurface_pars_f%pars_xyz(ind_p_rho_c,k,j,i)
     3806                   ENDDO
     3807
     3808                ENDIF
     3809             ENDDO
     3810          ENDDO
     3811       ENDIF
    24973812!
    24983813!--    Initial run actions
    24993814       IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
    25003815!
    2501 !--       First, for horizontal surfaces
    2502 !
    2503 !--       Set artifical values for ts and us so that r_a has its initial value
    2504 !--       for the first time step. Only for interior core domain, not for ghost points
     3816!--       First, initialize soil temperature and moisture.
     3817!--       According to the initialization for surface and soil parameters,
     3818!--       initialize soil moisture and temperature via a level approach. This
     3819!--       is to assure that all surface elements are initialized, even if
     3820!--       data provided from input file contains fill values at some locations.
     3821!--       Level 1, initialization via profiles given in parameter file
    25053822          DO  m = 1, surf_lsm_h%ns
    2506 
    2507              t_soil_h%var_2d(:,m)    = 0.0_wp
    2508              m_soil_h%var_2d(:,m)    = 0.0_wp
    2509              m_liq_h%var_1d(m)       = 0.0_wp
    2510 
    2511 !--          Map user settings of T and q for each soil layer
    2512 !--          (make sure that the soil moisture does not drop below the permanent
    2513 !--          wilting point) -> problems with devision by zero)
    2514              IF ( surf_lsm_h%water_surface(m) )  THEN
    2515 
    2516                 surf_lsm_h%z0(m)          = z0_water
    2517                 surf_lsm_h%z0h(m)         = z0h_water
    2518                 surf_lsm_h%z0q(m)         = z0h_water
    2519                 t_soil_h%var_2d(:,m)      = water_temperature
    2520                 surf_lsm_h%lambda_surface_s(m) = 1.0E10_wp
    2521                 surf_lsm_h%lambda_surface_u(m) = 1.0E10_wp               
    2522                 surf_lsm_h%c_surface(m)        = 0.0_wp
    2523                
    2524              ELSEIF ( surf_lsm_h%pavement_surface(m) )  THEN
    2525 
    2526                 surf_lsm_h%z0(m)          = z0_pavement
    2527                 surf_lsm_h%z0h(m)         = z0h_pavement
    2528                 surf_lsm_h%z0q(m)         = z0h_pavement   
     3823             IF ( surf_lsm_h%vegetation_surface(m)  .OR.                       &
     3824                  surf_lsm_h%pavement_surface(m) )  THEN
    25293825                DO  k = nzb_soil, nzt_soil 
    2530                    t_soil_h%var_2d(k,m)   = soil_temperature(k)
    2531                    m_soil_h%var_2d(k,m)   = soil_moisture(k)
     3826                   t_soil_h%var_2d(k,m) = soil_temperature(k)
     3827                   m_soil_h%var_2d(k,m) = soil_moisture(k)
    25323828                ENDDO
    25333829                t_soil_h%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1)
    2534                 surf_lsm_h%lambda_surface_s(m) = pavement_heat_conduct         &
    2535                                                  * ddz_soil(nzb_soil)          &
    2536                                                  * 2.0_wp
    2537                 surf_lsm_h%lambda_surface_u(m) = surf_lsm_h%lambda_surface_s(m)                                         
    2538                 surf_lsm_h%c_surface(m)        = pavement_heat_capacity        &
    2539                                                  * dz_soil(nzb_soil)           &
    2540                                                  * 0.25_wp
    2541 
    2542                 IF ( pavement_type /= 0 )  THEN
     3830             ENDIF
     3831          ENDDO
     3832          DO  l = 0, 3
     3833             DO  m = 1, surf_lsm_v(l)%ns
     3834                IF ( surf_lsm_v(l)%vegetation_surface(m)  .OR.                 &
     3835                     surf_lsm_v(l)%pavement_surface(m) )  THEN
    25433836                   DO  k = nzb_soil, nzt_soil 
    2544                       surf_lsm_h%lambda_h_def(k,m)      =  pavement_subsurface_pars_1(k,pavement_type)                       
    2545                       surf_lsm_h%rho_c_total_def(k,m)   =  pavement_subsurface_pars_2(k,pavement_type)
     3837                      t_soil_v(l)%var_2d(k,m) = soil_temperature(k)
     3838                      m_soil_v(l)%var_2d(k,m) = soil_moisture(k)
    25463839                   ENDDO
    2547                 ELSE
    2548                     surf_lsm_h%lambda_h_def(:,m)     = pavement_heat_conduct
    2549                     surf_lsm_h%rho_c_total_def(:,m)  = pavement_heat_capacity
     3840                   t_soil_v(l)%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1)
    25503841                ENDIF
    2551 
    2552              ELSEIF ( surf_lsm_h%vegetation_surface(m) )  THEN
    2553 
    2554                 surf_lsm_h%z0(m)          = z0_vegetation
    2555                 surf_lsm_h%z0h(m)         = z0h_vegetation
    2556                 surf_lsm_h%z0q(m)         = z0h_vegetation   
    2557                 DO  k = nzb_soil, nzt_soil                       
    2558                    t_soil_h%var_2d(k,m)   = soil_temperature(k)
    2559                    m_soil_h%var_2d(k,m)   = soil_moisture(k)       
    2560                 ENDDO
    2561                 t_soil_h%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1)
    2562                 surf_lsm_h%lambda_surface_s(m) = lambda_surface_stable
    2563                 surf_lsm_h%lambda_surface_u(m) = lambda_surface_unstable
    2564                 surf_lsm_h%c_surface(m)        = c_surface
    2565 
     3842             ENDDO
     3843          ENDDO
     3844!
     3845!--       Level 2, if soil moisture and/or temperature  are
     3846!--       provided from file, interpolate / extrapolate the provided data
     3847!--       onto the respective soil layers. Please note, both, zs as well as
     3848!--       init_3d%z_soil indicate a depth with positive values, so that no
     3849!--       distinction between atmosphere is required concerning interpolation.
     3850!--       Start with soil moisture
     3851          IF ( init_3d%from_file_msoil )  THEN
     3852
     3853             IF ( init_3d%lod_msoil == 1 )  THEN
     3854                DO  m = 1, surf_lsm_h%ns
     3855
     3856                   CALL netcdf_data_input_interpolate(                         &
     3857                                       m_soil_h%var_2d(nzb_soil:nzt_soil,m),   &
     3858                                       init_3d%msoil_init(:),                  &
     3859                                       zs(nzb_soil:nzt_soil), init_3d%z_soil,  &
     3860                                       nzb_soil, nzt_soil,                     &
     3861                                       nzb_soil, init_3d%nzs-1 )
     3862                ENDDO
     3863                DO  l = 0, 3
     3864                   DO  m = 1, surf_lsm_v(l)%ns
     3865
     3866                      CALL netcdf_data_input_interpolate(                      &
     3867                                       m_soil_v(l)%var_2d(nzb_soil:nzt_soil,m),&
     3868                                       init_3d%msoil_init(:),                  &
     3869                                       zs(nzb_soil:nzt_soil), init_3d%z_soil,  &
     3870                                       nzb_soil, nzt_soil,                     &
     3871                                       nzb_soil, init_3d%nzs-1 )
     3872                   ENDDO
     3873                ENDDO
     3874             ELSE
     3875
     3876                DO  m = 1, surf_lsm_h%ns
     3877                   i = surf_lsm_h%i(m)
     3878                   j = surf_lsm_h%j(m)
     3879
     3880                   IF ( init_3d%msoil(0,j,i) /= init_3d%fill_msoil )           &
     3881                      CALL netcdf_data_input_interpolate(                      &
     3882                                       m_soil_h%var_2d(nzb_soil:nzt_soil,m),   &
     3883                                       init_3d%msoil(:,j,i),                   &
     3884                                       zs(nzb_soil:nzt_soil), init_3d%z_soil,  &
     3885                                       nzb_soil, nzt_soil,                     &
     3886                                       nzb_soil, init_3d%nzs-1 )
     3887                ENDDO
     3888                DO  l = 0, 3
     3889                   DO  m = 1, surf_lsm_v(l)%ns
     3890                      i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,   &
     3891                                             surf_lsm_v(l)%building_covered(m) )
     3892                      j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,   &
     3893                                             surf_lsm_v(l)%building_covered(m) )
     3894
     3895                      IF ( init_3d%msoil(0,j,i) /= init_3d%fill_msoil )        &
     3896                         CALL netcdf_data_input_interpolate(                   &
     3897                                       m_soil_v(l)%var_2d(nzb_soil:nzt_soil,m),&
     3898                                       init_3d%msoil(:,j,i),                   &
     3899                                       zs(nzb_soil:nzt_soil), init_3d%z_soil,  &
     3900                                       nzb_soil, nzt_soil,                     &
     3901                                       nzb_soil, init_3d%nzs-1 )
     3902                   ENDDO
     3903                ENDDO
    25663904             ENDIF
    2567  
     3905
     3906          ENDIF
     3907!
     3908!--       Soil temperature
     3909          IF ( init_3d%from_file_tsoil )  THEN
     3910
     3911             IF ( init_3d%lod_tsoil == 1 )  THEN ! change to 1 if provided correctly by INIFOR
     3912                DO  m = 1, surf_lsm_h%ns
     3913
     3914                   CALL netcdf_data_input_interpolate(                         &
     3915                                       t_soil_h%var_2d(nzb_soil:nzt_soil,m),   &
     3916                                       init_3d%tsoil_init(:),                  &
     3917                                       zs(nzb_soil:nzt_soil), init_3d%z_soil,  &
     3918                                       nzb_soil, nzt_soil,                     &
     3919                                       nzb_soil, init_3d%nzs-1 )
     3920                   t_soil_h%var_2d(nzt_soil+1,m) = t_soil_h%var_2d(nzt_soil,m)
     3921                ENDDO
     3922                DO  l = 0, 3
     3923                   DO  m = 1, surf_lsm_v(l)%ns
     3924
     3925                      CALL netcdf_data_input_interpolate(                      &
     3926                                       t_soil_v(l)%var_2d(nzb_soil:nzt_soil,m),&
     3927                                       init_3d%tsoil_init(:),                  &
     3928                                       zs(nzb_soil:nzt_soil), init_3d%z_soil,  &
     3929                                       nzb_soil, nzt_soil,                     &
     3930                                       nzb_soil, init_3d%nzs-1 )
     3931                      t_soil_v(l)%var_2d(nzt_soil+1,m) =                       &
     3932                                                 t_soil_v(l)%var_2d(nzt_soil,m)
     3933                   ENDDO
     3934                ENDDO
     3935             ELSE
     3936
     3937                DO  m = 1, surf_lsm_h%ns
     3938                   i = surf_lsm_h%i(m)
     3939                   j = surf_lsm_h%j(m)
     3940
     3941                   IF ( init_3d%msoil(0,j,i) /= init_3d%fill_msoil )           &
     3942                      CALL netcdf_data_input_interpolate(                      &
     3943                                       t_soil_h%var_2d(nzb_soil:nzt_soil,m),   &
     3944                                       init_3d%tsoil(:,j,i),                   &
     3945                                       zs(nzb_soil:nzt_soil), init_3d%z_soil,  &
     3946                                       nzb_soil, nzt_soil,                     &
     3947                                       nzb_soil, init_3d%nzs-1 )
     3948                   t_soil_h%var_2d(nzt_soil+1,m) = t_soil_h%var_2d(nzt_soil,m)
     3949                ENDDO
     3950                DO  l = 0, 3
     3951                   DO  m = 1, surf_lsm_v(l)%ns
     3952                      i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,   &
     3953                                                surf_lsm_v(l)%building_covered(m) )
     3954                      j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,   &
     3955                                                surf_lsm_v(l)%building_covered(m) )
     3956
     3957                      IF ( init_3d%msoil(0,j,i) /= init_3d%fill_msoil )        &
     3958                         CALL netcdf_data_input_interpolate(                   &
     3959                                       t_soil_v(l)%var_2d(nzb_soil:nzt_soil,m),&
     3960                                       init_3d%tsoil(:,j,i),                   &
     3961                                       zs(nzb_soil:nzt_soil), init_3d%z_soil,  &
     3962                                       nzb_soil, nzt_soil,                     &
     3963                                       nzb_soil, init_3d%nzs-1 )
     3964                      t_soil_v(l)%var_2d(nzt_soil+1,m) =                       &
     3965                                                 t_soil_v(l)%var_2d(nzt_soil,m)
     3966                   ENDDO
     3967                ENDDO
     3968             ENDIF
     3969          ENDIF
     3970!
     3971!--       Further initialization
     3972          DO  m = 1, surf_lsm_h%ns
     3973
    25683974             i   = surf_lsm_h%i(m)           
    25693975             j   = surf_lsm_h%j(m)
     
    25753981             IF ( surf_lsm_h%lambda_surface_s(m) == 0.0_wp )  THEN
    25763982                t_surface_h%var_1d(:)    = t_soil_h%var_2d(nzb_soil,m)
    2577                 surf_lsm_h%pt_surface(:) = t_soil_h%var_2d(nzb_soil,m) / exn
     3983                surf_lsm_h%pt_surface(m) = t_soil_h%var_2d(nzb_soil,m) / exn
    25783984             ELSE
    2579                 t_surface_h%var_1d(:)    = pt_surface * exn
    2580                 surf_lsm_h%pt_surface(:) = pt_surface 
     3985                t_surface_h%var_1d(m)    = pt(k-1,j,i) * exn
     3986                surf_lsm_h%pt_surface(m) = pt(k-1,j,i) 
    25813987             ENDIF
    25823988             
    2583              IF ( cloud_physics  .OR.  cloud_droplets ) THEN
    2584                 pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
     3989             IF ( cloud_physics  .OR. cloud_droplets ) THEN
     3990                surf_lsm_h%pt1(m) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
    25853991             ELSE
    2586                 pt1 = pt(k,j,i)
    2587              ENDIF
     3992                surf_lsm_h%pt1(m) = pt(k,j,i)
     3993             ENDIF
     3994
    25883995
    25893996!
    25903997!--          Assure that r_a cannot be zero at model start
    2591              IF ( pt1 == surf_lsm_h%pt_surface(m) )  pt1 = pt1 + 1.0E-10_wp
     3998             IF ( surf_lsm_h%pt1(m) == surf_lsm_h%pt_surface(m) )              &
     3999                surf_lsm_h%pt1(m) = surf_lsm_h%pt1(m) + 1.0E-20_wp
    25924000
    25934001             surf_lsm_h%us(m)   = 0.1_wp
    2594              surf_lsm_h%ts(m)   = ( pt1 - surf_lsm_h%pt_surface(m) )           &
     4002             surf_lsm_h%ts(m)   = ( surf_lsm_h%pt1(m) - surf_lsm_h%pt_surface(m) )&
    25954003                                  / surf_lsm_h%r_a(m)
    25964004             surf_lsm_h%shf(m)  = - surf_lsm_h%us(m) * surf_lsm_h%ts(m)        &
    25974005                                  * rho_surface
    2598 
    25994006         ENDDO
    2600 
    26014007!
    26024008!--      Vertical surfaces
    26034009         DO  l = 0, 3
    26044010             DO  m = 1, surf_lsm_v(l)%ns
    2605              
    2606                 t_soil_v(l)%var_2d    = 0.0_wp
    2607                 m_soil_v(l)%var_2d    = 0.0_wp
    2608                 m_liq_v(l)%var_1d  = 0.0_wp
    2609 
    2610              
    2611 !--             Map user settings of T and q for each soil layer
    2612 !--             (make sure that the soil moisture does not drop below the permanent
    2613 !--             wilting point) -> problems with devision by zero)
    2614                 IF ( surf_lsm_v(l)%water_surface(m) )  THEN
    2615                    surf_lsm_v(l)%z0(m)          = z0_water
    2616                    surf_lsm_v(l)%z0h(m)         = z0h_water
    2617                    surf_lsm_v(l)%z0q(m)         = z0h_water                   
    2618                    t_soil_v(l)%var_2d(:,m)      = water_temperature
    2619                    surf_lsm_v(l)%lambda_surface_s(m) = 1.0E10_wp
    2620                    surf_lsm_v(l)%lambda_surface_u(m) = 1.0E10_wp             
    2621                    surf_lsm_v(l)%c_surface(m)        = 0.0_wp
    2622                    
    2623                 ELSEIF ( surf_lsm_v(l)%pavement_surface(m) )  THEN
    2624                    surf_lsm_v(l)%z0(m)          = z0_pavement
    2625                    surf_lsm_v(l)%z0h(m)         = z0h_pavement
    2626                    surf_lsm_v(l)%z0q(m)         = z0h_pavement 
    2627                    DO  k = nzb_soil, nzt_soil 
    2628                       t_soil_v(l)%var_2d(k,m)   = soil_temperature(k)
    2629                       m_soil_v(l)%var_2d(k,m)   = soil_moisture(k)
    2630                    ENDDO
    2631                    t_soil_v(l)%var_2d(nzt_soil+1,m)  = soil_temperature(nzt_soil+1)
    2632                    surf_lsm_v(l)%lambda_surface_s(m) = pavement_heat_conduct   &
    2633                                                        * ddz_soil(nzb_soil)    &
    2634                                                        * 2.0_wp
    2635                    surf_lsm_v(l)%lambda_surface_u(m) = surf_lsm_v(l)%lambda_surface_s(m)
    2636                    surf_lsm_v(l)%c_surface(m) = pavement_heat_capacity         &
    2637                                                 * dz_soil(nzb_soil)            &
    2638                                                 * 0.25_wp
    2639                 IF ( pavement_type /= 0 )  THEN
    2640                    DO  k = nzb_soil, nzt_soil 
    2641                       surf_lsm_v(l)%lambda_h_def(k,m)    =  pavement_subsurface_pars_1(k,pavement_type)                       
    2642                       surf_lsm_v(l)%rho_c_total_def(k,m) =  pavement_subsurface_pars_2(k,pavement_type) 
    2643                    ENDDO
    2644                 ELSE
    2645                     surf_lsm_v(l)%lambda_h_def(:,m)     = pavement_heat_conduct
    2646                     surf_lsm_v(l)%rho_c_total_def(:,m)  = pavement_heat_capacity
    2647                 ENDIF
    2648                        
    2649                    
    2650                 ELSEIF ( surf_lsm_v(l)%vegetation_surface(m) )  THEN
    2651 
    2652                    surf_lsm_v(l)%z0(m)          = z0_vegetation
    2653                    surf_lsm_v(l)%z0h(m)         = z0h_vegetation
    2654                    surf_lsm_v(l)%z0q(m)         = z0h_vegetation   
    2655                    DO  k = nzb_soil, nzt_soil                       
    2656                       t_soil_v(l)%var_2d(k,m)   = soil_temperature(k)
    2657                       m_soil_v(l)%var_2d(k,m)   = soil_moisture(k)       
    2658                    ENDDO
    2659                    t_soil_v(l)%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1)
    2660                    surf_lsm_v(l)%lambda_surface_s(m) = lambda_surface_stable
    2661                    surf_lsm_v(l)%lambda_surface_u(m) = lambda_surface_unstable
    2662                    surf_lsm_v(l)%c_surface(m)        = c_surface
    2663                 ENDIF       
    2664          
     4011                i   = surf_lsm_v(l)%i(m)           
     4012                j   = surf_lsm_v(l)%j(m)
     4013                k   = surf_lsm_v(l)%k(m)         
    26654014!
    26664015!--             Calculate surface temperature. In case of bare soil, the surface
     
    26684017!--             layer
    26694018                IF ( surf_lsm_v(l)%lambda_surface_s(m) == 0.0_wp )  THEN
    2670                    t_surface_v(l)%var_1d(:)      = t_soil_v(l)%var_2d(nzb_soil,m)
    2671                    surf_lsm_v(l)%pt_surface(:)   = t_soil_v(l)%var_2d(nzb_soil,m) / exn
     4019                   t_surface_v(l)%var_1d(m)      = t_soil_v(l)%var_2d(nzb_soil,m)
     4020                   surf_lsm_v(l)%pt_surface(m)   = t_soil_v(l)%var_2d(nzb_soil,m) / exn
    26724021                ELSE
    2673                    t_surface_v(l)%var_1d(:)      = pt_surface * exn
    2674                    surf_lsm_v(l)%pt_surface(:)   = pt_surface               
     4022                   j_off = surf_lsm_v(l)%joff
     4023                   i_off = surf_lsm_v(l)%ioff
     4024
     4025                   t_surface_v(l)%var_1d(m)      = pt(k,j+j_off,i+i_off) * exn
     4026                   surf_lsm_v(l)%pt_surface(m)   = pt(k,j+j_off,i+i_off)           
    26754027                ENDIF
    26764028
     4029
     4030                IF ( cloud_physics  .OR. cloud_droplets ) THEN
     4031                   surf_lsm_v(l)%pt1(m) = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
     4032                ELSE
     4033                   surf_lsm_v(l)%pt1(m) = pt(k,j,i)
     4034                ENDIF
     4035
     4036!
     4037!--             Assure that r_a cannot be zero at model start
     4038                IF ( surf_lsm_v(l)%pt1(m) == surf_lsm_v(l)%pt_surface(m) )     &
     4039                     surf_lsm_v(l)%pt1(m) = surf_lsm_v(l)%pt1(m) + 1.0E-20_wp
    26774040!
    26784041!--             Set artifical values for ts and us so that r_a has its initial value
    26794042!--             for the first time step. Only for interior core domain, not for ghost points
    2680 
    2681                 i   = surf_lsm_v(l)%i(m)           
    2682                 j   = surf_lsm_v(l)%j(m)
    2683                 k   = surf_lsm_v(l)%k(m)
    2684 
    2685                 IF ( cloud_physics  .OR.  cloud_droplets )  THEN
    2686                    pt1 = pt(k,j,i) + l_d_cp * pt_d_t(k) * ql(k,j,i)
    2687                 ELSE
    2688                    pt1 = pt(k,j,i)
    2689                 ENDIF
    2690 
    2691 !
    2692 !--             Assure that r_a cannot be zero at model start
    2693                 IF ( pt1 == surf_lsm_v(l)%pt_surface(m) )  pt1 = pt1 + 1.0E-10_wp
    2694 
    26954043                surf_lsm_v(l)%us(m)   = 0.1_wp
    2696                 surf_lsm_v(l)%ts(m)   = ( pt1 - surf_lsm_v(l)%pt_surface(m) ) / surf_lsm_v(l)%r_a(m)
    2697                 surf_lsm_v(l)%shf(m)  = - surf_lsm_v(l)%us(m) * surf_lsm_v(l)%ts(m) * rho_surface
     4044                surf_lsm_v(l)%ts(m)   = ( surf_lsm_v(l)%pt1(m) - surf_lsm_v(l)%pt_surface(m) ) /&
     4045                                          surf_lsm_v(l)%r_a(m)
     4046                surf_lsm_v(l)%shf(m)  = - surf_lsm_v(l)%us(m) *                &
     4047                                          surf_lsm_v(l)%ts(m) * rho_surface
     4048
    26984049             ENDDO
    26994050          ENDDO
     
    27154066!
    27164067!--          Set index offset of surface element, seen from atmospheric grid point
    2717              IF ( l == 0 )  THEN
    2718                 j_off = -1
    2719                 i_off = 0
    2720              ELSEIF ( l == 1 )  THEN
    2721                 j_off = 1
    2722                 i_off = 0
    2723              ELSEIF ( l == 2 )  THEN
    2724                 j_off = 0
    2725                 i_off = -1
    2726              ELSEIF ( l == 3 )  THEN
    2727                 j_off = 0
    2728                 i_off = 1
    2729              ENDIF
     4068             j_off = surf_lsm_v(l)%joff
     4069             i_off = surf_lsm_v(l)%ioff
     4070
    27304071             DO  m = 1, surf_lsm_v(l)%ns
    27314072                i   = surf_lsm_v(l)%i(m)           
     
    27384079       ENDIF
    27394080!
    2740 !--    Initialize root fraction
    2741 !--    Horizontal surfaces
     4081!--    Level 1 initialization of root distribution - provided by the user via
     4082!--    via namelist.
    27424083       DO  m = 1, surf_lsm_h%ns
    2743           i   = surf_lsm_h%i(m)           
    2744           j   = surf_lsm_h%j(m)
    2745 
    27464084          DO  k = nzb_soil, nzt_soil
    27474085             surf_lsm_h%root_fr(k,m) = root_fraction(k)
    27484086          ENDDO
    27494087       ENDDO
    2750 !
    2751 !--    Vertical surfaces
     4088
    27524089       DO  l = 0, 3
    27534090          DO  m = 1, surf_lsm_v(l)%ns
    2754              i   = surf_lsm_v(l)%i(m)           
    2755              j   = surf_lsm_v(l)%j(m)
    2756 
    27574091             DO  k = nzb_soil, nzt_soil
    27584092                surf_lsm_v(l)%root_fr(k,m) = root_fraction(k)
     
    27624096
    27634097!
     4098!--    Level 2 initialization of root distribution.
    27644099!--    When no root distribution is given by the user, use look-up table to prescribe
    2765 !--    the root fraction in the individual soil layers
     4100!--    the root fraction in the individual soil layers.
    27664101       IF ( ALL( root_fraction == 9999999.9_wp ) )  THEN
    2767 
    27684102!
    27694103!--       First, calculate the index bounds for integration
     
    28204154!--       Map calculated root fractions
    28214155          DO  m = 1, surf_lsm_h%ns
    2822              i   = surf_lsm_h%i(m)           
    2823              j   = surf_lsm_h%j(m)
    28244156             DO  k = nzb_soil, nzt_soil
    28254157                IF ( surf_lsm_h%pavement_surface(m)  .AND.                     &
    2826                      k <= nzt_pavement )  THEN
     4158                     k <= surf_lsm_h%nzt_pavement(m) )  THEN
    28274159                   surf_lsm_h%root_fr(k,m) = 0.0_wp
    28284160                ELSE
    28294161                   surf_lsm_h%root_fr(k,m) = root_fraction(k)
    28304162                ENDIF
    2831              ENDDO
    2832 
     4163
     4164             ENDDO
    28334165!
    28344166!--          Normalize so that the sum = 1. Only relevant when the root         
     
    28404172                ENDDO
    28414173             ENDIF
    2842 
    28434174          ENDDO
    28444175          DO  l = 0, 3
    28454176             DO  m = 1, surf_lsm_v(l)%ns
    2846                 i   = surf_lsm_v(l)%i(m)           
    2847                 j   = surf_lsm_v(l)%j(m)
    2848 
    28494177                DO  k = nzb_soil, nzt_soil
    28504178                   IF ( surf_lsm_v(l)%pavement_surface(m)  .AND.               &
    2851                         k <= nzt_pavement )  THEN
     4179                        k <= surf_lsm_h%nzt_pavement(m) )  THEN
    28524180                      surf_lsm_v(l)%root_fr(k,m) = 0.0_wp
    28534181                   ELSE
     
    28644192                   ENDDO
    28654193                ENDIF
    2866 
    28674194             ENDDO
    28684195           ENDDO
    28694196       ENDIF
    28704197!
    2871 !--    Map vegetation parameters to the respective 2D arrays
    2872        surf_lsm_h%r_canopy_min      = min_canopy_resistance
    2873        surf_lsm_h%lai               = leaf_area_index
    2874        surf_lsm_h%c_veg             = vegetation_coverage
    2875        surf_lsm_h%g_d               = canopy_resistance_coefficient
    2876        surf_lsm_h%f_sw_in           = f_shortwave_incoming
    2877 
    2878 !--    Vertical surfaces
    2879        DO  l = 0, 3
    2880 
    2881 !
    2882 !--       Map vegetation parameters to the respective 2D arrays
    2883           surf_lsm_v(l)%r_canopy_min      = min_canopy_resistance
    2884           surf_lsm_v(l)%lai               = leaf_area_index
    2885           surf_lsm_v(l)%c_veg             = vegetation_coverage
    2886           surf_lsm_v(l)%g_d               = canopy_resistance_coefficient
    2887           surf_lsm_v(l)%f_sw_in           = f_shortwave_incoming
    2888        ENDDO
     4198!--    Level 3 initialization of root distribution.
     4199!--    Take value from file
     4200       IF ( root_area_density_lsm_f%from_file )  THEN
     4201          DO  m = 1, surf_lsm_h%ns
     4202             IF ( surf_lsm_h%vegetation_surface(m) )  THEN
     4203                i = surf_lsm_h%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,            &
     4204                                             surf_lsm_v(l)%building_covered(m) )
     4205                j = surf_lsm_h%j(m) + MERGE( 0, surf_lsm_v(l)%joff,            &
     4206                                             surf_lsm_v(l)%building_covered(m) )
     4207                DO  k = nzb_soil, nzt_soil
     4208                   surf_lsm_h%root_fr(k,m) = root_area_density_lsm_f%var(k,j,i)
     4209                ENDDO
     4210
     4211             ENDIF
     4212          ENDDO
     4213
     4214          DO  l = 0, 3
     4215             DO  m = 1, surf_lsm_v(l)%ns
     4216                IF ( surf_lsm_v(l)%vegetation_surface(m) )  THEN
     4217                   i = surf_lsm_v(l)%i(m) + MERGE( 0, surf_lsm_v(l)%ioff,      &
     4218                                                   surf_lsm_v(l)%building_covered(m) )
     4219                   j = surf_lsm_v(l)%j(m) + MERGE( 0, surf_lsm_v(l)%joff,      &
     4220                                                   surf_lsm_v(l)%building_covered(m) )
     4221
     4222                   DO  k = nzb_soil, nzt_soil
     4223                      surf_lsm_v(l)%root_fr(k,m) = root_area_density_lsm_f%var(k,j,i)
     4224                   ENDDO
     4225
     4226                ENDIF
     4227             ENDDO
     4228          ENDDO
     4229
     4230       ENDIF
    28894231 
    28904232!
     
    29124254!--    Store initial profiles of t_soil and m_soil (assuming they are
    29134255!--    horizontally homogeneous on this PE)
    2914        hom(nzb_soil:nzt_soil,1,90,:)  = SPREAD( t_soil_h%var_2d(nzb_soil:nzt_soil,1),  &
    2915                                                 2, statistic_regions+1 )
    2916        hom(nzb_soil:nzt_soil,1,92,:)  = SPREAD( m_soil_h%var_2d(nzb_soil:nzt_soil,1),  &
    2917                                                 2, statistic_regions+1 )
     4256!--    DEACTIVATED FOR NOW - leads to error when number of locations with
     4257!--    soil model is zero on a PE.
     4258!        hom(nzb_soil:nzt_soil,1,90,:)  = SPREAD( t_soil_h%var_2d(nzb_soil:nzt_soil,1),  &
     4259!                                                 2, statistic_regions+1 )
     4260!        hom(nzb_soil:nzt_soil,1,92,:)  = SPREAD( m_soil_h%var_2d(nzb_soil:nzt_soil,1),  &
     4261!                                                 2, statistic_regions+1 )
     4262
     4263!
     4264!--    Finally, make some consistency checks.
     4265!--    Ceck for eck for illegal combination of LAI and vegetation coverage.
     4266       IF ( ANY( .NOT. surf_lsm_h%pavement_surface  .AND.                      &
     4267                 surf_lsm_h%lai == 0.0_wp  .AND.  surf_lsm_h%c_veg == 1.0_wp ) &
     4268          )  THEN
     4269          message_string = 'For non-pavement surfaces the combination ' //     &
     4270                           ' lai = 0.0 and c_veg = 1.0 is not allowed.'
     4271          CALL message( 'lsm_read_restart_data', 'PA0999', 2, 2, 0, 6, 0 )
     4272       ENDIF
     4273
     4274       DO  l = 0, 3
     4275          IF ( ANY( .NOT. surf_lsm_v(l)%pavement_surface  .AND.                &
     4276                    surf_lsm_v(l)%lai == 0.0_wp  .AND.                         &
     4277                    surf_lsm_v(l)%c_veg == 1.0_wp ) )  THEN
     4278             message_string = 'For non-pavement surfaces the combination ' //  &
     4279                              ' lai = 0.0 and c_veg = 1.0 is not allowed.'
     4280             CALL message( 'lsm_read_restart_data', 'PA0999', 2, 2, 0, 6, 0 )
     4281          ENDIF
     4282       ENDDO
     4283
    29184284
    29194285    END SUBROUTINE lsm_init
     
    29434309!
    29444310!--    Horizontal surfaces
    2945        ALLOCATE ( m_liq_h%var_1d(1:surf_lsm_h%ns)                        )
    29464311       ALLOCATE ( m_liq_h_p%var_1d(1:surf_lsm_h%ns)                      )
    29474312       ALLOCATE ( t_surface_h%var_1d(1:surf_lsm_h%ns)                    )
    29484313       ALLOCATE ( t_surface_h_p%var_1d(1:surf_lsm_h%ns)                  )
    2949        ALLOCATE ( m_soil_h%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns)     )
    29504314       ALLOCATE ( m_soil_h_p%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns)   )
    2951        ALLOCATE ( t_soil_h%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns)   )
    29524315       ALLOCATE ( t_soil_h_p%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) )
    29534316
     
    29634326          ALLOCATE ( t_soil_v(l)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(l)%ns)   )
    29644327          ALLOCATE ( t_soil_v_p(l)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(l)%ns) )
     4328       ENDDO
     4329!
     4330!--    Allocate soil temperature and moisture. As these variables might be
     4331!--    already allocated in case of restarts, check this.
     4332       IF ( .NOT. ALLOCATED( m_liq_h%var_1d ) )                                &
     4333          ALLOCATE ( m_liq_h%var_1d(1:surf_lsm_h%ns) )
     4334       IF ( .NOT. ALLOCATED( m_soil_h%var_2d ) )                               &
     4335          ALLOCATE ( m_soil_h%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
     4336       IF ( .NOT. ALLOCATED( t_soil_h%var_2d ) )                               &
     4337          ALLOCATE ( t_soil_h%var_2d(nzb_soil:nzt_soil,1:surf_lsm_h%ns) )
     4338
     4339       DO  l = 0, 3
     4340          IF ( .NOT. ALLOCATED( m_liq_v(l)%var_1d ) )                          &
     4341             ALLOCATE ( m_liq_v(l)%var_1d(1:surf_lsm_v(l)%ns) )
     4342          IF ( .NOT. ALLOCATED( m_soil_v(l)%var_2d ) )                         &
     4343             ALLOCATE ( m_soil_v(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
     4344          IF ( .NOT. ALLOCATED( t_soil_v(l)%var_2d ) )                         &
     4345             ALLOCATE ( t_soil_v(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns) )
    29654346       ENDDO
    29664347#else
     
    29884369       ENDDO
    29894370#endif
     4371!
     4372!--    Allocate array for heat flux in W/m2, required for radiation?
     4373!--    Consider to remove this array
     4374       ALLOCATE( surf_lsm_h%surfhf(1:surf_lsm_h%ns) )
     4375       DO  l = 0, 3
     4376          ALLOCATE( surf_lsm_v(l)%surfhf(1:surf_lsm_v(l)%ns) )
     4377       ENDDO
     4378
    29904379
    29914380!
     
    30044393          ALLOCATE ( tt_soil_v_m(l)%var_2d(nzb_soil:nzt_soil,1:surf_lsm_v(l)%ns)  )
    30054394       ENDDO
    3006 !
    3007 !--    Allocate skin-surface temperature
    3008        ALLOCATE ( surf_lsm_h%pt_surface(1:surf_lsm_h%ns) )
    3009        DO  l = 0, 3
    3010           ALLOCATE ( surf_lsm_v(l)%pt_surface(1:surf_lsm_v(l)%ns) )
    3011        ENDDO
    30124395
    30134396!
    30144397!--    Allocate 2D vegetation model arrays
    30154398!--    Horizontal surfaces
    3016        ALLOCATE ( surf_lsm_h%building_surface(1:surf_lsm_h%ns) )
    3017        ALLOCATE ( surf_lsm_h%c_liq(1:surf_lsm_h%ns)            )
    3018        ALLOCATE ( surf_lsm_h%c_surface(1:surf_lsm_h%ns)        )
    3019        ALLOCATE ( surf_lsm_h%c_veg(1:surf_lsm_h%ns)            )
    3020        ALLOCATE ( surf_lsm_h%f_sw_in(1:surf_lsm_h%ns)          )
    3021        ALLOCATE ( surf_lsm_h%ghf(1:surf_lsm_h%ns)              )
    3022        ALLOCATE ( surf_lsm_h%g_d(1:surf_lsm_h%ns)              )
    3023        ALLOCATE ( surf_lsm_h%lai(1:surf_lsm_h%ns)              )
    3024        ALLOCATE ( surf_lsm_h%lambda_surface_u(1:surf_lsm_h%ns) )
    3025        ALLOCATE ( surf_lsm_h%lambda_surface_s(1:surf_lsm_h%ns) )
     4399       ALLOCATE ( surf_lsm_h%building_surface(1:surf_lsm_h%ns)    )
     4400       ALLOCATE ( surf_lsm_h%c_liq(1:surf_lsm_h%ns)               )
     4401       ALLOCATE ( surf_lsm_h%c_surface(1:surf_lsm_h%ns)           )
     4402       ALLOCATE ( surf_lsm_h%c_veg(1:surf_lsm_h%ns)               )
     4403       ALLOCATE ( surf_lsm_h%f_sw_in(1:surf_lsm_h%ns)             )
     4404       ALLOCATE ( surf_lsm_h%ghf(1:surf_lsm_h%ns)                 )
     4405       ALLOCATE ( surf_lsm_h%g_d(1:surf_lsm_h%ns)                 )
     4406       ALLOCATE ( surf_lsm_h%lai(1:surf_lsm_h%ns)                 )
     4407       ALLOCATE ( surf_lsm_h%lambda_surface_u(1:surf_lsm_h%ns)    )
     4408       ALLOCATE ( surf_lsm_h%lambda_surface_s(1:surf_lsm_h%ns)    )
     4409       ALLOCATE ( surf_lsm_h%nzt_pavement(1:surf_lsm_h%ns)        )
     4410       ALLOCATE ( surf_lsm_h%pavement_surface(1:surf_lsm_h%ns)    )
     4411       ALLOCATE ( surf_lsm_h%qsws_soil(1:surf_lsm_h%ns)           )
     4412       ALLOCATE ( surf_lsm_h%qsws_liq(1:surf_lsm_h%ns)            )
     4413       ALLOCATE ( surf_lsm_h%qsws_veg(1:surf_lsm_h%ns)            )
     4414       ALLOCATE ( surf_lsm_h%rad_net_l(1:surf_lsm_h%ns)           )
     4415       ALLOCATE ( surf_lsm_h%r_a(1:surf_lsm_h%ns)                 )
     4416       ALLOCATE ( surf_lsm_h%r_canopy(1:surf_lsm_h%ns)            )
     4417       ALLOCATE ( surf_lsm_h%r_soil(1:surf_lsm_h%ns)              )
     4418       ALLOCATE ( surf_lsm_h%r_soil_min(1:surf_lsm_h%ns)          )
     4419       ALLOCATE ( surf_lsm_h%r_s(1:surf_lsm_h%ns)                 )
     4420       ALLOCATE ( surf_lsm_h%r_canopy_min(1:surf_lsm_h%ns)        )
    30264421       ALLOCATE ( surf_lsm_h%vegetation_surface(1:surf_lsm_h%ns)  )
    3027        ALLOCATE ( surf_lsm_h%pavement_surface(1:surf_lsm_h%ns) )
    3028        ALLOCATE ( surf_lsm_h%qsws_soil(1:surf_lsm_h%ns)        )
    3029        ALLOCATE ( surf_lsm_h%qsws_liq(1:surf_lsm_h%ns)         )
    3030        ALLOCATE ( surf_lsm_h%qsws_veg(1:surf_lsm_h%ns)         )
    3031        ALLOCATE ( surf_lsm_h%rad_net_l(1:surf_lsm_h%ns)        )
    3032        ALLOCATE ( surf_lsm_h%r_a(1:surf_lsm_h%ns)              )
    3033        ALLOCATE ( surf_lsm_h%r_canopy(1:surf_lsm_h%ns)         )
    3034        ALLOCATE ( surf_lsm_h%r_soil(1:surf_lsm_h%ns)           )
    3035        ALLOCATE ( surf_lsm_h%r_soil_min(1:surf_lsm_h%ns)       )
    3036        ALLOCATE ( surf_lsm_h%r_s(1:surf_lsm_h%ns)              )
    3037        ALLOCATE ( surf_lsm_h%r_canopy_min(1:surf_lsm_h%ns)     )
    3038        ALLOCATE ( surf_lsm_h%water_surface(1:surf_lsm_h%ns)    )
    3039 
    3040        surf_lsm_h%water_surface     = .FALSE.
    3041        surf_lsm_h%pavement_surface  = .FALSE.
     4422       ALLOCATE ( surf_lsm_h%water_surface(1:surf_lsm_h%ns)       )
     4423
     4424       surf_lsm_h%water_surface        = .FALSE.
     4425       surf_lsm_h%pavement_surface     = .FALSE.
    30424426       surf_lsm_h%vegetation_surface   = .FALSE.
    30434427!
    30444428!--    Vertical surfaces
    30454429       DO  l = 0, 3
    3046           ALLOCATE ( surf_lsm_v(l)%building_surface(1:surf_lsm_v(l)%ns) )
    3047           ALLOCATE ( surf_lsm_v(l)%c_liq(1:surf_lsm_v(l)%ns)            )
    3048           ALLOCATE ( surf_lsm_v(l)%c_surface(1:surf_lsm_v(l)%ns)        )
    3049           ALLOCATE ( surf_lsm_v(l)%c_veg(1:surf_lsm_v(l)%ns)            )
    3050           ALLOCATE ( surf_lsm_v(l)%f_sw_in(1:surf_lsm_v(l)%ns)          )
    3051           ALLOCATE ( surf_lsm_v(l)%ghf(1:surf_lsm_v(l)%ns)              )
    3052           ALLOCATE ( surf_lsm_v(l)%g_d(1:surf_lsm_v(l)%ns)              )
    3053           ALLOCATE ( surf_lsm_v(l)%lai(1:surf_lsm_v(l)%ns)              )
    3054           ALLOCATE ( surf_lsm_v(l)%lambda_surface_u(1:surf_lsm_v(l)%ns) )
    3055           ALLOCATE ( surf_lsm_v(l)%lambda_surface_s(1:surf_lsm_v(l)%ns) )
     4430          ALLOCATE ( surf_lsm_v(l)%building_surface(1:surf_lsm_v(l)%ns)    )
     4431          ALLOCATE ( surf_lsm_v(l)%c_liq(1:surf_lsm_v(l)%ns)               )
     4432          ALLOCATE ( surf_lsm_v(l)%c_surface(1:surf_lsm_v(l)%ns)           )
     4433          ALLOCATE ( surf_lsm_v(l)%c_veg(1:surf_lsm_v(l)%ns)               )
     4434          ALLOCATE ( surf_lsm_v(l)%f_sw_in(1:surf_lsm_v(l)%ns)             )
     4435          ALLOCATE ( surf_lsm_v(l)%ghf(1:surf_lsm_v(l)%ns)                 )
     4436          ALLOCATE ( surf_lsm_v(l)%g_d(1:surf_lsm_v(l)%ns)                 )
     4437          ALLOCATE ( surf_lsm_v(l)%lai(1:surf_lsm_v(l)%ns)                 )
     4438          ALLOCATE ( surf_lsm_v(l)%lambda_surface_u(1:surf_lsm_v(l)%ns)    )
     4439          ALLOCATE ( surf_lsm_v(l)%lambda_surface_s(1:surf_lsm_v(l)%ns)    )
     4440          ALLOCATE ( surf_lsm_v(l)%nzt_pavement(1:surf_lsm_v(l)%ns)        )
     4441          ALLOCATE ( surf_lsm_v(l)%pavement_surface(1:surf_lsm_v(l)%ns)    )
     4442          ALLOCATE ( surf_lsm_v(l)%qsws_soil(1:surf_lsm_v(l)%ns)           )
     4443          ALLOCATE ( surf_lsm_v(l)%qsws_liq(1:surf_lsm_v(l)%ns)            )
     4444          ALLOCATE ( surf_lsm_v(l)%qsws_veg(1:surf_lsm_v(l)%ns)            )
     4445          ALLOCATE ( surf_lsm_v(l)%rad_net_l(1:surf_lsm_v(l)%ns)           )
     4446          ALLOCATE ( surf_lsm_v(l)%r_a(1:surf_lsm_v(l)%ns)                 )
     4447          ALLOCATE ( surf_lsm_v(l)%r_canopy(1:surf_lsm_v(l)%ns)            )
     4448          ALLOCATE ( surf_lsm_v(l)%r_soil(1:surf_lsm_v(l)%ns)              )
     4449          ALLOCATE ( surf_lsm_v(l)%r_soil_min(1:surf_lsm_v(l)%ns)          )
     4450          ALLOCATE ( surf_lsm_v(l)%r_s(1:surf_lsm_v(l)%ns)                 )
     4451          ALLOCATE ( surf_lsm_v(l)%r_canopy_min(1:surf_lsm_v(l)%ns)        )
    30564452          ALLOCATE ( surf_lsm_v(l)%vegetation_surface(1:surf_lsm_v(l)%ns)  )
    3057           ALLOCATE ( surf_lsm_v(l)%pavement_surface(1:surf_lsm_v(l)%ns) )
    3058           ALLOCATE ( surf_lsm_v(l)%qsws_soil(1:surf_lsm_v(l)%ns)        )
    3059           ALLOCATE ( surf_lsm_v(l)%qsws_liq(1:surf_lsm_v(l)%ns)         )
    3060           ALLOCATE ( surf_lsm_v(l)%qsws_veg(1:surf_lsm_v(l)%ns)         )
    3061           ALLOCATE ( surf_lsm_v(l)%rad_net_l(1:surf_lsm_v(l)%ns)        )
    3062           ALLOCATE ( surf_lsm_v(l)%r_a(1:surf_lsm_v(l)%ns)              )
    3063           ALLOCATE ( surf_lsm_v(l)%r_canopy(1:surf_lsm_v(l)%ns)         )
    3064           ALLOCATE ( surf_lsm_v(l)%r_soil(1:surf_lsm_v(l)%ns)           )
    3065           ALLOCATE ( surf_lsm_v(l)%r_soil_min(1:surf_lsm_v(l)%ns)       )
    3066           ALLOCATE ( surf_lsm_v(l)%r_s(1:surf_lsm_v(l)%ns)              )
    3067           ALLOCATE ( surf_lsm_v(l)%r_canopy_min(1:surf_lsm_v(l)%ns)     )
    3068           ALLOCATE ( surf_lsm_v(l)%water_surface(1:surf_lsm_v(l)%ns)    )
     4453          ALLOCATE ( surf_lsm_v(l)%water_surface(1:surf_lsm_v(l)%ns)       )
    30694454
    30704455          surf_lsm_v(l)%water_surface     = .FALSE.
     
    30734458         
    30744459       ENDDO
     4460
    30754461   
    30764462#if ! defined( __nopointer )
     
    30884474       m_soil_v    => m_soil_v_1;    m_soil_v_p    => m_soil_v_2
    30894475       m_liq_v     => m_liq_v_1;     m_liq_v_p     => m_liq_v_2
     4476
    30904477#endif
    30914478
     
    31224509                                  residual_moisture, root_fraction,            &
    31234510                                  saturation_moisture, skip_time_do_lsm,       &
    3124                                   soil_moisture, soil_temperature, soil_type,  &
     4511                                  soil_moisture, soil_temperature,             &
     4512                                  soil_type,                                   &
    31254513                                  surface_type,                                &
    31264514                                  vegetation_coverage, vegetation_type,        &
     
    32044592       IF ( horizontal )  THEN
    32054593          surf           => surf_lsm_h
     4594
    32064595          surf_m_soil    => m_soil_h
    32074596          surf_m_soil_p  => m_soil_h_p
     
    32124601       ELSE
    32134602          surf           => surf_lsm_v(l)
     4603
    32144604          surf_m_soil    => m_soil_v(l)
    32154605          surf_m_soil_p  => m_soil_v_p(l)
     
    32254615             DO  k = nzb_soil, nzt_soil
    32264616
    3227                 IF ( surf%pavement_surface(m)  .AND.  k <= nzt_pavement )  THEN
    3228                    
     4617                IF ( surf%pavement_surface(m)  .AND.                           &
     4618                     k <= surf%nzt_pavement(m) )  THEN
     4619                   
    32294620                   surf%rho_c_total(k,m) = surf%rho_c_total_def(k,m)
    32304621                   lambda_temp(k)        = surf%lambda_h_def(k,m)
     
    32664657             tend(:) = 0.0_wp
    32674658
    3268              tend(nzb_soil) = ( 1.0_wp / surf%rho_c_total(nzb_soil,m) ) *      &
    3269                             ( surf%lambda_h(nzb_soil,m)                        &
    3270                             * ( surf_t_soil%var_2d(nzb_soil+1,m)               &
    3271                             - surf_t_soil%var_2d(nzb_soil,m) )                 &
    3272                             * ddz_soil_center(nzb_soil) + surf%ghf(m) )        &
    3273                             * ddz_soil(nzb_soil)
     4659             tend(nzb_soil) = ( 1.0_wp / surf%rho_c_total(nzb_soil,m) ) *            &
     4660                    ( surf%lambda_h(nzb_soil,m) * ( surf_t_soil%var_2d(nzb_soil+1,m) &
     4661                      - surf_t_soil%var_2d(nzb_soil,m) ) * ddz_soil_center(nzb_soil) &
     4662                      + surf%ghf(m) ) * ddz_soil(nzb_soil)
    32744663
    32754664             DO  k = nzb_soil+1, nzt_soil
    3276                 tend(k) = ( 1.0_wp / surf%rho_c_total(k,m) )                      &
    3277                         * (   surf%lambda_h(k,m)                                  &
    3278                         * ( surf_t_soil%var_2d(k+1,m) - surf_t_soil%var_2d(k,m) ) &
    3279                         * ddz_soil_center(k) - surf%lambda_h(k-1,m)               &
    3280                         * ( surf_t_soil%var_2d(k,m) - surf_t_soil%var_2d(k-1,m) ) &
    3281                         * ddz_soil_center(k-1) ) * ddz_soil(k)
     4665                tend(k) = ( 1.0_wp / surf%rho_c_total(k,m) )                   &
     4666                          * (   surf%lambda_h(k,m)                             &
     4667                     * ( surf_t_soil%var_2d(k+1,m) - surf_t_soil%var_2d(k,m) ) &
     4668                     * ddz_soil_center(k)                                      &
     4669                     - surf%lambda_h(k-1,m)                                    &
     4670                     * ( surf_t_soil%var_2d(k,m) - surf_t_soil%var_2d(k-1,m) ) &
     4671                     * ddz_soil_center(k-1)                                    &
     4672                            ) * ddz_soil(k)
    32824673
    32834674             ENDDO
    32844675
    32854676             surf_t_soil_p%var_2d(nzb_soil:nzt_soil,m) =                       &
    3286                                   surf_t_soil%var_2d(nzb_soil:nzt_soil,m)      &
    3287                                   + dt_3d * ( tsc(2)                           &
    3288                                   * tend(nzb_soil:nzt_soil) + tsc(3)           &
    3289                                   * surf_tt_soil_m%var_2d(nzb_soil:nzt_soil,m) )
     4677                                       surf_t_soil%var_2d(nzb_soil:nzt_soil,m) &
     4678                                               + dt_3d * ( tsc(2)              &
     4679                                               * tend(nzb_soil:nzt_soil)       &
     4680                                               + tsc(3)                        &
     4681                                               * surf_tt_soil_m%var_2d(nzb_soil:nzt_soil,m) )
    32904682
    32914683!
     
    32994691                         intermediate_timestep_count_max )  THEN
    33004692                   DO  k = nzb_soil, nzt_soil
    3301                       surf_tt_soil_m%var_2d(k,m) = -9.5625_wp * tend(k)        &
    3302                                                    + 5.3125_wp                 &
    3303                                                    * surf_tt_soil_m%var_2d(k,m)
     4693                      surf_tt_soil_m%var_2d(k,m) = -9.5625_wp * tend(k) +      &
     4694                                                    5.3125_wp *                &
     4695                                                     surf_tt_soil_m%var_2d(k,m)
    33044696                   ENDDO
    33054697                ENDIF
     
    33124704!--             In order to prevent water tranport through paved surfaces,
    33134705!--             conductivity and diffusivity are set to zero
    3314                 IF ( surf%pavement_surface(m)  .AND.  k <= nzt_pavement )  THEN
    3315 
     4706                IF ( surf%pavement_surface(m)  .AND.                           &
     4707                     k <= surf%nzt_pavement(m) )  THEN
    33164708                   lambda_temp(k) = 0.0_wp
    33174709                   gamma_temp(k)  = 0.0_wp
    3318 
     4710   
    33194711                ELSE
    3320 
     4712   
    33214713!
    33224714!--                Calculate soil diffusivity at the center of the soil layers
     
    33324724                   h_vg = ( ( ( surf%m_res(k,m) - surf%m_sat(k,m) ) /          &
    33334725                              ( surf%m_res(k,m) -                              &
    3334                                 MAX(surf_m_soil%var_2d(k,m), surf%m_wilt(k,m)) &
     4726                                MAX( surf_m_soil%var_2d(k,m), surf%m_wilt(k,m) )&
    33354727                              )                                                &
    3336                             )**( surf%n_vg(k,m) / ( surf%n_vg(k,m) - 1.0_wp )  &
     4728                            )**(                                               &
     4729                          surf%n_vg(k,m) / ( surf%n_vg(k,m) - 1.0_wp )         &
    33374730                               ) - 1.0_wp                                      &
    33384731                          )**( 1.0_wp / surf%n_vg(k,m) ) / surf%alpha_vg(k,m)
     
    33404733                   gamma_temp(k) = surf%gamma_w_sat(k,m) * ( ( ( 1.0_wp +      &
    33414734                          ( surf%alpha_vg(k,m) * h_vg )**surf%n_vg(k,m)        &
    3342                           )**( 1.0_wp - 1.0_wp / surf%n_vg(k,m) )              &
    3343                           - ( surf%alpha_vg(k,m) * h_vg )**( surf%n_vg(k,m)    &
    3344                           - 1.0_wp) )**2 ) / ( ( 1.0_wp + ( surf%alpha_vg(k,m) &
    3345                           * h_vg )**surf%n_vg(k,m) )**(( 1.0_wp  - 1.0_wp      &
    3346                               / surf%n_vg(k,m) ) * ( surf%l_vg(k,m) + 2.0_wp) ))
     4735                                                                  )**(         &
     4736                              1.0_wp - 1.0_wp / surf%n_vg(k,m)) - (            &
     4737                          surf%alpha_vg(k,m) * h_vg )**( surf%n_vg(k,m)        &
     4738                              - 1.0_wp) )**2 )                                 &
     4739                              / ( ( 1.0_wp + ( surf%alpha_vg(k,m) * h_vg       &
     4740                              )**surf%n_vg(k,m) )**( ( 1.0_wp  - 1.0_wp        &
     4741                              / surf%n_vg(k,m) ) *                             &
     4742                              ( surf%l_vg(k,m) + 2.0_wp) ) )
     4743
    33474744                ENDIF
    33484745
     
    33814778!--             would cause accumulation of (non-existing) water in the lowest
    33824779!--             soil layer
    3383                 IF ( conserve_water_content .AND. surf_m_soil%var_2d(nzt_soil,m) /= 0.0_wp )  THEN
     4780                IF ( conserve_water_content .AND.                              &
     4781                     surf_m_soil%var_2d(nzt_soil,m) /= 0.0_wp )  THEN
     4782
    33844783                   surf%gamma_w(nzt_soil,m) = 0.0_wp
    33854784                ELSE
     
    34044803                    IF ( surf_m_soil%var_2d(k,m) > surf%m_wilt(k,m) )  THEN
    34054804                       m_total = m_total + surf%root_fr(k,m)                   &
    3406                                  * surf_m_soil%var_2d(k,m)
     4805                                         * surf_m_soil%var_2d(k,m)
    34074806                    ENDIF
    34084807                ENDDO 
     
    34114810                      IF ( surf_m_soil%var_2d(k,m) > surf%m_wilt(k,m) )  THEN
    34124811                         root_extr(k) = surf%root_fr(k,m)                      &
    3413                                         * surf_m_soil%var_2d(k,m) / m_total
     4812                                      * surf_m_soil%var_2d(k,m) / m_total
    34144813                      ELSE
    34154814                         root_extr(k) = 0.0_wp
     
    34284827                            + surf%qsws_soil(m) ) * drho_l_lv )                &
    34294828                            * ddz_soil(nzb_soil)
     4829
    34304830
    34314831                DO  k = nzb_soil+1, nzt_soil-1
     
    34514851                surf_m_soil_p%var_2d(nzb_soil:nzt_soil,m) =                    &
    34524852                                       surf_m_soil%var_2d(nzb_soil:nzt_soil,m) &
    3453                                        + dt_3d * ( tsc(2) * tend(:)            &
    3454                                        + tsc(3) * surf_tm_soil_m%var_2d(:,m) )   
     4853                                         + dt_3d * ( tsc(2) * tend(:)          &
     4854                                         + tsc(3) * surf_tm_soil_m%var_2d(:,m) )   
     4855   
    34554856!
    34564857!--             Account for dry soils (find a better solution here!)
     
    34734874                                                    * surf_tm_soil_m%var_2d(k,m)
    34744875                      ENDDO
     4876
    34754877                   ENDIF
    34764878                ENDIF
     
    35264928                m_liq_h   => m_liq_h_1;     m_liq_h_p      => m_liq_h_2
    35274929             ENDIF
     4930
    35284931!
    35294932!--          Vertical surfaces
     
    35334936                m_soil_v  => m_soil_v_1;    m_soil_v_p     => m_soil_v_2
    35344937                m_liq_v   => m_liq_v_1;     m_liq_v_p      => m_liq_v_2
     4938
    35354939             ENDIF
    35364940
     
    35454949                m_soil_h  => m_soil_h_2;    m_soil_h_p     => m_soil_h_1
    35464950                m_liq_h   => m_liq_h_2;     m_liq_h_p      => m_liq_h_1
     4951
    35474952             ENDIF
    35484953!
     
    38915296                ENDDO
    38925297             ENDDO
     5298
    38935299!
    38945300!--
     
    43275733!> Write restart data for land surface model
    43285734!------------------------------------------------------------------------------!
    4329  SUBROUTINE lsm_last_actions
     5735 SUBROUTINE lsm_write_restart_data
    43305736 
    43315737
     
    43355741
    43365742    IMPLICIT NONE
     5743
     5744    CHARACTER (LEN=1) ::  dum    !< dummy to create correct string for creating variable string
     5745    INTEGER(iwp)      ::  l      !< index variable for surface orientation
    43375746
    43385747    IF ( write_binary )  THEN
     
    43525761          WRITE ( 14 )  'lai_av              ';  WRITE ( 14 )  lai_av
    43535762       ENDIF
    4354        WRITE ( 14 )     'm_liq               ';  WRITE ( 14 )  m_liq_h%var_1d
    43555763       IF ( ALLOCATED( m_liq_av ) )  THEN
    43565764          WRITE ( 14 )  'm_liq_av            ';  WRITE ( 14 )  m_liq_av
    43575765       ENDIF
    4358        WRITE ( 14 )     'm_soil              ';  WRITE ( 14 )  m_soil_h%var_2d
    43595766       IF ( ALLOCATED( m_soil_av ) )  THEN
    43605767          WRITE ( 14 )  'm_soil_av           ';  WRITE ( 14 )  m_soil_av
     
    43695776          WRITE ( 14 )  'qsws_veg_av         ';  WRITE ( 14 )  qsws_veg_av
    43705777       ENDIF
    4371        WRITE ( 14 )     't_soil              ';  WRITE ( 14 )  t_soil_h%var_2d
    43725778       IF ( ALLOCATED( t_soil_av ) )  THEN
    43735779          WRITE ( 14 )  't_soil_av           ';  WRITE ( 14 )  t_soil_av
    43745780       ENDIF
    43755781
     5782
     5783       WRITE ( 14 ) 'ns_h_on_file_lsm    '
     5784       WRITE ( 14 ) surf_lsm_h%ns
     5785       WRITE ( 14 ) 'ns_v_on_file_lsm    '
     5786       WRITE ( 14 ) surf_lsm_v(0:3)%ns
     5787       
     5788       WRITE ( 14 ) 'lsm_start_index_h   '
     5789       WRITE ( 14 ) surf_lsm_h%start_index
     5790       WRITE ( 14 ) 'lsm_end_index_h     '
     5791       WRITE ( 14 ) surf_lsm_h%end_index
     5792       WRITE ( 14 ) 't_soil_h            '
     5793       WRITE ( 14 ) t_soil_h%var_2d
     5794       
     5795       DO  l = 0, 3
     5796          WRITE ( 14 ) 'lsm_start_index_v   '
     5797          WRITE ( 14 ) surf_lsm_v(l)%start_index
     5798          WRITE ( 14 ) 'lsm_end_index_v     '
     5799          WRITE ( 14 ) surf_lsm_v(l)%end_index
     5800          WRITE( dum, '(I1)')  l         
     5801          WRITE ( 14 ) 't_soil_v(' // dum // ')         '
     5802          WRITE ( 14 ) t_soil_v(l)%var_2d         
     5803       ENDDO
     5804
     5805       WRITE ( 14 ) 'lsm_start_index_h   '
     5806       WRITE ( 14 ) surf_lsm_h%start_index
     5807       WRITE ( 14 ) 'lsm_end_index_h     '
     5808       WRITE ( 14 ) surf_lsm_h%end_index
     5809       WRITE ( 14 ) 'm_soil_h            '
     5810       WRITE ( 14 ) m_soil_h%var_2d
     5811       
     5812       DO  l = 0, 3
     5813          WRITE ( 14 ) 'lsm_start_index_v   '
     5814          WRITE ( 14 ) surf_lsm_v(l)%start_index
     5815          WRITE ( 14 ) 'lsm_end_index_v     '
     5816          WRITE ( 14 ) surf_lsm_v(l)%end_index
     5817          WRITE( dum, '(I1)')  l         
     5818          WRITE ( 14 ) 'm_soil_v(' // dum // ')         '
     5819          WRITE ( 14 ) m_soil_v(l)%var_2d         
     5820       ENDDO
     5821
     5822
     5823       WRITE ( 14 ) 'lsm_start_index_h   '
     5824       WRITE ( 14 ) surf_lsm_h%start_index
     5825       WRITE ( 14 ) 'lsm_end_index_h     '
     5826       WRITE ( 14 ) surf_lsm_h%end_index
     5827       WRITE ( 14 ) 'm_liq_h             '
     5828       WRITE ( 14 ) m_liq_h%var_1d
     5829       
     5830       DO  l = 0, 3
     5831          WRITE ( 14 ) 'lsm_start_index_v   '
     5832          WRITE ( 14 ) surf_lsm_v(l)%start_index
     5833          WRITE ( 14 ) 'lsm_end_index_v     '
     5834          WRITE ( 14 ) surf_lsm_v(l)%end_index
     5835          WRITE( dum, '(I1)')  l         
     5836          WRITE ( 14 ) 'm_liq_v(' // dum // ')          '
     5837          WRITE ( 14 ) m_liq_v(l)%var_1d         
     5838       ENDDO
     5839
    43765840       WRITE ( 14 )  '*** end lsm ***     '
    43775841
    43785842    ENDIF
    43795843
    4380  END SUBROUTINE lsm_last_actions
     5844 END SUBROUTINE lsm_write_restart_data
    43815845
    43825846
     
    43995863    CHARACTER (LEN=20) :: field_char   !<
    44005864
    4401     INTEGER(iwp) ::  i               !<
    4402     INTEGER(iwp) ::  k               !<
    4403     INTEGER(iwp) ::  nxlc            !<
    4404     INTEGER(iwp) ::  nxlf            !<
    4405     INTEGER(iwp) ::  nxl_on_file     !<
    4406     INTEGER(iwp) ::  nxrc            !<
    4407     INTEGER(iwp) ::  nxrf            !<
    4408     INTEGER(iwp) ::  nxr_on_file     !<
    4409     INTEGER(iwp) ::  nync            !<
    4410     INTEGER(iwp) ::  nynf            !<
    4411     INTEGER(iwp) ::  nyn_on_file     !<
    4412     INTEGER(iwp) ::  nysc            !<
    4413     INTEGER(iwp) ::  nysf            !<
    4414     INTEGER(iwp) ::  nys_on_file     !<
    4415     INTEGER(iwp) ::  overlap_count   !<
     5865    INTEGER(iwp) ::  i                 !<
     5866    INTEGER(iwp) ::  k                 !<
     5867    INTEGER(iwp) ::  l                 !< running index surface orientation
     5868    INTEGER(iwp) ::  ns_h_on_file_lsm  !< number of horizontal surface elements (natural type) on file
     5869    INTEGER(iwp) ::  nxlc              !<
     5870    INTEGER(iwp) ::  nxlf              !<
     5871    INTEGER(iwp) ::  nxl_on_file       !<
     5872    INTEGER(iwp) ::  nxrc              !<
     5873    INTEGER(iwp) ::  nxrf              !<
     5874    INTEGER(iwp) ::  nxr_on_file       !<
     5875    INTEGER(iwp) ::  nync              !<
     5876    INTEGER(iwp) ::  nynf              !<
     5877    INTEGER(iwp) ::  nyn_on_file       !<
     5878    INTEGER(iwp) ::  nysc              !<
     5879    INTEGER(iwp) ::  nysf              !<
     5880    INTEGER(iwp) ::  nys_on_file       !<
     5881    INTEGER(iwp) ::  overlap_count     !<
     5882
     5883    INTEGER(iwp) ::  ns_v_on_file_lsm(0:3) !< number of vertical surface elements (natural type) on file
    44165884
    44175885    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  nxlfa       !<
     
    44225890    INTEGER(iwp), DIMENSION(numprocs_previous_run,1000) ::  offset_ya   !<
    44235891
     5892    INTEGER(iwp), DIMENSION(nys_on_file:nyn_on_file,nxl_on_file:nxr_on_file) ::  start_index_on_file
     5893    INTEGER(iwp), DIMENSION(nys_on_file:nyn_on_file,nxl_on_file:nxr_on_file) ::  end_index_on_file
     5894
    44245895    REAL(wp),                                                                  &
    44255896       DIMENSION(nys_on_file-nbgp:nyn_on_file+nbgp,nxl_on_file-nbgp:nxr_on_file+nbgp) :: &
     
    44345905          tmp_3d2   !<
    44355906
    4436     REAL(wp),                                                                  &
    4437        DIMENSION(1:surf_lsm_h%ns) ::                                           &
    4438           tmp_walltype_1d   !<
    4439 
    4440     REAL(wp),                                                                  &
    4441        DIMENSION(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) ::                       &
    4442           tmp_walltype_2d   !<
    4443 
    4444     REAL(wp),                                                                  &
    4445        DIMENSION(nzb_soil:nzt_soil,1:surf_lsm_h%ns) ::                         &
    4446           tmp_walltype_2d2  !<
    4447 
     5907    TYPE(surf_type_lsm) :: tmp_walltype_h_1d   !< temporary 1D array containing the respective surface variable stored on file, horizontal surfaces
     5908    TYPE(surf_type_lsm) :: tmp_walltype_h_2d   !< temporary 2D array containing the respective surface variable stored on file, horizontal surfaces
     5909    TYPE(surf_type_lsm) :: tmp_walltype_h_2d2  !< temporary 2D array containing the respective surface variable stored on file, horizontal surfaces
     5910
     5911    TYPE(surf_type_lsm), DIMENSION(0:3) :: tmp_walltype_v_1d   !< temporary 1D array containing the respective surface variable stored on file, vertical surfaces
     5912    TYPE(surf_type_lsm), DIMENSION(0:3) :: tmp_walltype_v_2d   !< temporary 2D array containing the respective surface variable stored on file, vertical surfaces
     5913    TYPE(surf_type_lsm), DIMENSION(0:3) :: tmp_walltype_v_2d2  !< temporary 2D array containing the respective surface variable stored on file, vertical surfaces
    44485914
    44495915   IF ( initializing_actions == 'read_restart_data' )  THEN
     5916      READ ( 13 )  field_char
     5917
     5918!
     5919!--   At first, determine the number of surface elements (with certain orientation) on file
     5920      IF ( TRIM( field_char ) /= 'ns_h_on_file_lsm' )  THEN
     5921!
     5922!--      Add a proper error message
     5923      ENDIF
     5924      READ ( 13 ) ns_h_on_file_lsm
     5925
     5926      READ ( 13 )  field_char
     5927      IF ( TRIM( field_char ) /= 'ns_v_on_file_lsm' )  THEN
     5928!
     5929!--      Add a proper error message
     5930      ENDIF
     5931      READ ( 13 ) ns_v_on_file_lsm
     5932
     5933!
     5934!--   Allocate temporary arrays to store surface data
     5935      ALLOCATE( tmp_walltype_h_1d%var_1d(1:ns_h_on_file_lsm)                     )
     5936      ALLOCATE( tmp_walltype_h_2d%var_2d(nzb_soil:nzt_soil+1,1:ns_h_on_file_lsm) )
     5937      ALLOCATE( tmp_walltype_h_2d2%var_2d(nzb_soil:nzt_soil,1:ns_h_on_file_lsm)  )
     5938
     5939      DO  l = 0, 3
     5940         ALLOCATE( tmp_walltype_v_1d(l)%var_1d(1:ns_v_on_file_lsm(l))                     )
     5941         ALLOCATE( tmp_walltype_v_2d(l)%var_2d(nzb_soil:nzt_soil+1,1:ns_v_on_file_lsm(l)) )
     5942         ALLOCATE( tmp_walltype_v_2d2(l)%var_2d(nzb_soil:nzt_soil,1:ns_v_on_file_lsm(l))  )
     5943      ENDDO
     5944     
    44505945      READ ( 13 )  field_char
    44515946
     
    44635958            nync = nynfa(i,k) + offset_ya(i,k)
    44645959
     5960
    44655961            SELECT CASE ( TRIM( field_char ) )
    44665962
     
    44955991                   ENDIF
    44965992                   IF ( k == 1 )  READ ( 13 )  tmp_2d
    4497                    ghf_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =        &
     5993                   ghf_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =           &
    44985994                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    4499 
    4500                 CASE ( 'm_liq' )
    4501                    IF ( k == 1 )  READ ( 13 )  tmp_walltype_1d !tmp_2d
    4502                    m_liq_h%var_1d(1:surf_lsm_h%ns)  =                       &
    4503                                  tmp_walltype_1d(1:surf_lsm_h%ns)
    45045995
    45055996                CASE ( 'lai_av' )
     
    45186009                   m_liq_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =         &
    45196010                                  tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    4520 
    4521                 CASE ( 'm_soil' )
    4522                    IF ( k == 1 )  READ ( 13 )  tmp_walltype_2d2(:,:)
    4523                    m_soil_h%var_2d(:,1:surf_lsm_h%ns) =                        &
    4524                           tmp_walltype_2d2(:,1:surf_lsm_h%ns)   
    45256011
    45266012                CASE ( 'm_soil_av' )
     
    45566042                                          tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    45576043
    4558                 CASE ( 't_soil' )
    4559                    IF ( k == 1 )  READ ( 13 )  tmp_walltype_2d(:,:)
    4560                    t_soil_h%var_2d(:,1:surf_lsm_h%ns) =                        &
    4561                          tmp_walltype_2d(:,1:surf_lsm_h%ns)     
    4562 
    45636044                CASE ( 't_soil_av' )
    45646045                   IF ( .NOT. ALLOCATED( t_soil_av ) )  THEN
     
    45696050                                    tmp_3d2(:,nysf-nbgp:nynf+nbgp,             &
    45706051                                    nxlf-nbgp:nxrf+nbgp)
     6052
     6053                CASE ( 'lsm_start_index_h', 'lsm_start_index_v'  )   
     6054                   IF ( k == 1 )                                               &
     6055                      READ ( 13 )  start_index_on_file
     6056                     
     6057                CASE ( 'lsm_end_index_h', 'lsm_end_index_v' )   
     6058                   IF ( k == 1 )                                               &
     6059                      READ ( 13 )  end_index_on_file
     6060               
     6061                CASE ( 't_soil_h' )
     6062                 
     6063                   IF ( k == 1 )  THEN
     6064                      IF ( .NOT.  ALLOCATED( t_soil_h%var_2d ) )               &
     6065                         ALLOCATE( t_soil_h%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) )
     6066                      READ ( 13 )  tmp_walltype_h_2d%var_2d
     6067                   ENDIF
     6068                   CALL surface_restore_elements(                              &
     6069                                              t_soil_h%var_2d,                 &
     6070                                              tmp_walltype_h_2d%var_2d,        &
     6071                                              surf_lsm_h%start_index,          &
     6072                                              start_index_on_file,             &
     6073                                              end_index_on_file,               &
     6074                                              nxlc, nysc,                      &
     6075                                              nxlf, nxrf, nysf, nynf,          &
     6076                                              nys_on_file, nyn_on_file,        &
     6077                                              nxl_on_file,nxr_on_file )
     6078
     6079                CASE ( 't_soil_v(0)' )
     6080                 
     6081                   IF ( k == 1 )  THEN
     6082                      IF ( .NOT.  ALLOCATED( t_soil_v(0)%var_2d ) )            &
     6083                         ALLOCATE( t_soil_v(0)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(0)%ns) )
     6084                      READ ( 13 )  tmp_walltype_v_2d(0)%var_2d
     6085                   ENDIF
     6086                   CALL surface_restore_elements(                              &
     6087                                           t_soil_v(0)%var_2d,                 &
     6088                                           tmp_walltype_v_2d(0)%var_2d,        &
     6089                                           surf_lsm_v(0)%start_index,          &
     6090                                           start_index_on_file,                &
     6091                                           end_index_on_file,                  &
     6092                                           nxlc, nysc,                         &
     6093                                           nxlf, nxrf, nysf, nynf,             &
     6094                                           nys_on_file, nyn_on_file,           &
     6095                                           nxl_on_file,nxr_on_file )
     6096
     6097                CASE ( 't_soil_v(1)' )
     6098                 
     6099                   IF ( k == 1 )  THEN
     6100                      IF ( .NOT.  ALLOCATED( t_soil_v(1)%var_2d ) )            &
     6101                         ALLOCATE( t_soil_v(1)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(1)%ns) )
     6102                      READ ( 13 )  tmp_walltype_v_2d(1)%var_2d
     6103                   ENDIF
     6104                   CALL surface_restore_elements(                              &
     6105                                           t_soil_v(1)%var_2d,                 &
     6106                                           tmp_walltype_v_2d(1)%var_2d,        &
     6107                                           surf_lsm_v(1)%start_index,          &
     6108                                           start_index_on_file,                &
     6109                                           end_index_on_file,                  &
     6110                                           nxlc, nysc,                         &
     6111                                           nxlf, nxrf, nysf, nynf,             &
     6112                                           nys_on_file, nyn_on_file,           &
     6113                                           nxl_on_file,nxr_on_file )
     6114
     6115                CASE ( 't_soil_v(2)' )
     6116                 
     6117                   IF ( k == 1 )  THEN
     6118                      IF ( .NOT.  ALLOCATED( t_soil_v(2)%var_2d ) )            &
     6119                         ALLOCATE( t_soil_v(2)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(2)%ns) )
     6120                      READ ( 13 )  tmp_walltype_v_2d(2)%var_2d
     6121                   ENDIF
     6122                   CALL surface_restore_elements(                              &
     6123                                           t_soil_v(2)%var_2d,                 &
     6124                                           tmp_walltype_v_2d(2)%var_2d,        &
     6125                                           surf_lsm_v(2)%start_index,          &
     6126                                           start_index_on_file,                &
     6127                                           end_index_on_file,                  &
     6128                                           nxlc, nysc,                         &
     6129                                           nxlf, nxrf, nysf, nynf,             &
     6130                                           nys_on_file, nyn_on_file,           &
     6131                                           nxl_on_file,nxr_on_file )
     6132
     6133                CASE ( 't_soil_v(3)' )
     6134                 
     6135                   IF ( k == 1 )  THEN
     6136                      IF ( .NOT.  ALLOCATED( t_soil_v(3)%var_2d ) )            &
     6137                         ALLOCATE( t_soil_v(1)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(3)%ns) )
     6138                      READ ( 13 )  tmp_walltype_v_2d(3)%var_2d
     6139                   ENDIF
     6140                   CALL surface_restore_elements(                              &
     6141                                           t_soil_v(3)%var_2d,                 &
     6142                                           tmp_walltype_v_2d(3)%var_2d,        &
     6143                                           surf_lsm_v(3)%start_index,          &
     6144                                           start_index_on_file,                &
     6145                                           end_index_on_file,                  &
     6146                                           nxlc, nysc,                         &
     6147                                           nxlf, nxrf, nysf, nynf,             &
     6148                                           nys_on_file, nyn_on_file,           &
     6149                                           nxl_on_file,nxr_on_file )
     6150
     6151                CASE ( 'm_soil_h' )
     6152                 
     6153                   IF ( k == 1 )  THEN
     6154                      IF ( .NOT.  ALLOCATED( m_soil_h%var_2d ) )               &
     6155                         ALLOCATE( m_soil_h%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_h%ns) )
     6156                      READ ( 13 )  tmp_walltype_h_2d2%var_2d
     6157                   ENDIF
     6158                   CALL surface_restore_elements(                              &
     6159                                             m_soil_h%var_2d,                  &
     6160                                             tmp_walltype_h_2d2%var_2d,        &
     6161                                             surf_lsm_h%start_index,           &
     6162                                             start_index_on_file,              &
     6163                                             end_index_on_file,                &
     6164                                             nxlc, nysc,                       &
     6165                                             nxlf, nxrf, nysf, nynf,           &
     6166                                             nys_on_file, nyn_on_file,         &
     6167                                             nxl_on_file,nxr_on_file )
     6168
     6169                CASE ( 'm_soil_v(0)' )
     6170                 
     6171                   IF ( k == 1 )  THEN
     6172                      IF ( .NOT.  ALLOCATED( m_soil_v(0)%var_2d ) )            &
     6173                         ALLOCATE( m_soil_v(0)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(0)%ns) )
     6174                      READ ( 13 )  tmp_walltype_v_2d2(0)%var_2d
     6175                   ENDIF
     6176                   CALL surface_restore_elements(                              &
     6177                                          m_soil_v(0)%var_2d,                  &
     6178                                          tmp_walltype_v_2d2(0)%var_2d,        &
     6179                                          surf_lsm_v(0)%start_index,           &
     6180                                          start_index_on_file,                 &
     6181                                          end_index_on_file,                   &
     6182                                          nxlc, nysc,                          &
     6183                                          nxlf, nxrf, nysf, nynf,              &
     6184                                          nys_on_file, nyn_on_file,            &
     6185                                          nxl_on_file,nxr_on_file )
     6186
     6187                CASE ( 'm_soil_v(1)' )
     6188                 
     6189                   IF ( k == 1 )  THEN
     6190                      IF ( .NOT.  ALLOCATED( m_soil_v(1)%var_2d ) )            &
     6191                         ALLOCATE( m_soil_v(1)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(1)%ns) )
     6192                      READ ( 13 )  tmp_walltype_v_2d2(1)%var_2d
     6193                   ENDIF
     6194                   CALL surface_restore_elements(                              &
     6195                                          m_soil_v(1)%var_2d,                  &
     6196                                          tmp_walltype_v_2d2(1)%var_2d,        &
     6197                                          surf_lsm_v(1)%start_index,           &
     6198                                          start_index_on_file,                 &
     6199                                          end_index_on_file,                   &
     6200                                          nxlc, nysc,                          &
     6201                                          nxlf, nxrf, nysf, nynf,              &
     6202                                          nys_on_file, nyn_on_file,            &
     6203                                          nxl_on_file,nxr_on_file )
     6204
     6205
     6206                CASE ( 'm_soil_v(2)' )
     6207                 
     6208                   IF ( k == 1 )  THEN
     6209                      IF ( .NOT.  ALLOCATED( m_soil_v(2)%var_2d ) )            &
     6210                         ALLOCATE( m_soil_v(2)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(2)%ns) )
     6211                      READ ( 13 )  tmp_walltype_v_2d2(2)%var_2d
     6212                   ENDIF
     6213                   CALL surface_restore_elements(                              &
     6214                                          m_soil_v(2)%var_2d,                  &
     6215                                          tmp_walltype_v_2d2(2)%var_2d,        &
     6216                                          surf_lsm_v(2)%start_index,           &
     6217                                          start_index_on_file,                 &
     6218                                          end_index_on_file,                   &
     6219                                          nxlc, nysc,                          &
     6220                                          nxlf, nxrf, nysf, nynf,              &
     6221                                          nys_on_file, nyn_on_file,            &
     6222                                          nxl_on_file,nxr_on_file )
     6223
     6224
     6225                CASE ( 'm_soil_v(3)' )
     6226                 
     6227                   IF ( k == 1 )  THEN
     6228                      IF ( .NOT.  ALLOCATED( m_soil_v(3)%var_2d ) )            &
     6229                         ALLOCATE( m_soil_v(1)%var_2d(nzb_soil:nzt_soil+1,1:surf_lsm_v(3)%ns) )
     6230                      READ ( 13 )  tmp_walltype_v_2d2(3)%var_2d
     6231                   ENDIF
     6232                   CALL surface_restore_elements(                              &
     6233                                          m_soil_v(3)%var_2d,                  &
     6234                                          tmp_walltype_v_2d2(3)%var_2d,        &
     6235                                          surf_lsm_v(3)%start_index,           &
     6236                                          start_index_on_file,                 &
     6237                                          end_index_on_file,                   &
     6238                                          nxlc, nysc,                          &
     6239                                          nxlf, nxrf, nysf, nynf,              &
     6240                                          nys_on_file, nyn_on_file,            &
     6241                                          nxl_on_file,nxr_on_file )
     6242
     6243
     6244                CASE ( 'm_liq_h' )
     6245                 
     6246                   IF ( k == 1 )  THEN
     6247                      IF ( .NOT.  ALLOCATED( m_liq_h%var_1d ) )                &
     6248                         ALLOCATE( m_liq_h%var_1d(1:surf_lsm_h%ns) )
     6249                      READ ( 13 )  tmp_walltype_h_1d%var_1d
     6250                   ENDIF
     6251                   CALL surface_restore_elements(                              &
     6252                                              m_liq_h%var_1d,                  &
     6253                                              tmp_walltype_h_1d%var_1d,        &
     6254                                              surf_lsm_h%start_index,          &
     6255                                              start_index_on_file,             &
     6256                                              end_index_on_file,               &
     6257                                              nxlc, nysc,                      &
     6258                                              nxlf, nxrf, nysf, nynf,          &
     6259                                              nys_on_file, nyn_on_file,        &
     6260                                              nxl_on_file,nxr_on_file )
     6261
     6262
     6263                CASE ( 'm_liq_v(0)' )
     6264                 
     6265                   IF ( k == 1 )  THEN
     6266                      IF ( .NOT.  ALLOCATED( m_liq_v(0)%var_1d ) )             &
     6267                         ALLOCATE( m_liq_v(0)%var_1d(1:surf_lsm_v(0)%ns) )
     6268                      READ ( 13 )  tmp_walltype_v_1d(0)%var_1d
     6269                   ENDIF
     6270                   CALL surface_restore_elements(                              &
     6271                                           m_liq_v(0)%var_1d,                  &
     6272                                           tmp_walltype_v_1d(0)%var_1d,        &
     6273                                           surf_lsm_v(0)%start_index,          &
     6274                                           start_index_on_file,                &
     6275                                           end_index_on_file,                  &
     6276                                           nxlc, nysc,                         &
     6277                                           nxlf, nxrf, nysf, nynf,             &
     6278                                           nys_on_file, nyn_on_file,           &
     6279                                           nxl_on_file,nxr_on_file )
     6280
     6281
     6282                CASE ( 'm_liq_v(1)' )
     6283                 
     6284                   IF ( k == 1 )  THEN
     6285                      IF ( .NOT.  ALLOCATED( m_liq_v(1)%var_1d ) )             &
     6286                         ALLOCATE( m_liq_v(1)%var_1d(1:surf_lsm_v(1)%ns) )
     6287                      READ ( 13 )  tmp_walltype_v_1d(1)%var_1d
     6288                   ENDIF
     6289                   CALL surface_restore_elements(                              &
     6290                                           m_liq_v(1)%var_1d,                  &
     6291                                           tmp_walltype_v_1d(1)%var_1d,        &
     6292                                           surf_lsm_v(1)%start_index,          &
     6293                                           start_index_on_file,                &
     6294                                           end_index_on_file,                  &
     6295                                           nxlc, nysc,                         &
     6296                                           nxlf, nxrf, nysf, nynf,             &
     6297                                           nys_on_file, nyn_on_file,           &
     6298                                           nxl_on_file,nxr_on_file )
     6299
     6300
     6301                CASE ( 'm_liq_v(2)' )
     6302                 
     6303                   IF ( k == 1 )  THEN
     6304                      IF ( .NOT.  ALLOCATED( m_liq_v(2)%var_1d ) )             &
     6305                         ALLOCATE( m_liq_v(2)%var_1d(1:surf_lsm_v(2)%ns) )
     6306                      READ ( 13 )  tmp_walltype_v_1d(2)%var_1d
     6307                   ENDIF
     6308                   CALL surface_restore_elements(                              &
     6309                                           m_liq_v(2)%var_1d,                  &
     6310                                           tmp_walltype_v_1d(2)%var_1d,        &
     6311                                           surf_lsm_v(2)%start_index,          &
     6312                                           start_index_on_file,                &
     6313                                           end_index_on_file,                  &
     6314                                           nxlc, nysc,                         &
     6315                                           nxlf, nxrf, nysf, nynf,             &
     6316                                           nys_on_file, nyn_on_file,           &
     6317                                           nxl_on_file,nxr_on_file )
     6318
     6319                CASE ( 'm_liq_v(3)' )
     6320                 
     6321                   IF ( k == 1 )  THEN
     6322                      IF ( .NOT.  ALLOCATED( m_liq_v(3)%var_1d ) )             &
     6323                         ALLOCATE( m_liq_v(3)%var_1d(1:surf_lsm_v(3)%ns) )
     6324                      READ ( 13 )  tmp_walltype_v_1d(3)%var_1d
     6325                   ENDIF
     6326                   CALL surface_restore_elements(                              &
     6327                                           m_liq_v(3)%var_1d,                  &
     6328                                           tmp_walltype_v_1d(3)%var_1d,        &
     6329                                           surf_lsm_v(3)%start_index,          &
     6330                                           start_index_on_file,                &
     6331                                           end_index_on_file,                  &
     6332                                           nxlc, nysc,                         &
     6333                                           nxlf, nxrf, nysf, nynf,             &
     6334                                           nys_on_file, nyn_on_file,           &
     6335                                           nxl_on_file,nxr_on_file )
    45716336
    45726337
     
    46456410!--          Charnock relation is used.
    46466411             surf_lsm_h%z0(m)  = ( 0.11_wp * molecular_viscosity /             &
    4647                                  surf_lsm_h%us(m) )  &
     6412                                 surf_lsm_h%us(m) )                            &
    46486413                               + ( alpha_ch * surf_lsm_h%us(m)**2 / g )
    46496414
Note: See TracChangeset for help on using the changeset viewer.