Index: palm/trunk/SOURCE/CURRENT_MODIFICATIONS
===================================================================
--- palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 137)
+++ palm/trunk/SOURCE/CURRENT_MODIFICATIONS	(revision 138)
@@ -1,4 +1,19 @@
 New:
 ---
+
+Plant canopy model of Watanabe (2004,BLM 112,307-341) added. 
+It can be switched on by the inipar parameter plant_canopy. 
+The inipar parameter canopy_mode can be used to prescribe a 
+plant canopy type. The default case is a homogeneous plant 
+canopy. Heterogeneous distributions of the leaf area 
+density and the canopy drag coefficient can be defined in the 
+new routine user_init_plant_canopy (user_interface). 
+The inipar parameters lad_surface, lad_vertical_gradient and 
+lad_vertical_gradient_level can be used in order to 
+prescribe the vertical profile of leaf area density. The 
+inipar parameter drag_coefficient determines the canopy 
+drag coefficient.
+Finally, the inipar parameter pch_index determines the 
+index of the upper boundary of the plant canopy.  
 
 Allow new case bc_uv_t = 'dirichlet_0' for channel flow.
@@ -15,5 +30,5 @@
 advec_particles after the end of this loop.
 
-advec_particles, check_parameters, data_output_dvrp, header, init_3d_model, init_grid, init_particles, init_pegrid, modules, package_parin, read_var_list, user_interface, write_var_list
+advec_particles, check_parameters, data_output_dvrp, header, init_3d_model, init_grid, init_particles, init_pegrid, modules, package_parin, parin, plant_canopy_model, read_var_list, read_3d_binary, user_interface, write_var_list, write_3d_binary
 
 
Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 137)
+++ palm/trunk/SOURCE/Makefile	(revision 138)
@@ -4,4 +4,6 @@
 # Actual revisions:
 # -----------------
+# +plant_canopy_model
+#
 # +surface_coupler
 #
@@ -59,5 +61,6 @@
 	local_flush.f90 local_getenv.f90 local_stop.f90 local_system.f90 local_tremain.f90 \
 	local_tremain_ini.f90 modules.f90 netcdf.f90 package_parin.f90 \
-	palm.f90 parin.f90 particle_boundary_conds.f90 poisfft.f90 \
+	palm.f90 parin.f90 particle_boundary_conds.f90 \
+        plant_canopy_model.f90 poisfft.f90 \
 	poisfft_hybrid.f90 poismg.f90 prandtl_fluxes.f90 pres.f90 print_1d.f90 \
 	production_e.f90 prognostic_equations.f90 random_function.f90 \
@@ -89,5 +92,6 @@
 	interaction_droplets_ptq.o local_flush.o local_getenv.o local_stop.o \
 	local_system.o local_tremain.o local_tremain_ini.o modules.o netcdf.o \
-	package_parin.o palm.o parin.o particle_boundary_conds.o poisfft.o \
+	package_parin.o palm.o parin.o particle_boundary_conds.o \
+        plant_canopy_model.o poisfft.o \
 	poisfft_hybrid.o poismg.o prandtl_fluxes.o pres.o print_1d.o \
 	production_e.o prognostic_equations.o random_function.o random_gauss.o \
@@ -201,4 +205,5 @@
 parin.o: modules.o
 particle_boundary_conds.o: modules.o
+plant_canopy_model.o: modules.o
 poisfft.o: modules.o fft_xy.o
 poisfft_hybrid.o: modules.o fft_xy.o
@@ -212,5 +217,6 @@
 	buoyancy.o calc_precipitation.o calc_radiation.o coriolis.o \
 	diffusion_e.o diffusion_s.o diffusion_u.o diffusion_v.o diffusion_w.o \
-	eqn_state_seawater.o impact_of_latent_heat.o production_e.o \
+	eqn_state_seawater.o impact_of_latent_heat.o \
+        plant_canopy_model.o production_e.o \
 	user_interface.o
 random_gauss.o: random_function.o
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 137)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 138)
@@ -4,4 +4,5 @@
 ! Actual revisions:
 ! -----------------
+! Plant canopy added
 ! Allow new case bc_uv_t = 'dirichlet_0' for channel flow.
 ! Multigrid solver allows topography, checking of dt_sort_particles
@@ -592,4 +593,10 @@
     ENDIF
 
+    IF ( plant_canopy .AND. ( drag_coefficient == 0.0 ) ) THEN
+       IF ( myid == 0 )  PRINT*, '+++ check_parameters: plant_canopy = TRUE', &
+                                 'requires a non-zero drag_coefficient'
+       CALL local_stop
+    ENDIF 
+
 !
 !-- In case of no model continuation run, check initialising parameters and
@@ -605,4 +612,5 @@
        IF ( ocean )           sa_init = sa_surface
        IF ( passive_scalar )  q_init  = s_surface
+       IF ( plant_canopy )    lad = 0.0
 
 !
@@ -934,6 +942,51 @@
        ENDIF
 
