Index: palm/trunk/SOURCE/plant_canopy_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/plant_canopy_model.f90	(revision 1818)
+++ palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 1826)
@@ -1,3 +1,3 @@
-!> @file plant_canopy_model.f90
+!> @file plant_canopy_model_mod.f90
 !--------------------------------------------------------------------------------!
 ! This file is part of PALM.
@@ -19,5 +19,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Further modularization
 !
 ! Former revisions:
@@ -35,10 +35,10 @@
 ! 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),
+!   initialization of the canopy model (pcm_init),
 !   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
+!   and preliminary velocities in subroutine pcm_tendency,
+!   some redundant MPI communication removed in subroutine pcm_init
 !   (was previously in init_3d_model),
 !   unnecessary 3d-arrays lad_u, lad_v, lad_w removed - lad information on the
@@ -73,7 +73,7 @@
 ! ------------
 !> 1) Initialization of the canopy model, e.g. construction of leaf area density 
-!> profile (subroutine init_plant_canopy).
+!> profile (subroutine pcm_init).
 !> 2) Calculation of sinks and sources of momentum, heat and scalar concentration 
-!> due to canopy elements (subroutine plant_canopy_model).
+!> due to canopy elements (subroutine pcm_tendency).
 !------------------------------------------------------------------------------!
  MODULE plant_canopy_model_mod
@@ -134,27 +134,211 @@
 
     PRIVATE
-    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
-       MODULE PROCEDURE plant_canopy_model
-       MODULE PROCEDURE plant_canopy_model_ij
-    END INTERFACE plant_canopy_model
-
-
+  
+!
+!-- Public functions
+    PUBLIC pcm_check_parameters, pcm_header, pcm_init, pcm_parin, pcm_tendency
+
+!
+!-- Public variables and constants
+    PUBLIC canopy_mode, cthf, dt_plant_canopy, plant_canopy
+
+
+    INTERFACE pcm_check_parameters
+       MODULE PROCEDURE pcm_check_parameters
+    END INTERFACE pcm_check_parameters     
+    
+     INTERFACE pcm_header
+       MODULE PROCEDURE pcm_header
+    END INTERFACE pcm_header        
+    
+    INTERFACE pcm_init
+       MODULE PROCEDURE pcm_init
+    END INTERFACE pcm_init
+
+    INTERFACE pcm_parin
+       MODULE PROCEDURE pcm_parin
+    END INTERFACE pcm_parin  
+    
+    INTERFACE pcm_tendency
+       MODULE PROCEDURE pcm_tendency
+       MODULE PROCEDURE pcm_tendency_ij
+    END INTERFACE pcm_tendency
 
 
  CONTAINS
 
