Index: palm/trunk/SOURCE/plant_canopy_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 2008)
+++ palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 2011)
@@ -20,5 +20,7 @@
 ! Current revisions:
 ! -----------------
-! 
+! Renamed canopy_heat_flux to pc_heating_rate, since the original meaning/
+! calculation of the quantity has changed, related to the urban surface model
+! and similar future applications. 
 ! 
 ! Former revisions:
@@ -146,5 +148,5 @@
     
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::                                 &
-       canopy_heat_flux                                    !< canopy heat flux
+       pc_heating_rate                                    !< plant canopy heating rate
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  cum_lai_hf !< cumulative lai for heatflux calc.
     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  lad_s      !< lad on scalar-grid
@@ -162,5 +164,5 @@
 !
 !-- Public variables and constants
-    PUBLIC canopy_heat_flux, canopy_mode, cthf, dt_plant_canopy, lad, lad_s,   &
+    PUBLIC pc_heating_rate, canopy_mode, cthf, dt_plant_canopy, lad, lad_s,   &
            pch_index, plant_canopy
            
@@ -503,5 +505,5 @@
        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) )
+                    pc_heating_rate(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
        ENDIF
 
@@ -548,8 +550,18 @@
 
 !
-!--    Initialization of the canopy heat source distribution
+!--    Initialization of the canopy heat source distribution due to heating
+!--    of the canopy layers by incoming solar radiation, in case that a non-zero
+!--    value is set for the canopy top heat flux (cthf), which equals the
+!--    available net radiation at canopy top.
+!--    The heat source distribution is calculated by a decaying exponential 
+!--    function of the downward cumulative leaf area index (cum_lai_hf), 
+!--    assuming that the foliage inside the plant canopy is heated by solar
+!--    radiation penetrating the canopy layers according to the distribution
+!--    of net radiation as suggested by Brown & Covey (1966; Agric. Meteorol. 3, 
+!--    73â96). This approach has been applied e.g. by Shaw & Schumann (1992;
+!--    Bound.-Layer Meteorol. 61, 47â64)
        IF ( cthf /= 0.0_wp )  THEN
 !
-!--       Piecewise calculation of the leaf area index by vertical
+!--       Piecewise calculation of the cumulative leaf area index by vertical
 !--       integration of the leaf area density
           cum_lai_hf(:,:,:) = 0.0_wp
@@ -581,23 +593,25 @@
 
 !
-!--       Calculation of the upward kinematic vertical heat flux within the 
-!--       canopy
+!--       Calculation of the heating rate (K/s) within the different layers of 
+!--       the plant canopy
           DO  i = nxlg, nxrg
              DO  j = nysg, nyng
-                DO  k = 0, pch_index
-                   canopy_heat_flux(k,j,i) = cthf *                            &
-                                             exp( -ext_coef * cum_lai_hf(k,j,i) )
-                ENDDO
-             ENDDO
-          ENDDO
-
-!
-!--       In areas covered with canopy, the surface heat flux is set to
-!--       the surface value of the above calculated in-canopy heat flux
-!--       distribution
-          DO  i = nxlg,nxrg
-             DO  j = nysg, nyng
-                IF ( canopy_heat_flux(0,j,i) /= cthf )  THEN
-                   shf(j,i) = canopy_heat_flux(0,j,i)
+!
+!--             Calculation only necessary in areas covered with canopy
+                IF ( cum_lai_hf(0,j,i) /= 0.0_wp )  THEN
+!--             
+!--                In areas with canopy the surface value of the canopy heat 
+!--                flux distribution overrides the surface heat flux (shf) 
+                   shf(j,i) = cthf * exp( -ext_coef * cum_lai_hf(0,j,i) )
+!
+!--                Within the different canopy layers the plant-canopy heating
+!--                rate (pc_heating_rate) is calculated as the vertical 
+!--                divergence of the canopy heat fluxes at the top and bottom
+!--                of the respective layer
+                   DO  k = 1, pch_index
+                      pc_heating_rate(k,j,i) = cthf *                         &
+                                                ( exp(-ext_coef*cum_lai_hf(k,j,i)) -  &
+                                                  exp(-ext_coef*cum_lai_hf(k-1,j,i)) ) / dzw(k)
+                   ENDDO
                 ENDIF
              ENDDO
@@ -689,14 +703,5 @@
         INTEGER(iwp)                            :: i, j, dtype, nzp, nzpltop, nzpl, kk
         REAL(wp), DIMENSION(:), ALLOCATABLE     :: col
-       
-        lad_s = 0.0_wp
-!         cdc = 0.0_wp
-!         if ( passive_scalar )  then
-!            lsc = 0.0_wp
-!            lsec = 0.0_wp
-!         endif
-
-        WRITE(9,*) 'Reading PLANT_CANOPY_DATA_3D', nzt
-        FLUSH(9)
+
         OPEN(152, file='PLANT_CANOPY_DATA_3D', access='SEQUENTIAL', &
                 action='READ', status='OLD', form='FORMATTED', err=515)
@@ -709,6 +714,4 @@
             READ(152, *, err=516, end=517) dtype, i, j, col(:)
             IF ( i < nxlg .or. i > nxrg .or. j < nysg .or. j > nyng ) CYCLE
-            WRITE(9,*) 'Read ', i,j,nzb_s_inner(j,i),col(:)
-            FLUSH(9)
 
             SELECT CASE (dtype)
@@ -974,7 +977,5 @@
                    DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
                       kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
-                      tend(k,j,i) = tend(k,j,i) +                              &
-                                       ( canopy_heat_flux(kk,j,i) -            &
-                                         canopy_heat_flux(kk-1,j,i) ) / dzw(k)
+                      tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i)
                    ENDDO
                 ENDDO
@@ -1281,7 +1282,5 @@
              DO  k = nzb_s_inner(j,i)+1, nzb_s_inner(j,i)+pch_index
                 kk = k - nzb_s_inner(j,i)  !- lad arrays are defined flat
-                tend(k,j,i) = tend(k,j,i) +                                    &
-                                 ( canopy_heat_flux(kk,j,i) -                  &
-                                   canopy_heat_flux(kk-1,j,i) ) / dzw(k)
+                tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i)
              ENDDO
 
