Index: /palm/trunk/SOURCE/check_parameters.f90
===================================================================
--- /palm/trunk/SOURCE/check_parameters.f90	(revision 2208)
+++ /palm/trunk/SOURCE/check_parameters.f90	(revision 2209)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! -----------------
-! 
+! Check for plant canopy model output
 ! 
 ! Former revisions:
@@ -476,5 +476,5 @@
     USE pegrid
     USE plant_canopy_model_mod,                                                &
-        ONLY:  pcm_check_parameters, plant_canopy
+        ONLY:  pcm_check_data_output, pcm_check_parameters, plant_canopy
         
     USE pmc_interface,                                                         &
@@ -3213,4 +3213,10 @@
              ENDIF
 
+!
+!--          Block of plant canopy model outputs
+             IF ( unit == 'illegal' .AND. plant_canopy .AND. var(1:4) == 'pcm_' ) THEN
+                 CALL pcm_check_data_output( var, unit )
+             ENDIF
+             
              IF ( unit == 'illegal' )  THEN
                 unit = ''
Index: /palm/trunk/SOURCE/data_output_3d.f90
===================================================================
--- /palm/trunk/SOURCE/data_output_3d.f90	(revision 2208)
+++ /palm/trunk/SOURCE/data_output_3d.f90	(revision 2209)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Added plant canopy model output
 ! 
 ! Former revisions:
@@ -188,4 +188,7 @@
     USE pegrid
 
+    USE plant_canopy_model_mod,                                                &
+        ONLY:  pcm_data_output_3d, plant_canopy
+        
     USE radiation_model_mod,                                                   &
         ONLY:  radiation, radiation_data_output_3d
@@ -620,4 +623,12 @@
              IF ( .NOT. found  .AND.  radiation )  THEN
                 CALL radiation_data_output_3d( av, do3d(av,if), found,         &
+                                               local_pf )
+                resorted = .TRUE.
+             ENDIF
+
+!
+!--          Plant canopy model output
+             IF ( .NOT. found  .AND.  plant_canopy )  THEN
+                CALL pcm_data_output_3d( av, do3d(av,if), found,         &
                                                local_pf )
                 resorted = .TRUE.
Index: /palm/trunk/SOURCE/netcdf_interface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 2208)
+++ /palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 2209)
@@ -20,5 +20,5 @@
 ! Current revisions:
 ! ------------------
-! 
+! Added support for plant canopy model output
 ! 
 ! Former revisions:
@@ -461,4 +461,7 @@
         ONLY:  maximum_number_of_particles, number_of_particle_groups
 
+    USE plant_canopy_model_mod,                                                &
+        ONLY:  pcm_define_netcdf_grid, plant_canopy
+
     USE profil_parameter,                                                      &
         ONLY:  crmax, cross_profiles, dopr_index, profile_columns, profile_rows
@@ -839,4 +842,11 @@
                                                    grid_x, grid_y, grid_z )
                    ENDIF
+                   
+!
+!--                Check for plant canopy quantities
+                   IF ( plant_canopy )  THEN
+                      CALL pcm_define_netcdf_grid( domask(mid,av,i), found,    &
+                                                   grid_x, grid_y, grid_z )
+                   ENDIF
 
 !
@@ -1361,4 +1371,11 @@
                    ENDIF
 
+!
+!--                Check for plant canopy quantities
+                   IF ( plant_canopy )  THEN
+                      CALL pcm_define_netcdf_grid( do3d(av,i), found, grid_x,  &
+                                                   grid_y, grid_z )
+                   ENDIF
+                   
 !
 !--                Check for radiation quantities
Index: /palm/trunk/SOURCE/plant_canopy_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 2208)
+++ /palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 2209)
@@ -20,5 +20,6 @@
 ! Current revisions:
 ! -----------------
-! 
+! Added 3d output of leaf area density (pcm_lad) and canopy 
+! heat rate (pcm_heatrate)
 ! 
 ! Former revisions:
@@ -166,5 +167,6 @@
 !
 !-- Public functions
