Index: palm/trunk/SOURCE/Makefile
===================================================================
--- palm/trunk/SOURCE/Makefile	(revision 1483)
+++ palm/trunk/SOURCE/Makefile	(revision 1484)
@@ -20,5 +20,6 @@
 # Current revisions:
 # ------------------
-# 
+# plant_canopy_model-dependency added for check_parameters, header, init_3d_model,
+# package_parin, read_var_list, user_init_plant_canopy, write_var_list
 # 
 # Former revisions:
@@ -288,5 +289,5 @@
 check_for_restart.o: modules.o mod_kinds.o
 check_open.o: modules.o mod_kinds.o mod_particle_attributes.o
-check_parameters.o: modules.o mod_kinds.o subsidence.o
+check_parameters.o: modules.o mod_kinds.o subsidence.o plant_canopy_model.o
 close_file.o: modules.o mod_kinds.o
 compute_vpt.o: modules.o mod_kinds.o
@@ -318,10 +319,10 @@
 flow_statistics.o: modules.o cpulog.o mod_kinds.o
 global_min_max.o: modules.o mod_kinds.o
-header.o: modules.o cpulog.o mod_kinds.o subsidence.o
+header.o: modules.o cpulog.o mod_kinds.o plant_canopy_model.o subsidence.o
 impact_of_latent_heat.o: modules.o mod_kinds.o
 inflow_turbulence.o: modules.o cpulog.o mod_kinds.o
 init_1d_model.o: modules.o mod_kinds.o
 init_3d_model.o: modules.o cpulog.o mod_kinds.o random_function.o advec_ws.o \
-	ls_forcing.o lpm_init.o random_generator_parallel.o
+	ls_forcing.o lpm_init.o plant_canopy_model.o random_generator_parallel.o
 init_advec.o: modules.o mod_kinds.o
 init_cloud_physics.o: modules.o mod_kinds.o
@@ -379,5 +380,5 @@
 netcdf.o: modules.o mod_kinds.o
 nudging.o: modules.o cpulog.o mod_kinds.o
-package_parin.o: modules.o mod_kinds.o
+package_parin.o: modules.o mod_kinds.o plant_canopy_model.o
 palm.o: modules.o cpulog.o ls_forcing.o mod_kinds.o nudging.o
 parin.o: modules.o cpulog.o mod_kinds.o progress_bar.o
@@ -400,5 +401,5 @@
 random_generator_parallel.o: mod_kinds.o
 read_3d_binary.o: modules.o cpulog.o mod_kinds.o random_function.o random_generator_parallel.o
-read_var_list.o: modules.o mod_kinds.o
+read_var_list.o: modules.o mod_kinds.o plant_canopy_model.o
 run_control.o: modules.o cpulog.o mod_kinds.o
 set_slicer_attributes_dvrp.o: modules.o mod_kinds.o
@@ -434,5 +435,5 @@
 user_init_3d_model.o: modules.o mod_kinds.o user_module.o
 user_init_grid.o: modules.o mod_kinds.o user_module.o
-user_init_plant_canopy.o: modules.o mod_kinds.o user_module.o
+user_init_plant_canopy.o: modules.o mod_kinds.o user_module.o plant_canopy_model.o
 user_last_actions.o: modules.o mod_kinds.o user_module.o
 user_lpm_advec.o: modules.o mod_kinds.o user_module.o
@@ -446,3 +447,3 @@
 wall_fluxes.o: modules.o mod_kinds.o
 write_3d_binary.o: modules.o cpulog.o mod_kinds.o random_function.o random_generator_parallel.o
-write_var_list.o: modules.o mod_kinds.o
+write_var_list.o: modules.o mod_kinds.o plant_canopy_model.o
Index: palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- palm/trunk/SOURCE/check_parameters.f90	(revision 1483)
+++ palm/trunk/SOURCE/check_parameters.f90	(revision 1484)
@@ -20,5 +20,12 @@
 ! Current revisions:
 ! -----------------
-! 
+! Changes due to new module structure of the plant canopy model:
+!   module plant_canopy_model_mod added,
+!   checks regarding usage of new method for leaf area density profile 
+!   construction added,
+!   lad-profile construction moved to new subroutine init_plant_canopy within 
+!   the module plant_canopy_model_mod,
+!   drag_coefficient renamed to canopy_drag_coeff.
+! Missing KIND-attribute for REAL constant added
 ! 
 ! Former revisions:
@@ -246,4 +253,5 @@
     USE particle_attributes
     USE pegrid
+    USE plant_canopy_model_mod
     USE profil_parameter
     USE spectrum
@@ -891,15 +899,44 @@
     ENDIF
 
-    IF ( plant_canopy .AND. ( drag_coefficient == 0.0_wp ) ) THEN
-       message_string = 'plant_canopy = .TRUE. requires a non-zero drag ' // &
-                        'coefficient & given value is drag_coefficient = 0.0'
-       CALL message( 'check_parameters', 'PA0041', 1, 2, 0, 6, 0 )
-    ENDIF 
-
-    IF ( plant_canopy  .AND.  cloud_physics  .AND.  icloud_scheme == 0 ) THEN
-       message_string = 'plant_canopy = .TRUE. requires cloud_scheme /=' //  &
-                        ' seifert_beheng'
-       CALL message( 'check_parameters', 'PA0360', 1, 2, 0, 6, 0 )
-    ENDIF 
+    IF ( plant_canopy )  THEN
+    
+       IF ( canopy_drag_coeff == 0.0_wp )  THEN
+          message_string = 'plant_canopy = .TRUE. requires a non-zero drag '// &
+                           'coefficient & given value is canopy_drag_coeff = 0.0'
+          CALL message( 'check_parameters', 'PA0041', 1, 2, 0, 6, 0 )
+       ENDIF
+    
+       IF ( ( alpha_lad /= 9999999.9_wp  .AND.  beta_lad == 9999999.9_wp )  .OR.&
+              beta_lad /= 9999999.9_wp   .AND.  alpha_lad == 9999999.9_wp )  THEN
+          message_string = 'using the beta function for the construction ' //  &
+                           'of the leaf area density profile requires '    //  &
+                           'both alpha_lad and beta_lad to be /= 9999999.9'
+          CALL message( 'check_parameters', 'PA0118', 1, 2, 0, 6, 0 )
+       ENDIF
+    
+       IF ( calc_beta_lad_profile  .AND.  lai_beta == 0.0_wp )  THEN
+          message_string = 'using the beta function for the construction ' //  &
+                           'of the leaf area density profile requires '    //  &
+                           'a non-zero lai_beta, but given value is '      //  &
+                           'lai_beta = 0.0'
+          CALL message( 'check_parameters', 'PA0119', 1, 2, 0, 6, 0 )
+       ENDIF
+
+       IF ( calc_beta_lad_profile  .AND.  lad_surface /= 0.0_wp )  THEN
+          message_string = 'simultaneous setting of alpha_lad /= 9999999.9' // &
+                           'and lad_surface /= 0.0 is not possible, '       // &
+                           'use either vertical gradients or the beta '     // &
+                           'function for the construction of the leaf area '// &
+                           'density profile'
+          CALL message( 'check_parameters', 'PA0120', 1, 2, 0, 6, 0 )
+       ENDIF 
+
+       IF ( cloud_physics  .AND.  icloud_scheme == 0 )  THEN
+          message_string = 'plant_canopy = .TRUE. requires cloud_scheme /=' // &
+                           ' seifert_beheng'
+          CALL message( 'check_parameters', 'PA0360', 1, 2, 0, 6, 0 )
+       ENDIF 
+
+    ENDIF
 
     IF ( .NOT. ( loop_optimization == 'cache'  .OR.                            &
@@ -924,5 +961,4 @@
        IF ( ocean )           sa_init = sa_surface
        IF ( passive_scalar )  q_init  = s_surface
-       IF ( plant_canopy )    lad = 0.0_wp
 
 !
@@ -1264,47 +1300,4 @@
        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_wp
-
-          IF ( .NOT. ocean ) THEN
-
-             lad(0) = lad_surface
- 
-             lad_vertical_gradient_level_ind(1) = 0
-             DO k = 1, pch_index
-                IF ( i < 11 ) THEN
-                   IF ( lad_vertical_gradient_level(i) < zu(k) .AND.  &
-                        lad_vertical_gradient_level(i) >= 0.0_wp ) THEN
-                      gradient = lad_vertical_gradient(i)
-                      lad_vertical_gradient_level_ind(i) = k - 1
-                      i = i + 1
-                   ENDIF
-                ENDIF
-                IF ( gradient /= 0.0_wp ) THEN
-                   IF ( k /= 1 ) THEN
-                      lad(k) = lad(k-1) + dzu(k) * gradient
-                   ELSE
-                      lad(k) = lad_surface + 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_wp ) THEN
-             lad_vertical_gradient_level(1) = 0.0_wp
-          ENDIF
-
-       ENDIF
          
     ENDIF
@@ -1631,5 +1624,5 @@
        ENDIF
 
-       IF ( top_heatflux /= 0.0  .AND.  top_heatflux /= 9999999.9_wp ) &
+       IF ( top_heatflux /= 0.0_wp  .AND.  top_heatflux /= 9999999.9_wp ) &
        THEN
           message_string = 'heatflux must not be set for pure neutral flow'
@@ -1828,5 +1821,5 @@
        rayleigh_damping_factor = 0.0_wp
     ELSE
-       IF ( rayleigh_damping_factor < 0.0 .OR. rayleigh_damping_factor > 1.0_wp ) &
+       IF ( rayleigh_damping_factor < 0.0_wp .OR. rayleigh_damping_factor > 1.0_wp ) &
        THEN
           WRITE( message_string, * )  'rayleigh_damping_factor = ', &
Index: palm/trunk/SOURCE/header.f90
===================================================================
--- palm/trunk/SOURCE/header.f90	(revision 1483)
+++ palm/trunk/SOURCE/header.f90	(revision 1484)
@@ -20,5 +20,12 @@
 ! Current revisions:
 ! -----------------
-! 
+! Changes due to new module structure of the plant canopy model:
+!   module plant_canopy_model_mod and output for new canopy model parameters 
+!   (alpha_lad, beta_lad, lai_beta,...) added,
+!   drag_coefficient, leaf_surface_concentration and scalar_exchange_coefficient
+!   renamed to canopy_drag_coeff, leaf_surface_conc and leaf_scalar_exch_coeff,
+!   learde renamed leaf_area_density.
+! Bugfix: DO-WHILE-loop for lad header information additionally restricted 
+! by maximum number of gradient levels (currently 10)
 !
 ! Former revisions:
@@ -165,5 +172,5 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  lad, pt_init, qsws, q_init, sa_init, shf, ug, vg, w_subs, zu
+        ONLY:  pt_init, qsws, q_init, sa_init, shf, ug, vg, w_subs, zu
         
     USE control_parameters
@@ -205,4 +212,11 @@
         
     USE pegrid
+
+    USE plant_canopy_model_mod,                                                &
+        ONLY:  alpha_lad, beta_lad, calc_beta_lad_profile, canopy_drag_coeff,  &
+               canopy_mode, cthf, lad, lad_surface, lad_vertical_gradient,     &
+               lad_vertical_gradient_level, lad_vertical_gradient_level_ind,   &
+               lai_beta, leaf_scalar_exch_coeff, leaf_surface_conc, pch_index, &
+               plant_canopy
     
     USE spectrum,                                                              &
@@ -242,5 +256,5 @@
     CHARACTER (LEN=86) ::  coordinates         !:
     CHARACTER (LEN=86) ::  gradients           !:
-    CHARACTER (LEN=86) ::  learde              !:
+    CHARACTER (LEN=86) ::  leaf_area_density   !:
     CHARACTER (LEN=86) ::  slices              !:
     CHARACTER (LEN=86) ::  temperatures        !:
@@ -270,8 +284,10 @@
     INTEGER(iwp) ::  io        !:
     INTEGER(iwp) ::  j         !:
+    INTEGER(iwp) ::  k         !:
     INTEGER(iwp) ::  l         !:
     INTEGER(iwp) ::  ll        !:
     INTEGER(iwp) ::  mpi_type  !:
     
+    REAL(wp) ::  canopy_height                    !: canopy height (in m)
     REAL(wp) ::  cpuseconds_per_simulated_second  !:
 
@@ -760,44 +776,71 @@
 
     IF ( plant_canopy )  THEN
-
-       WRITE ( io, 280 ) canopy_mode, pch_index, drag_coefficient
+    
+       canopy_height = pch_index * dz
+
+       WRITE ( io, 280 )  canopy_mode, canopy_height, pch_index,               &
+                          canopy_drag_coeff
        IF ( passive_scalar )  THEN
-          WRITE ( io, 281 ) scalar_exchange_coefficient,   &
-                            leaf_surface_concentration
+          WRITE ( io, 281 )  leaf_scalar_exch_coeff,                           &
+                             leaf_surface_conc
        ENDIF
 
 !
 !--    Heat flux at the top of vegetation
-       WRITE ( io, 282 ) cthf
-
-!
-!--    Leaf area density profile
-!--    Building output strings, starting with surface value
-       WRITE ( learde, '(F6.4)' )  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, 283 )  TRIM( coordinates ), TRIM( learde ), &
-                          TRIM( gradients ), TRIM( slices )
-
-    ENDIF
+       WRITE ( io, 282 )  cthf
+
+!
+!--    Leaf area density profile, calculated either from given vertical 
+!--    gradients or from beta probability density function.
+       IF (  .NOT.  calc_beta_lad_profile )  THEN
+
+!--       Building output strings, starting with surface value
+          WRITE ( leaf_area_density, '(F6.4)' )  lad_surface
+          gradients = '------'
+          slices = '     0'
+          coordinates = '   0.0'
+          i = 1
+          DO  WHILE ( i < 11  .AND.  lad_vertical_gradient_level_ind(i) /= -9999 )
+
+             WRITE (coor_chr,'(F7.2)')  lad(lad_vertical_gradient_level_ind(i))
+             leaf_area_density = TRIM( leaf_area_density ) // ' ' // 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, 283 )  TRIM( coordinates ), TRIM( leaf_area_density ),              &
+                             TRIM( gradients ), TRIM( slices )
+
+       ELSE
+       
+          WRITE ( leaf_area_density, '(F6.4)' )  lad_surface
+          coordinates = '   0.0'
+          
+          DO  k = 1, pch_index
+
+             WRITE (coor_chr,'(F7.2)')  lad(k)
+             leaf_area_density = TRIM( leaf_area_density ) // ' ' // TRIM( coor_chr )
+ 
+             WRITE (coor_chr,'(F7.1)')  zu(k)
+             coordinates = TRIM( coordinates ) // ' '  // TRIM( coor_chr )
+
+          ENDDO       
+
+          WRITE ( io, 284 ) TRIM( coordinates ), TRIM( leaf_area_density ), alpha_lad,    &
+                            beta_lad, lai_beta
+
+       ENDIF  
+
+    ENDIF
+
 
 !