-    ENDIF
-
+!
+!--    If required compute the profile of leaf area density used in the plant canopy model 
+       IF ( plant_canopy ) THEN
+       
+          i = 1
+          gradient = 0.0
+
+          IF ( .NOT. ocean ) THEN
+ 
+             lad_vertical_gradient_level_ind(1) = 0
+             DO k = 1, pch_index
+                IF ( lad_vertical_gradient_level(i) < zu(k) .AND.  &
+                     lad_vertical_gradient_level(i) >= 0.0 ) THEN
+                   gradient = lad_vertical_gradient(i)
+                   lad_vertical_gradient_level_ind(i) = k - 1
+                   i = i + 1
+                   IF ( i > 10 ) THEN
+                      IF ( myid == 0 ) THEN
+                         PRINT*, '+++ user_init_3d_model: upper bound 10 of array',  &
+                                 ' "lad_vertical_gradient_level" exceeded'
+                      ENDIF
+                      CALL local_stop
+                   ENDIF
+                ENDIF
+                IF ( gradient /= 0.0 ) THEN
+                   IF ( k /= 1 ) THEN
+                      lad(k) = lad(k-1) + dzu(k) * gradient
+                   ELSE
+                      lad(k) = lad_surface + 0.5 * dzu(k) *gradient
+                   ENDIF
+                ELSE
+                   lad(k) = lad(k-1)
+                ENDIF
+             ENDDO
+
+          ENDIF
+
+!
+!--       In case of no given leaf area density gradients, choose a vanishing gradient
+          IF ( lad_vertical_gradient_level(1) == -9999999.9 ) THEN
+             lad_vertical_gradient_level(1) = 0.0
+          ENDIF
+
+       ENDIF
+         
+    ENDIF
+             
 !
 !-- Compute Coriolis parameter
Index: palm/trunk/SOURCE/header.f90
===================================================================
--- palm/trunk/SOURCE/header.f90	(revision 137)
+++ palm/trunk/SOURCE/header.f90	(revision 138)
@@ -78,5 +78,5 @@
     CHARACTER (LEN=70) ::  char1, char2, coordinates, gradients, dopr_chr, &
                            do2d_xy, do2d_xz, do2d_yz, do3d_chr, &
-                           run_classification, slices, temperatures, &
+                           run_classification, slices, temperatures, learde, &
                            ugcomponent, vgcomponent
     CHARACTER (LEN=85) ::  roben, runten
@@ -364,4 +364,38 @@
 
     END SELECT
+
+    IF ( plant_canopy ) THEN
+
+       WRITE ( io, 280 ) canopy_mode, pch_index, drag_coefficient
+
+!
+!--    Leaf area density profile
+!--    Building output strings, starting with surface value
+       WRITE ( learde, '(F6.2)' )  lad_surface
+       gradients = '------'
+       slices = '     0'
+       coordinates = '   0.0'
+       i = 1
+       DO  WHILE ( lad_vertical_gradient_level_ind(i) /= -9999 )
+
+          WRITE (coor_chr,'(F7.2)')  lad(lad_vertical_gradient_level_ind(i))
+          learde = TRIM( learde ) // ' ' // TRIM( coor_chr )
+
+          WRITE (coor_chr,'(F7.2)')  lad_vertical_gradient(i)
+          gradients = TRIM( gradients ) // ' ' // TRIM( coor_chr )
+
+          WRITE (coor_chr,'(I7)')  lad_vertical_gradient_level_ind(i)
+          slices = TRIM( slices ) // ' ' // TRIM( coor_chr )
+
+          WRITE (coor_chr,'(F7.1)')  lad_vertical_gradient_level(i)
+          coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
+
+          i = i + 1
+       ENDDO
+
+       WRITE ( io, 281 )  TRIM( coordinates ), TRIM( learde ), &
+                          TRIM( gradients ), TRIM( slices )
+
+    ENDIF
 
 !
@@ -1333,4 +1367,15 @@
               ' Horizontal index bounds (l/r/s/n): ',I4,' / ',I4,' / ',I4, &
                 ' / ',I4)
+280 FORMAT (//' Vegetation canopy (drag) model:'/ &
+              ' ------------------------------'// &
+              ' Canopy mode: ', A / &
+              ' Canopy top: ',I4 / &
+              ' Leaf drag coefficient: ',F6.2 /)
+281 FORMAT (/ ' Characteristic levels of the leaf area density:'// &
+              ' Height:              ',A,'  m'/ &
+              ' Leaf area density:   ',A,'  m**2/m**3'/ &
+              ' Gradient:            ',A,'  m**2/m**4'/ &
+              ' Gridpoint:           ',A)
+               
 300 FORMAT (//' Boundary conditions:'/ &
              ' -------------------'// &
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 137)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 138)
@@ -235,4 +235,14 @@
 
 !
