Index: /palm/trunk/SOURCE/plant_canopy_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 4215)
+++ /palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 4216)
@@ -27,4 +27,7 @@
 ! -----------------
 ! $Id$
+! Bugfixes in 3d data output
+! 
+! 4205 2019-08-30 13:25:00Z suehring
 ! Missing working precision + bugfix in calculation of wind speed
 ! 
@@ -75,4 +78,5 @@
 !
 ! @todo - precalculate constant terms in pcm_calc_transpiration_rate
+! @todo - unify variable names (pcm_, pc_, ...)
 !------------------------------------------------------------------------------!
  MODULE plant_canopy_model_mod
@@ -452,5 +456,5 @@
 !> Subroutine for averaging 3D data
 !------------------------------------------------------------------------------!
-SUBROUTINE pcm_3d_data_averaging( mode, variable )
+ SUBROUTINE pcm_3d_data_averaging( mode, variable )
 
 
@@ -602,5 +606,5 @@
     ENDIF
 
-END SUBROUTINE pcm_3d_data_averaging
+ END SUBROUTINE pcm_3d_data_averaging
 
 !------------------------------------------------------------------------------!
@@ -623,17 +627,18 @@
     IMPLICIT NONE
 
-    CHARACTER (LEN=*) ::  variable !< 
-
-    INTEGER(iwp) ::  av     !< 
-    INTEGER(iwp) ::  i      !< 
-    INTEGER(iwp) ::  j      !< 
-    INTEGER(iwp) ::  k      !< 
+    CHARACTER (LEN=*) ::  variable !< treated variable
+
+    INTEGER(iwp) ::  av     !< flag indicating instantaneous or averaged data output
+    INTEGER(iwp) ::  i      !< grid index x-direction
+    INTEGER(iwp) ::  j      !< grid index y-direction
+    INTEGER(iwp) ::  k      !< grid index z-direction
+    INTEGER(iwp) ::  kk     !< grid index z-direction relative to canopy arrays
     INTEGER(iwp) ::  nzb_do !< lower limit of the data output (usually 0)
     INTEGER(iwp) ::  nzt_do !< vertical upper limit of the data output (usually nz_do3d)
 
-    LOGICAL      ::  found !< 
-
-    REAL(wp)     ::  fill_value
-    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< 
+    LOGICAL      ::  found  !< flag indicating if variable is found
+
+    REAL(wp)     ::  fill_value !< fill value
+    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< data output array
 
 
@@ -643,14 +648,18 @@
 
     SELECT CASE ( TRIM( variable ) )
-
+!
+!--    Note, to save memory arrays for heating are allocated from 0:pch_index.
+!--    Thus, output must be relative to these array indices. Further, check
+!--    whether the output is within the vertical output range,
+!--    i.e. nzb_do:nzt_do.
        CASE ( 'pcm_heatrate' )
           IF ( av == 0 )  THEN
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   IF ( pch_index_ji(j,i) /= 0 )  THEN
-                      DO  k = nzb_do, nzt_do
-                         local_pf(i,j,k) = pc_heating_rate(k,j,i)
-                      ENDDO
-                   ENDIF
+                   DO  k = MAX( topo_top_ind(j,i,0)+1, nzb_do ),               &
+                           MIN( topo_top_ind(j,i,0) + pch_index_ji(j,i), nzt_do )
+                      kk = k - topo_top_ind(j,i,0)
+                      local_pf(i,j,k) = pc_heating_rate(kk,j,i)
+                   ENDDO
                 ENDDO
              ENDDO
@@ -658,6 +667,8 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = pcm_heatrate_av(k,j,i)
+                   DO  k = MAX( topo_top_ind(j,i,0)+1, nzb_do ),               &
+                           MIN( topo_top_ind(j,i,0) + pch_index_ji(j,i), nzt_do )
+                      kk = k - topo_top_ind(j,i,0)
+                      local_pf(i,j,k) = pcm_heatrate_av(kk,j,i)
                    ENDDO
                 ENDDO
@@ -669,9 +680,9 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   IF ( pch_index_ji(j,i) /= 0 )  THEN
-                      DO  k = nzb_do, nzt_do
-                         local_pf(i,j,k) = pc_latent_rate(k,j,i)
-                      ENDDO
-                   ENDIF
+                   DO  k = MAX( topo_top_ind(j,i,0)+1, nzb_do ),               &
+                           MIN( topo_top_ind(j,i,0) + pch_index_ji(j,i), nzt_do )
+                      kk = k - topo_top_ind(j,i,0)
+                      local_pf(i,j,k) = pc_latent_rate(kk,j,i)
+                   ENDDO
                 ENDDO
              ENDDO
