Index: palm/trunk/SOURCE/plant_canopy_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 3885)
+++ palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 4127)
@@ -27,4 +27,10 @@
 ! -----------------
 ! $Id$
+! Output of 3D plant canopy variables changed. It is now relative to the local
+! terrain rather than located at the acutal vertical level in the model. This 
+! way, the vertical dimension of the output can be significantly reduced. 
+! (merge from branch resler) 
+! 
+! 3885 2019-04-11 11:29:34Z kanani
 ! Changes related to global restructuring of location messages and introduction 
 ! of additional debug messages
@@ -311,9 +317,13 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE ::  pre_lad        !< preliminary lad
     
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  cum_lai_hf             !< cumulative lai for heatflux calc.
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  lad_s                  !< lad on scalar-grid
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pc_heating_rate        !< plant canopy heating rate
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pc_transpiration_rate  !< plant canopy transpiration rate
-    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pc_latent_rate         !< plant canopy latent heating rate
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  cum_lai_hf               !< cumulative lai for heatflux calc.
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  lad_s                    !< lad on scalar-grid
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pc_heating_rate          !< plant canopy heating rate
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pc_transpiration_rate    !< plant canopy transpiration rate
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pc_latent_rate           !< plant canopy latent heating rate
+
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pcm_heatrate_av          !< array for averaging plant canopy sensible heating rate
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pcm_latentrate_av        !< array for averaging plant canopy latent heating rate
+    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  pcm_transpirationrate_av !< array for averaging plant canopy transpiration rate
 
     SAVE
@@ -325,5 +335,6 @@
 !-- Public functions
     PUBLIC pcm_calc_transpiration_rate, pcm_check_data_output,                &
-           pcm_check_parameters, pcm_data_output_3d, pcm_define_netcdf_grid,  &
+           pcm_check_parameters, pcm_3d_data_averaging,                       &
+           pcm_data_output_3d, pcm_define_netcdf_grid,                        &
            pcm_header, pcm_init, pcm_parin, pcm_tendency
 
@@ -345,4 +356,8 @@
        MODULE PROCEDURE pcm_check_parameters
     END INTERFACE pcm_check_parameters
+
+    INTERFACE pcm_3d_data_averaging
+       MODULE PROCEDURE pcm_3d_data_averaging
+    END INTERFACE pcm_3d_data_averaging
 
     INTERFACE pcm_data_output_3d
@@ -618,5 +633,166 @@
 ! Description:
 ! ------------
