Index: /palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2327)
+++ /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 2328)
@@ -25,4 +25,12 @@
 ! -----------------
 ! $Id$
+! Revised skin layer concept.
+! Bugfix for runs with pavement surface and humidity
+! Revised some standard values in vegetation_pars
+! Added emissivity and default albedo_type as variable to tables
+! Changed default surface type to vegetation
+! Revised input of soil layer configuration
+! 
+! 2307 2017-07-07 11:32:10Z suehring
 ! Bugfix, variable names corrected
 ! 
@@ -257,6 +265,6 @@
 
     USE radiation_model_mod,                                                   &
-        ONLY:  force_radiation_call, rad_net, rad_sw_in, rad_lw_out,           &
-               rad_lw_out_change_0, radiation_scheme,                          &
+        ONLY:  albedo, albedo_type, emissivity, force_radiation_call, rad_net,         &
+               rad_sw_in, rad_lw_out, rad_lw_out_change_0, radiation_scheme,   &
                unscheduled_radiation_calls
         
@@ -288,15 +296,15 @@
 
 
-    REAL(wp), DIMENSION(0:7), PARAMETER  :: zs_default =                       & ! default soil layer configuration
-                                            (/ 0.005_wp, 0.02_wp, 0.04_wp,     &
-                                               0.07_wp, 0.15_wp, 0.28_wp,      &
-                                               1.00_wp, 2.89_wp /)
+    REAL(wp), DIMENSION(0:7), PARAMETER  :: dz_soil_default =                  & ! default soil layer configuration
+                                            (/ 0.01_wp, 0.02_wp, 0.04_wp,      &
+                                               0.07_wp, 0.15_wp, 0.21_wp,      &
+                                               0.72_wp, 1.89_wp/)
 
 
 !
 !-- LSM variables
-    CHARACTER(10) :: surface_type = 'netcdf'  !< general classification. Allowed are:
-                                             !< 'vegetation', 'pavement', ('building'), 
-                                             !< 'water', and 'netcdf'
+    CHARACTER(10) :: surface_type = 'vegetation'  !< general classification. Allowed are:
+                                                  !< 'vegetation', 'pavement', ('building'), 
+                                                  !< 'water', and 'netcdf'
 
 
@@ -362,9 +370,9 @@
                 
                 
-    REAL(wp), DIMENSION(:), ALLOCATABLE  :: ddz_soil,         & !< 1/dz_soil
-                                            ddz_soil_layer,   & !< 1/dz_soil_layer
-                                            dz_soil,          & !< soil grid spacing (center-center)
-                                            dz_soil_layer,    & !< soil grid spacing (edge-edge)
-                                            root_extr           !< root extraction
+    REAL(wp), DIMENSION(:), ALLOCATABLE  :: ddz_soil_center, & !< 1/dz_soil_center
+                                            ddz_soil,        & !< 1/dz_soil
+                                            dz_soil_center,  & !< soil grid spacing (center-center)
+                                            zs,              & !< depth of the temperature/moisute levels
+                                            root_extr          !< root extraction
 
 
@@ -373,5 +381,5 @@
                                    soil_moisture = 0.0_wp,                     & !< NAMELIST soil moisture content (m3/m3)
                                    soil_temperature = 300.0_wp,                & !< NAMELIST soil temperature (K) +1
-                                   zs = 9999999.9_wp,                          & !< (NAMELIST) soil layer depths (center)
+                                   dz_soil  = 9999999.9_wp,                    & !< (NAMELIST) soil layer depths (spacing)
                                    zs_layer = 9999999.9_wp                       !< soil layer depths (edge)
                                  
@@ -536,10 +544,11 @@
 !
 !-- Water classes