+!-- 3D-arrays for the leaf area density and the canopy drag coefficient
+    IF ( plant_canopy ) THEN
+       ALLOCATE ( lad_s(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
+                  lad_u(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
+                  lad_v(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1),  &
+                  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) )
+    ENDIF
+
+!
 !-- 4D-array for storing the Rif-values at vertical walls
     IF ( topography /= 'flat' )  THEN
@@ -923,4 +933,54 @@
                                                     'problem'
        CALL local_stop
+    ENDIF
+
+!
+!-- Initialization of the leaf area density
+    IF ( plant_canopy ) THEN
+ 
+       SELECT CASE ( TRIM( canopy_mode ) )
+
+          CASE( 'block' )
+
+             DO  i = nxl-1, nxr+1
+                DO  j = nys-1, nyn+1
+                   lad_s(:,j,i) = lad(:)
+                   cdc(:,j,i)   = drag_coefficient
+                ENDDO
+             ENDDO
+
+          CASE DEFAULT
+
+!
+!--          The DEFAULT case is reached either if the parameter 
+!--          canopy mode contains a wrong character string or if the 
+!--          user has coded a special case in the user interface. 
+!--          There, the subroutine user_init_plant_canopy checks 
+!--          which of these two conditions applies.
+             CALL user_init_plant_canopy
+ 
+          END SELECT
+
+       CALL exchange_horiz( lad_s )
+       CALL exchange_horiz( cdc )
+
+       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) )
+             ENDDO
+             DO  k = nzb, nzt
+                lad_w(k,j,i) = 0.5 * ( lad_s(k+1,j,i) + lad_s(k,j,i) )
+             ENDDO
+          ENDDO
+       ENDDO
+
+       lad_w(nzt+1,:,:) = lad_w(nzt,:,:)
+
+       CALL exchange_horiz( lad_u )
+       CALL exchange_horiz( lad_v )
+       CALL exchange_horiz( lad_w )
+  
     ENDIF
 
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 137)
+++ palm/trunk/SOURCE/modules.f90	(revision 138)
@@ -5,4 +5,8 @@
 ! Actual revisions:
 ! -----------------
+! +drag_coefficient, pch_index, lad_surface, lad_vertical_gradient, 
+! lad_vertical_gradient_level, plant_canopy, lad, lad_s, lad_u, lad_v, 
+! lad_w, cdc, lad_vertical_gradient_level_ind, canopy_mode
+!
 ! +dt_sort_particles, ngp_2dh_s_inner, time_sort_particles, flags, 
 ! wall_flags_1..10, wall_humidityflux(0:4), wall_qflux(0:4), 
@@ -107,5 +111,5 @@
 
     REAL, DIMENSION(:), ALLOCATABLE ::                                         &
-          ddzu, dd2zu, dzu, ddzw, dzw, hyp, km_damp_x, km_damp_y, l_grid,      &
+          ddzu, dd2zu, dzu, ddzw, dzw, hyp, km_damp_x, km_damp_y, lad, l_grid, &
           pt_init, q_init, rdf, sa_init, ug, u_init, u_nzb_p1_for_vfc, vg,     &
           v_init, v_nzb_p1_for_vfc, zu, zw
@@ -126,6 +130,7 @@
 
     REAL, DIMENSION(:,:,:), ALLOCATABLE ::                                     &
-          d, diss, 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
+          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
 
     REAL, DIMENSION(:,:,:), ALLOCATABLE, TARGET ::                             &
@@ -259,4 +264,5 @@
                             bc_sa_t = 'neumann', &
                             bc_uv_b = 'dirichlet', bc_uv_t = 'dirichlet', &
+                            canopy_mode = 'block', &
                             coupling_mode = 'uncoupled', &
                             coupling_mode_remote = 'uncoupled', &
@@ -307,4 +313,6 @@
                 do2d_xz_time_count(0:1), do2d_yz_time_count(0:1), &
                 do3d_no(0:1) = 0, do3d_time_count(0:1), &
+                lad_vertical_gradient_level_ind(10) = -9999, &
+                pch_index = 0, &
                 pt_vertical_gradient_level_ind(10) = -9999, &
                 q_vertical_gradient_level_ind(10) = -9999, &
@@ -340,4 +348,5 @@
                 outflow_l = .FALSE., outflow_n = .FALSE., outflow_r = .FALSE., &
                 outflow_s = .FALSE., passive_scalar = .FALSE., &
+                plant_canopy = .FALSE., &
                 prandtl_layer = .TRUE., precipitation = .FALSE., &
                 profil_output = .FALSE., radiation = .FALSE., &
@@ -365,4 +374,5 @@
              disturbance_level_b = -9999999.9, &
              disturbance_level_t = -9999999.9, &
+             drag_coefficient = 0.0, &
              dt = -1.0, dt_averaging_input = 0.0, &
              dt_averaging_input_pr = 9999999.9, dt_coupling = 9999999.9, &
@@ -379,5 +389,5 @@
              end_time = 0.0, &
              f = 0.0, fs = 0.0, g = 9.81, kappa = 0.4, km_constant = -1.0, &