@@ -1836,7 +1879,7 @@
               ' ------------------------------'// &
               ' Canopy mode: ', A / &
-              ' Canopy top: ',I4 / &
+              ' Canopy height: ',F6.2,'m (',I4,' grid points)' / &
               ' Leaf drag coefficient: ',F6.2 /)
-281 FORMAT (/ ' Scalar_exchange_coefficient: ',F6.2 / &
+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')
@@ -1846,4 +1889,10 @@
               ' Gradient:            ',A,'  m**2/m**4'/ &
               ' Gridpoint:           ',A)
+284 FORMAT (//' Characteristic levels of the leaf area density and coefficients:'// &
+              ' Height:              ',A,'  m'/ &
+              ' Leaf area density:   ',A,'  m**2/m**3'/ &
+              ' Coefficient alpha: ',F6.2 / &
+              ' Coefficient beta: ',F6.2 / &
+              ' Leaf area index: ',F6.2,'  m**2/m**2' /)
                
 300 FORMAT (//' Boundary conditions:'/ &
Index: palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- palm/trunk/SOURCE/init_3d_model.f90	(revision 1483)
+++ palm/trunk/SOURCE/init_3d_model.f90	(revision 1484)
@@ -20,5 +20,8 @@
 ! Current revisions:
 ! ------------------
-! 
+! Changes due to new module structure of the plant canopy model:
+!   canopy-related initialization (e.g. lad and canopy_heat_flux) moved to new 
+!   subroutine init_plant_canopy within the module plant_canopy_model_mod,
+!   call of subroutine init_plant_canopy added.
 ! 
 ! Former revisions:
@@ -217,4 +220,7 @@
     USE pegrid
     
+    USE plant_canopy_model_mod,                                                &
+        ONLY:  init_plant_canopy, plant_canopy
+    
     USE random_function_mod 
     
@@ -517,25 +523,4 @@
     
 !
-!-- 3D-arrays for the leaf area density and the canopy drag coefficient
-    IF ( plant_canopy ) THEN
-       ALLOCATE ( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
-                  lad_u(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
-                  lad_v(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
-                  lad_w(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
-                  cdc(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-
-       IF ( passive_scalar ) THEN
-          ALLOCATE ( sls(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
-                     sec(nzb:nzt+1,nysg:nyng,nxlg:nxrg) ) 
-       ENDIF
-
-       IF ( cthf /= 0.0_wp ) THEN
-          ALLOCATE ( lai(nzb:nzt+1,nysg:nyng,nxlg:nxrg),  &
-                     canopy_heat_flux(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-       ENDIF
-
-    ENDIF
-
-!
 !-- 4D-array for storing the Rif-values at vertical walls
     IF ( topography /= 'flat' )  THEN
@@ -1571,113 +1556,6 @@
 
 !
-!-- Initialization of the leaf area density
-    IF ( plant_canopy )  THEN
- 
-       SELECT CASE ( TRIM( canopy_mode ) )
-
-          CASE( 'block' )
-
-             DO  i = nxlg, nxrg
-                DO  j = nysg, nyng
-                   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
-
-          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, nbgp )
-       CALL exchange_horiz( cdc, nbgp )
-
-       IF ( passive_scalar )  THEN
-          CALL exchange_horiz( sls, nbgp )
-          CALL exchange_horiz( sec, nbgp )
-       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 
-                IF ( lad_s(k,j,i) > 0.0_wp )  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
-                lad_w(k,j,i) = 0.5_wp * ( lad_s(k+1,j,i) + lad_s(k,j,i) )
-             ENDDO
-          ENDDO
-       ENDDO
-
-       lad_w(pch_index,:,:) = 0.0_wp
-       lad_w(nzt+1,:,:)     = lad_w(nzt,:,:)
-
-       CALL exchange_horiz( lad_u, nbgp )
-       CALL exchange_horiz( lad_v, nbgp )
-       CALL exchange_horiz( lad_w, nbgp )
-
-!
-!--    Initialisation of the canopy heat source distribution
-       IF ( cthf /= 0.0_wp )  THEN
-!
-!--       Piecewise evaluation of the leaf area index by 
-!--       integration of the leaf area density
-          lai(:,:,:) = 0.0_wp
-          DO  i = nxlg, nxrg
-             DO  j = nysg, nyng
-                DO  k = pch_index-1, 0, -1
-                   lai(k,j,i) = lai(k+1,j,i) +                   &
-                                ( 0.5_wp * ( lad_w(k+1,j,i) +    &
-                                          lad_s(k+1,j,i) ) *     &
-                                  ( zw(k+1) - zu(k+1) ) )  +     &
-                                ( 0.5_wp * ( 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 = nxlg, nxrg
-             DO  j = nysg, nyng
-                DO  k = 0, pch_index
-                   canopy_heat_flux(k,j,i) = cthf *                    &
-                                             exp( -0.6_wp * 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,:,:)
-
-       ENDIF
-
-    ENDIF
+!-- If required, initialize quantities needed for the plant canopy model
+    IF ( plant_canopy )  CALL init_plant_canopy
 
 !
Index: palm/trunk/SOURCE/modules.f90
===================================================================
--- palm/trunk/SOURCE/modules.f90	(revision 1483)
+++ palm/trunk/SOURCE/modules.f90	(revision 1484)
@@ -20,5 +20,7 @@
 ! Current revisions:
 ! ------------------
-! 
+! Changes due to new module structure of the plant canopy model:
+!   canopy-model related parameters/variables moved to module 
+!   plant_canopy_model_mod
 ! 
 ! Former revisions:
@@ -38,5 +40,5 @@
 ! 1429 2014-07-15 12:53:45Z knoop
 ! +ensemble_member_nr to prepare the random_generator for ensemble runs
-! 
+!
 ! 1382 2014-04-30 12:15:41Z boeske
 ! Renamed variables which store large scale forcing tendencies
@@ -294,5 +296,5 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE ::                                     &
           c_u_m, c_u_m_l, c_v_m, c_v_m_l, c_w_m, c_w_m_l, ddzu, ddzu_pres,     &
-          dd2zu, dzu, ddzw, dzw, hyp, inflow_damping_factor, lad, l_grid,      &
+          dd2zu, dzu, ddzw, dzw, hyp, inflow_damping_factor, l_grid,           &
           nc_1d, nr_1d, ptdf_x, ptdf_y, p_surf, pt_surf, pt_1d, pt_init,       &
           qsws_surf, q_1d, q_init, q_surf, qc_1d, qr_1d, rdf, rdf_sc,          & 
@@ -313,9 +315,9 @@
 
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
-          canopy_heat_flux, cdc, d, de_dx, de_dy, de_dz, diss, diss_l_e,       &
+          d, de_dx, de_dy, de_dz, diss, diss_l_e,                              &
           diss_l_nr, diss_l_pt, diss_l_q, diss_l_qr, diss_l_sa, diss_l_u,      &
           diss_l_v, diss_l_w, flux_l_e, flux_l_nr, flux_l_pt, flux_l_q,        &
-          flux_l_qr, flux_l_sa, flux_l_u, flux_l_v, flux_l_w, kh, km, lad_s,   &
-          lad_u, lad_v, lad_w, lai, l_wall, p_loc, sec, sls, tend, tric,       &
+          flux_l_qr, flux_l_sa, flux_l_u, flux_l_v, flux_l_w, kh, km,          &
+          l_wall, p_loc, tend, tric,                                           &
           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
@@ -533,5 +535,4 @@
                              bc_sa_t = 'neumann', &
                              bc_uv_b = 'dirichlet', bc_uv_t = 'dirichlet', &
-                             canopy_mode = 'block', &
                              cloud_scheme = 'seifert_beheng', &
                              coupling_mode = 'uncoupled', &
@@ -586,5 +587,5 @@
                      nr_timesteps_this_run = 0, &
                      nsor = 20, nsor_ini = 100, n_sor, normalizing_region = 0, &
-                     nz_do3d = -9999, pch_index = 0, prt_time_count = 0, &
+                     nz_do3d = -9999, prt_time_count = 0, &
                      recycling_plane, runnr = 0, &
                      skip_do_avs = 0, subdomain_size, terminate_coupled = 0, &
@@ -596,5 +597,4 @@
                      do3d_no(0:1) = 0, do3d_time_count(0:1), &
                      domask_no(max_masks,0:1) = 0, domask_time_count(max_masks,0:1),&
-                     lad_vertical_gradient_level_ind(10) = -9999, &
                      mask_size(max_masks,3) = -1, mask_size_l(max_masks,3) = -1, &
                      mask_start_l(max_masks,3) = -1, &
@@ -652,5 +652,5 @@
                 outflow_l = .FALSE., outflow_n = .FALSE., outflow_r = .FALSE., &
                 outflow_s = .FALSE., passive_scalar = .FALSE., &
-                plant_canopy = .FALSE., prandtl_layer = .TRUE., &
+                prandtl_layer = .TRUE., &
                 precipitation = .FALSE., &
                 radiation = .FALSE., random_heatflux = .FALSE., &
@@ -682,10 +682,10 @@
                  canyon_width_x = 9999999.9_wp, canyon_width_y = 9999999.9_wp, &
                  canyon_wall_left = 9999999.9_wp, canyon_wall_south = 9999999.9_wp, &
-                 cthf = 0.0_wp, cfl_factor = -1.0_wp, cos_alpha_surface, &
+                 cfl_factor = -1.0_wp, cos_alpha_surface, &
                  coupling_start_time = 0.0_wp, disturbance_amplitude = 0.25_wp, &
                  disturbance_energy_limit = 0.01_wp, &
                  disturbance_level_b = -9999999.9_wp, &
                  disturbance_level_t = -9999999.9_wp, &
-                 dp_level_b = 0.0_wp, drag_coefficient = 0.0_wp, &
+                 dp_level_b = 0.0_wp, &
                  dt = -1.0_wp, dt_averaging_input = 0.0_wp, &
                  dt_averaging_input_pr = 9999999.9_wp, dt_coupling = 9999999.9_wp, &
@@ -703,5 +703,4 @@
                  f = 0.0_wp, fs = 0.0_wp, g = 9.81_wp, inflow_damping_height = 9999999.9_wp, &
                  inflow_damping_width = 9999999.9_wp, kappa = 0.4_wp, km_constant = -1.0_wp,&
-                 lad_surface = 0.0_wp, leaf_surface_concentration = 0.0_wp, &
                  mask_scale_x = 1.0_wp, mask_scale_y = 1.0_wp, mask_scale_z = 1.0_wp, &
                  maximum_cpu_time_allowed = 0.0_wp,  &
@@ -719,5 +718,5 @@
                  restart_time = 9999999.9_wp, rho_reference, rho_surface, &
                  rif_max = 1.0_wp, rif_min = -5.0_wp, roughness_length = 0.1_wp, &
-                 sa_surface = 35.0_wp, scalar_exchange_coefficient = 0.0_wp, &
+                 sa_surface = 35.0_wp, &
                  simulated_time = 0.0_wp, simulated_time_at_begin, sin_alpha_surface, &
                  skip_time_data_output = 0.0_wp, skip_time_data_output_av = 9999999.9_wp,&
@@ -745,6 +744,5 @@
     REAL(wp) ::  do2d_xy_last_time(0:1) = -1.0_wp, do2d_xz_last_time(0:1) = -1.0_wp, &
                  do2d_yz_last_time(0:1) = -1.0_wp, dpdxy(1:2) = 0.0_wp, &
-                 dt_domask(max_masks) = 9999999.9_wp, lad_vertical_gradient(10) = 0.0_wp,&
-                 lad_vertical_gradient_level(10) = -9999999.9_wp, &
+                 dt_domask(max_masks) = 9999999.9_wp, &
                  mask_scale(3), &
                  pt_vertical_gradient(10) = 0.0_wp, &
Index: palm/trunk/SOURCE/package_parin.f90
===================================================================
--- palm/trunk/SOURCE/package_parin.f90	(revision 1483)
+++ palm/trunk/SOURCE/package_parin.f90	(revision 1484)
@@ -20,5 +20,11 @@
 ! Current revisions:
 ! -----------------
-! 
+! Changes due to new module structure of the plant canopy model:
+!   module plant_canopy_model_mod added,
+!   new package/namelist canopy_par added, i.e. the canopy model is no longer 
+!   steered over the inipar-namelist,
+!   drag_coefficient, leaf_surface_concentration and scalar_exchange_coefficient
+!   renamed to canopy_drag_coeff, leaf_surface_conc and leaf_scalar_exch_coeff.
+! Changed statement tags in CONTINUE-statement
 ! 
 ! Former revisions:
@@ -103,4 +109,11 @@
                write_particle_statistics
 
+    USE plant_canopy_model_mod,                                                &
+        ONLY:  alpha_lad, beta_lad, calc_beta_lad_profile, canopy_drag_coeff,  &
+               canopy_mode, cthf, lad_surface,                                 &
+               lad_vertical_gradient, lad_vertical_gradient_level, lai_beta,   &
+               leaf_scalar_exch_coeff, leaf_surface_conc, pch_index,           &
+               plant_canopy
+
     USE spectrum,                                                              &
         ONLY:  comp_spectra_level, data_output_sp, plot_spectra_level,         &
@@ -110,4 +123,13 @@
 
     CHARACTER (LEN=80) ::  line  !:
+
+    NAMELIST /canopy_par/         alpha_lad, beta_lad, canopy_drag_coeff,      &
+                                  canopy_mode, cthf,                           &
+                                  lad_surface,                                 & 
+                                  lad_vertical_gradient,                       &
+                                  lad_vertical_gradient_level,                 &
+                                  lai_beta,                                    &
+                                  leaf_scalar_exch_coeff,                      &
+                                  leaf_surface_conc, pch_index
 
     NAMELIST /dvrp_graphics_par/  clip_dvrp_l, clip_dvrp_n, clip_dvrp_r,       &
@@ -159,9 +181,35 @@
     line = ' '
 
+!
+!-- Try to find canopy package
+    REWIND ( 11 )
+    line = ' '
+    DO   WHILE ( INDEX( line, '&canopy_par' ) == 0 )
+       READ ( 11, '(A)', END=10 )  line
+    ENDDO
+    BACKSPACE ( 11 )
+
+!
+!-- Read user-defined namelist
+    READ ( 11, canopy_par )
+
+!
+!-- Set flag that indicates that canopy model is switched on
+    plant_canopy = .TRUE.
+
+!
+!-- Set flag that indicates that the lad-profile shall be calculated by using
+!-- a beta probability density function
+    IF ( alpha_lad /= 9999999.9_wp  .AND.  beta_lad /= 9999999.9_wp )          &
+       calc_beta_lad_profile = .TRUE.
+
+ 10 CONTINUE
+
+
 #if defined( __dvrp_graphics )
     REWIND ( 11 )
     line = ' '
     DO   WHILE ( INDEX( line, '&dvrp_graphics_par' ) == 0 )
-       READ ( 11, '(A)', END=10 )  line
+       READ ( 11, '(A)', END=20 )  line
     ENDDO
     BACKSPACE ( 11 )
@@ -171,5 +219,5 @@
     READ ( 11, dvrp_graphics_par )
 
- 10 CONTINUE
+ 20 CONTINUE
 #endif
 
@@ -179,5 +227,5 @@
     line = ' '
     DO   WHILE ( INDEX( line, '&particles_par' ) == 0 )
-       READ ( 11, '(A)', END=20 )  line
+       READ ( 11, '(A)', END=30 )  line
     ENDDO
     BACKSPACE ( 11 )
@@ -191,5 +239,5 @@
     particle_advection = .TRUE.
 
- 20 CONTINUE
+ 30 CONTINUE
 
 
@@ -198,5 +246,5 @@
     line = ' '
     DO   WHILE ( INDEX( line, '&spectra_par' ) == 0 )
-       READ ( 11, '(A)', END=30 )  line
+       READ ( 11, '(A)', END=40 )  line
     ENDDO
     BACKSPACE ( 11 )
@@ -211,5 +259,5 @@
     IF ( dt_dosp == 9999999.9_wp )  dt_dosp = dt_data_output
 
- 30 CONTINUE
+ 40 CONTINUE
 #endif
 
Index: palm/trunk/SOURCE/parin.f90
===================================================================
--- palm/trunk/SOURCE/parin.f90	(revision 1483)
+++ palm/trunk/SOURCE/parin.f90	(revision 1484)
@@ -19,6 +19,8 @@
 !
 ! Current revisions:
-! ------------------
-! 
+! -----------------
+! Changes due to new module structure of the plant canopy model:
+!   canopy-model related parameters moved to new package canopy_par in 
+!   subroutine package_parin
 ! 
 ! Former revisions:
@@ -151,5 +153,5 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  lad, pt_init, q_init, ref_state, sa_init, ug, u_init, v_init,   &
+        ONLY:  pt_init, q_init, ref_state, sa_init, ug, u_init, v_init,        &
                vg
 
@@ -165,10 +167,10 @@
                building_length_y, building_wall_left, building_wall_south,     &
                call_microphysics_at_all_substeps, call_psolver_at_all_substeps,&
-               canopy_mode, canyon_height,                                     &
+               canyon_height,                                                  &
                canyon_width_x, canyon_width_y, canyon_wall_left,               &
                canyon_wall_south, cfl_factor,                                  &
                cloud_droplets, cloud_physics, cloud_scheme,                    &
                conserve_volume_flow, conserve_volume_flow_mode,                &
-               coupling_start_time, create_disturbances, cthf, cycle_mg,       &
+               coupling_start_time, create_disturbances, cycle_mg,             &
                data_output, data_output_masks,                                 &
                data_output_pr, data_output_2d_on_each_pe,                      &
@@ -176,5 +178,5 @@
                disturbance_level_b, disturbance_level_t, dissipation_1d,       &
                do2d_at_begin, do3d_at_begin,                                   &
-               dp_external, dp_level_b, dp_smooth, dpdxy, drag_coefficient,    &
+               dp_external, dp_level_b, dp_smooth, dpdxy,                      &
                drizzle, dt, dt_averaging_input, dt_averaging_input_pr,         &
                dt_coupling, dt_data_output, dt_data_output_av, dt_disturb,     &
@@ -187,7 +189,7 @@
                inflow_damping_width, inflow_disturbance_begin,                 &
                inflow_disturbance_end,  initializing_actions, io_blocks,       &
-               io_group, km_constant, lad_surface, lad_vertical_gradient,      &
-               lad_vertical_gradient_level, large_scale_forcing,               &
-               large_scale_subsidence, leaf_surface_concentration,             &
+               io_group, km_constant,                                          &
+               large_scale_forcing,                                            &
+               large_scale_subsidence,                                         &
                loop_optimization, masking_method, mask_scale_x, mask_scale_y,  &
                mask_scale_z, mask_x, mask_y, mask_z, mask_x_loop,              &
@@ -197,6 +199,6 @@
                momentum_advec, netcdf_data_format, netcdf_precision, neutral,  &
                ngsrb, normalizing_region, nsor, nsor_ini, nudging, ocean,      &
-               omega, omega_sor, passive_scalar, pch_index, phi, nz_do3d,      &
-               plant_canopy, prandtl_layer, prandtl_number, precipitation,     &
+               omega, omega_sor, passive_scalar, phi, nz_do3d,                 &
+               prandtl_layer, prandtl_number, precipitation,                   &
                precipitation_amount_interval, psolver, pt_damping_factor,      &
                pt_damping_width, pt_reference, pt_surface,                     &
@@ -211,5 +213,5 @@
                run_identifier, sa_surface, sa_vertical_gradient,               &
                sa_vertical_gradient_level, scalar_advec,                       &
-               scalar_exchange_coefficient, scalar_rayleigh_damping,           &
+               scalar_rayleigh_damping,                                        &
                section_xy, section_xz, section_yz, skip_time_data_output,      &
                skip_time_data_output_av, skip_time_dopr, skip_time_do2d_xy,    &
@@ -270,13 +272,13 @@
              building_length_y, building_wall_left, building_wall_south,       &
              call_psolver_at_all_substeps, call_microphysics_at_all_substeps,  &
-             canopy_mode, canyon_height,                                       &
+             canyon_height,                                                    &
              canyon_width_x, canyon_width_y, canyon_wall_left,                 &
              canyon_wall_south, c_sedimentation, cfl_factor, cloud_droplets,   &
              cloud_physics, cloud_scheme, collective_wait,                     &
              conserve_volume_flow,                                             &
-             conserve_volume_flow_mode, coupling_start_time, cthf,             &
+             conserve_volume_flow_mode, coupling_start_time,                   &
              curvature_solution_effects, cycle_mg, damp_level_1d,              &
              dissipation_1d,                                                   &
-             dp_external, dp_level_b, dp_smooth, dpdxy, drag_coefficient,      &
+             dp_external, dp_level_b, dp_smooth, dpdxy,                        &
              drizzle, dt, dt_pr_1d, dt_run_control_1d, dx, dy, dz, dz_max,     & 
              dz_stretch_factor, dz_stretch_level, end_time_1d,                 &
@@ -285,13 +287,12 @@
              inflow_damping_height, inflow_damping_width,                      &
              inflow_disturbance_begin, inflow_disturbance_end,                 &
-             initializing_actions, km_constant, lad_surface,                   &
-             lad_vertical_gradient, lad_vertical_gradient_level,               &
+             initializing_actions, km_constant,                                &
              large_scale_forcing, large_scale_subsidence,                      &
-             leaf_surface_concentration, limiter_sedimentation,                &
+             limiter_sedimentation, &
              loop_optimization, masking_method, mg_cycles,                     &
              mg_switch_to_pe0_level, mixing_length_1d, momentum_advec,         &
              nc_const, netcdf_precision, neutral, ngsrb,                       &
              nsor, nsor_ini, nudging, nx, ny, nz, ocean, omega, omega_sor,     &
-             passive_scalar, pch_index, phi, plant_canopy, prandtl_layer,      &
+             passive_scalar, phi, prandtl_layer, &
              prandtl_number, precipitation, psolver, pt_damping_factor,        &
              pt_damping_width, pt_reference, pt_surface,                       &
@@ -304,5 +305,5 @@
              rif_max, rif_min, roughness_length, sa_surface,                   &
              sa_vertical_gradient, sa_vertical_gradient_level, scalar_advec,   &
-             scalar_exchange_coefficient, scalar_rayleigh_damping,             &
+             scalar_rayleigh_damping,                                          &
              statistic_regions, subs_vertical_gradient,                        &
              subs_vertical_gradient_level, surface_heatflux, surface_pressure, &
@@ -478,5 +479,5 @@
 !--          ATTENTION: in case of changes to the following statement please
 !--                  also check the allocate statement in routine read_var_list
-             ALLOCATE( lad(0:nz+1), pt_init(0:nz+1), q_init(0:nz+1),           &
+             ALLOCATE( pt_init(0:nz+1), q_init(0:nz+1),                       &
                        ref_state(0:nz+1), sa_init(0:nz+1), ug(0:nz+1),         &
                        u_init(0:nz+1), v_init(0:nz+1), vg(0:nz+1),             &
Index: palm/trunk/SOURCE/plant_canopy_model.f90
===================================================================
--- palm/trunk/SOURCE/plant_canopy_model.f90	(revision 1483)
+++ palm/trunk/SOURCE/plant_canopy_model.f90	(revision 1484)
@@ -20,5 +20,22 @@
 ! Current revisions:
 ! -----------------
-! 
+! Changes due to new module structure of the plant canopy model:
+!   module plant_canopy_model_mod now contains a subroutine for the 
+!   initialization of the canopy model (init_plant_canopy),
+!   limitation of the canopy drag (previously accounted for by calculation of 
+!   a limiting canopy timestep for the determination of the maximum LES timestep
+!   in subroutine timestep) is now realized by the calculation of pre-tendencies
+!   and preliminary velocities in subroutine plant_canopy_model,
+!   some redundant MPI communication removed in subroutine init_plant_canopy
+!   (was previously in init_3d_model),
+!   unnecessary 3d-arrays lad_u, lad_v, lad_w removed - lad information on the
+!   respective grid is now provided only by lad_s (e.g. in the calculation of 
+!   the tendency terms or of cum_lai_hf),
+!   drag_coefficient, lai, leaf_surface_concentration, 
+!   scalar_exchange_coefficient, sec and sls renamed to canopy_drag_coeff, 
+!   cum_lai_hf, leaf_surface_conc, leaf_scalar_exch_coeff, lsec and lsc,
+!   respectively,
+!   unnecessary 3d-arrays cdc, lsc and lsec now defined as single-value constants,
+!   USE-statements and ONLY-lists modified accordingly
 ! 
 ! Former revisions:
@@ -46,10 +63,75 @@
 ! Description:
 ! ------------
-! Evaluation of sinks and sources of momentum, heat and scalar concentration 
-! due to canopy elements
+! 1) Initialization of the canopy model, e.g. construction of leaf area density 
+! profile (subroutine init_plant_canopy).
+! 2) Calculation of sinks and sources of momentum, heat and scalar concentration 
+! due to canopy elements (subroutine plant_canopy_model).
 !------------------------------------------------------------------------------!
+    USE arrays_3d,                                                             &
+        ONLY:  dzu, dzw, e, q, shf, tend, u, v, w, zu, zw 
+
+    USE indices,                                                               &
+        ONLY:  nbgp, nxl, nxlg, nxlu, nxr, nxrg, nyn, nyng, nys, nysg, nysv,   &
+               nz, nzb, nzb_s_inner, nzb_u_inner, nzb_v_inner, nzb_w_inner, nzt
+
+    USE kinds
+
+
+    IMPLICIT NONE
+
+
+    CHARACTER (LEN=20)   ::  canopy_mode = 'block' !: canopy coverage
+
+    INTEGER(iwp) ::  pch_index = 0                 !: plant canopy height/top index
+    INTEGER(iwp) ::                                                            &
+       lad_vertical_gradient_level_ind(10) = -9999 !: lad-profile levels (index)
+
+    LOGICAL ::  calc_beta_lad_profile = .FALSE. !: switch for calc. of lad from beta func.
+    LOGICAL ::  plant_canopy = .FALSE.          !: switch for use of canopy model
+
+    REAL(wp) ::  alpha_lad = 9999999.9_wp   !: coefficient for lad calculation
+    REAL(wp) ::  beta_lad = 9999999.9_wp    !: coefficient for lad calculation
+    REAL(wp) ::  canopy_drag_coeff = 0.0_wp !: canopy drag coefficient (parameter)
+    REAL(wp) ::  cdc = 0.0_wp               !: canopy drag coeff. (abbreviation used in equations)
+    REAL(wp) ::  cthf = 0.0_wp              !: canopy top heat flux
+    REAL(wp) ::  dt_plant_canopy = 0.0_wp   !: timestep account. for canopy drag
+    REAL(wp) ::  lad_surface = 0.0_wp       !: lad surface value
+    REAL(wp) ::  lai_beta = 0.0_wp          !: leaf area index (lai) for lad calc.
+    REAL(wp) ::                                                                &
+       leaf_scalar_exch_coeff = 0.0_wp      !: canopy scalar exchange coeff.
+    REAL(wp) ::                                                                &
+       leaf_surface_conc = 0.0_wp           !: leaf surface concentration
+    REAL(wp) ::  lsec = 0.0_wp              !: leaf scalar exchange coeff.
+    REAL(wp) ::  lsc = 0.0_wp               !: leaf surface concentration
+
+    REAL(wp) ::                                                                &
+       lad_vertical_gradient(10) = 0.0_wp              !: lad gradient
+    REAL(wp) ::                                                                &
+       lad_vertical_gradient_level(10) = -9999999.9_wp !: lad-prof. levels (in m)
+
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  lad            !: leaf area density
+    REAL(wp), DIMENSION(:), ALLOCATABLE ::  pre_lad        !: preliminary lad
+    
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
+       canopy_heat_flux                                    !: canopy heat flux
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  cum_lai_hf !: cumulative lai for heatflux calc.
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  lad_s      !: lad on scalar-grid
+
+
+    SAVE
+
 
     PRIVATE
-    PUBLIC plant_canopy_model
+    PUBLIC alpha_lad, beta_lad, calc_beta_lad_profile, canopy_drag_coeff,      &
+           canopy_mode, cdc, cthf, dt_plant_canopy, init_plant_canopy, lad,    &
+           lad_s, lad_surface, lad_vertical_gradient,                          &
+           lad_vertical_gradient_level, lad_vertical_gradient_level_ind,       &
+           lai_beta, leaf_scalar_exch_coeff, leaf_surface_conc, lsc, lsec,     &
+           pch_index, plant_canopy, plant_canopy_model
+
+
+    INTERFACE init_plant_canopy
+       MODULE PROCEDURE init_plant_canopy
+    END INTERFACE init_plant_canopy
 
     INTERFACE plant_canopy_model
@@ -58,22 +140,255 @@
     END INTERFACE plant_canopy_model
 
+
+
+
  CONTAINS
 
 
 !------------------------------------------------------------------------------!
-! Call for all grid points
+! Description:
+! ------------
+!-- Initialization of the plant canopy model
+!------------------------------------------------------------------------------!
+    SUBROUTINE init_plant_canopy
+    
+
+       USE control_parameters,                                                 &
+           ONLY: dz, ocean, passive_scalar
+
+
+       IMPLICIT NONE
+
+       INTEGER(iwp) ::  i !: running index
+       INTEGER(iwp) ::  j !: running index
+       INTEGER(iwp) ::  k !: running index
+
+       REAL(wp) ::  int_bpdf      !: vertical integral for lad-profile construction
+       REAL(wp) ::  dzh           !: vertical grid spacing in units of canopy height
+       REAL(wp) ::  gradient      !: gradient for lad-profile construction
+       REAL(wp) ::  canopy_height !: canopy height for lad-profile construction
+
+!
+!--    Allocate one-dimensional arrays for the computation of the 
+!--    leaf area density (lad) profile
+       ALLOCATE( lad(0:nz+1), pre_lad(0:nz+1) )
+       lad = 0.0_wp
+       pre_lad = 0.0_wp
+
+!
+!--    Compute the profile of leaf area density used in the plant
+!--    canopy model. The profile can either be constructed from
+!--    prescribed vertical gradients of the leaf area density or by
+!--    using a beta probability density function (see e.g. Markkanen et al., 
+!--    2003: Boundary-Layer Meteorology, 106, 437-459) 
+       IF (  .NOT.  calc_beta_lad_profile )  THEN   
+
+!
+!--       Use vertical gradients for lad-profile construction    
+          i = 1
+          gradient = 0.0_wp
+
+          IF (  .NOT.  ocean )  THEN
+
+             lad(0) = lad_surface
+             lad_vertical_gradient_level_ind(1) = 0
+ 
+             DO k = 1, pch_index
+                IF ( i < 11 )  THEN
+                   IF ( lad_vertical_gradient_level(i) < zu(k)  .AND.          &
+                        lad_vertical_gradient_level(i) >= 0.0_wp )  THEN
+                      gradient = lad_vertical_gradient(i)
+                      lad_vertical_gradient_level_ind(i) = k - 1
+                      i = i + 1
+                   ENDIF
+                ENDIF
+                IF ( gradient /= 0.0_wp )  THEN
+                   IF ( k /= 1 )  THEN
+                      lad(k) = lad(k-1) + dzu(k) * gradient
+                   ELSE
+                      lad(k) = lad_surface + 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. This information is used for the HEADER and the RUN_CONTROL
+!--       file.
+          IF ( lad_vertical_gradient_level(1) == -9999999.9_wp )  THEN
+             lad_vertical_gradient_level(1) = 0.0_wp
+          ENDIF
+
+       ELSE
+
+! 
+!--       Use beta function for lad-profile construction
+          int_bpdf = 0.0_wp
+          canopy_height = pch_index * dz
+
+          DO k = nzb, pch_index
+             int_bpdf = int_bpdf +                                             &
+                           ( ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) ) * &
+                           ( ( 1.0_wp - ( zw(k) / canopy_height ) )**( beta_lad-1.0_wp ) ) *   &
+                           ( ( zw(k+1)-zw(k) ) / canopy_height ) )
+          ENDDO
+
+!
+!--       Preliminary lad profile (defined on w-grid)
+          DO k = nzb, pch_index
+             pre_lad(k) =  lai_beta *                                              &
+                           (   ( ( zw(k) / canopy_height )**( alpha_lad-1.0_wp ) ) *          &
+                               ( ( 1.0_wp - ( zw(k) / canopy_height ) )**( beta_lad-1.0_wp ) ) /   &
+                               int_bpdf                                            &
+                           ) / canopy_height
+          ENDDO
+
+!
+!--       Final lad profile (defined on scalar-grid level, since most prognostic
+!--       quantities are defined there, hence, less interpolation is required 
+!--       when calculating the canopy tendencies)
+          lad(0) = pre_lad(0)
+          DO k = nzb+1, pch_index
+             lad(k) = 0.5 * ( pre_lad(k-1) + pre_lad(k) )
+          ENDDO          
+
+       ENDIF
+
+!
+!--    Allocate 3D-array for the leaf area density (lad_s). In case of a 
+!--    prescribed canopy-top heat flux (cthf), allocate 3D-arrays for
+!--    the cumulative leaf area index (cum_lai_hf) and the canopy heat flux.
+       ALLOCATE( lad_s(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+
+       IF ( cthf /= 0.0_wp )  THEN
+          ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                 &
+                    canopy_heat_flux(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+       ENDIF
+
+!
+!--    Initialize canopy parameters cdc (canopy drag coefficient), 
+!--    lsec (leaf scalar exchange coefficient), lsc (leaf surface concentration)
+!--    with the prescribed values
+       cdc = canopy_drag_coeff
+       lsec = leaf_scalar_exch_coeff
+       lsc = leaf_surface_conc
+
+!
+!--    Initialization of the canopy coverage in the model domain:
+!--    Setting the parameter canopy_mode = 'block' initializes a canopy, which
+!--    fully covers the domain surface
+       SELECT CASE ( TRIM( canopy_mode ) )
+
+          CASE( 'block' )
+
+             DO  i = nxlg, nxrg
+                DO  j = nysg, nyng
+                   lad_s(:,j,i) = lad(:)
+                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
+
+!
+!--    Initialization of the canopy heat source distribution
+       IF ( cthf /= 0.0_wp )  THEN
+!
+!--       Piecewise calculation of the leaf area index by vertical
+!--       integration of the leaf area density
+          cum_lai_hf(:,:,:) = 0.0_wp
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
+                DO  k = pch_index-1, 0, -1
+                   IF ( k == pch_index-1 )  THEN
+                      cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) +                &
+                         ( 0.5_wp * lad_s(k+1,j,i) *                           &
+                           ( zw(k+1) - zu(k+1) ) )  +                          &
+                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) +              &
+                                                 lad_s(k,j,i) ) +              &
+                                      lad_s(k+1,j,i) ) *                       &
+                           ( zu(k+1) - zw(k) ) ) 
+                   ELSE
+                      cum_lai_hf(k,j,i) = cum_lai_hf(k+1,j,i) +                &
+                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+2,j,i) +              &
+                                                 lad_s(k+1,j,i) ) +            &
+                                      lad_s(k+1,j,i) ) *                       &
+                           ( zw(k+1) - zu(k+1) ) )  +                          &
+                         ( 0.5_wp * ( 0.5_wp * ( lad_s(k+1,j,i) +              &
+                                                 lad_s(k,j,i) ) +              &
+                                      lad_s(k+1,j,i) ) *                       &
+                           ( zu(k+1) - zw(k) ) )
+                   ENDIF
+                ENDDO
+             ENDDO
+          ENDDO
+
+!
+!--       Calculation of the upward kinematic vertical heat flux within the 
+!--       canopy
+          DO  i = nxlg, nxrg
+             DO  j = nysg, nyng
+                DO  k = 0, pch_index
+                   canopy_heat_flux(k,j,i) = cthf *                            &
+                                             exp( -0.6_wp * cum_lai_hf(k,j,i) )
+                ENDDO
+             ENDDO
+          ENDDO
+
+!
+!--       The surface heat flux is set to the surface value of the calculated
+!--       in-canopy heat flux distribution  
+          shf(:,:) = canopy_heat_flux(0,:,:)
+
+       ENDIF
+
+
+
+    END SUBROUTINE init_plant_canopy
+
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!-- Calculation of the tendency terms, accounting for the effect of the plant
+!-- canopy on momentum and scalar quantities.
+!--
+!-- The canopy is located where the leaf area density lad_s(k,j,i) > 0.0 
+!-- (defined on scalar grid), as initialized in subroutine init_plant_canopy. 
+!-- The lad on the w-grid is vertically interpolated from the surrounding
+!-- lad_s. The upper boundary of the canopy is defined on the w-grid at 
+!-- k = pch_index. Here, the lad is zero.
+!--
+!-- The canopy drag must be limited (previously accounted for by calculation of 
+!-- a limiting canopy timestep for the determination of the maximum LES timestep
+!-- in subroutine timestep), since it is physically impossible that the canopy
+!-- drag alone can locally change the sign of a velocity component. This 
+!-- limitation is realized by calculating preliminary tendencies and velocities.
+!-- It is subsequently checked if the preliminary new velocity has a different
+!-- sign than the current velocity. If so, the tendency is limited in a way that
+!-- the velocity can at maximum be reduced to zero by the canopy drag. 
+!--
+!--
+!-- Call for all grid points
 !------------------------------------------------------------------------------!
     SUBROUTINE plant_canopy_model( component )
 
-       USE arrays_3d,                                                          &
-           ONLY:  canopy_heat_flux, cdc, dzw, e, lad_s, lad_u, lad_v, lad_w,   &
-                  q, sec, sls, tend, u, v, w
 
        USE control_parameters,                                                 &
-           ONLY:  pch_index, message_string
-
-       USE indices,                                                            &
-           ONLY:  nxl, nxlu, nxr, nys, nysv, nyn, nzb_s_inner, nzb_u_inner,    &
-                  nzb_v_inner, nzb_w_inner
+           ONLY:  dt_3d, message_string
 
        USE kinds
@@ -81,11 +396,21 @@
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  component  !:
-       INTEGER(iwp) ::  i          !:
-       INTEGER(iwp) ::  j          !:
-       INTEGER(iwp) ::  k          !:
+       INTEGER(iwp) ::  component !: prognostic variable (u,v,w,pt,q,e)
+       INTEGER(iwp) ::  i         !: running index
+       INTEGER(iwp) ::  j         !: running index
+       INTEGER(iwp) ::  k         !: running index
+
+       REAL(wp) ::  ddt_3d    !: inverse of the LES timestep (dt_3d)
+       REAL(wp) ::  lad_local !: local lad value
+       REAL(wp) ::  pre_tend  !: preliminary tendency
+       REAL(wp) ::  pre_u     !: preliminary u-value 
+       REAL(wp) ::  pre_v     !: preliminary v-value 
+       REAL(wp) ::  pre_w     !: preliminary w-value 
+
+
+       ddt_3d = 1.0_wp / dt_3d
  
 !
-!--    Compute drag for the three velocity components and the SGS-TKE
+!--    Compute drag for the three velocity components and the SGS-TKE:
        SELECT CASE ( component )
 
@@ -96,18 +421,54 @@
                 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_wp )**2     +  &
-                                          ( ( w(k-1,j,i-1)    +  &
-                                              w(k-1,j,i)      +  &
-                                              w(k,j,i-1)      +  &
-                                              w(k,j,i) )         &
-                                            / 4.0_wp )**2 )   *  &
-                                    u(k,j,i)
+
+!
+!--                   In order to create sharp boundaries of the plant canopy, 
+!--                   the lad on the u-grid at index (k,j,i) is equal to 
+!--                   lad_s(k,j,i), rather than being interpolated from the 
+!--                   surrounding lad_s, because this would yield smaller lad 
+!--                   at the canopy boundaries than inside of the canopy.
+!--                   For the same reason, the lad at the rightmost(i+1)canopy 
+!--                   boundary on the u-grid equals lad_s(k,j,i).
+                      lad_local = lad_s(k,j,i)
+                      IF ( lad_local == 0.0_wp  .AND.                          &
+                           lad_s(k,j,i-1) > 0.0_wp )  THEN
+                         lad_local = lad_s(k,j,i-1)
+                      ENDIF
+
+                      pre_tend = 0.0_wp
+                      pre_u = 0.0_wp
+!
+!--                   Calculate preliminary value (pre_tend) of the tendency
+                      pre_tend = - cdc *                                       &
+                                   lad_local *                                 &
+                                   SQRT( u(k,j,i)**2 +                         &
+                                         ( 0.25_wp * ( v(k,j,i-1) +            &
+                                                       v(k,j,i)   +            &
+                                                       v(k,j+1,i) +            &
+                                                       v(k,j+1,i-1) )          &
+                                         )**2 +                                &
+                                         ( 0.25_wp * ( w(k-1,j,i-1) +          &
+                                                       w(k-1,j,i)   +          &
+                                                       w(k,j,i-1)   +          &
+                                                       w(k,j,i) )              &
+                                         )**2                                  &
+                                       ) *                                     &
+                                   u(k,j,i)
+
+!
+!--                   Calculate preliminary new velocity, based on pre_tend
+                      pre_u = u(k,j,i) + dt_3d * pre_tend
+!
+!--                   Compare sign of old velocity and new preliminary velocity,
+!--                   and in case the signs are different, limit the tendency
+                      IF ( SIGN(pre_u,u(k,j,i)) /= pre_u )  THEN
+                         pre_tend = - u(k,j,i) * ddt_3d
+                      ELSE
+                         pre_tend = pre_tend
+                      ENDIF
+!
+!--                   Calculate final tendency
+                      tend(k,j,i) = tend(k,j,i) + pre_tend
+
                    ENDDO
                 ENDDO