@@ -679,6 +690,8 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = pcm_latentrate_av(k,j,i)
+                   DO  k = MAX( topo_top_ind(j,i,0)+1, nzb_do ),               &
+                           MIN( topo_top_ind(j,i,0) + pch_index_ji(j,i), nzt_do )
+                      kk = k - topo_top_ind(j,i,0)
+                      local_pf(i,j,k) = pcm_latentrate_av(kk,j,i)
                    ENDDO
                 ENDDO
@@ -690,9 +703,9 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   IF ( pch_index_ji(j,i) /= 0 )  THEN
-                      DO  k = nzb_do, nzt_do
-                         local_pf(i,j,k) = pc_transpiration_rate(k,j,i)
-                      ENDDO
-                   ENDIF
+                   DO  k = MAX( topo_top_ind(j,i,0)+1, nzb_do ),               &
+                           MIN( topo_top_ind(j,i,0) + pch_index_ji(j,i), nzt_do )
+                      kk = k - topo_top_ind(j,i,0)
+                      local_pf(i,j,k) = pc_transpiration_rate(kk,j,i)
+                   ENDDO
                 ENDDO
              ENDDO
@@ -700,6 +713,8 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb_do, nzt_do
-                      local_pf(i,j,k) = pcm_transpirationrate_av(k,j,i)
+                   DO  k = MAX( topo_top_ind(j,i,0)+1, nzb_do ),               &
+                           MIN( topo_top_ind(j,i,0) + pch_index_ji(j,i), nzt_do )
+                      kk = k - topo_top_ind(j,i,0)
+                      local_pf(i,j,k) = pcm_transpirationrate_av(kk,j,i)
                    ENDDO
                 ENDDO
@@ -711,12 +726,25 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   IF ( pch_index_ji(j,i) /= 0 )  THEN
-                      DO  k = nzb_do, nzt_do
-                         IF ( pc_latent_rate(k,j,i) /= 0._wp ) THEN
-                            local_pf(i,j,k) = pc_heating_rate(k,j,i) / &
-                                              pc_latent_rate(k,j,i)
-                         ENDIF
-                      ENDDO
-                   ENDIF
+                   DO  k = MAX( topo_top_ind(j,i,0)+1, nzb_do ),               &
+                           MIN( topo_top_ind(j,i,0) + pch_index_ji(j,i), nzt_do )
+                      kk = k - topo_top_ind(j,i,0)
+                      IF ( pc_latent_rate(kk,j,i) /= 0.0_wp ) THEN
+                         local_pf(i,j,k) = pc_heating_rate(kk,j,i) /           &
+                                           pc_latent_rate(kk,j,i)
+                      ENDIF
+                   ENDDO
+                ENDDO
+             ENDDO
+          ELSE
+             DO  i = nxl, nxr
+                DO  j = nys, nyn
+                   DO  k = MAX( topo_top_ind(j,i,0)+1, nzb_do ),               &
+                           MIN( topo_top_ind(j,i,0) + pch_index_ji(j,i), nzt_do )
+                      kk = k - topo_top_ind(j,i,0)
+                      IF ( pcm_latentrate_av(kk,j,i) /= 0.0_wp ) THEN
+                         local_pf(i,j,k) = pcm_heatrate_av(kk,j,i) /           &
+                                           pcm_latentrate_av(kk,j,i)
+                      ENDIF
+                   ENDDO
                 ENDDO
              ENDDO
@@ -727,9 +755,9 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   IF ( pch_index_ji(j,i) /= 0 )  THEN
-                      DO  k = nzb_do, nzt_do
-                         local_pf(i,j,k) = lad_s(k,j,i)
-                      ENDDO
-                   ENDIF
+                   DO  k = MAX( topo_top_ind(j,i,0)+1, nzb_do ),               &
+                           MIN( topo_top_ind(j,i,0) + pch_index_ji(j,i), nzt_do )
+                      kk = k - topo_top_ind(j,i,0)
+                      local_pf(i,j,k) = lad_s(kk,j,i)
+                   ENDDO
                 ENDDO
              ENDDO