-             km_damp_max = -1.0, long_filter_factor = 0.0, &
+             km_damp_max = -1.0, lad_surface = 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, &
@@ -419,5 +429,8 @@
 
     REAL ::  do2d_xy_last_time(0:1) = -1.0, do2d_xz_last_time(0:1) = -1.0, &
-             do2d_yz_last_time(0:1) = -1.0, pt_vertical_gradient(10) = 0.0, &
+             do2d_yz_last_time(0:1) = -1.0, &
+             lad_vertical_gradient(10) = 0.0, &
+             lad_vertical_gradient_level(10) = -9999999.9, &
+             pt_vertical_gradient(10) = 0.0, &
              pt_vertical_gradient_level(10) = -9999999.9, &
              q_vertical_gradient(10) = 0.0, &
Index: palm/trunk/SOURCE/parin.f90
===================================================================
--- palm/trunk/SOURCE/parin.f90	(revision 137)
+++ palm/trunk/SOURCE/parin.f90	(revision 138)
@@ -4,5 +4,9 @@
 ! Actual revisions:
 ! -----------------
-! 
+! +canopy_mode
+! +drag_coefficient, lad_surface, lad_vertical_gradient, 
+! +lad_vertical_gradient_level, pch_index, plant_canopy 
+! +allocation of leaf area density field
+!
 !
 ! Former revisions:
@@ -70,21 +74,22 @@
                        bc_ns, bc_p_b, bc_p_t, bc_pt_b, bc_pt_t, bc_q_b, &
              bc_q_t,bc_s_b, bc_s_t, bc_sa_t, bc_uv_b, bc_uv_t, &
-             bottom_salinityflux, building_height, &
-             building_length_x, building_length_y, building_wall_left, &
-             building_wall_south, cloud_droplets, cloud_physics, &
+             bottom_salinityflux, building_height, building_length_x, &
+             building_length_y, building_wall_left, &
+             building_wall_south, canopy_mode, cloud_droplets, cloud_physics, &
              conserve_volume_flow, cut_spline_overshoot, damp_level_1d, &
-             dissipation_1d, dt, dt_pr_1d, dt_run_control_1d, dx, dy, dz, &
-             dz_max, dz_stretch_factor, dz_stretch_level, e_init, e_min, &
-             end_time_1d, &
-             fft_method, galilei_transformation, grid_matching, humidity, &
+             dissipation_1d, drag_coefficient, dt, dt_pr_1d, &
+             dt_run_control_1d, dx, dy, dz, dz_max, dz_stretch_factor, & 
+             dz_stretch_level, e_init, e_min, end_time_1d, fft_method, &
+             galilei_transformation, grid_matching, humidity, &
              inflow_disturbance_begin, inflow_disturbance_end, &
-             initializing_actions, km_constant, km_damp_max, &
+             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, npex, npey, nsor_ini, nx, ny, &
              nz, ocean, omega, outflow_damping_width, overshoot_limit_e, &
              overshoot_limit_pt, overshoot_limit_u, overshoot_limit_v, &
-             overshoot_limit_w, passive_scalar, phi, prandtl_layer, &
-             precipitation, pt_reference, pt_surface, &
-             pt_surface_initial_change, pt_vertical_gradient, &
+             overshoot_limit_w, passive_scalar, pch_index, phi, &  
+             plant_canopy, prandtl_layer, precipitation, pt_reference, &
+             pt_surface, pt_surface_initial_change, pt_vertical_gradient, &
              pt_vertical_gradient_level, q_surface, q_surface_initial_change, &
              q_vertical_gradient, q_vertical_gradient_level, radiation, &
@@ -242,4 +247,9 @@
        hom = 0.0
 
+       IF ( plant_canopy ) THEN
+          ALLOCATE( lad(0:nz+1) )
+          lad(:) = 0.0
+       ENDIF
+
     ENDIF
 