@@ -120,18 +481,54 @@
                 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_wp )**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_wp )**2 )   *  &
-                                    v(k,j,i)  
+
+!
+!--                   In order to create sharp boundaries of the plant canopy, 
+!--                   the lad on the v-grid at index (k,j,i) is equal to 
+!--                   lad_s(k,j,i), rather than being interpolated from the 
+!--                   surrounding lad_s, because this would yield smaller lad 
+!--                   at the canopy boundaries than inside of the canopy.
+!--                   For the same reason, the lad at the northmost(j+1) canopy 
+!--                   boundary on the v-grid equals lad_s(k,j,i).
+                      lad_local = lad_s(k,j,i)
+                      IF ( lad_local == 0.0_wp  .AND.                          &
+                           lad_s(k,j-1,i) > 0.0_wp )  THEN
+                         lad_local = lad_s(k,j-1,i)
+                      ENDIF
+
+                      pre_tend = 0.0_wp
+                      pre_v = 0.0_wp
+!
+!--                   Calculate preliminary value (pre_tend) of the tendency
+                      pre_tend = - cdc *                                       &
+                                   lad_local *                                 &
+                                   SQRT( ( 0.25_wp * ( u(k,j-1,i)   +          &
+                                                       u(k,j-1,i+1) +          &
+                                                       u(k,j,i)     +          &
+                                                       u(k,j,i+1) )            &
+                                         )**2 +                                &
+                                         v(k,j,i)**2 +                         &
+                                         ( 0.25_wp * ( w(k-1,j-1,i) +          &
+                                                       w(k-1,j,i)   +          &
+                                                       w(k,j-1,i)   +          &
+                                                       w(k,j,i) )              &
+                                         )**2                                  &
+                                       ) *                                     &
+                                   v(k,j,i)
+
+!
+!--                   Calculate preliminary new velocity, based on pre_tend
+                      pre_v = v(k,j,i) + dt_3d * pre_tend
+!
+!--                   Compare sign of old velocity and new preliminary velocity,
+!--                   and in case the signs are different, limit the tendency
+                      IF ( SIGN(pre_v,v(k,j,i)) /= pre_v )  THEN
+                         pre_tend = - v(k,j,i) * ddt_3d
+                      ELSE
+                         pre_tend = pre_tend
+                      ENDIF
+!
+!--                   Calculate final tendency
+                      tend(k,j,i) = tend(k,j,i) + pre_tend
+
                    ENDDO
                 ENDDO