-
+ 
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Check parameters routine for plant canopy model
+!------------------------------------------------------------------------------!
+    SUBROUTINE pcm_check_parameters
+
+       USE control_parameters,                                                 &
+           ONLY: cloud_physics, message_string, microphysics_seifert
+                 
+    
+       IMPLICIT NONE
+
+    
+       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.  microphysics_seifert )  THEN
+          message_string = 'plant_canopy = .TRUE. requires cloud_scheme /=' // &
+                          ' seifert_beheng'
+          CALL message( 'check_parameters', 'PA0360', 1, 2, 0, 6, 0 )
+       ENDIF
+
+ 
+    END SUBROUTINE pcm_check_parameters 
+ 
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Header output for plant canopy model
+!------------------------------------------------------------------------------!
+    SUBROUTINE pcm_header ( io )
+
+       USE control_parameters,                                                 &
+           ONLY: dz, passive_scalar
+
+
+       IMPLICIT NONE
+ 
+       CHARACTER (LEN=10) ::  coor_chr            !<
+
+       CHARACTER (LEN=86) ::  coordinates         !<
+       CHARACTER (LEN=86) ::  gradients           !<
+       CHARACTER (LEN=86) ::  leaf_area_density   !<
+       CHARACTER (LEN=86) ::  slices              !<
+ 
+       INTEGER(iwp) :: i                !< 
+       INTEGER(iwp),  INTENT(IN) ::  io !< Unit of the output file
+       INTEGER(iwp) :: k                !<       
+   
+       REAL(wp) ::  canopy_height       !< canopy height (in m)
+       
+       canopy_height = pch_index * dz
+
+       WRITE ( io, 1 )  canopy_mode, canopy_height, pch_index,                 &
+                          canopy_drag_coeff
+       IF ( passive_scalar )  THEN
+          WRITE ( io, 2 )  leaf_scalar_exch_coeff,                             &
+                             leaf_surface_conc
+       ENDIF
+
+!
+!--    Heat flux at the top of vegetation
+       WRITE ( io, 3 )  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, '(F7.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, 4 )  TRIM( coordinates ), TRIM( leaf_area_density ),     &
+                             TRIM( gradients ), TRIM( slices )
+
+       ELSE
+       
+          WRITE ( leaf_area_density, '(F7.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, 5 ) TRIM( coordinates ), TRIM( leaf_area_density ),      &
+                          alpha_lad, beta_lad, lai_beta
+
+       ENDIF  
+
+1 FORMAT (//' Vegetation canopy (drag) model:'/                                &
+              ' ------------------------------'//                              &
+              ' Canopy mode: ', A /                                            &
+              ' Canopy height: ',F6.2,'m (',I4,' grid points)' /               &
+              ' Leaf drag coefficient: ',F6.2 /)
+2 FORMAT (/ ' Scalar exchange coefficient: ',F6.2 /                            &
+              ' Scalar concentration at leaf surfaces in kg/m**3: ',F6.2 /)
+3 FORMAT (' Predefined constant heatflux at the top of the vegetation: ',F6.2, &
+          ' K m/s')
+4 FORMAT (/ ' Characteristic levels of the leaf area density:'//               &
+              ' Height:              ',A,'  m'/                                &
+              ' Leaf area density:   ',A,'  m**2/m**3'/                        &
+              ' Gradient:            ',A,'  m**2/m**4'/                        &
+              ' Gridpoint:           ',A)
+5 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' /)   
+       
+    END SUBROUTINE pcm_header
+ 
+ 
 !------------------------------------------------------------------------------!
 ! Description:
@@ -162,5 +346,5 @@
 !> Initialization of the plant canopy model
 !------------------------------------------------------------------------------!
-    SUBROUTINE init_plant_canopy
+    SUBROUTINE pcm_init
     
 
@@ -187,4 +371,12 @@
        pre_lad = 0.0_wp
 
+!
+!--    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 )  THEN
+          calc_beta_lad_profile = .TRUE.
+       ENDIF
+       
+       
 !
 !--    Compute the profile of leaf area density used in the plant
@@ -244,7 +436,8 @@
           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 ) )
+                      ( ( ( 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
 
@@ -252,9 +445,9 @@
 !--       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
+             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
 
@@ -376,5 +569,5 @@
 
 
-    END SUBROUTINE init_plant_canopy
+    END SUBROUTINE pcm_init
 
 
@@ -387,5 +580,5 @@
 !>
 !> 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. 
+!> (defined on scalar grid), as initialized in subroutine pcm_init. 
 !> 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 
@@ -404,5 +597,5 @@
 !> Call for all grid points
 !------------------------------------------------------------------------------!
-    SUBROUTINE plant_canopy_model( component )
+    SUBROUTINE pcm_tendency( component )
 
 
@@ -672,9 +865,56 @@
 
              WRITE( message_string, * ) 'wrong component: ', component
-             CALL message( 'plant_canopy_model', 'PA0279', 1, 2, 0, 6, 0 ) 
+             CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 ) 
 
        END SELECT
 
-    END SUBROUTINE plant_canopy_model
+    END SUBROUTINE pcm_tendency
+
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Parin for &canopy_par for plant canopy model
+!------------------------------------------------------------------------------!
+    SUBROUTINE pcm_parin
+
+
+       IMPLICIT NONE
+
+       CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file 
+       
+        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
+       
+       line = ' '
+       
+!
+!--    Try to find radiation model 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 the radiation model is switched on
+       plant_canopy = .TRUE.
+
+ 10    CONTINUE
+       
+
+    END SUBROUTINE pcm_parin    
+    
 
 
@@ -686,5 +926,5 @@
 !>
 !> 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. 
+!> (defined on scalar grid), as initialized in subroutine pcm_init. 
 !> 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 
@@ -703,5 +943,5 @@
 !> Call for grid point i,j
 !------------------------------------------------------------------------------!
-    SUBROUTINE plant_canopy_model_ij( i, j, component )
+    SUBROUTINE pcm_tendency_ij( i, j, component )
 
 
@@ -944,9 +1184,9 @@
 
           WRITE( message_string, * ) 'wrong component: ', component
-          CALL message( 'plant_canopy_model', 'PA0279', 1, 2, 0, 6, 0 ) 
+          CALL message( 'pcm_tendency', 'PA0279', 1, 2, 0, 6, 0 ) 
 
        END SELECT
 
-    END SUBROUTINE plant_canopy_model_ij
+    END SUBROUTINE pcm_tendency_ij
 
  END MODULE plant_canopy_model_mod