-    CHARACTER(12), DIMENSION(0:4), PARAMETER :: water_type_name = (/ &
+    CHARACTER(12), DIMENSION(0:5), PARAMETER :: water_type_name = (/ &
                                    'user defined',                   & ! 0 
                                    'lake        ',                   & ! 1
                                    'river       ',                   & ! 2
                                    'ocean       ',                   & ! 3
-                                   'pond        '                    & ! 4
+                                   'pond        ',                   & ! 4
+                                   'fountain    '                    & ! 5
                                                                   /)                                                                  
                                                                  
@@ -549,25 +558,25 @@
 !
 !-- Land surface parameters
-!-- r_canopy_min,     lai,   c_veg,     g_d         z0,         z0h, lambda_s_s, lambda_s_u, f_sw_in,  c_surface
-    REAL(wp), DIMENSION(0:9,1:18), PARAMETER :: vegetation_pars = RESHAPE( (/ &
-          0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, 0.5E-2_wp,   0.5E-4_wp,     0.0_wp,     0.0_wp, 0.00_wp, 0.00_wp, & !  1
-        180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,   0.25_wp,  0.25E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & !  2
-        110.0_wp, 2.00_wp, 0.85_wp, 0.00_wp,   0.20_wp,  0.20E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & !  3
-        500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,   2.00_wp,     2.00_wp,   280.0_wp,   196.0_wp, 0.03_wp, 0.00_wp, & !  4
-        500.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,   2.00_wp,     2.00_wp,   280.0_wp,   196.0_wp, 0.03_wp, 0.00_wp, & !  5
-        175.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,   2.00_wp,     2.00_wp,   280.0_wp,   196.0_wp, 0.03_wp, 0.00_wp, & !  6
-        240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp,   2.00_wp,     2.00_wp,   280.0_wp,   196.0_wp, 0.03_wp, 0.00_wp, & !  7
-        100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp,   0.47_wp,  0.47E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & !  8
-        250.0_wp, 0.05_wp, 0.00_wp, 0.00_wp,  0.013_wp, 0.013E-2_wp,   196.0_wp,   196.0_wp, 0.00_wp, 0.00_wp, & !  9
-         80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp,  0.034_wp, 0.034E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & ! 10
-        180.0_wp, 3.00_wp, 0.90_wp, 0.00_wp,    0.5_wp,  0.50E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & ! 11
-        150.0_wp, 0.50_wp, 0.10_wp, 0.00_wp,   0.17_wp,  0.17E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & ! 12
-          0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, 1.3E-3_wp,   1.3E-4_wp,   812.0_wp,   812.0_wp, 0.00_wp, 0.00_wp, & ! 13
-        240.0_wp, 4.00_wp, 0.60_wp, 0.00_wp,   0.83_wp,  0.83E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & ! 14
-        225.0_wp, 3.00_wp, 0.50_wp, 0.00_wp,   0.10_wp,  0.10E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & ! 15
-        225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp,   0.25_wp,  0.25E-2_wp,   140.0_wp,   140.0_wp, 0.05_wp, 0.00_wp, & ! 16
-        250.0_wp, 5.00_wp, 0.90_wp, 0.03_wp,   2.00_wp,     2.00_wp,   280.0_wp,   196.0_wp, 0.03_wp, 0.00_wp, & ! 17
-        175.0_wp, 2.50_wp, 0.90_wp, 0.03_wp,   1.10_wp,     1.10_wp,   280.0_wp,   196.0_wp, 0.03_wp, 0.00_wp  & ! 18
-                                                               /), (/ 10, 18 /) )
+!-- r_canopy_min,     lai,   c_veg,     g_d         z0,         z0h, lambda_s_s, lambda_s_u, f_sw_in,  c_surface, albedo_type, emissivity
+    REAL(wp), DIMENSION(0:11,1:18), PARAMETER :: vegetation_pars = RESHAPE( (/ &
+          0.0_wp, 0.00_wp, 1.00_wp, 0.00_wp,  0.005_wp,   0.5E-4_wp,     0.0_wp,    0.0_wp, 0.00_wp, 0.00_wp,  0.0_wp, 0.94_wp, & !  1
+        180.0_wp, 3.00_wp, 1.00_wp, 0.00_wp,   0.10_wp,    0.001_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp,  2.0_wp, 0.95_wp, & !  2
+        110.0_wp, 2.00_wp, 1.00_wp, 0.00_wp,   0.03_wp,   0.3E-4_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp,  2.0_wp, 0.95_wp, & !  3
+        500.0_wp, 5.00_wp, 1.00_wp, 0.03_wp,   2.00_wp,     2.00_wp,    20.0_wp,   15.0_wp, 0.03_wp, 0.00_wp,  5.0_wp, 0.97_wp, & !  4
+        500.0_wp, 5.00_wp, 1.00_wp, 0.03_wp,   2.00_wp,     2.00_wp,    20.0_wp,   15.0_wp, 0.03_wp, 0.00_wp,  6.0_wp, 0.97_wp, & !  5
+        175.0_wp, 5.00_wp, 1.00_wp, 0.03_wp,   2.00_wp,     2.00_wp,    20.0_wp,   15.0_wp, 0.03_wp, 0.00_wp,  8.0_wp, 0.97_wp, & !  6
+        240.0_wp, 6.00_wp, 0.99_wp, 0.13_wp,   2.00_wp,     2.00_wp,    20.0_wp,   15.0_wp, 0.03_wp, 0.00_wp,  9.0_wp, 0.97_wp, & !  7
+        100.0_wp, 2.00_wp, 0.70_wp, 0.00_wp,   0.47_wp,  0.47E-2_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp,  8.0_wp, 0.97_wp, & !  8
+        250.0_wp, 0.05_wp, 0.00_wp, 0.00_wp,  0.013_wp, 0.013E-2_wp,    15.0_wp,   15.0_wp, 0.00_wp, 0.00_wp,  3.0_wp, 0.94_wp, & !  9
+         80.0_wp, 1.00_wp, 0.50_wp, 0.00_wp,  0.034_wp, 0.034E-2_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp, 11.0_wp, 0.97_wp, & ! 10
+        180.0_wp, 3.00_wp, 1.00_wp, 0.00_wp,    0.5_wp,  0.50E-2_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp, 13.0_wp, 0.97_wp, & ! 11
+        150.0_wp, 0.50_wp, 0.10_wp, 0.00_wp,   0.17_wp,  0.17E-2_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp,  2.0_wp, 0.97_wp, & ! 12
+          0.0_wp, 0.00_wp, 0.00_wp, 0.00_wp, 1.3E-3_wp,   1.3E-4_wp,    58.0_wp,   58.0_wp, 0.00_wp, 0.00_wp, 11.0_wp, 0.97_wp, & ! 13
+        240.0_wp, 4.00_wp, 0.60_wp, 0.00_wp,   0.83_wp,  0.83E-2_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp,  4.0_wp, 0.97_wp, & ! 14
+        225.0_wp, 3.00_wp, 0.50_wp, 0.00_wp,   0.10_wp,  0.10E-2_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp,  4.0_wp, 0.97_wp, & ! 15
+        225.0_wp, 1.50_wp, 0.50_wp, 0.00_wp,   0.25_wp,  0.25E-2_wp,    10.0_wp,   10.0_wp, 0.05_wp, 0.00_wp,  4.0_wp, 0.97_wp, & ! 16
+        250.0_wp, 5.00_wp, 1.00_wp, 0.03_wp,   2.00_wp,     2.00_wp,    20.0_wp,   15.0_wp, 0.03_wp, 0.00_wp,  7.0_wp, 0.97_wp, & ! 17
+        175.0_wp, 2.50_wp, 1.00_wp, 0.03_wp,   1.10_wp,     1.10_wp,    20.0_wp,   15.0_wp, 0.03_wp, 0.00_wp,  8.0_wp, 0.97_wp  & ! 18
+                                                               /), (/ 12, 18 /) )
 
                                     
@@ -631,24 +640,25 @@
 !
 !-- TO BE FILLED
-!-- Pavement parameters  depth,        z0,       z0h, lambda_h,      rho_c  
-    REAL(wp), DIMENSION(0:4,1:7), PARAMETER :: pavement_pars = RESHAPE( (/ &
-                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, & ! 1
-                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, & ! 2
-                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, & ! 3                                  
-                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, & ! 4
-                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, & ! 5
-                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, & ! 6
-                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp  & ! 7
-                                                                 /), (/ 5, 7 /) )                             
+!-- Pavement parameters  depth,        z0,       z0h, lambda_h,      rho_c, albedo_type, emissivity  
+    REAL(wp), DIMENSION(0:5,1:7), PARAMETER :: pavement_pars = RESHAPE( (/ &
+                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, 17.0_wp, 0.97_wp, & ! 1
+                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, 18.0_wp, 0.94_wp,  & ! 2
+                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, 19.0_wp, 0.98_wp,  & ! 3                                  
+                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, 20.0_wp, 0.93_wp,  & ! 4
+                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, 21.0_wp, 0.97_wp,  & ! 5
+                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, 22.0_wp, 0.97_wp,  & ! 6
+                      0.050_wp, 1.0E-4_wp, 1.0E-5_wp,  1.00_wp,  1.94E6_wp, 23.0_wp, 0.97_wp  & ! 7
+                                                                 /), (/ 6, 7 /) )                             
 
 !
 !-- TO BE FILLED
-!-- Water parameters                    temperature,     z0,      z0h
-    REAL(wp), DIMENSION(0:2,1:4), PARAMETER :: water_pars = RESHAPE( (/ &
-                                           283.0_wp, 0.01_wp, 0.001_wp, & ! 1
-                                           283.0_wp, 0.01_wp, 0.001_wp, & ! 2
-                                           283.0_wp, 0.01_wp, 0.001_wp, & ! 3
-                                           283.0_wp, 10.01_wp, 0.001_wp & ! 4
-                                                                /), (/ 3, 4 /) )                                                                     
+!-- Water parameters                    temperature,     z0,      z0h, albedo_type, emissivity,
+    REAL(wp), DIMENSION(0:4,1:5), PARAMETER :: water_pars = RESHAPE( (/ &
+                                           283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp, & ! 1
+                                           283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp, & ! 2
+                                           283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp, & ! 3
+                                           283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp, & ! 4
+                                           283.0_wp, 0.01_wp, 0.001_wp, 1.0_wp, 0.99_wp  & ! 5
+                                                                /), (/ 5, 5 /) )                                                                     
                                                                                                                                       
     SAVE
@@ -1099,6 +1109,4 @@
     ENDIF
     
-    
-     
     IF ( TRIM( surface_type ) == 'pavement' )  THEN
 
@@ -1144,4 +1152,12 @@
     ENDIF
     
+!
+!-- Temporary message as long as NetCDF input is not available
+    IF ( TRIM( surface_type ) == 'netcdf' )  THEN
+             message_string = 'surface_type = netcdf '//                       &
+                              'is not supported at the moment'
+             CALL message( 'check_parameters', 'PA0465', 1, 2, 0, 6, 0 )
+    ENDIF
+
     IF ( soil_type == 0 )  THEN
 
@@ -1209,13 +1225,13 @@
     nzb_soil = 0
     nzt_soil = -1
-    IF ( ALL( zs == 9999999.9_wp ) )  THEN
+    IF ( ALL( dz_soil == 9999999.9_wp ) )  THEN
        nzt_soil = 7
-       zs(nzb_soil:nzt_soil) = zs_default
+       dz_soil(nzb_soil:nzt_soil) = dz_soil_default
     ELSE
        DO k = 0, 19
-          IF ( zs(k) /= 9999999.9_wp )  THEN
+          IF ( dz_soil(k) /= 9999999.9_wp )  THEN
              nzt_soil = nzt_soil + 1
              IF ( root_fraction(k) == 9999999.9_wp )  THEN
-                message_string = 'manual setting of zs '//                     &
+                message_string = 'manual setting of dz_soil '//                &
                                  'requires adequate setting of root_fraction'//&
                                  '/= 9999999.9 ' //                            &
@@ -1292,4 +1308,5 @@
                 f_qsws_liq,  & !< factor for qsws_liq
                 f_shf,       & !< factor for shf
+                lambda_soil, & !< Thermal conductivity of the uppermost soil layer (W/m2/K)
                 lambda_surface, & !< Current value of lambda_surface (W/m2/K)
                 m_liq_max,   & !< maxmimum value of the liq. water reservoir
@@ -1370,14 +1387,21 @@
 !
 !--    Define heat conductivity between surface and soil depending on surface
-!--    type. For vegetation, a skin layer parameterization is used. For bare
-!--    soil and pavements, no skin layer is applied. In these cases, the
-!--    temperature is assumed to be constant between the surface and the first
-!--    soil layer. The heat conductivity is then derived from the soil/
-!--    pavement properties. For water surfaces, the conductivity is already set
-!--    to 1E10.
+!--    type. For vegetation, a skin layer parameterization is used. The new
+!--    parameterization uses a combination of two conductivities: a constant
+!--    conductivity for the skin layer, and a conductivity according to the
+!--    uppermost soil layer. For bare soil and pavements, no skin layer is 
+!--    applied. In these cases, the temperature is assumed to be constant 
+!--    between the surface and the first soil layer. The heat conductivity is 
+!--    then derived from the soil/pavement properties. 
+!--    For water surfaces, the conductivity is already set to 1E10.
 !--    Moreover, the heat capacity is set. For bare soil the heat capacity is
-!--    the capacity of the uppermost soil layer
-       IF ( surf%lambda_surface_s(m) == 0.0_wp )  THEN 
-       
+!--    the capacity of the uppermost soil layer, for pavement it is that of
+!--    the material involved.
+
+!
+!--    for vegetation type surfaces, the thermal conductivity of the soil is
+!--    needed
+       IF ( surf%vegetation_surface(m) )  THEN
+
           lambda_h_sat = lambda_h_sm**(1.0_wp - surf%m_sat(nzb_soil,m)) *      &
                          lambda_h_water ** surf_m_soil%var_2d(nzb_soil,m)
@@ -1386,23 +1410,34 @@
                                                      surf%m_sat(nzb_soil,m) ) )                   
                          