@@ -143,19 +540,41 @@
              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_wp )**2     +  &
-                                          ( ( v(k,j,i)        +  &
-                                              v(k,j+1,i)      +  &
-                                              v(k+1,j,i)      +  &
-                                              v(k+1,j+1,i) )     &
-                                            / 4.0_wp )**2     +  &
-                                              w(k,j,i)**2 )   *  &
-                                    w(k,j,i)
+                   DO  k = nzb_w_inner(j,i)+1, pch_index-1
+
+                      pre_tend = 0.0_wp
+                      pre_w = 0.0_wp
+!
+!--                   Calculate preliminary value (pre_tend) of the tendency
+                      pre_tend = - cdc *                                       &
+                                   (0.5_wp *                                   &
+                                      ( lad_s(k+1,j,i) + lad_s(k,j,i) )) *     &
+                                   SQRT( ( 0.25_wp * ( u(k,j,i)   +            &
+                                                       u(k,j,i+1) +            &
+                                                       u(k+1,j,i) +            &
+                                                       u(k+1,j,i+1) )          &
+                                         )**2 +                                &
+                                         ( 0.25_wp * ( v(k,j,i)   +            &
+                                                       v(k,j+1,i) +            &
+                                                       v(k+1,j,i) +            &
+                                                       v(k+1,j+1,i) )          &
+                                         )**2 +                                &
+                                         w(k,j,i)**2                           &
+                                       ) *                                     &
+                                   w(k,j,i)
+!
+!--                   Calculate preliminary new velocity, based on pre_tend
+                      pre_w = w(k,j,i) + dt_3d * pre_tend
+!
+!--                   Compare sign of old velocity and new preliminary velocity,
+!--                   and in case the signs are different, limit the tendency
+                      IF ( SIGN(pre_w,w(k,j,i)) /= pre_w )  THEN
+                         pre_tend = - w(k,j,i) * ddt_3d
+                      ELSE
+                         pre_tend = pre_tend
+                      ENDIF
+!
+!--                   Calculate final tendency
+                      tend(k,j,i) = tend(k,j,i) + pre_tend
+
                    ENDDO
                 ENDDO
