Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 152)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 153)
@@ -4,4 +4,5 @@
 ! Actual revisions:
 ! -----------------
+! Leaf area density (LAD) explicitly set to its surface value at k=0
 ! Case of reading data for recycling included in initializing_actions,
 ! check of turbulent_inflow and calculation of recycling_plane
@@ -955,4 +956,6 @@
 
           IF ( .NOT. ocean ) THEN
+
+             lad(0) = lad_surface
  
              lad_vertical_gradient_level_ind(1) = 0
Index: palm/trunk/SOURCE/header.f90
===================================================================
--- palm/trunk/SOURCE/header.f90	(revision 152)
+++ palm/trunk/SOURCE/header.f90	(revision 153)
@@ -371,4 +371,12 @@
 
        WRITE ( io, 280 ) canopy_mode, pch_index, drag_coefficient
+       IF ( passive_scalar ) THEN
+          WRITE ( io, 281 ) scalar_exchange_coefficient,   &
+                            leaf_surface_concentration
+       ENDIF
+
+!
+!--    Heat flux at the top of vegetation
+       WRITE ( io, 282 ) cthf
 
 !
@@ -397,5 +405,5 @@
        ENDDO
 
-       WRITE ( io, 281 )  TRIM( coordinates ), TRIM( learde ), &
+       WRITE ( io, 283 )  TRIM( coordinates ), TRIM( learde ), &
                           TRIM( gradients ), TRIM( slices )
 
@@ -1379,5 +1387,8 @@
               ' Canopy top: ',I4 / &
               ' Leaf drag coefficient: ',F6.2 /)
-281 FORMAT (/ ' Characteristic levels of the leaf area density:'// &
+281 FORMAT (/ ' Scalar_exchange_coefficient: ',F6.2 / &
+              ' Scalar concentration at leaf surfaces in kg/m**3: ',F6.2 /)
+282 FORMAT (' Predefined constant heatflux at the top of the vegetation: ',F6.2,' K m/s')
+283 FORMAT (/ ' Characteristic levels of the leaf area density:'// &
               ' Height:              ',A,'  m'/ &
               ' Leaf area density:   ',A,'  m**2/m**3'/ &
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 152)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 153)
@@ -7,4 +7,6 @@
 ! Actual revisions:
 ! -----------------
+! Further amendments to and modifications in the initialisation of the plant 
+! canopy model 
 ! Allocation of hom_sum moved to parin, initialization of spectrum_x|y directly
 ! after allocating theses arrays,
@@ -90,4 +92,6 @@
     INTEGER, DIMENSION(:,:), ALLOCATABLE ::  ngp_2dh_outer_l,  &
          ngp_2dh_s_inner_l
+
+    REAL ::  a, b
 
     REAL, DIMENSION(1:2) ::  volume_flow_area_l, volume_flow_initial_l