-          lambda_surface = (ke * (lambda_h_sat - lambda_h_dry) + lambda_h_dry )&
-                           * ddz_soil_layer(nzb_soil) * 2.0_wp
-          
-          c_surface_tmp = (rho_c_soil * ( 1.0_wp - surf%m_sat(nzb_soil,m) )    &
-                           + rho_c_water * surf_m_soil%var_2d(nzb_soil,m) )    &
-                           * dz_soil_layer(nzb_soil) * 0.5_wp               
-
-       ELSEIF ( surf_t_surface%var_1d(m) >= surf_t_soil%var_2d(nzb_soil,m)) &
-       THEN
-
+          lambda_soil = (ke * (lambda_h_sat - lambda_h_dry) + lambda_h_dry )   &
+                           * ddz_soil(nzb_soil) * 2.0_wp
+
+!
+!--       When bare soil is set without a thermal conductivity (no skin layer),
+!--       a heat capacity is that of the soil layer, otherwise it is a
+!--       combination of the conductivities from the skin and the soil layer
+          IF ( surf%lambda_surface_s(m) == 0.0_wp )  THEN 
+            surf%c_surface(m) = (rho_c_soil * (1.0_wp - surf%m_sat(nzb_soil,m))&
+                              + rho_c_water * surf_m_soil%var_2d(nzb_soil,m) ) &
+                              * dz_soil(nzb_soil) * 0.5_wp    
+            lambda_surface = lambda_soil
+
+          ELSE IF ( surf_t_surface%var_1d(m) >= surf_t_soil%var_2d(nzb_soil,m))&
+          THEN
+             lambda_surface = surf%lambda_surface_s(m) * lambda_soil           &
+                              / ( surf%lambda_surface_s(m) + lambda_soil )
+          ELSE
+
+             lambda_surface = surf%lambda_surface_u(m) * lambda_soil           &
+                              / ( surf%lambda_surface_u(m) + lambda_soil )
+          ENDIF
+       ELSE
           lambda_surface = surf%lambda_surface_s(m)