@@ -168,8 +587,7 @@
                 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)
+                      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
@@ -182,16 +600,18 @@
                 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_wp )**2     +       &
-                                          ( ( v(k,j,i)        +       &
-                                              v(k,j+1,i) )            &
-                                            / 2.0_wp )**2     +       &
-                                          ( ( w(k-1,j,i)      +       & 
-                                              w(k,j,i) )              &
-                                            / 2.0_wp )**2 )   *       &
-                                    ( q(k,j,i) - sls(k,j,i) )
+                      tend(k,j,i) = tend(k,j,i) -                              &
+                                       lsec *                                  &
+                                       lad_s(k,j,i) *                          &
+                                       SQRT( ( 0.5_wp * ( u(k,j,i) +           &
+                                                          u(k,j,i+1) )         &
+                                             )**2 +                            &
+                                             ( 0.5_wp * ( v(k,j,i) +           &
+                                                          v(k,j+1,i) )         &
+                                             )**2 +                            &
+                                             ( 0.5_wp * ( w(k-1,j,i) +         & 
+                                                          w(k,j,i) )           &
+                                             )**2                              &
+                                           ) *                                 &
+                                       ( q(k,j,i) - lsc )
                    ENDDO
                 ENDDO
@@ -204,20 +624,23 @@
                 DO  j = nys, nyn
                    DO  k = nzb_s_inner(j,i)+1, pch_index