@@ -248,4 +252,15 @@
                   lad_w(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
                   cdc(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+
+       IF ( passive_scalar ) THEN
+          ALLOCATE ( sls(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),   &
+                     sec(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) ) 
+       ENDIF
+
+       IF ( cthf /= 0.0 ) THEN
+          ALLOCATE ( lai(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),   &
+                     canopy_heat_flux(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
+       ENDIF
+
     ENDIF
 
@@ -1078,4 +1093,8 @@
                    lad_s(:,j,i) = lad(:)
                    cdc(:,j,i)   = drag_coefficient
+                   IF ( passive_scalar ) THEN
+                      sls(:,j,i) = leaf_surface_concentration
+                      sec(:,j,i) = scalar_exchange_coefficient
+                   ENDIF
                 ENDDO
              ENDDO
@@ -1096,9 +1115,25 @@
        CALL exchange_horiz( cdc )
 
+       IF ( passive_scalar ) THEN
+          CALL exchange_horiz( sls )
+          CALL exchange_horiz( sec )
+       ENDIF
+
+!
+!--    Sharp boundaries of the plant canopy in horizontal directions
+!--    In vertical direction the interpolation is retained, as the leaf
+!--    area density is initialised by prescribing a vertical profile 
+!--    consisting of piecewise linear segments. The upper boundary 
+!--    of the plant canopy is now defined by lad_w(pch_index,:,:) = 0.0.
+
        DO  i = nxl, nxr
           DO  j = nys, nyn
              DO  k = nzb, nzt+1 
-                lad_u(k,j,i) = 0.5 * ( lad_s(k,j,i-1) + lad_s(k,j,i) )
-                lad_v(k,j,i) = 0.5 * ( lad_s(k,j-1,i) + lad_s(k,j,i) )
+                IF ( lad_s(k,j,i) > 0.0 ) THEN
+                   lad_u(k,j,i)   = lad_s(k,j,i) 
+                   lad_u(k,j,i+1) = lad_s(k,j,i)
+                   lad_v(k,j,i)   = lad_s(k,j,i)
+                   lad_v(k,j+1,i) = lad_s(k,j,i)
+                ENDIF
              ENDDO
              DO  k = nzb, nzt
@@ -1108,10 +1143,53 @@
        ENDDO
 
-       lad_w(nzt+1,:,:) = lad_w(nzt,:,:)
+       lad_w(pch_index,:,:) = 0.0
+       lad_w(nzt+1,:,:)     = lad_w(nzt,:,:)
 
        CALL exchange_horiz( lad_u )
        CALL exchange_horiz( lad_v )
        CALL exchange_horiz( lad_w )
-  
+
+!
+!--    Initialisation of the canopy heat source distribution
+       IF ( cthf /= 0.0 ) THEN
+!
+!--       Piecewise evaluation of the leaf area index by 
+!--       integration of the leaf area density
+          lai(:,:,:) = 0.0
+          DO  i = nxl-1, nxr+1
+             DO  j = nys-1, nyn+1
+                DO  k = pch_index-1, 0, -1
+                   lai(k,j,i) = lai(k+1,j,i) +                   &
+                                ( 0.5 * ( lad_w(k+1,j,i) +       &
+                                          lad_s(k+1,j,i) ) *     &
+                                  ( zw(k+1) - zu(k+1) ) )  +     &
+                                ( 0.5 * ( lad_w(k,j,i)   +       &
+                                          lad_s(k+1,j,i) ) *     &
+                                  ( zu(k+1) - zw(k) ) )
+                ENDDO
+             ENDDO
+          ENDDO
+
+!
+!--       Evaluation of the upward kinematic vertical heat flux within the 
+!--       canopy
+          DO  i = nxl-1, nxr+1
+             DO  j = nys-1, nyn+1
+                DO  k = 0, pch_index
+                   canopy_heat_flux(k,j,i) = cthf *                    &
+                                             exp( -0.6 * lai(k,j,i) )
+                ENDDO
+             ENDDO
+          ENDDO
+
+!
+!--       The near surface heat flux is derived from the heat flux 
+!--       distribution within the canopy
+          shf(:,:) = canopy_heat_flux(0,:,:)
+
+          IF ( ASSOCIATED( shf_m ) ) shf_m = shf
+
+       ENDIF
+
     ENDIF
 
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 152)
+++ palm/trunk/SOURCE/modules.f90	(revision 153)
@@ -5,4 +5,6 @@
 ! Actual revisions:
 ! -----------------
+! +canopy_heat_flux, cthf, lai,
+! +leaf_surface_concentration, scalar_exchange_coefficient, sec, sls
 ! +hor_index_bounds, hor_index_bounds_previous_run, id_inflow,
 ! inflow_damping_*, mean_inflow_profiles, numprocs_previous_run, nx_on_file,
@@ -135,7 +137,7 @@
 
     REAL, DIMENSION(:,:,:), ALLOCATABLE ::                                     &
-          cdc, d, diss, lad_s, lad_u, lad_v, lad_w, l_wall, tend,              &
-          u_m_l, u_m_n, u_m_r, u_m_s, v_m_l, v_m_n, v_m_r, v_m_s, w_m_l,       &
-          w_m_n, w_m_r, w_m_s
+          canopy_heat_flux, cdc, d, diss, lad_s, lad_u, lad_v, lad_w, lai,     &
+          l_wall, sec, sls, tend, u_m_l, u_m_n, u_m_r, u_m_s, v_m_l, v_m_n,    &
+          v_m_r, v_m_s, w_m_l, w_m_n, w_m_r, w_m_s
 
     REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                             &
@@ -375,5 +377,5 @@
              building_height = 50.0, building_length_x = 50.0, &
              building_length_y = 50.0, building_wall_left = 9999999.9, &
-             building_wall_south = 9999999.9, cfl_factor = -1.0, &
+             building_wall_south = 9999999.9, cthf = 0.0, cfl_factor = -1.0, &
              cos_alpha_surface, disturbance_amplitude = 0.25, &
              disturbance_energy_limit = 0.01, &
@@ -396,5 +398,6 @@
              f = 0.0, fs = 0.0, g = 9.81, inflow_damping_height = 9999999.9, &
              inflow_damping_width = 9999999.9, kappa = 0.4, km_constant = -1.0,&
-             km_damp_max = -1.0, lad_surface = 0.0, long_filter_factor = 0.0, &
+             km_damp_max = -1.0, lad_surface = 0.0,  &
+             leaf_surface_concentration = 0.0, long_filter_factor = 0.0, &
              maximum_cpu_time_allowed = 0.0, molecular_viscosity = 1.461E-5, &
              old_dt = 1.0E-10, omega = 7.29212E-5, omega_sor = 1.8, &
@@ -412,4 +415,5 @@
              rif_max = 1.0, &
              rif_min = -5.0, roughness_length = 0.1, sa_surface = 35.0, &
+             scalar_exchange_coefficient = 0.0, &
              simulated_time = 0.0, simulated_time_at_begin, sin_alpha_surface, &
              skip_time_data_output = 0.0, skip_time_data_output_av = 9999999.9,&
Index: palm/trunk/SOURCE/parin.f90
===================================================================
--- palm/trunk/SOURCE/parin.f90	(revision 152)
+++ palm/trunk/SOURCE/parin.f90	(revision 153)
@@ -4,4 +4,6 @@
 ! Actual revisions:
 ! -----------------
+! +cthf
+! +leaf_surface_concentration, scalar_exchange_coefficient
 ! +inflow_damping_height, inflow_damping_width, recycling_width,
 ! turbulent_inflow in inipar
@@ -82,5 +84,5 @@
              building_length_y, building_wall_left, &
              building_wall_south, canopy_mode, cloud_droplets, cloud_physics, &
-             conserve_volume_flow, cut_spline_overshoot, damp_level_1d, &
+             conserve_volume_flow, cthf, cut_spline_overshoot, damp_level_1d, &
              dissipation_1d, drag_coefficient, dt, dt_pr_1d, &
              dt_run_control_1d, dx, dy, dz, dz_max, dz_stretch_factor, & 
@@ -91,6 +93,7 @@
              initializing_actions, km_constant, km_damp_max, lad_surface, &
              lad_vertical_gradient, lad_vertical_gradient_level, &
-             long_filter_factor, loop_optimization, mixing_length_1d, &
-             momentum_advec, netcdf_precision, nsor_ini, nx, ny, &
+             leaf_surface_concentration, long_filter_factor, &
+             loop_optimization, mixing_length_1d, momentum_advec, &
+             netcdf_precision, nsor_ini, nx, ny, &
              nz, ocean, omega, outflow_damping_width, overshoot_limit_e, &
              overshoot_limit_pt, overshoot_limit_u, overshoot_limit_v, &
@@ -102,5 +105,6 @@
              random_generator, random_heatflux, recycling_width, rif_max, &
              rif_min, roughness_length, sa_surface, sa_vertical_gradient, &
-             sa_vertical_gradient_level, scalar_advec, statistic_regions, &
+             sa_vertical_gradient_level, scalar_advec, &
+             scalar_exchange_coefficient, statistic_regions, &
              surface_heatflux, surface_pressure, surface_scalarflux, &
              surface_waterflux, s_surface, s_surface_initial_change, &
Index: palm/trunk/SOURCE/plant_canopy_model.f90
===================================================================
--- palm/trunk/SOURCE/plant_canopy_model.f90	(revision 152)
+++ palm/trunk/SOURCE/plant_canopy_model.f90	(revision 153)
@@ -4,4 +4,10 @@
 ! Actual revisions:
 ! -----------------
+! Consideration of heat sources within the forest canopy which represent the
+! rate of heat input into the air from the forest leaves
+!
+! Evaluation of sinks and sources for scalar concentration due to canopy 
+! elements
+!
 ! Bugfix: remove IF statement in plant_canopy_model_ij
 !
@@ -15,5 +21,6 @@
 ! Description:
 ! ------------
-! Evaluation of the drag due to vegetation
+! Evaluation of sinks and sources of momentum, heat and scalar concentration 
+! due to canopy elements
 !------------------------------------------------------------------------------!
 
@@ -59,5 +66,5 @@
                                               v(k,j,i)        +  &
                                               v(k,j+1,i)      +  &
-                                              v(k,j+1,i+1) )     &
+                                              v(k,j+1,i-1) )     &
                                             / 4.0 )**2        +  &
                                           ( ( w(k-1,j,i-1)    +  &
@@ -120,6 +127,42 @@
 
 !
+!--       potential temperature
+          CASE ( 4 )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_s_inner(j,i)+1, pch_index
+                      tend(k,j,i) = tend(k,j,i) +                     &
+                                    ( canopy_heat_flux(k,j,i) -     &
+                                      canopy_heat_flux(k-1,j,i) ) / &
+                                      dzw(k)
+                   ENDDO
+                ENDDO
+             ENDDO
+
+!
+!--       scalar concentration
+          CASE ( 5 )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_s_inner(j,i)+1, pch_index
+                      tend(k,j,i) = tend(k,j,i) -                     &
+                                    sec(k,j,i) * lad_s(k,j,i) *       &
+                                    SQRT( ( ( u(k,j,i)        +       &
+                                              u(k,j,i+1) )            &
+                                            / 2.0 )**2        +       &
+                                          ( ( v(k,j,i)        +       &
+                                              v(k,j+1,i) )            &
+                                            / 2.0 )**2        +       &
+                                          ( ( w(k-1,j,i)      +       & 
+                                              w(k,j,i) )              &
+                                            / 2.0 )**2 )      *       &
+                                    ( q(k,j,i) - sls(k,j,i) )
+                   ENDDO
+                ENDDO
+             ENDDO
+
+!
 !--       sgs-tke
-          CASE ( 4 )
+          CASE ( 6 )
              DO  i = nxl, nxr
                 DO  j = nys, nyn
@@ -180,5 +223,5 @@
                                         v(k,j,i)    +        &
                                         v(k,j+1,i)  +        &
-                                        v(k,j+1,i+1) )       &
+                                        v(k,j+1,i-1) )       &
                                       / 4.0 )**2    +        &
                                     ( ( w(k-1,j,i-1) +       &
@@ -232,6 +275,35 @@
 
 !
+!--       potential temperature
+          CASE ( 4 )
+             DO  k = nzb_s_inner(j,i)+1, pch_index
+                tend(k,j,i) = tend(k,j,i) +                     &
+                              ( canopy_heat_flux(k,j,i) -     &
+                                canopy_heat_flux(k-1,j,i) ) / &
+                                dzw(k)
+             ENDDO
+
+
+!
+!--       scalar concentration 
+          CASE ( 5 )
+             DO  k = nzb_s_inner(j,i)+1, pch_index
+                tend(k,j,i) = tend(k,j,i) -                     &
+                              sec(k,j,i) * lad_s(k,j,i) *       &
+                              SQRT( ( ( u(k,j,i)        +       &
+                                        u(k,j,i+1) )            &
+                                      / 2.0 )**2        +       &
+                                    ( ( v(k,j,i)        +       &
+                                        v(k,j+1,i) )            &
+                                      / 2.0 )**2        +       &
+                                    ( ( w(k-1,j,i)      +       &
+                                        w(k,j,i) )              &
+                                      / 2.0 )**2 )      *       &
+                              ( q(k,j,i) - sls(k,j,i) )
+             ENDDO   
+
+!
 !--       sgs-tke
-       CASE ( 4 )
+       CASE ( 6 )
           DO  k = nzb_s_inner(j,i)+1, pch_index   
              tend(k,j,i) = tend(k,j,i) -                     &
Index: palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 152)
+++ palm/trunk/SOURCE/prognostic_equations.f90	(revision 153)
@@ -4,5 +4,8 @@
 ! Actual revisions:
 ! -----------------
-! 
+! add call of plant_canopy_model in the prognostic equation for
+! the potential temperature
+! add call of plant_canopy_model in the prognostic equation for
+! the concentration of a passive scalar
 !
 ! Former revisions:
@@ -422,4 +425,11 @@
              CALL impact_of_latent_heat( i, j )
           ENDIF
+
+!
+!--       Consideration of heat sources within the plant canopy
+          IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN
+             CALL plant_canopy_model( i, j, 4 )
+          ENDIF
+
           CALL user_actions( i, j, 'pt-tendency' )
 
@@ -612,4 +622,9 @@
                 CALL calc_precipitation( i, j )
              ENDIF
+
+!
+!--          Sink or source of scalar concentration due to canopy elements
+             IF ( plant_canopy ) CALL plant_canopy_model( i, j, 5 )
+
              CALL user_actions( i, j, 'q-tendency' )
 
@@ -754,5 +769,5 @@
 !
 !--          Additional sink term for flows through plant canopies
-             IF ( plant_canopy )  CALL plant_canopy_model( i, j, 4 )         
+             IF ( plant_canopy )  CALL plant_canopy_model( i, j, 6 )         
 
              CALL user_actions( i, j, 'e-tendency' )
@@ -1034,4 +1049,11 @@
              CALL impact_of_latent_heat( i, j )
           ENDIF
+
+!
+!--       Consideration of heat sources within the plant canopy
+          IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN
+             CALL plant_canopy_model( i, j, 4 )
+          ENDIF
+
           CALL user_actions( i, j, 'pt-tendency' )
 
@@ -1142,4 +1164,10 @@
                 CALL calc_precipitation( i, j )
              ENDIF
+
+!
+!--          Sink or source of scalar concentration due to canopy elements
+             IF ( plant_canopy ) CALL plant_canopy_model( i, j, 5 )
+
+
              CALL user_actions( i, j, 'q-tendency' )
 
@@ -1219,5 +1247,5 @@
 !
 !--          Additional sink term for flows through plant canopies
-             IF ( plant_canopy )  CALL plant_canopy_model( i, j, 4 ) 
+             IF ( plant_canopy )  CALL plant_canopy_model( i, j, 6 ) 
 
              CALL user_actions( i, j, 'e-tendency' )
@@ -1582,4 +1610,12 @@
        CALL impact_of_latent_heat
     ENDIF
+
+!
+!-- Consideration of heat sources within the plant canopy
+    IF ( plant_canopy .AND. ( cthf /= 0.0 ) ) THEN
+       CALL plant_canopy_model( 4 )
+    ENDIF
+
+
     CALL user_actions( 'pt-tendency' )
 
@@ -1782,4 +1818,9 @@
           CALL calc_precipitation
        ENDIF
+
+!
+!--    Sink or source of scalar concentration due to canopy elements
+       IF ( plant_canopy ) CALL plant_canopy_model( 5 )
+
        CALL user_actions( 'q-tendency' )
 
@@ -1919,5 +1960,5 @@
 !
 !--    Additional sink term for flows through plant canopies
-       IF ( plant_canopy )  CALL plant_canopy_model( 4 )
+       IF ( plant_canopy )  CALL plant_canopy_model( 6 )
        CALL user_actions( 'e-tendency' )
 
Index: palm/trunk/SOURCE/read_var_list.f90
===================================================================
--- palm/trunk/SOURCE/read_var_list.f90	(revision 152)
+++ palm/trunk/SOURCE/read_var_list.f90	(revision 153)
@@ -4,4 +4,5 @@
 ! Actual revisions:
 ! -----------------
+! +cthf, leaf_surface_concentration, scalar_exchange_coefficient
 ! +numprocs_previous_run, hor_index_bounds_previous_run, inflow_damping_factor,
 ! inflow_damping_height, inflow_damping_width, mean_inflow_profiles,
@@ -233,4 +234,6 @@
           CASE ( 'conserve_volume_flow' )
              READ ( 13 )  conserve_volume_flow
+          CASE ( 'cthf' )
+             READ ( 13 )  cthf
           CASE ( 'current_timestep_number' )
              READ ( 13 )  current_timestep_number
@@ -312,4 +315,6 @@
           CASE ( 'last_dt_change' )
              READ ( 13 )  last_dt_change
+          CASE ( 'leaf_surface_concentration' )
+             READ ( 13 )  leaf_surface_concentration
           CASE ( 'long_filter_factor' )
              READ ( 13 )  long_filter_factor
@@ -417,4 +422,6 @@
           CASE ( 'scalar_advec' )
              READ ( 13 )  scalar_advec
+          CASE ( 'scalar_exchange_coefficient' )
+             READ ( 13 )  scalar_exchange_coefficient
           CASE ( 'simulated_time' )
              READ ( 13 )  simulated_time
Index: palm/trunk/SOURCE/write_var_list.f90
===================================================================
--- palm/trunk/SOURCE/write_var_list.f90	(revision 152)
+++ palm/trunk/SOURCE/write_var_list.f90	(revision 153)
@@ -4,4 +4,5 @@
 ! Actual revisions:
 ! -----------------
+! +cthf, leaf_surface_concentration, scalar_exchange_coefficient
 ! +numprocs, hor_index_bounds, inflow_damping_height, inflow_damping_width,
 ! mean_inflow_profiles, recycling_width, turbulent_inflow,
@@ -155,4 +156,6 @@
     WRITE ( 14 )  'current_timestep_number       '
     WRITE ( 14 )  current_timestep_number
+    WRITE ( 14 )  'cthf                          '
+    WRITE ( 14 )  cthf
     WRITE ( 14 )  'cut_spline_overshoot          '
     WRITE ( 14 )  cut_spline_overshoot
@@ -231,4 +234,6 @@
     WRITE ( 14 )  'last_dt_change                '
     WRITE ( 14 )  last_dt_change
+    WRITE ( 14 )  'leaf_surface_concentration    '
+    WRITE ( 14 )  leaf_surface_concentration
     WRITE ( 14 )  'long_filter_factor            '
     WRITE ( 14 )  long_filter_factor
@@ -333,4 +338,6 @@
     WRITE ( 14 )  'scalar_advec                  '
     WRITE ( 14 )  scalar_advec
+    WRITE ( 14 )  'scalar_exchange_coefficient   '
+    WRITE ( 14 )  scalar_exchange_coefficient
     WRITE ( 14 )  'simulated_time                '
     WRITE ( 14 )  simulated_time