-    PUBLIC pcm_check_parameters, pcm_header, pcm_init, pcm_parin, pcm_tendency
+    PUBLIC pcm_check_data_output, pcm_check_parameters, pcm_data_output_3d,    &
+           pcm_define_netcdf_grid, pcm_header, pcm_init, pcm_parin, pcm_tendency
 
 !
@@ -174,8 +176,19 @@
            
 
-
+    INTERFACE pcm_check_data_output
+       MODULE PROCEDURE pcm_check_data_output
+    END INTERFACE pcm_check_data_output
+    
     INTERFACE pcm_check_parameters
        MODULE PROCEDURE pcm_check_parameters
-    END INTERFACE pcm_check_parameters     
+    END INTERFACE pcm_check_parameters
+
+    INTERFACE pcm_data_output_3d
+       MODULE PROCEDURE pcm_data_output_3d
+    END INTERFACE pcm_data_output_3d
+
+    INTERFACE pcm_define_netcdf_grid
+       MODULE PROCEDURE pcm_define_netcdf_grid
+    END INTERFACE pcm_define_netcdf_grid
     
      INTERFACE pcm_header
@@ -203,4 +216,39 @@
  CONTAINS
 
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Check data output for plant canopy model
+!------------------------------------------------------------------------------!
+ SUBROUTINE pcm_check_data_output( var, unit )
+ 
+ 
+    USE control_parameters,                                                 &
+        ONLY:  data_output, message_string
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*) ::  unit  !< 
+    CHARACTER (LEN=*) ::  var   !<
+
+
+    SELECT CASE ( TRIM( var ) )
+
+       CASE ( 'pcm_heatrate' )
+          unit = 'K s-1'
+    
+       CASE ( 'pcm_lad' )
+          unit = 'm2 m-3'
+
+
+       CASE DEFAULT
+          unit = 'illegal'
+
+    END SELECT
+
+
+ END SUBROUTINE pcm_check_data_output
+ 
  
 !------------------------------------------------------------------------------!