-                      tend(k,j,i) = tend(k,j,i) -                     &
-                                    2.0_wp * cdc(k,j,i) * lad_s(k,j,i) * &
-                                    SQRT( ( ( u(k,j,i)              + &
-                                              u(k,j,i+1) )            &
-                                            / 2.0_wp )**2           + &
-                                          ( ( v(k,j,i)              + &
-                                              v(k,j+1,i) )            &
-                                            / 2.0_wp )**2           + &
-                                          ( ( w(k,j,i)              + &
-                                              w(k+1,j,i) )            &
-                                            / 2.0_wp )**2 )         * &
-                                    e(k,j,i)
+                      tend(k,j,i) = tend(k,j,i) -                              &
+                                       2.0_wp * cdc *                          &
+                                       lad_s(k,j,i) *                          &
+                                       SQRT( ( 0.5_wp * ( u(k,j,i) +           &
+                                                          u(k,j,i+1) )         &
+                                             )**2 +                            &
+                                             ( 0.5_wp * ( v(k,j,i) +           &
+                                                          v(k,j+1,i) )         &
+                                             )**2 +                            &
+                                             ( 0.5_wp * ( w(k,j,i) +           &
+                                                          w(k+1,j,i) )         &
+                                             )**2                              &
+                                           ) *                                 &
+                                       e(k,j,i)
                    ENDDO
                 ENDDO
              ENDDO 
-                        
+
+
           CASE DEFAULT
 
@@ -231,18 +654,32 @@
 
 !------------------------------------------------------------------------------!
-! Call for grid point i,j
+! Description:
+! ------------
+!-- Calculation of the tendency terms, accounting for the effect of the plant
+!-- canopy on momentum and scalar quantities.
+!--
+!-- The canopy is located where the leaf area density lad_s(k,j,i) > 0.0 
+!-- (defined on scalar grid), as initialized in subroutine init_plant_canopy. 
+!-- The lad on the w-grid is vertically interpolated from the surrounding
+!-- lad_s. The upper boundary of the canopy is defined on the w-grid at 
+!-- k = pch_index. Here, the lad is zero.
+!--
+!-- The canopy drag must be limited (previously accounted for by calculation of 
+!-- a limiting canopy timestep for the determination of the maximum LES timestep
+!-- in subroutine timestep), since it is physically impossible that the canopy
+!-- drag alone can locally change the sign of a velocity component. This 
+!-- limitation is realized by calculating preliminary tendencies and velocities.
+!-- It is subsequently checked if the preliminary new velocity has a different
+!-- sign than the current velocity. If so, the tendency is limited in a way that
+!-- the velocity can at maximum be reduced to zero by the canopy drag. 
+!--
+!--
+!-- Call for grid point i,j
 !------------------------------------------------------------------------------!
     SUBROUTINE plant_canopy_model_ij( i, j, component )
 
-       USE arrays_3d,                                                          &
-           ONLY:  canopy_heat_flux, cdc, dzw, e, lad_s, lad_u, lad_v, lad_w,   &
-                  q, sec, sls, tend, u, v, w
 
        USE control_parameters,                                                 &
-           ONLY:  pch_index, message_string
-
-       USE indices,                                                            &
-           ONLY:  nxl, nxlu, nxr, nys, nysv, nyn, nzb_s_inner, nzb_u_inner,    &
-                  nzb_v_inner, nzb_w_inner
+           ONLY:  dt_3d, message_string
 
        USE kinds
@@ -250,73 +687,175 @@
        IMPLICIT NONE
 
-       INTEGER(iwp) ::  component  !:
-       INTEGER(iwp) ::  i          !:
-       INTEGER(iwp) ::  j          !:
-       INTEGER(iwp) ::  k          !:
-
-!
-!--    Compute drag for the three velocity components
+       INTEGER(iwp) ::  component !: prognostic variable (u,v,w,pt,q,e)
+       INTEGER(iwp) ::  i         !: running index
+       INTEGER(iwp) ::  j         !: running index
+       INTEGER(iwp) ::  k         !: running index
+
+       REAL(wp) ::  ddt_3d    !: inverse of the LES timestep (dt_3d)
+       REAL(wp) ::  lad_local !: local lad value
+       REAL(wp) ::  pre_tend  !: preliminary tendency
+       REAL(wp) ::  pre_u     !: preliminary u-value 
+       REAL(wp) ::  pre_v     !: preliminary v-value 
+       REAL(wp) ::  pre_w     !: preliminary w-value 
+
+
+       ddt_3d = 1.0_wp / dt_3d
+
+!
+!--    Compute drag for the three velocity components and the SGS-TKE
        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_wp )**2 +        &
-                                    ( ( w(k-1,j,i-1) +       &
-                                        w(k-1,j,i)   +       &
-                                        w(k,j,i-1)   +       &
-                                        w(k,j,i) )           &
-                                      / 4.0_wp )**2 ) *      &
-                              u(k,j,i)
-          ENDDO
+          CASE ( 1 )
+             DO  k = nzb_u_inner(j,i)+1, pch_index
+
+!
+!--             In order to create sharp boundaries of the plant canopy, 
+!--             the lad on the u-grid at index (k,j,i) is equal to lad_s(k,j,i),
+!--             rather than being interpolated from the surrounding lad_s, 
+!--             because this would yield smaller lad at the canopy boundaries 
+!--             than inside of the canopy.
+!--             For the same reason, the lad at the rightmost(i+1)canopy 
+!--             boundary on the u-grid equals lad_s(k,j,i).
+                lad_local = lad_s(k,j,i)
+                IF ( lad_local == 0.0_wp  .AND.                                &
+                     lad_s(k,j,i-1) > 0.0_wp )  THEN
+                   lad_local = lad_s(k,j,i-1)
+                ENDIF
+
+                pre_tend = 0.0_wp
+                pre_u = 0.0_wp
+!
+!--             Calculate preliminary value (pre_tend) of the tendency
+                pre_tend = - cdc *                                             &
+                             lad_local *                                       &   
+                             SQRT( u(k,j,i)**2 +                               &
+                                   ( 0.25_wp * ( v(k,j,i-1)  +                 &
+                                                 v(k,j,i)    +                 &
+                                                 v(k,j+1,i)  +                 &
+                                                 v(k,j+1,i-1) )                &
+                                   )**2 +                                      &
+                                   ( 0.25_wp * ( w(k-1,j,i-1) +                &
+                                                 w(k-1,j,i)   +                &
+                                                 w(k,j,i-1)   +                &
+                                                 w(k,j,i) )                    &
+                                   )**2                                        &
+                                 ) *                                           &
+                             u(k,j,i)
+
+!
+!--             Calculate preliminary new velocity, based on pre_tend
+                pre_u = u(k,j,i) + dt_3d * pre_tend
+!
+!--             Compare sign of old velocity and new preliminary velocity,
+!--             and in case the signs are different, limit the tendency
+                IF ( SIGN(pre_u,u(k,j,i)) /= pre_u )  THEN
+                   pre_tend = - u(k,j,i) * ddt_3d
+                ELSE
+                   pre_tend = pre_tend
+                ENDIF
+!
+!--             Calculate final tendency
+                tend(k,j,i) = tend(k,j,i) + pre_tend
+             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_wp )**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_wp )**2 ) *      &
-                              v(k,j,i)
-          ENDDO
+          CASE ( 2 )
+             DO  k = nzb_v_inner(j,i)+1, pch_index
+
+!
+!--             In order to create sharp boundaries of the plant canopy, 
+!--             the lad on the v-grid at index (k,j,i) is equal to lad_s(k,j,i),
+!--             rather than being interpolated from the surrounding lad_s, 
+!--             because this would yield smaller lad at the canopy boundaries 
+!--             than inside of the canopy.
+!--             For the same reason, the lad at the northmost(j+1)canopy 
+!--             boundary on the v-grid equals lad_s(k,j,i).
+                lad_local = lad_s(k,j,i)
+                IF ( lad_local == 0.0_wp  .AND.                                &
+                     lad_s(k,j-1,i) > 0.0_wp )  THEN
+                   lad_local = lad_s(k,j-1,i)
+                ENDIF
+
+                pre_tend = 0.0_wp
+                pre_v = 0.0_wp
+!
+!--             Calculate preliminary value (pre_tend) of the tendency
+                pre_tend = - cdc *                                             &
+                             lad_local *                                       &
+                             SQRT( ( 0.25_wp * ( u(k,j-1,i)   +                &
+                                                 u(k,j-1,i+1) +                &
+                                                 u(k,j,i)     +                &
+                                                 u(k,j,i+1) )                  &
+                                   )**2 +                                      &
+                                   v(k,j,i)**2 +                               &
+                                   ( 0.25_wp * ( w(k-1,j-1,i) +                &
+                                                 w(k-1,j,i)   +                &
+                                                 w(k,j-1,i)   +                &
+                                                 w(k,j,i) )                    &
+                                   )**2                                        &
+                                 ) *                                           &
+                             v(k,j,i)
+
+!
+!--             Calculate preliminary new velocity, based on pre_tend
+                pre_v = v(k,j,i) + dt_3d * pre_tend
+!
+!--             Compare sign of old velocity and new preliminary velocity,
+!--             and in case the signs are different, limit the tendency
+                IF ( SIGN(pre_v,v(k,j,i)) /= pre_v )  THEN
+                   pre_tend = - v(k,j,i) * ddt_3d
+                ELSE
+                   pre_tend = pre_tend
+                ENDIF
+!
+!--             Calculate final tendency
+                tend(k,j,i) = tend(k,j,i) + pre_tend
+             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_wp )**2 +        &
-                                    ( ( v(k,j,i)    +        &
-                                        v(k,j+1,i)  +        &
-                                        v(k+1,j,i)  +        &
-                                        v(k+1,j+1,i) )       &
-                                      / 4.0_wp )**2 +        &
-                                        w(k,j,i)**2 ) *      &
-                              w(k,j,i)
-   
-          ENDDO
+          CASE ( 3 )
+             DO  k = nzb_w_inner(j,i)+1, pch_index-1
+
+                pre_tend = 0.0_wp
+                pre_w = 0.0_wp
+!
+!--             Calculate preliminary value (pre_tend) of the tendency
+                pre_tend = - cdc *                                             &
+                             (0.5_wp *                                         &
+                                ( lad_s(k+1,j,i) + lad_s(k,j,i) )) *           &
+                             SQRT( ( 0.25_wp * ( u(k,j,i)    +                 &  
+                                                 u(k,j,i+1)  +                 &
+                                                 u(k+1,j,i)  +                 &
+                                                 u(k+1,j,i+1) )                &
+                                   )**2 +                                      &
+                                   ( 0.25_wp * ( v(k,j,i)    +                 &
+                                                 v(k,j+1,i)  +                 &
+                                                 v(k+1,j,i)  +                 &
+                                                 v(k+1,j+1,i) )                &
+                                   )**2 +                                      &
+                                   w(k,j,i)**2                                 &
+                                 ) *                                           &
+                             w(k,j,i)
+!
+!--             Calculate preliminary new velocity, based on pre_tend
+                pre_w = w(k,j,i) + dt_3d * pre_tend
+!
+!--             Compare sign of old velocity and new preliminary velocity,
+!--             and in case the signs are different, limit the tendency
+                IF ( SIGN(pre_w,w(k,j,i)) /= pre_w )  THEN
+                   pre_tend = - w(k,j,i) * ddt_3d
+                ELSE
+                   pre_tend = pre_tend
+                ENDIF
+!
+!--             Calculate final tendency
+                tend(k,j,i) = tend(k,j,i) + pre_tend
+             ENDDO
 
 !