Index: palm/trunk/SOURCE/plant_canopy_model.f90
===================================================================
--- palm/trunk/SOURCE/plant_canopy_model.f90	(revision 138)
+++ palm/trunk/SOURCE/plant_canopy_model.f90	(revision 138)
@@ -0,0 +1,263 @@
+ MODULE plant_canopy_model_mod
+
+!------------------------------------------------------------------------------!
+! Actual revisions:
+! -----------------
+! Initial revision
+!
+! Former revisions:
+! -----------------
+! $Id$
+!
+! Description:
+! ------------
+! Evaluation of the drag due to vegetation
+!------------------------------------------------------------------------------!
+
+    PRIVATE
+    PUBLIC plant_canopy_model
+
+    INTERFACE plant_canopy_model
+       MODULE PROCEDURE plant_canopy_model
+       MODULE PROCEDURE plant_canopy_model_ij
+    END INTERFACE plant_canopy_model
+
+ CONTAINS
+
+
+!------------------------------------------------------------------------------!
+! Call for all grid points
+!------------------------------------------------------------------------------!
+    SUBROUTINE plant_canopy_model( component )
+
+       USE arrays_3d
+       USE control_parameters
+       USE indices
+       USE pegrid
+
+       IMPLICIT NONE
+
+       INTEGER ::  component, i, j, k
+ 
+!
+!--    Compute drag for the three velocity components and the SGS-TKE
+       SELECT CASE ( component )
+
+!
+!--       u-component
+          CASE ( 1 )
+             DO  i = nxlu, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_u_inner(j,i)+1, pch_index
+                      tend(k,j,i) = tend(k,j,i) -                &
+                                    cdc(k,j,i) * lad_u(k,j,i) *  &
+                                    SQRT(     u(k,j,i)**2     +  &
+                                          ( ( v(k,j,i-1)      +  &
+                                              v(k,j,i)        +  &
+                                              v(k,j+1,i)      +  &
+                                              v(k,j+1,i+1) )     &
+                                            / 4.0 )**2        +  &
+                                          ( ( w(k-1,j,i-1)    +  &
+                                              w(k-1,j,i)      +  &
+                                              w(k,j,i-1)      +  &
+                                              w(k,j,i) )         &
+                                            / 4.0 )**2 )      *  &
+                                    u(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+
+!
+!--       v-component
+          CASE ( 2 )
+             DO  i = nxl, nxr
+                DO  j = nysv, nyn
+                   DO  k = nzb_v_inner(j,i)+1, pch_index
+                      tend(k,j,i) = tend(k,j,i) -                &
+                                    cdc(k,j,i) * lad_v(k,j,i) *  &
+                                    SQRT( ( ( u(k,j-1,i)      +  &
+                                              u(k,j-1,i+1)    +  &
+                                              u(k,j,i)        +  &
+                                              u(k,j,i+1) )       &
+                                            / 4.0 )**2        +  &
+                                              v(k,j,i)**2     +  &
+                                          ( ( w(k-1,j-1,i)    +  &
+                                              w(k-1,j,i)      +  &
+                                              w(k,j-1,i)      +  &
+                                              w(k,j,i) )         &
+                                            / 4.0 )**2 )      *  &
+                                    v(k,j,i)  
+                   ENDDO
+                ENDDO
+             ENDDO
+
+!
+!--       w-component
+          CASE ( 3 )
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = nzb_w_inner(j,i)+1, pch_index
+                      tend(k,j,i) = tend(k,j,i) -                & 
+                                    cdc(k,j,i) * lad_w(k,j,i) *  &
+                                    SQRT( ( ( u(k,j,i)        +  &
+                                              u(k,j,i+1)      +  &
+                                              u(k+1,j,i)      +  &
+                                              u(k+1,j,i+1) )     &
+                                            / 4.0 )**2        +  &
+                                          ( ( v(k,j,i)        +  &
+                                              v(k,j+1,i)      +  &
+                                              v(k+1,j,i)      +  &
+                                              v(k+1,j+1,i) )     &
+                                            / 4.0 )**2        +  &
+                                              w(k,j,i)**2 )   *  &
+                                    w(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO
+
+!
+!--       sgs-tke
+          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) -                     &
+                                    2.0 * cdc(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,j,i)              + &
+                                              w(k+1,j,i) )            &
+                                            / 2.0 )**2 )            * &
+                                    e(k,j,i)
+                   ENDDO
+                ENDDO
+             ENDDO 
+                        
+          CASE DEFAULT
+
+             IF ( myid == 0 )  PRINT*,'+++ pcm:  wrong component: ', &
+                                      component
+             CALL local_stop
+
+       END SELECT
+
+    END SUBROUTINE plant_canopy_model
+
+
+!------------------------------------------------------------------------------!
+! Call for grid point i,j
+!------------------------------------------------------------------------------!
+    SUBROUTINE plant_canopy_model_ij( i, j, component )
+
+       USE arrays_3d
+       USE control_parameters
+       USE indices
+       USE pegrid
+
+       IMPLICIT NONE
+
+       INTEGER ::  component, i, j, k
+
+       IF ( drag_coefficient /= 0.0 ) THEN
+
+!
+!--       Compute drag for the three velocity components 
+          SELECT CASE ( component )
+
+!
+!--          u-component
+             CASE ( 1 )
+                DO  k = nzb_u_inner(j,i)+1, pch_index
+                   tend(k,j,i) = tend(k,j,i) -                  &
+                                 cdc(k,j,i) * lad_u(k,j,i) *    &   
+                                 SQRT(     u(k,j,i)**2 +        &
+                                       ( ( v(k,j,i-1)  +        &
+                                           v(k,j,i)    +        &
+                                           v(k,j+1,i)  +        &
+                                           v(k,j+1,i+1) )       &
+                                         / 4.0 )**2    +        &
+                                       ( ( w(k-1,j,i-1) +       &
+                                           w(k-1,j,i)   +       &
+                                           w(k,j,i-1)   +       &
+                                           w(k,j,i) )           &
+                                         / 4.0 )**2 ) *         &
+                                 u(k,j,i)
+                ENDDO
+
+!
+!--          v-component
+             CASE ( 2 )
+                DO  k = nzb_v_inner(j,i)+1, pch_index
+                   tend(k,j,i) = tend(k,j,i) -                  &
+                                 cdc(k,j,i) * lad_v(k,j,i) *    &
+                                 SQRT( ( ( u(k,j-1,i)   +       &
+                                           u(k,j-1,i+1) +       &
+                                           u(k,j,i)     +       &
+                                           u(k,j,i+1) )         &
+                                         / 4.0 )**2     +       &
+                                           v(k,j,i)**2  +       &
+                                       ( ( w(k-1,j-1,i) +       &
+                                           w(k-1,j,i)   +       &
+                                           w(k,j-1,i)   +       &
+                                           w(k,j,i) )           &
+                                         / 4.0 )**2 ) *         &
+                                 v(k,j,i)
+                ENDDO
+
+!
+!--          w-component
+             CASE ( 3 )
+                DO  k = nzb_w_inner(j,i)+1, pch_index
+                   tend(k,j,i) = tend(k,j,i) -                  &
+                                 cdc(k,j,i) * lad_w(k,j,i) *    & 
+                                 SQRT( ( ( u(k,j,i)    +        &  
+                                           u(k,j,i+1)  +        &
+                                           u(k+1,j,i)  +        &
+                                           u(k+1,j,i+1) )       &
+                                         / 4.0 )**2    +        &
+                                       ( ( v(k,j,i)    +        &
+                                           v(k,j+1,i)  +        &
+                                           v(k+1,j,i)  +        &
+                                           v(k+1,j+1,i) )       &
+                                         / 4.0 )**2    +        &
+                                           w(k,j,i)**2 ) *      &
+                                 w(k,j,i)
+   
+                ENDDO
+
+!
+!--          sgs-tke
+             CASE ( 4 )
+                DO  k = nzb_s_inner(j,i)+1, pch_index   
+                   tend(k,j,i) = tend(k,j,i) -                     &
+                                 2.0 * cdc(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,j,i)           +    &
+                                           w(k+1,j,i) )            &
+                                         / 2.0 )**2 )         *    &
+                                 e(k,j,i)
+
+                ENDDO
+
+             CASE DEFAULT
+
+                IF ( myid == 0 )  PRINT*,'+++ pcm:  wrong component: ', &
+                                         component
+                CALL local_stop
+
+          END SELECT
+
+       ENDIF
+
+    END SUBROUTINE plant_canopy_model_ij
+
+ END MODULE plant_canopy_model_mod
Index: palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 137)
+++ palm/trunk/SOURCE/prognostic_equations.f90	(revision 138)
@@ -4,4 +4,5 @@
 ! Actual revisions:
 ! -----------------