-          c_surface_tmp = surf%c_surface(m) 
-          
-       ELSE
-
-          lambda_surface = surf%lambda_surface_u(m)
-          c_surface_tmp = surf%c_surface(m) 
-          
-       ENDIF
+       ENDIF
+
+!
+!--    Set heat capacity of the skin/surface. It is ususally zero when a skin
+!--    layer is used, and non-zero otherwise.
+       c_surface_tmp = surf%c_surface(m) 
 
 !
@@ -1686,5 +1721,5 @@
                           * surf_t_surface%var_1d(m) - dq_s_dt *               &
                             surf_t_surface_p%var_1d(m) ) /                     &
-                            surf%qsws(m) - surf%r_a(m)
+                            (surf%qsws(m) + 1.0E-20)  - surf%r_a(m)
        ENDIF
 
@@ -2069,32 +2104,76 @@
        ENDDO
 
+
+
 !
 !--    Calculate grid spacings. Temperature and moisture are defined at
 !--    the center of the soil layers, whereas gradients/fluxes are 
 !--    defined at the edges (_layer)
+       zs(nzb_soil) = 0.5_wp * dz_soil(nzb_soil)
+       zs_layer(nzb_soil) = dz_soil(nzb_soil)
+
+       DO  k = nzb_soil+1, nzt_soil
+          zs_layer(k) = zs_layer(k-1) + dz_soil(k)
+          zs(k) = (zs_layer(k) +  zs_layer(k-1)) * 0.5_wp
+
+       ENDDO
+
+       dz_soil(nzt_soil+1) = zs_layer(nzt_soil) + dz_soil(nzt_soil)
+       zs(nzt_soil+1) = zs_layer(nzt_soil) + 0.5_wp * dz_soil(nzt_soil)
+ 
        DO  k = nzb_soil, nzt_soil-1