@@ -324,8 +863,7 @@
           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)
+                tend(k,j,i) = tend(k,j,i) +                                    &
+                                 ( canopy_heat_flux(k,j,i) -                   &
+                                   canopy_heat_flux(k-1,j,i) ) / dzw(k)
              ENDDO
 
@@ -335,35 +873,39 @@
           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_wp )**2     +       &
-                                    ( ( v(k,j,i)        +       &
-                                        v(k,j+1,i) )            &
-                                      / 2.0_wp )**2     +       &
-                                    ( ( w(k-1,j,i)      +       &
-                                        w(k,j,i) )              &
-                                      / 2.0_wp )**2 )   *       &
-                              ( q(k,j,i) - sls(k,j,i) )
+                tend(k,j,i) = tend(k,j,i) -                                    &
+                                 lsec *                                        &
+                                 lad_s(k,j,i) *                                &
+                                 SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
+                                                    u(k,j,i+1) )               &
+                                       )**2  +                                 &
+                                       ( 0.5_wp * ( v(k,j,i) +                 &
+                                                    v(k,j+1,i) )               &
+                                       )**2 +                                  &
+                                       ( 0.5_wp * ( w(k-1,j,i) +               &
+                                                    w(k,j,i) )                 &
+                                       )**2                                    &
+                                     ) *                                       &
+                                 ( q(k,j,i) - lsc )
              ENDDO   
 
 !
 !--       sgs-tke
-       CASE ( 6 )
-          DO  k = nzb_s_inner(j,i)+1, pch_index   
-             tend(k,j,i) = tend(k,j,i) -                        &
-                              2.0_wp * cdc(k,j,i) * lad_s(k,j,i) * &
-                              SQRT( ( ( u(k,j,i)           +    &
-                                        u(k,j,i+1) )            &
-                                      / 2.0_wp )**2        +    &  
-                                    ( ( v(k,j,i)           +    &
-                                        v(k,j+1,i) )            &
-                                      / 2.0_wp )**2        +    &
-                                    ( ( w(k,j,i)           +    &
-                                        w(k+1,j,i) )            &
-                                      / 2.0_wp )**2 )      *    &
-                              e(k,j,i)
-          ENDDO
+          CASE ( 6 )
+             DO  k = nzb_s_inner(j,i)+1, pch_index   
+                tend(k,j,i) = tend(k,j,i) -                                    &
+                                 2.0_wp * cdc *                                &
+                                 lad_s(k,j,i) *                                &
+                                 SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
+                                                    u(k,j,i+1) )               &
+                                       )**2 +                                  &  
+                                       ( 0.5_wp * ( v(k,j,i) +                 &
+                                                    v(k,j+1,i) )               &
+                                       )**2 +                                  &
+                                       ( 0.5_wp * ( w(k,j,i) +                 &
+                                                    w(k+1,j,i) )               &
+                                       )**2                                    &
+                                     ) *                                       &
+                                 e(k,j,i)
+             ENDDO
 
        CASE DEFAULT
Index: palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- palm/trunk/SOURCE/prognostic_equations.f90	(revision 1483)
+++ palm/trunk/SOURCE/prognostic_equations.f90	(revision 1484)
@@ -20,5 +20,8 @@
 ! Current revisions:
 ! ------------------
-! 
+! Changes due to new module structure of the plant canopy model:
+!   parameters cthf and plant_canopy moved to module plant_canopy_model_mod.
+! Removed double-listing of use_upstream_for_tke in ONLY-list of module
+! control_parameters
 ! 
 ! Former revisions:
@@ -163,15 +166,15 @@
     USE control_parameters,                                                    &
         ONLY:  call_microphysics_at_all_substeps, cloud_physics,               &
-               constant_diffusion, cthf, dp_external,                          &
+               constant_diffusion, dp_external,                                &
                dp_level_ind_b, dp_smooth_factor, dpdxy, dt_3d, humidity,       &
                icloud_scheme, inflow_l, intermediate_timestep_count,           &
                intermediate_timestep_count_max, large_scale_forcing,           &
                large_scale_subsidence, neutral, nudging, ocean, outflow_l,     &
-               outflow_s, passive_scalar, plant_canopy, precipitation,         &
+               outflow_s, passive_scalar, precipitation,                       &
                prho_reference, prho_reference, prho_reference, pt_reference,   &
                pt_reference, pt_reference, radiation, scalar_advec,            &
                scalar_advec, simulated_time, sloping_surface, timestep_scheme, &
                tsc, use_subsidence_tendencies, use_upstream_for_tke,           &
-               use_upstream_for_tke, use_upstream_for_tke, wall_heatflux,      &
+               wall_heatflux,                                                  &
                wall_nrflux, wall_qflux, wall_qflux, wall_qflux, wall_qrflux,   &
                wall_salinityflux, ws_scheme_mom, ws_scheme_sca
@@ -257,5 +260,5 @@
 
     USE plant_canopy_model_mod,                                                &
-        ONLY:  plant_canopy_model
+        ONLY:  cthf, plant_canopy, plant_canopy_model
 
     USE production_e_mod,                                                      &
Index: palm/trunk/SOURCE/read_var_list.f90
===================================================================
--- palm/trunk/SOURCE/read_var_list.f90	(revision 1483)
+++ palm/trunk/SOURCE/read_var_list.f90	(revision 1484)
@@ -20,5 +20,9 @@
 ! Current revisions:
 ! ------------------
-!
+! Changes in the course of the canopy-model modularization:
+!   parameters alpha_lad, beta_lad, lai_beta added,
+!   module plant_canopy_model_mod added,
+!   drag_coefficient, leaf_surface_concentration and scalar_exchange_coefficient
+!   renamed to canopy_drag_coeff, leaf_surface_conc and leaf_scalar_exch_coeff
 !
 ! Former revisions:
@@ -117,5 +121,5 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  inflow_damping_factor, lad, mean_inflow_profiles, pt_init,      &
+        ONLY:  inflow_damping_factor, mean_inflow_profiles, pt_init,           &
                q_init, ref_state, sa_init, u_init, ug, v_init, vg
 
@@ -141,4 +145,11 @@
 
     USE pegrid
+
+    USE plant_canopy_model_mod,                                                &
+        ONLY:  alpha_lad, beta_lad, canopy_drag_coeff, canopy_mode, cthf, lad, &
+               lad_surface, lad_vertical_gradient, lad_vertical_gradient_level,&
+               lad_vertical_gradient_level_ind,                                &
+               lai_beta, leaf_scalar_exch_coeff, leaf_surface_conc, pch_index, &
+               plant_canopy
 
     USE statistics,                                                            &
@@ -242,4 +253,6 @@
           CASE ( 'advected_distance_y' )
              READ ( 13 )  advected_distance_y
+          CASE ( 'alpha_lad' )
+             READ ( 13 )  alpha_lad
           CASE ( 'alpha_surface' )
              READ ( 13 )  alpha_surface
@@ -282,4 +295,6 @@
           CASE ( 'bc_uv_t' )
              READ ( 13 )  bc_uv_t
+          CASE ( 'beta_lad' )
+             READ ( 13 )  beta_lad
           CASE ( 'bottom_salinityflux' )
              READ ( 13 )  bottom_salinityflux
@@ -296,4 +311,6 @@
           CASE ( 'call_psolver_at_all_substeps' )
              READ ( 13 )  call_psolver_at_all_substeps
+          CASE ( 'canopy_drag_coeff' )
+             READ ( 13 )  canopy_drag_coeff
           CASE ( 'canopy_mode' )
              READ ( 13 )  canopy_mode
@@ -354,6 +371,4 @@
           CASE ( 'dpdxy' )
              READ ( 13 )  dpdxy
-          CASE ( 'drag_coefficient' )
-             READ ( 13 )  drag_coefficient
           CASE ( 'drizzle' )
              READ ( 13 )  drizzle
@@ -419,10 +434,14 @@
           CASE ( 'lad_vertical_gradient_level_in' )
              READ ( 13 )  lad_vertical_gradient_level_ind
+          CASE ( 'lai_beta' )
+             READ ( 13 )  lai_beta
           CASE ( 'large_scale_forcing' )
              READ ( 13 )  large_scale_forcing
            CASE ( 'large_scale_subsidence' )
              READ ( 13 )  large_scale_subsidence
-          CASE ( 'leaf_surface_concentration' )
-             READ ( 13 )  leaf_surface_concentration
+          CASE ( 'leaf_scalar_exch_coeff' )
+             READ ( 13 )  leaf_scalar_exch_coeff
+          CASE ( 'leaf_surface_conc' )
+             READ ( 13 )  leaf_surface_conc
           CASE ( 'limiter_sedimentation' )
              READ ( 13 )  limiter_sedimentation
@@ -558,6 +577,4 @@
           CASE ( 'scalar_advec' )
              READ ( 13 )  scalar_advec
-          CASE ( 'scalar_exchange_coefficient' )
-             READ ( 13 )  scalar_exchange_coefficient
           CASE ( 'simulated_time' )
              READ ( 13 )  simulated_time
@@ -731,5 +748,5 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  inflow_damping_factor, lad, mean_inflow_profiles, pt_init,      &
+        ONLY:  inflow_damping_factor, mean_inflow_profiles, pt_init,           &
                q_init, ref_state, sa_init, u_init, ug, v_init, vg
 
Index: palm/trunk/SOURCE/timestep.f90
===================================================================
--- palm/trunk/SOURCE/timestep.f90	(revision 1483)
+++ palm/trunk/SOURCE/timestep.f90	(revision 1484)
@@ -20,5 +20,9 @@
 ! Current revisions:
 ! ------------------