@@ -259,4 +307,110 @@
  
 
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Subroutine defining 3D output variables
+!------------------------------------------------------------------------------!
+ SUBROUTINE pcm_data_output_3d( av, variable, found, local_pf )
+ 
+    USE control_parameters,                                                    &
+        ONLY :  nz_do3d
+ 
+    USE indices
+
+    USE kinds
+
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*) ::  variable !< 
+
+    INTEGER(iwp) ::  av    !< 
+    INTEGER(iwp) ::  i     !< 
+    INTEGER(iwp) ::  j     !< 
+    INTEGER(iwp) ::  k     !< 
+
+    LOGICAL      ::  found !< 
+
+    REAL(sp), DIMENSION(nxlg:nxrg,nysg:nyng,nzb:nz_do3d) ::  local_pf !< 
+
+
+    found = .TRUE.
+
+
+    SELECT CASE ( TRIM( variable ) )
+
+      CASE ( 'pcm_heatrate' )
+         IF ( av == 0 )  THEN
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb_s_inner(j,i), nz_do3d
+                     local_pf(i,j,k) = pc_heating_rate(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+    
+    
+      CASE ( 'pcm_lad' )
+
+         IF ( av == 0 )  THEN
+            DO  i = nxlg, nxrg
+               DO  j = nysg, nyng
+                  DO  k = nzb_s_inner(j,i), nz_do3d
+                     local_pf(i,j,k) = lad_s(k,j,i)
+                  ENDDO
+               ENDDO
+            ENDDO
+         ENDIF
+                  
+         
+       CASE DEFAULT
+          found = .FALSE.
+
+    END SELECT
+
+
+ END SUBROUTINE pcm_data_output_3d 
+         
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Subroutine defining appropriate grid for netcdf variables.
+!> It is called from subroutine netcdf.
+!------------------------------------------------------------------------------!
+ SUBROUTINE pcm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
+    
+     IMPLICIT NONE
+
+     CHARACTER (LEN=*), INTENT(IN)  ::  var         !< 
+     LOGICAL, INTENT(OUT)           ::  found       !< 
+     CHARACTER (LEN=*), INTENT(OUT) ::  grid_x      !< 
+     CHARACTER (LEN=*), INTENT(OUT) ::  grid_y      !< 
+     CHARACTER (LEN=*), INTENT(OUT) ::  grid_z      !< 
+
+     found  = .TRUE.
+
+!
+!--  Check for the grid
+     SELECT CASE ( TRIM( var ) )
+
+        CASE ( 'pcm_heatrate', 'pcm_lad' )
+           grid_x = 'x'
+           grid_y = 'y'
+           grid_z = 'zu'
+
+        CASE DEFAULT
+           found  = .FALSE.
+           grid_x = 'none'
+           grid_y = 'none'
+           grid_z = 'none'
+     END SELECT
+
+ END SUBROUTINE pcm_define_netcdf_grid
+ 
+ 
 !------------------------------------------------------------------------------!
 ! Description:
Index: /palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/urban_surface_mod.f90	(revision 2208)
+++ /palm/trunk/SOURCE/urban_surface_mod.f90	(revision 2209)
@@ -21,5 +21,7 @@
 ! Current revisions:
 ! ------------------
-! 
+! cpp switch __mpi3 removed,
+! minor formatting,
+! small bugfix for division by zero (Krc)
 ! 
 ! Former revisions:
@@ -1489,5 +1491,5 @@
         CALL location_message( '    calculation of SVF and CSF', .TRUE. )
 
-#if defined( __mpi3 )
+
 !--     precalculate face areas for different face directions using normal vector
         DO d = 0, 9
@@ -1964,5 +1966,5 @@
         CALL message( 'init_urban_surface', 'PA0502', 2, 2, 0, 6, 0 )
        
-#endif
+
     END SUBROUTINE usm_calc_svf
 
@@ -2487,44 +2489,59 @@
 !------------------------------------------------------------------------------!
     PURE SUBROUTINE usm_find_boundary_face(origin, uvect, bdycross)
-        IMPLICIT NONE
-        REAL(wp), DIMENSION(3), INTENT(in)      :: origin    !< ray origin
-        REAL(wp), DIMENSION(3), INTENT(in)      :: uvect     !< ray unit vector
-        INTEGER(iwp), DIMENSION(4), INTENT(out) :: bdycross  !< found boundary crossing (d, z, y, x)
-        REAL(wp), DIMENSION(3)                  :: crossdist !< crossing distance
-        INTEGER(iwp), DIMENSION(3)              :: bdyd      !< boundary direction
-        REAL(wp)                                :: bdydim    !< 
-        REAL(wp)                                :: dist      !< 
-        INTEGER(iwp)                            :: seldim    !< found fist crossing index
-        INTEGER(iwp)                            :: d         !< 
-
-        bdydim = nzut + .5_wp !< top boundary
-        bdyd(1) = isky
-        crossdist(1) = (bdydim - origin(1)) / uvect(1)
-
-        IF ( uvect(2) >= 0._wp )  THEN
-            bdydim = ny + .5_wp !< north global boundary
-            bdyd(2) = inorthb
-        ELSE
-            bdydim = -.5_wp !< south global boundary
-            bdyd(2) = isouthb
-        ENDIF
-        crossdist(2) = (bdydim - origin(2)) / uvect(2)
-
-        IF ( uvect(3) >= 0._wp )  THEN
-            bdydim = nx + .5_wp !< east global boundary
-            bdyd(3) = ieastb
-        ELSE
-            bdydim = -.5_wp !< west global boundary
-            bdyd(3) = iwestb
-        ENDIF
-        crossdist(3) = (bdydim - origin(3)) / uvect(3)
-
-        seldim = minloc(crossdist, 1)
-        dist = crossdist(seldim)
-        d = bdyd(seldim)
-
-        bdycross(1) = d
-        bdycross(2:4) = NINT( origin(:) + uvect(:)*dist &
-                        + .5_wp * (/ kdir(d), jdir(d), idir(d) /) )
+    
+       IMPLICIT NONE
+       
+       INTEGER(iwp) ::  d       !<
+       INTEGER(iwp) ::  seldim  !< found fist crossing index
+
+       INTEGER(iwp), DIMENSION(3)              ::  bdyd      !< boundary direction       
+       INTEGER(iwp), DIMENSION(4), INTENT(out) ::  bdycross  !< found boundary crossing (d, z, y, x)
+       
+       REAL(wp)                                ::  bdydim  !< 
+       REAL(wp)                                ::  dist    !<
+       
+       REAL(wp), DIMENSION(3)             ::  crossdist  !< crossing distance
+       REAL(wp), DIMENSION(3), INTENT(in) ::  origin     !< ray origin
+       REAL(wp), DIMENSION(3), INTENT(in) ::  uvect      !< ray unit vector
+ 
+
+       bdydim       = nzut + .5_wp  !< top boundary
+       bdyd(1)      = isky
+       crossdist(1) = ( bdydim - origin(1) ) / uvect(1)  !< subroutine called only when uvect(1)>0
+
+       IF ( uvect(2) == 0._wp )  THEN
+          crossdist(2) = huge(1._wp)
+       ELSE
+          IF ( uvect(2) >= 0._wp )  THEN
+             bdydim  = ny + .5_wp  !< north global boundary
+             bdyd(2) = inorthb
+          ELSE
+             bdydim  = -.5_wp  !< south global boundary
+             bdyd(2) = isouthb
+          ENDIF
+          crossdist(2) = ( bdydim - origin(2) ) / uvect(2)
+       ENDIF
+
+       IF ( uvect(3) == 0._wp )  THEN
+          crossdist(3) = huge(1._wp)
+       ELSE
+          IF ( uvect(3) >= 0._wp )  THEN
+             bdydim  = nx + .5_wp  !< east global boundary
+             bdyd(3) = ieastb
+          ELSE
+             bdydim  = -.5_wp  !< west global boundary
+             bdyd(3) = iwestb
+          ENDIF
+          crossdist(3) = ( bdydim - origin(3) ) / uvect(3)
+       ENDIF
+
+       seldim = minloc(crossdist, 1)
+       dist   = crossdist(seldim)
+       d      = bdyd(seldim)
+
+       bdycross(1)   = d
+       bdycross(2:4) = NINT( origin(:) + uvect(:) * dist &
+                                       + .5_wp * (/ kdir(d), jdir(d), idir(d) /) )
+                        
     END SUBROUTINE
 
@@ -2864,15 +2881,4 @@
        urban_surface = .TRUE.
        
-!
-!--    Check whether pre-processor (cpp) option "__mpi3" is set. It is required
-!--    for the full functionality of the USM. "__mpi3" directive is implemented,
-!--    because some compilers cannot handle MPI-3 operations, hence, these parts
-!--    of code shall only be compiled if explicitly enabled.
-#if ! defined ( __mpi3 )
-          message_string = 'urban surface model requires compilation of ' //   & 
-                           'PALM with pre-processor directive -D__mpi3'
-          CALL message( 'usm_parin', 'PA0503', 1, 2, 0, 6, 0 )
-#endif
-
 
  10    CONTINUE
@@ -3220,5 +3226,5 @@
         REAL(wp), PARAMETER                    :: grow_factor = 1.5_wp !< factor of expansion of grow arrays
 
-#if defined( __mpi3 )
+
 !--     Maximum number of gridboxes visited equals to maximum number of boundaries crossed in each dimension plus one. That's also
 !--     the maximum number of plant canopy boxes written. We grow the acsf array accordingly using exponential factor.
@@ -3363,9 +3369,6 @@
         
         visible = .TRUE.
-        
-#else
-        visible      = .FALSE.                          !Set variables to avoid compiler warnimngs
-        transparency = 0.0
-#endif
+
+        
     END SUBROUTINE usm_raytrace
     
@@ -3936,5 +3939,5 @@
         REAL(wp)                              :: acoef              !< actual coefficient of diurnal profile of anthropogenic heat
 
-#if defined( __mpi3 )
+
         dxdir = (/dz,dy,dy,dx,dx/)
         
@@ -4147,5 +4150,5 @@
        ENDIF
 
-#endif
+
     END SUBROUTINE usm_surface_energy_balance
 