-          dz_soil(k) = zs(k+1) - zs(k)
+          dz_soil_center(k) = zs(k+1) - zs(k)
           
-          IF ( dz_soil(k) == 0.0_wp )  THEN
+          IF ( dz_soil_center(k) == 0.0_wp )  THEN
              message_string = 'invalid soil layer configuration found ' //     &
-                              '(dz_soil(k) = 0.0)'
+                              '(dz_soil_center(k) = 0.0)'
              CALL message( 'lsm_read_restart_data', 'PA0140', 1, 2, 0, 6, 0 )
-          ENDIF
-          
+          ENDIF  
        ENDDO
-
-       dz_soil_layer(nzb_soil) = 2.0_wp * zs(nzb_soil)
-       zs_layer(nzb_soil)      = 2.0_wp * zs(nzb_soil)
   
-       DO  k = nzb_soil+1, nzt_soil
-          dz_soil_layer(k) = ( zs(k) - zs_layer(k-1) ) * 2.0_wp
-          zs_layer(k) = zs_layer(k-1) + dz_soil_layer(k)
-       ENDDO
-       dz_soil_layer(nzt_soil+1) = 0.5_wp * dz_soil_layer(nzt_soil)
-       dz_soil(nzt_soil) = zs_layer(k-1) + dz_soil_layer(k) - zs(nzt_soil)
+       dz_soil_center(nzt_soil) = zs_layer(k-1) + dz_soil(k) - zs(nzt_soil)
+
+       IF ( myid == 0 )  THEN
+          PRINT*, "zs:", zs
+          PRINT*, "zs_layer:", zs_layer
+          PRINT*, "dz_soil:", dz_soil
+          PRINT*, "dz_soil_center:", dz_soil_center
+       ENDIF
 
        
-       ddz_soil       = 1.0_wp / dz_soil
-       ddz_soil_layer = 1.0_wp / dz_soil_layer
+       ddz_soil_center       = 1.0_wp / dz_soil_center
+       ddz_soil = 1.0_wp / dz_soil
+
+
+
+
+
+
+!        DO  k = nzb_soil, nzt_soil-1
+!           dz_soil_center(k) = zs(k+1) - zs(k)
+!           
+!           IF ( dz_soil_center(k) == 0.0_wp )  THEN
+!              message_string = 'invalid soil layer configuration found ' //     &
+!                               '(dz_soil_center(k) = 0.0)'
+!              CALL message( 'lsm_read_restart_data', 'PA0140', 1, 2, 0, 6, 0 )
+!           ENDIF
+!           
+!        ENDDO
+! 
+!        dz_soil(nzb_soil) = 2.0_wp * zs(nzb_soil)
+!        zs_layer(nzb_soil)      = 2.0_wp * zs(nzb_soil)
+!   
+!        DO  k = nzb_soil+1, nzt_soil
+!           dz_soil(k) = ( zs(k) - zs_layer(k-1) ) * 2.0_wp
+!           zs_layer(k) = zs_layer(k-1) + dz_soil(k)
+!        ENDDO
+!        dz_soil(nzt_soil+1) = 0.5_wp * dz_soil(nzt_soil)
+!        dz_soil_center(nzt_soil) = zs_layer(k-1) + dz_soil(k) - zs(nzt_soil)
+! 
+!        
+!        ddz_soil_center       = 1.0_wp / dz_soil_center
+!        ddz_soil = 1.0_wp / dz_soil
+
+
 
 
@@ -2208,19 +2287,10 @@
           ENDIF   
           
-!
-!--       The heat conductivity between canopy and first soil layer must
-!--       depend on the depth of the soil layer. The value in the lookup
-!--       table refers to the first soil temperature at zs = -0.005 m. Lambda 
-!--       is thus interpolated to the actual soil layer depth used.
           IF ( lambda_surface_stable == 9999999.9_wp )  THEN
-             lambda_surface_stable = vegetation_pars(6,vegetation_type)        &
-                                     * 0.005_wp * ddz_soil_layer(nzb_soil)     &
-                                     * 2.0_wp    
+             lambda_surface_stable = vegetation_pars(6,vegetation_type) 
           ENDIF
 
           IF ( lambda_surface_unstable == 9999999.9_wp )  THEN
-             lambda_surface_unstable = vegetation_pars(7,vegetation_type)      &
-                                       * 0.005_wp  * ddz_soil_layer(nzb_soil)  &
-                                       * 2.0_wp             
+             lambda_surface_unstable = vegetation_pars(7,vegetation_type)           
           ENDIF
 
@@ -2232,5 +2302,13 @@
              c_surface = vegetation_pars(9,vegetation_type)        
           ENDIF
-              
+ 
+          IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
+             albedo_type = INT(vegetation_pars(10,vegetation_type))       
+          ENDIF
+   
+          IF ( emissivity == 9999999.9_wp )  THEN
+             emissivity = vegetation_pars(11,vegetation_type)     
+          ENDIF
+          
        ENDIF
 
@@ -2242,10 +2320,18 @@
 
           IF ( z0_water == 9999999.9_wp )  THEN
-             z0_water = water_pars(1,water_type)       
+             z0_water = water_pars(1,water_type)
           ENDIF       
 
           IF ( z0h_water == 9999999.9_wp )  THEN
-             z0h_water = water_pars(2,water_type)       
+             z0h_water = water_pars(2,water_type)    
           ENDIF  