+! add call of subroutines that evaluate the canopy drag terms
 ! add wall_*flux to parameter list of calls of diffusion_s
 !
@@ -80,4 +81,5 @@
     USE diffusion_w_mod
     USE impact_of_latent_heat_mod
+    USE plant_canopy_model_mod
     USE production_e_mod
     USE user_actions_mod
@@ -160,4 +162,8 @@
           CALL coriolis( i, j, 1 )
           IF ( sloping_surface )  CALL buoyancy( i, j, pt, pt_reference, 1, 4 )
+
+!
+!--       Drag by plant canopy
+          IF ( plant_canopy )  CALL plant_canopy_model( i, j, 1 )
           CALL user_actions( i, j, 'u-tendency' )
 
@@ -227,4 +233,9 @@
           ENDIF
           CALL coriolis( i, j, 2 )
+
+!
+!--       Drag by plant canopy
+          IF ( plant_canopy )  CALL plant_canopy_model( i, j, 2 )     
+
           CALL user_actions( i, j, 'v-tendency' )
 
@@ -303,4 +314,9 @@
              ENDIF
           ENDIF
+
+!
+!--       Drag by plant canopy
+          IF ( plant_canopy )  CALL plant_canopy_model( i, j, 3 )
+
           CALL user_actions( i, j, 'w-tendency' )
 
@@ -732,4 +748,9 @@
              ENDIF
              CALL production_e( i, j )
+
+!
+!--          Additional sink term for flows through plant canopies
+             IF ( plant_canopy )  CALL plant_canopy_model( i, j, 4 )         
+
              CALL user_actions( i, j, 'e-tendency' )
 
@@ -833,4 +854,9 @@
              IF ( sloping_surface )  CALL buoyancy( i, j, pt, pt_reference, 1, &
                                                     4 )
+
+!
+!--          Drag by plant canopy
+             IF ( plant_canopy )  CALL plant_canopy_model( i, j, 1 )
+
              CALL user_actions( i, j, 'u-tendency' )
 