-!> Subroutine defining 3D output variables
+!> Subroutine for averaging 3D data
+!------------------------------------------------------------------------------!
+SUBROUTINE pcm_3d_data_averaging( mode, variable )
+
+
+    USE control_parameters
+
+    USE indices
+
+    USE kinds
+
+    IMPLICIT NONE
+
+    CHARACTER (LEN=*) ::  mode    !<
+    CHARACTER (LEN=*) :: variable !<
+
+    INTEGER(iwp) ::  i            !<
+    INTEGER(iwp) ::  j            !<
+    INTEGER(iwp) ::  k            !<
+
+
+    IF ( mode == 'allocate' )  THEN
+
+       SELECT CASE ( TRIM( variable ) )
+
+          CASE ( 'pcm_heatrate' )
+             IF ( .NOT. ALLOCATED( pcm_heatrate_av ) )  THEN
+                ALLOCATE( pcm_heatrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
+             ENDIF
+             pcm_heatrate_av = 0.0_wp
+
+
+          CASE ( 'pcm_latentrate' )
+             IF ( .NOT. ALLOCATED( pcm_latentrate_av ) )  THEN
+                ALLOCATE( pcm_latentrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
+             ENDIF
+             pcm_latentrate_av = 0.0_wp
+
+
+          CASE ( 'pcm_transpirationrate' )
+             IF ( .NOT. ALLOCATED( pcm_transpirationrate_av ) )  THEN
+                ALLOCATE( pcm_transpirationrate_av(0:pch_index,nysg:nyng,nxlg:nxrg) )
+             ENDIF
+             pcm_transpirationrate_av = 0.0_wp
+
+          CASE DEFAULT
+             CONTINUE
+
+       END SELECT
+
+    ELSEIF ( mode == 'sum' )  THEN
+
+       SELECT CASE ( TRIM( variable ) )
+
+          CASE ( 'pcm_heatrate' )
+             IF ( ALLOCATED( pcm_heatrate_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      IF ( pch_index_ji(j,i) /= 0 )  THEN
+                         DO  k = 0, pch_index_ji(j,i)
+                            pcm_heatrate_av(k,j,i) = pcm_heatrate_av(k,j,i) + pc_heating_rate(k,j,i)
+                         ENDDO
+                      ENDIF
+                   ENDDO
+                ENDDO
+             ENDIF
+
+
+          CASE ( 'pcm_latentrate' )
+             IF ( ALLOCATED( pcm_latentrate_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      IF ( pch_index_ji(j,i) /= 0 )  THEN
+                         DO  k = 0, pch_index_ji(j,i)
+                            pcm_latentrate_av(k,j,i) = pcm_latentrate_av(k,j,i) + pc_latent_rate(k,j,i)
+                         ENDDO
+                      ENDIF
+                   ENDDO
+                ENDDO
+             ENDIF
+
+
+          CASE ( 'pcm_transpirationrate' )
+             IF ( ALLOCATED( pcm_transpirationrate_av ) ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+                      IF ( pch_index_ji(j,i) /= 0 )  THEN
+                         DO  k = 0, pch_index_ji(j,i)
+                            pcm_transpirationrate_av(k,j,i) = pcm_transpirationrate_av(k,j,i) + pc_transpiration_rate(k,j,i)
+                         ENDDO
+                      ENDIF
+                   ENDDO
+                ENDDO
+             ENDIF
+
+          CASE DEFAULT
+             CONTINUE
+
+       END SELECT
+
+    ELSEIF ( mode == 'average' )  THEN
+
+       SELECT CASE ( TRIM( variable ) )
+
+          CASE ( 'pcm_heatrate' )
+             IF ( ALLOCATED( pcm_heatrate_av ) ) THEN
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      IF ( pch_index_ji(j,i) /= 0 )  THEN
+                         DO  k = 0, pch_index_ji(j,i)
+                            pcm_heatrate_av(k,j,i) = pcm_heatrate_av(k,j,i)                 &
+                                                     / REAL( average_count_3d, KIND=wp )
+                         ENDDO
+                      ENDIF
+                   ENDDO
+                ENDDO
+             ENDIF
+
+
+          CASE ( 'pcm_latentrate' )
+             IF ( ALLOCATED( pcm_latentrate_av ) ) THEN
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      IF ( pch_index_ji(j,i) /= 0 )  THEN
+                         DO  k = 0, pch_index_ji(j,i)
+                            pcm_latentrate_av(k,j,i) = pcm_latentrate_av(k,j,i)              &
+                                                       / REAL( average_count_3d, KIND=wp )
+                         ENDDO
+                      ENDIF
+                   ENDDO
+                ENDDO
+             ENDIF
+
+
+          CASE ( 'pcm_transpirationrate' )
+             IF ( ALLOCATED( pcm_transpirationrate_av ) ) THEN
+                DO  i = nxlg, nxrg
+                   DO  j = nysg, nyng
+                      IF ( pch_index_ji(j,i) /= 0 )  THEN
+                         DO  k = 0, pch_index_ji(j,i)
+                            pcm_transpirationrate_av(k,j,i) = pcm_transpirationrate_av(k,j,i)  &
+                                                              / REAL( average_count_3d, KIND=wp )
+                         ENDDO
+                      ENDIF
+                   ENDDO
+                ENDDO
+             ENDIF
+
+       END SELECT
+
+    ENDIF
+
+END SUBROUTINE pcm_3d_data_averaging
+
+!------------------------------------------------------------------------------!
+!
+! Description:
+! ------------
+!> Subroutine defining 3D output variables. 
+!> Note, 3D plant-canopy output has it's own vertical output dimension, meaning
+!> that 3D output is relative to the model surface now rather than at the actual
+!> grid point where the plant canopy is located. 
 !------------------------------------------------------------------------------!
  SUBROUTINE pcm_data_output_3d( av, variable, found, local_pf, fill_value,     &
@@ -636,5 +812,4 @@
     INTEGER(iwp) ::  j      !< 
     INTEGER(iwp) ::  k      !< 
-    INTEGER(iwp) ::  k_topo !< topography top index
     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)
@@ -652,78 +827,96 @@
     SELECT CASE ( TRIM( variable ) )
 
-      CASE ( 'pcm_heatrate' )
-         IF ( av == 0 )  THEN
-            DO  i = nxl, nxr
-               DO  j = nys, nyn
-                  IF ( pch_index_ji(j,i) /= 0 )  THEN
-                     k_topo = get_topography_top_index_ji( j, i, 's' )
-                     DO  k = k_topo, k_topo + pch_index_ji(j,i)
-                        local_pf(i,j,k) = pc_heating_rate(k-k_topo,j,i)
-                     ENDDO
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-   
+       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
+                ENDDO
+             ENDDO
+          ELSE
+             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)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE ( 'pcm_latentrate' )
+          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_latent_rate(k,j,i)
+                      ENDDO
+                   ENDIF
+                ENDDO
+             ENDDO
+          ELSE
+             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)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
+
        CASE ( 'pcm_transpirationrate' )
-         IF ( av == 0 )  THEN
-            DO  i = nxl, nxr
-               DO  j = nys, nyn
-                  IF ( pch_index_ji(j,i) /= 0 )  THEN
-                     k_topo = get_topography_top_index_ji( j, i, 's' )
-                     DO  k = k_topo, k_topo + pch_index_ji(j,i)
-                        local_pf(i,j,k) = pc_transpiration_rate(k-k_topo,j,i)
-                     ENDDO
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-       CASE ( 'pcm_latentrate' )
-         IF ( av == 0 )  THEN
-            DO  i = nxl, nxr
-               DO  j = nys, nyn
-                  IF ( pch_index_ji(j,i) /= 0 )  THEN
-                     k_topo = get_topography_top_index_ji( j, i, 's' )
-                     DO  k = k_topo, k_topo + pch_index_ji(j,i)
-                        local_pf(i,j,k) = pc_latent_rate(k-k_topo,j,i)
-                     ENDDO
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
+          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_transpiration_rate(k,j,i)
+                      ENDDO
+                   ENDIF
+                ENDDO
+             ENDDO
+          ELSE
+             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)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDIF
 
        CASE ( 'pcm_bowenratio' )
-         IF ( av == 0 )  THEN
-            DO  i = nxl, nxr
-               DO  j = nys, nyn
-                  IF ( pch_index_ji(j,i) /= 0 )  THEN
-                     k_topo = get_topography_top_index_ji( j, i, 's' )
-                     DO  k = k_topo, k_topo + pch_index_ji(j,i)
-                        IF ( pc_latent_rate(k-k_topo,j,i) /= 0._wp ) THEN
-                           local_pf(i,j,k) = pc_heating_rate(k-k_topo,j,i) / &
-                                             pc_latent_rate(k-k_topo,j,i)
-                        ENDIF
-                     ENDDO
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-
-      CASE ( 'pcm_lad' )
-         IF ( av == 0 )  THEN
-            DO  i = nxl, nxr
-               DO  j = nys, nyn
-                  IF ( pch_index_ji(j,i) /= 0 )  THEN
-                     k_topo = get_topography_top_index_ji( j, i, 's' )
-                     DO  k = k_topo, k_topo + pch_index_ji(j,i)
-                        local_pf(i,j,k) = lad_s(k-k_topo,j,i)
-                     ENDDO
-                  ENDIF
-               ENDDO
-            ENDDO
-         ENDIF
-                  
-         
+          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
+                         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
+                ENDDO
+             ENDDO
+          ENDIF
+
+       CASE ( 'pcm_lad' )
+          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) = lad_s(k,j,i)
+                      ENDDO
+                   ENDIF
+                ENDDO
+             ENDDO
+          ENDIF
+
        CASE DEFAULT
           found = .FALSE.
@@ -760,5 +953,5 @@
            grid_x = 'x'
            grid_y = 'y'
-           grid_z = 'zu'
+           grid_z = 'zpc'
 
         CASE DEFAULT