+
+          IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
+             albedo_type = INT(water_pars(3,water_type))        
+          ENDIF
+   
+          IF ( emissivity == 9999999.9_wp )  THEN
+             emissivity = water_pars(4,water_type)        
+          ENDIF
 
        ENDIF 
@@ -2256,5 +2342,5 @@
              pavement_depth = pavement_pars(0,pavement_type)
           ELSE
-             pavement_depth = MAX( zs(nzb_soil), pavement_depth )
+             pavement_depth = MAX( 0.5_wp * dz_soil(nzb_soil), pavement_depth )
           ENDIF
 
@@ -2275,4 +2361,12 @@
           ENDIF   
     
+          IF ( albedo_type == 9999999  .AND.  albedo == 9999999.9_wp )  THEN
+             albedo_type = INT(pavement_pars(5,pavement_type))        
+          ENDIF
+   
+          IF ( emissivity == 9999999.9_wp )  THEN
+             emissivity = pavement_pars(6,pavement_type)        
+          ENDIF
+
        ENDIF
 !
@@ -2340,9 +2434,9 @@
                 t_soil_h%var_2d(nzt_soil+1,m) = soil_temperature(nzt_soil+1)
                 surf_lsm_h%lambda_surface_s(m) = pavement_heat_conduct         &
-                                                 * ddz_soil_layer(nzb_soil)    &
+                                                 * ddz_soil(nzb_soil)    &
                                                  * 2.0_wp
                 surf_lsm_h%lambda_surface_u(m) = surf_lsm_h%lambda_surface_s(m)                                         
                 surf_lsm_h%c_surface(m)        = pavement_heat_capacity        &
-                                                 * dz_soil_layer(nzb_soil)     &
+                                                 * dz_soil(nzb_soil)     &
                                                  * 0.25_wp
                 surf_lsm_h%lambda_h_def(m)    = pavement_heat_conduct                                     
@@ -2428,9 +2522,9 @@
                    t_soil_v(l)%var_2d(nzt_soil+1,m)  = soil_temperature(nzt_soil+1)
                    surf_lsm_v(l)%lambda_surface_s(m) = pavement_heat_conduct      &
-                                                       * ddz_soil_layer(nzb_soil) &
+                                                       * ddz_soil(nzb_soil) &
                                                        * 2.0_wp
                    surf_lsm_v(l)%lambda_surface_u(m) = surf_lsm_v(l)%lambda_surface_s(m)
                    surf_lsm_v(l)%c_surface(m) = pavement_heat_capacity         &
-                                                * dz_soil_layer(nzb_soil)      &
+                                                * dz_soil(nzb_soil)      &
                                                 * 0.25_wp
                    surf_lsm_v(l)%lambda_h_def(m)   = pavement_heat_conduct                                     
@@ -2642,10 +2736,10 @@
 !
 !--    Allocate global 1D arrays
-       ALLOCATE ( ddz_soil(nzb_soil:nzt_soil) )
-       ALLOCATE ( ddz_soil_layer(nzb_soil:nzt_soil+1) )
-       ALLOCATE ( dz_soil(nzb_soil:nzt_soil) )
-       ALLOCATE ( dz_soil_layer(nzb_soil:nzt_soil+1) )
+       ALLOCATE ( ddz_soil_center(nzb_soil:nzt_soil) )
+       ALLOCATE ( ddz_soil(nzb_soil:nzt_soil+1) )
+       ALLOCATE ( dz_soil_center(nzb_soil:nzt_soil) )
        ALLOCATE ( root_extr(nzb_soil:nzt_soil) )
-  
+       ALLOCATE ( zs(nzb_soil:nzt_soil+1) )
+
        root_extr = 0.0_wp 
        
@@ -2829,4 +2923,5 @@
                                   constant_roughness,                          &
                                   conserve_water_content,                      &
+                                  dz_soil,                                     &
                                   f_shortwave_incoming, field_capacity,        & 
                                   aero_resist_kray, hydraulic_conductivity,    &
@@ -2843,5 +2938,5 @@
                                   vegetation_coverage, vegetation_type,        &
                                   water_temperature, water_type,               &
-                                  wilting_point, zs, z0_vegetation,            &
+                                  wilting_point, z0_vegetation,                &
                                   z0h_vegetation, z0q_vegetation, z0_water,    &
                                   z0h_water, z0q_water, z0_pavement,           &
@@ -2979,7 +3074,7 @@
                      .AND.  zs(k+1) > pavement_depth )                         &
                 THEN
-                   surf%lambda_h(k,m) = ( pavement_depth - zs(k) ) * ddz_soil(k+1)   &
+                   surf%lambda_h(k,m) = ( pavement_depth - zs(k) ) * ddz_soil_center(k+1)   &
                                       * lambda_temp(k)                               &
-                                      + ( zs(k+1) - pavement_depth ) * ddz_soil(k+1) &
+                                      + ( zs(k+1) - pavement_depth ) * ddz_soil_center(k+1) &
                                       * lambda_temp(k+1)
                 ELSE
@@ -2997,6 +3092,6 @@
              tend(nzb_soil) = ( 1.0_wp / surf%rho_c_total(nzb_soil,m) ) *            &
                     ( surf%lambda_h(nzb_soil,m) * ( surf_t_soil%var_2d(nzb_soil+1,m) &
-                      - surf_t_soil%var_2d(nzb_soil,m) ) * ddz_soil(nzb_soil)        &
-                      + surf%ghf(m) ) * ddz_soil_layer(nzb_soil)
+                      - surf_t_soil%var_2d(nzb_soil,m) ) * ddz_soil_center(nzb_soil)        &
+                      + surf%ghf(m) ) * ddz_soil(nzb_soil)
 
              DO  k = nzb_soil+1, nzt_soil