@@ -882,4 +908,9 @@
              ENDIF
              CALL coriolis( i, j, 2 )
+
+!
+!--          Drag by plant canopy
+             IF ( plant_canopy )  CALL plant_canopy_model( i, j, 2 )        
+
              CALL user_actions( i, j, 'v-tendency' )
 
@@ -938,4 +969,9 @@
              ENDIF
           ENDIF
+
+!
+!--       Drag by plant canopy
+          IF ( plant_canopy )  CALL plant_canopy_model( i, j, 3 )
+
           CALL user_actions( i, j, 'w-tendency' )
 
@@ -1177,4 +1213,9 @@
              ENDIF
              CALL production_e( i, j )
+
+!
+!--          Additional sink term for flows through plant canopies
+             IF ( plant_canopy )  CALL plant_canopy_model( i, j, 4 ) 
+
              CALL user_actions( i, j, 'e-tendency' )
 
@@ -1269,4 +1310,9 @@
     CALL coriolis( 1 )
     IF ( sloping_surface )  CALL buoyancy( pt, pt_reference, 1, 4 )
+
+!
+!-- Drag by plant canopy
+    IF ( plant_canopy )  CALL plant_canopy_model( 1 )
+
     CALL user_actions( 'u-tendency' )
 
@@ -1340,4 +1386,8 @@
     ENDIF
     CALL coriolis( 2 )
+
+!
+!-- Drag by plant canopy
+    IF ( plant_canopy )  CALL plant_canopy_model( 2 )
     CALL user_actions( 'v-tendency' )
 
@@ -1420,4 +1470,9 @@
        ENDIF
     ENDIF
+
+!
+!-- Drag by plant canopy
+    IF ( plant_canopy )  CALL plant_canopy_model( 3 )
+
     CALL user_actions( 'w-tendency' )
 
@@ -1858,4 +1913,8 @@
        ENDIF
        CALL production_e
+
+!
+!--    Additional sink term for flows through plant canopies
+       IF ( plant_canopy )  CALL plant_canopy_model( 4 )
        CALL user_actions( 'e-tendency' )
 
Index: palm/trunk/SOURCE/read_var_list.f90
===================================================================
--- palm/trunk/SOURCE/read_var_list.f90	(revision 137)
+++ palm/trunk/SOURCE/read_var_list.f90	(revision 138)
@@ -4,4 +4,9 @@
 ! Actual revisions:
 ! -----------------
+! +canopy_mode
+! +drag_coefficient, lad, lad_surface, lad_vertical_gradient, 
+! lad_vertical_gradient_level, lad_vertical_gradient_level_ind,
+! pch_index, plant_canopy
+!
 ! +time_sort_particles
 !
@@ -186,4 +191,6 @@
           CASE ( 'building_wall_south' )
              READ ( 13 )  building_wall_south
+          CASE ( 'canopy_mode' )
+             READ ( 13 )  canopy_mode
           CASE ( 'cloud_droplets' )
              READ ( 13 )  cloud_droplets
@@ -200,4 +207,6 @@
           CASE ( 'dissipation_1d' )
              READ ( 13 )  dissipation_1d
+          CASE ( 'drag_coefficient' )
+             READ ( 13 )  drag_coefficient
           CASE ( 'dt_fixed' )
              READ ( 13 )  dt_fixed
@@ -244,4 +253,14 @@
           CASE ( 'km_damp_max' )
              READ ( 13 )  km_damp_max
+          CASE ( 'lad' )
+             READ ( 13 )  lad
+          CASE ( 'lad_surface' )
+             READ ( 13 )  lad_surface
+          CASE ( 'lad_vertical_gradient' )
+             READ ( 13 )  lad_vertical_gradient
+          CASE ( 'lad_vertical_gradient_level' )
+             READ ( 13 )  lad_vertical_gradient_level
+          CASE ( 'lad_vertical_gradient_level_in' )
+             READ ( 13 )  lad_vertical_gradient_level_ind
           CASE ( 'last_dt_change' )
              READ ( 13 )  last_dt_change
@@ -288,6 +307,10 @@
           CASE ( 'passive_scalar' )
              READ ( 13 )  passive_scalar
+          CASE ( 'pch_index' )
+             READ ( 13 )  pch_index
           CASE ( 'phi' )
              READ ( 13 )  phi
+          CASE ( 'plant_canopy' )
+             READ ( 13 )  plant_canopy
           CASE ( 'prandtl_layer' )
              READ ( 13 )  prandtl_layer
Index: palm/trunk/SOURCE/user_interface.f90
===================================================================
--- palm/trunk/SOURCE/user_interface.f90	(revision 137)
+++ palm/trunk/SOURCE/user_interface.f90	(revision 138)
@@ -4,4 +4,6 @@
 ! Actual revisions:
 ! -----------------
+! new subroutine user_init_plant_canopy
+!
 ! new subroutine user_data_output_dvrp
 ! +argument gls in user_init_grid