-! 
+! Changes due to new module structure of the plant canopy model:
+!   calculations and parameters related to the plant canopy model removed
+!   (the limitation of the canopy drag, i.e. that the canopy drag itself should
+!   not change the sign of the velocity components, is now assured for in the
+!   calculation of the canopy tendency terms in subroutine plant_canopy_model) 
 ! 
 ! Former revisions:
@@ -76,5 +80,5 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  cdc, dzu, dzw, kh, km, lad_u, lad_v, lad_w, u, v, w
+        ONLY:  dzu, dzw, kh, km, u, v, w
 
     USE cloud_parameters,                                                      &
@@ -83,5 +87,5 @@
     USE control_parameters,                                                    &
         ONLY:  cfl_factor, coupling_mode, dt_3d, dt_fixed, dt_max,             &
-               galilei_transformation, old_dt, plant_canopy, message_string,   &
+               galilei_transformation, old_dt, message_string,                 &
                stop_dt, terminate_coupled, terminate_coupled_remote,           &
                timestep_reason, u_gtrans, use_ug_for_galilei_tr, v_gtrans
@@ -115,9 +119,4 @@
     REAL(wp) ::  dt_diff           !: 
     REAL(wp) ::  dt_diff_l         !: 
-    REAL(wp) ::  dt_plant_canopy   !: 
-    REAL(wp) ::  dt_plant_canopy_l !: 
-    REAL(wp) ::  dt_plant_canopy_u !: 
-    REAL(wp) ::  dt_plant_canopy_v !: 
-    REAL(wp) ::  dt_plant_canopy_w !: 
     REAL(wp) ::  dt_u              !: 
     REAL(wp) ::  dt_u_l            !: 
@@ -413,94 +412,8 @@
 
 !
-!--    Additional timestep criterion with plant canopies: 
-!--    it is not allowed to extract more than the available momentum
-       IF ( plant_canopy ) THEN
-
-          dt_plant_canopy_l = 0.0_wp
-          DO  i = nxl, nxr
-             DO  j = nys, nyn
-                DO k = nzb+1, nzt
-                   dt_plant_canopy_u = 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_wp )**2     +  &
-                                             ( ( w(k-1,j,i-1)    +  &
-                                                 w(k-1,j,i)      +  &
-                                                 w(k,j,i-1)      +  &
-                                                 w(k,j,i) )         &
-                                                 / 4.0_wp )**2 ) 
-                   IF ( dt_plant_canopy_u > dt_plant_canopy_l ) THEN
-                      dt_plant_canopy_l = dt_plant_canopy_u  
-                   ENDIF
-                   dt_plant_canopy_v = 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_wp )**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_wp )**2 ) 
-                   IF ( dt_plant_canopy_v > dt_plant_canopy_l ) THEN
-                      dt_plant_canopy_l = dt_plant_canopy_v
-                   ENDIF                   
-                   dt_plant_canopy_w = 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_wp )**2     +  &
-                                             ( ( v(k,j,i)        +  &
-                                                 v(k,j+1,i)      +  &
-                                                 v(k+1,j,i)      +  &
-                                                 v(k+1,j+1,i) )     &
-                                               / 4.0_wp )**2        +  &
-                                                 w(k,j,i)**2 )     
-                   IF ( dt_plant_canopy_w > dt_plant_canopy_l ) THEN
-                      dt_plant_canopy_l = dt_plant_canopy_w
-                   ENDIF
-                ENDDO
-             ENDDO
-          ENDDO 
-
-          IF ( dt_plant_canopy_l > 0.0_wp ) THEN
-!
-!--          Invert dt_plant_canopy_l and apply a security timestep factor 0.1
-             dt_plant_canopy_l = 0.1_wp / dt_plant_canopy_l
-          ELSE
-!
-!--          In case of inhomogeneous plant canopy, some processors may have no
-!--          canopy at all. Then use dt_max as dummy instead.
-             dt_plant_canopy_l = dt_max
-          ENDIF
-
-!
-!--       Determine the global minumum
-#if defined( __parallel )
-          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
-          CALL MPI_ALLREDUCE( dt_plant_canopy_l, dt_plant_canopy, 1, MPI_REAL, &
-                              MPI_MIN, comm2d, ierr )
-#else
-          dt_plant_canopy = dt_plant_canopy_l
-#endif
-
-       ELSE
-!
-!--       Use dt_diff as dummy value to avoid extra IF branches further below
-          dt_plant_canopy = dt_diff
-
-       ENDIF
-
-!
 !--    The time step is the minimum of the 3-4 components and the diffusion time
 !--    step minus a reduction (cfl_factor) to be on the safe side.
 !--    The time step must not exceed the maximum allowed value.
-       dt_3d = cfl_factor * MIN( dt_diff, dt_plant_canopy, dt_u, dt_v, dt_w,   &
+       dt_3d = cfl_factor * MIN( dt_diff, dt_u, dt_v, dt_w,   &
                                  dt_precipitation )
        dt_3d = MIN( dt_3d, dt_max )
@@ -508,8 +421,6 @@
 !
 !--    Remember the restricting time step criterion for later output.
-       IF ( MIN( dt_u, dt_v, dt_w ) < MIN( dt_diff, dt_plant_canopy ) )  THEN
+       IF ( MIN( dt_u, dt_v, dt_w ) < dt_diff )  THEN
           timestep_reason = 'A'
-       ELSEIF ( dt_plant_canopy < dt_diff )  THEN
-          timestep_reason = 'C'
        ELSE
           timestep_reason = 'D'
@@ -528,5 +439,4 @@
                '&dt_w            = ', dt_w, ' s',                             &
                '&dt_diff         = ', dt_diff, ' s',                          &
-               '&dt_plant_canopy = ', dt_plant_canopy, ' s',                  &
                '&u_max           = ', u_max, ' m/s   k=', u_max_ijk(1),       &
                '  j=', u_max_ijk(2), '  i=', u_max_ijk(3),                    &
Index: palm/trunk/SOURCE/user_init_plant_canopy.f90
===================================================================
--- palm/trunk/SOURCE/user_init_plant_canopy.f90	(revision 1483)
+++ palm/trunk/SOURCE/user_init_plant_canopy.f90	(revision 1484)
@@ -20,5 +20,8 @@
 ! Current revisions:
 ! -----------------
-! 
+! Changes in the course of the canopy-model modularization:
+!   module plant_canopy_model_mod added,
+!   definition of array cdc (canopy drag coefficient) removed, since it is now
+!   defined purely as a single constant value (see module plant_canopy_model)
 ! 
 ! Former revisions:
@@ -57,4 +60,6 @@
     
     USE kinds
+
+    USE plant_canopy_model_mod
     
     USE user
@@ -62,6 +67,6 @@
     IMPLICIT NONE
 
-    INTEGER(iwp) :: i   !: 
-    INTEGER(iwp) :: j   !: 
+    INTEGER(iwp) :: i   !: running index
+    INTEGER(iwp) :: j   !: running index
 
 !
@@ -69,5 +74,5 @@
 
 !
-!-- Set the 3D-arrays lad_s and cdc for user defined canopies
+!-- Set the 3D-array lad_s for user defined canopies
     SELECT CASE ( TRIM( canopy_mode ) )
 
@@ -78,19 +83,20 @@
        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
+!--       Here the user can define his own forest topography. 
+!--       The following lines contain an example, where the plant canopy extends
+!--       only over the second half of the model domain along x.
+!--       Attention: DO-loops have to include the ghost points (nxlg-nxrg, 
+!--       nysg-nyng), because no exchange of ghost point information is intended,
+!--       in order to minimize communication between CPUs
 !          DO  i = nxlg, nxrg
 !             IF ( i >= INT(nx+1/2) ) THEN
 !                DO  j = nysg, nyng
 !                   lad_s(:,j,i) = lad(:)
-!                   cdc(:,j,i)   = drag_coefficient
 !                ENDDO
 !             ELSE
 !                lad_s(:,:,i) = 0.0_wp
-!                cdc(:,:,i)   = 0.0_wp
 !             ENDIF
 !          ENDDO 
+!
 !--       After definition, please
 !--       remove the following three lines!
Index: palm/trunk/SOURCE/write_var_list.f90
===================================================================
--- palm/trunk/SOURCE/write_var_list.f90	(revision 1483)
+++ palm/trunk/SOURCE/write_var_list.f90	(revision 1484)
@@ -20,5 +20,9 @@
 ! Current revisions:
 ! -----------------
-! 
+! Changes in the course of the canopy-model modularization:
+!   parameters alpha_lad, beta_lad, lai_beta added,
+!   module plant_canopy_model_mod added,
+!   drag_coefficient, leaf_surface_concentration and scalar_exchange_coefficient
+!   renamed to canopy_drag_coeff, leaf_surface_conc and leaf_scalar_exch_coeff
 !
 ! Former revisions:
@@ -104,5 +108,5 @@
 
     USE arrays_3d,                                                             &
-        ONLY:  inflow_damping_factor, lad, mean_inflow_profiles, pt_init,      &
+        ONLY:  inflow_damping_factor, mean_inflow_profiles, pt_init,           &
                q_init, ref_state, sa_init, u_init, ug, v_init, vg
     
@@ -128,4 +132,11 @@
     
     USE pegrid
+
+    USE plant_canopy_model_mod,                                                &
+        ONLY:  alpha_lad, beta_lad, canopy_drag_coeff, canopy_mode, cthf, lad, &
+               lad_surface, lad_vertical_gradient, lad_vertical_gradient_level,&
+               lad_vertical_gradient_level_ind,                                &
+               lai_beta, leaf_scalar_exch_coeff, leaf_surface_conc, pch_index, &
+               plant_canopy
     
     USE statistics,                                                            &
@@ -163,4 +174,6 @@
     WRITE ( 14 )  'advected_distance_y           '
     WRITE ( 14 )  advected_distance_y
+    WRITE ( 14 )  'alpha_lad                     '
+    WRITE ( 14 )  alpha_lad
     WRITE ( 14 )  'alpha_surface                 '
     WRITE ( 14 )  alpha_surface
@@ -203,4 +216,6 @@
     WRITE ( 14 )  'bc_uv_t                       '
     WRITE ( 14 )  bc_uv_t
+    WRITE ( 14 )  'beta_lad                      '
+    WRITE ( 14 )  beta_lad
     WRITE ( 14 )  'bottom_salinityflux           '
     WRITE ( 14 )  bottom_salinityflux
@@ -217,4 +232,6 @@
     WRITE ( 14 )  'call_psolver_at_all_substeps  '
     WRITE ( 14 )  call_psolver_at_all_substeps
+    WRITE ( 14 )  'canopy_drag_coeff              '
+    WRITE ( 14 )  canopy_drag_coeff
     WRITE ( 14 )  'canopy_mode                   '
     WRITE ( 14 )  canopy_mode
@@ -275,6 +292,4 @@
     WRITE ( 14 )  'dpdxy                         '
     WRITE ( 14 )  dpdxy
-    WRITE ( 14 )  'drag_coefficient              '
-    WRITE ( 14 )  drag_coefficient
     WRITE ( 14 )  'drizzle                       '
     WRITE ( 14 )  drizzle
@@ -339,10 +354,14 @@
     WRITE ( 14 )  'lad_vertical_gradient_level_in'
     WRITE ( 14 )  lad_vertical_gradient_level_ind
+    WRITE ( 14 )  'lai_beta                      '
+    WRITE ( 14 )  lai_beta
     WRITE ( 14 )  'large_scale_forcing           '
     WRITE ( 14 )  large_scale_forcing
     WRITE ( 14 )  'large_scale_subsidence        '
     WRITE ( 14 )  large_scale_subsidence
-    WRITE ( 14 )  'leaf_surface_concentration    '
-    WRITE ( 14 )  leaf_surface_concentration
+    WRITE ( 14 )  'leaf_scalar_exch_coeff   '
+    WRITE ( 14 )  leaf_scalar_exch_coeff
+    WRITE ( 14 )  'leaf_surface_conc    '
+    WRITE ( 14 )  leaf_surface_conc
     WRITE ( 14 )  'limiter_sedimentation         '
     WRITE ( 14 )  limiter_sedimentation
@@ -475,6 +494,4 @@
     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