@@ -3004,9 +3099,9 @@
                           * (   surf%lambda_h(k,m)                                      &
                               * ( surf_t_soil%var_2d(k+1,m) - surf_t_soil%var_2d(k,m) ) &
-                              * ddz_soil(k)                                             &
+                              * ddz_soil_center(k)                                             &
                               - surf%lambda_h(k-1,m)                                    &
                               * ( surf_t_soil%var_2d(k,m) - surf_t_soil%var_2d(k-1,m) ) &
-                              * ddz_soil(k-1)                                           &
-                            ) * ddz_soil_layer(k)
+                              * ddz_soil_center(k-1)                                           &
+                            ) * ddz_soil(k)
 
              ENDDO
@@ -3143,19 +3238,19 @@
                 tend(nzb_soil) = ( surf%lambda_w(nzb_soil,m) *   (       &
                          surf_m_soil%var_2d(nzb_soil+1,m) - surf_m_soil%var_2d(nzb_soil,m) )           &
-                         * ddz_soil(nzb_soil) - surf%gamma_w(nzb_soil,m) - &
+                         * ddz_soil_center(nzb_soil) - surf%gamma_w(nzb_soil,m) - &
                                                      (                         &
                             root_extr(nzb_soil) * surf%qsws_veg(m)    &
                             + surf%qsws_soil(m) ) * drho_l_lv )       &
-                            * ddz_soil_layer(nzb_soil)
+                            * ddz_soil(nzb_soil)
 
                 DO  k = nzb_soil+1, nzt_soil-1
                    tend(k) = ( surf%lambda_w(k,m) * ( surf_m_soil%var_2d(k+1,m)      &
-                             - surf_m_soil%var_2d(k,m) ) * ddz_soil(k)                   &
+                             - surf_m_soil%var_2d(k,m) ) * ddz_soil_center(k)                   &
                              - surf%gamma_w(k,m)                         &
                              - surf%lambda_w(k-1,m) * ( surf_m_soil%var_2d(k,m) -    &
-                             surf_m_soil%var_2d(k-1,m)) * ddz_soil(k-1)                      &
+                             surf_m_soil%var_2d(k-1,m)) * ddz_soil_center(k-1)                      &
                              + surf%gamma_w(k-1,m) - (root_extr(k)       &
                              * surf%qsws_veg(m) * drho_l_lv)          &
-                             ) * ddz_soil_layer(k)
+                             ) * ddz_soil(k)
                 ENDDO
                 tend(nzt_soil) = ( - surf%gamma_w(nzt_soil,m)            &
@@ -3163,9 +3258,9 @@
                                      * ( surf_m_soil%var_2d(nzt_soil,m)                    &
                                      - surf_m_soil%var_2d(nzt_soil-1,m))                   &
-                                     * ddz_soil(nzt_soil-1)                      &
+                                     * ddz_soil_center(nzt_soil-1)                      &
                                      + surf%gamma_w(nzt_soil-1,m) - (    &
                                        root_extr(nzt_soil)                     &
                                      * surf%qsws_veg(m) * drho_l_lv ) &
-                                  ) * ddz_soil_layer(nzt_soil)             
+                                  ) * ddz_soil(nzt_soil)             
 
                 surf_m_soil_p%var_2d(nzb_soil:nzt_soil,m) = surf_m_soil%var_2d(nzb_soil:nzt_soil,m)    &
Index: /palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 2327)
+++ /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 2328)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! Emissivity can now be set individually for each pixel.
+! Albedo type can be inferred from land surface model.
+! Added default albedo type for bare soil
+! 
+! 2318 2017-07-20 17:27:44Z suehring
 ! Get topography top index via Function call 
 ! 
@@ -211,5 +216,5 @@
 !
 !-- Predefined Land surface classes (albedo_type) after Briegleb (1992)