@@ -293,4 +295,73 @@
 
  END SUBROUTINE user_init_grid
+
+
+
+ SUBROUTINE user_init_plant_canopy
+
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+! Initialisation of the leaf area density array (for scalar grid points) and 
+! the array of the canopy drag coefficient, if the user has not chosen any 
+! of the default cases 
+!------------------------------------------------------------------------------!
+
+    USE arrays_3d
+    USE control_parameters
+    USE indices
+    USE user
+
+    IMPLICIT NONE
+
+    INTEGER :: i, j 
+
+!
+!-- Here the user-defined grid initializing actions follow:
+
+!
+!-- Set the 3D-arrays lad_s and cdc for user defined canopies
+    SELECT CASE ( TRIM( canopy_mode ) )
+
+       CASE ( 'block' )
+!
+!--       Not allowed here since this is the standard case used in init_3d_model.
+
+       CASE ( 'user_defined_canopy_1' )
+!
+!--       Here the user can define his own topography. 
+!--       The following lines contain an example in that the 
+!--       plant canopy extends only over the second half of the 
+!--       model domain along x
+!          DO  i = nxl-1, nxr+1
+!             IF ( i >= INT(nx+1/2) ) THEN
+!                DO  j = nys-1, nyn+1
+!                   lad_s(:,j,i) = lad(:)
+!                   cdc(:,j,i)   = drag_coefficient
+!                ENDDO
+!             ELSE
+!                lad_s(:,:,i) = 0.0
+!                cdc(:,:,i)   = 0.0
+!             ENDIF
+!          ENDDO 
+!--       After definition, please
+!--       remove the following three lines!
+          PRINT*, '+++ user_init_plant_canopy: canopy_mode "', &
+               canopy_mode, '" not available yet'
+
+       CASE DEFAULT
+!
+!--       The DEFAULT case is reached if the parameter canopy_mode contains a
+!--       wrong character string that is neither recognized in init_3d_model nor 
+!--       here in user_init_plant_canopy.
+          PRINT*, '+++ user_init_plant_canopy: unknown canopy_mode "', &
+               canopy_mode, '"'
+          CALL local_stop
+
+    END SELECT
+
+
+ END SUBROUTINE user_init_plant_canopy
 
 
Index: palm/trunk/SOURCE/write_var_list.f90
===================================================================
--- palm/trunk/SOURCE/write_var_list.f90	(revision 137)
+++ palm/trunk/SOURCE/write_var_list.f90	(revision 138)
@@ -4,4 +4,9 @@
 ! Actual revisions:
 ! -----------------
+! +canopy_mode
+! +drag_coefficient, lad, lad_surface, lad_vertical_gradient, 
+! lad_vertical_gradient_level, lad_vertical_gradient_level_ind, 
+! pch_index, plant_canopy
+!
 ! +time_sort_particles
 !
@@ -132,4 +137,6 @@
     WRITE ( 14 )  'building_wall_south           '
     WRITE ( 14 )  building_wall_south
+    WRITE ( 14 )  'canopy_mode                   '
+    WRITE ( 14 )  canopy_mode
     WRITE ( 14 )  'cloud_droplets                '
     WRITE ( 14 )  cloud_droplets
@@ -146,4 +153,6 @@
     WRITE ( 14 )  'dissipation_1d                '
     WRITE ( 14 )  dissipation_1d
+    WRITE ( 14 )  'drag_coefficient              '
+    WRITE ( 14 )  drag_coefficient
     WRITE ( 14 )  'dt_fixed                      '
     WRITE ( 14 )  dt_fixed
@@ -190,4 +199,14 @@
     WRITE ( 14 )  'km_damp_max                   '
     WRITE ( 14 )  km_damp_max
+    WRITE ( 14 )  'lad                           '
+    WRITE ( 14 )  lad
+    WRITE ( 14 )  'lad_surface                   '
+    WRITE ( 14 )  lad_surface
+    WRITE ( 14 )  'lad_vertical_gradient         '
+    WRITE ( 14 )  lad_vertical_gradient
+    WRITE ( 14 )  'lad_vertical_gradient_level   '
+    WRITE ( 14 )  lad_vertical_gradient_level
+    WRITE ( 14 )  'lad_vertical_gradient_level_in'
+    WRITE ( 14 )  lad_vertical_gradient_level_ind
     WRITE ( 14 )  'last_dt_change                '
     WRITE ( 14 )  last_dt_change
@@ -234,6 +253,10 @@
     WRITE ( 14 )  'passive_scalar                '
     WRITE ( 14 )  passive_scalar
+    WRITE ( 14 )  'pch_index                     '
+    WRITE ( 14 )  pch_index
     WRITE ( 14 )  'phi                           '
     WRITE ( 14 )  phi
+    WRITE ( 14 )  'plant_canopy                  '
+    WRITE ( 14 )  plant_canopy
     WRITE ( 14 )  'prandtl_layer                 '
     WRITE ( 14 )  prandtl_layer