-    CHARACTER(37), DIMENSION(0:17), PARAMETER :: albedo_type_name = (/      &
+    CHARACTER(37), DIMENSION(0:18), PARAMETER :: albedo_type_name = (/      &
                                    'user defined                         ', & !  0 
                                    'ocean                                ', & !  1
@@ -229,11 +234,12 @@
                                    'sea ice                              ', & ! 15 
                                    'snow                                 ', & ! 16
-                                   'pavement/roads                       '  & ! 17
+                                   'pavement/roads                       ', & ! 17
+                                   'bare soil                            '  & ! 18
                                                          /)
 
-    INTEGER(iwp) :: albedo_type  = 5,    & !< Albedo surface type (default: short grassland)
-                    day,                 & !< current day of the year
-                    day_init     = 172,  & !< day of the year at model start (21/06)
-                    dots_rad     = 0       !< starting index for timeseries output 
+    INTEGER(iwp) :: albedo_type  = 9999999, & !< Albedo surface type
+                    day,                    & !< current day of the year
+                    day_init     = 172,     & !< day of the year at model start (21/06)
+                    dots_rad     = 0          !< starting index for timeseries output 
 
     LOGICAL ::  unscheduled_radiation_calls = .TRUE., & !< flag parameter indicating whether additional calls of the radiation code are allowed
@@ -261,5 +267,5 @@
                 decl_3,                          & !< declination coef. 3
                 dt_radiation = 0.0_wp,           & !< radiation model timestep
-                emissivity = 0.98_wp,            & !< NAMELIST surface emissivity
+                emissivity = 9999999.9_wp,       & !< NAMELIST surface emissivity
                 lambda = 0.0_wp,                 & !< longitude in degrees
                 lon = 0.0_wp,                    & !< longitude in radians
@@ -278,4 +284,5 @@
     REAL(wp), DIMENSION(:,:), ALLOCATABLE :: &
                 alpha,                       & !< surface broadband albedo (used for clear-sky scheme)
+                emis,                        & !< surface broadband emissitivity
                 rad_lw_out_change_0,         & !< change in LW out due to change in surface temperature
                 rad_net,                     & !< net radiation at the surface
@@ -285,5 +292,5 @@
 !-- Land surface albedos for solar zenith angle of 60Â° after Briegleb (1992)     
 !-- (shortwave, longwave, broadband):   sw,      lw,      bb,
-    REAL(wp), DIMENSION(0:2,1:17), PARAMETER :: albedo_pars = RESHAPE( (/& 
+    REAL(wp), DIMENSION(0:2,1:18), PARAMETER :: albedo_pars = RESHAPE( (/& 
                                    0.06_wp, 0.06_wp, 0.06_wp,            & !  1
                                    0.09_wp, 0.28_wp, 0.19_wp,            & !  2
@@ -302,6 +309,7 @@
                                    0.90_wp, 0.65_wp, 0.77_wp,            & ! 15
                                    0.95_wp, 0.70_wp, 0.82_wp,            & ! 16
-                                   0.08_wp, 0.08_wp, 0.08_wp             & ! 17
-                                 /), (/ 3, 17 /) )
+                                   0.08_wp, 0.08_wp, 0.08_wp,            & ! 17
+                                   0.17_wp, 0.17_wp, 0.17_wp             & ! 18
+                                 /), (/ 3, 18 /) )
 
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE, TARGET :: &
@@ -510,5 +518,5 @@
 !
 !-- Public variables and constants / NEEDS SORTING
-    PUBLIC decl_1, decl_2, decl_3, dots_rad, dt_radiation, force_radiation_call,&
+    PUBLIC albedo, albedo_type, decl_1, decl_2, decl_3, dots_rad, dt_radiation, emissivity, force_radiation_call,&
            lat, lon, rad_net, rad_net_av, radiation, radiation_scheme, rad_lw_in,        &
            rad_lw_in_av, rad_lw_out, rad_lw_out_av, rad_lw_out_change_0,       &
@@ -888,4 +896,11 @@
 
 !
+!--    Allocate array for storing emissivity
+       IF ( .NOT. ALLOCATED ( emis ) )  THEN
+          ALLOCATE ( emis(nysg:nyng,nxlg:nxrg) )
+          emis = emissivity
+       ENDIF
+
+!
 !--    Allocate array for storing the surface net radiation
        IF ( .NOT. ALLOCATED ( rad_net ) )  THEN
@@ -956,8 +971,9 @@
 !
 !--       Overwrite albedo if manually set in parameter file
-          IF ( albedo_type /= 0 .AND. albedo == 9999999.9_wp )  THEN
+          IF ( albedo_type /= 0  .AND.  albedo_type /= 9999999 .AND. albedo == 9999999.9_wp )  THEN
              albedo = albedo_pars(2,albedo_type)
           ENDIF
-    
+!
+!--       Write albedo to 2d array alpha to allow surface heterogeneities    
           alpha = albedo
   
@@ -1194,5 +1210,5 @@
              rad_sw_in(0,j,i)  = solar_constant * sky_trans * zenith(0)
              rad_sw_out(0,j,i) = alpha(j,i) * rad_sw_in(0,j,i)
-             rad_lw_out(0,j,i) = emissivity * sigma_sb * (pt(k,j,i) * exn)**4
+             rad_lw_out(0,j,i) = emis(j,i) * sigma_sb * (pt(k,j,i) * exn)**4
 
              IF ( cloud_physics )  THEN
@@ -1207,5 +1223,5 @@
 
 
-             rad_lw_out_change_0(j,i) = 3.0_wp * sigma_sb * emissivity         &
+             rad_lw_out_change_0(j,i) = 3.0_wp * sigma_sb * emis(j,i)         &
                                         * (pt(k,j,i) * exn) ** 3
 
@@ -1254,5 +1270,5 @@
              ENDIF
 
-             rad_lw_out(0,j,i) = emissivity * sigma_sb * (pt(k,j,i) * exn)**4
+             rad_lw_out(0,j,i) = emis(j,i) * sigma_sb * (pt(k,j,i) * exn)**4
 
              rad_sw_in(0,j,i) = ( rad_net(j,i) - rad_lw_in(0,j,i)              &
@@ -1533,4 +1549,8 @@
              rrtm_tsfc = pt(nzb,j,i) * (surface_pressure / 1000.0_wp )**0.286_wp
 
+!
+!--          Set surface emissivity
+             rrtm_emis = emis(j,i)
+
              IF ( lw_radiation )  THEN
                CALL rrtmg_lw( 1, nzt_rad      , rrtm_icld    , rrtm_idrv      ,&
@@ -1723,4 +1743,11 @@
                  rrtm_aldir(0,:,:) = aldif
                  rrtm_asdir(0,:,:) = asdif
+
+!
+!--        Bare soil
+           ELSEIF ( albedo_type == 18 )  THEN
+                 rrtm_aldir(0,:,:) = aldif
+                 rrtm_asdir(0,:,:) = asdif
+
 !
 !--        Land surfaces
