Index: /palm/trunk/SOURCE/Makefile
===================================================================
--- /palm/trunk/SOURCE/Makefile	(revision 4126)
+++ /palm/trunk/SOURCE/Makefile	(revision 4127)
@@ -25,4 +25,8 @@
 # -----------------
 # $Id$
+# Add dependency of data_output_3d on plant_canopy_model_mod 
+# (merge from branch resler)
+#
+# 4106 2019-07-19 08:54:42Z gronemeier
 # Remove dependency on pmc_interface for boundary_conds
 # 
@@ -975,4 +979,5 @@
 	modules.o \
 	netcdf_interface_mod.o \
+	plant_canopy_model_mod.o \
 	radiation_model_mod.o \
 	urban_surface_mod.o
Index: /palm/trunk/SOURCE/biometeorology_mod.f90
===================================================================
--- /palm/trunk/SOURCE/biometeorology_mod.f90	(revision 4126)
+++ /palm/trunk/SOURCE/biometeorology_mod.f90	(revision 4127)
@@ -27,4 +27,7 @@
 ! -----------------
 ! $Id$
+! Output for bio_mrt added (merge from branch resler)
+! 
+! 4126 2019-07-30 11:09:11Z gronemeier
 ! renamed vitd3_exposure_av into vitd3_dose,
 ! renamed uvem_calc_exposure into bio_calculate_uv_exposure
@@ -226,4 +229,5 @@
     LOGICAL ::  average_trigger_utci   = .FALSE.  !< update averaged input on call to bio_utci?
     LOGICAL ::  average_trigger_pet    = .FALSE.  !< update averaged input on call to bio_pet?
+    LOGICAL ::  average_trigger_mrt    = .FALSE.  !< update averaged input on call to bio_pet?
     LOGICAL ::  do_calculate_perct     = .FALSE.  !< Turn index PT (instant. input) on or off
     LOGICAL ::  do_calculate_perct_av  = .FALSE.  !< Turn index PT (averaged input) on or off
@@ -232,4 +236,5 @@
     LOGICAL ::  do_calculate_utci      = .FALSE.  !< Turn index UTCI (instant. input) on or off
     LOGICAL ::  do_calculate_utci_av   = .FALSE.  !< Turn index UTCI (averaged input) on or off
+    LOGICAL ::  do_calculate_mrt2d     = .FALSE.  !< Turn index MRT 2D (averaged or inst) on or off
 
 !
@@ -432,5 +437,5 @@
 
 
-          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' )
+          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*', 'bio_mrt*' )
 
 !
@@ -447,5 +452,6 @@
              IF ( .NOT. average_trigger_perct  .AND.                           &
                   .NOT. average_trigger_utci   .AND.                           &
-                  .NOT. average_trigger_pet )  THEN
+                  .NOT. average_trigger_pet    .AND.                           &
+                  .NOT. average_trigger_mrt )  THEN
 !
 !--             Memorize the first index called to control averaging
@@ -458,4 +464,7 @@
                 IF ( TRIM( variable ) == 'bio_pet*' )  THEN
                     average_trigger_pet = .TRUE.
+                ENDIF
+                IF ( TRIM( variable ) == 'bio_mrt*' )  THEN
+                    average_trigger_mrt = .TRUE.
                 ENDIF
              ENDIF
@@ -499,5 +508,5 @@
              ENDIF
 
-          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' )
+          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*', 'bio_mrt*' )
 !
 !--          Only continue if the current index is the one to trigger the input 
@@ -509,4 +518,6 @@
              IF ( average_trigger_pet    .AND.  TRIM( variable ) /=            &
                 'bio_pet*')    RETURN
+             IF ( average_trigger_mrt    .AND.  TRIM( variable ) /=            &
+                'bio_mrt*')    RETURN
 !
 !--          Now memorize which of the input grids are not averaged by other 
@@ -637,5 +648,5 @@
              ENDIF
 
-          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*' )
+          CASE ( 'bio_perct*', 'bio_utci*', 'bio_pet*', 'bio_mrt*' )
 !
 !--          Only continue if update index, see above
@@ -646,4 +657,6 @@
              IF ( average_trigger_pet   .AND.                                  &
                 TRIM( variable ) /= 'bio_pet*' )  RETURN
+             IF ( average_trigger_mrt   .AND.                                  &
+                TRIM( variable ) /= 'bio_mrt*' )  RETURN
 
              IF ( ALLOCATED( pt_av )  .AND.  do_average_theta )  THEN
@@ -748,5 +761,5 @@
 !--    derived in a single step based on priorly averaged arrays (see
 !--    bio_calculate_thermal_index_maps).
-       CASE ( 'bio_mrt' )
+       CASE ( 'bio_mrt', 'bio_mrt*' )
           unit = 'degree_C'
           thermal_comfort = .TRUE.  !< enable thermal_comfort if user forgot to do so
@@ -755,4 +768,7 @@
              tmrt_grid = REAL( bio_fill_value, KIND = wp )
           ENDIF
+          IF ( TRIM( var ) == 'bio_mrt*' )  THEN
+             do_calculate_mrt2d = .TRUE.
+          END IF
 
        CASE ( 'bio_perct*' )
@@ -961,4 +977,21 @@
            ENDDO
 
+        CASE ( 'bio_mrt*_xy' )        ! 2d-array
+           grid = 'zu1'
+           two_d = .TRUE.
+           IF ( av == 0 )  THEN
+              DO  i = nxl, nxr
+                 DO  j = nys, nyn
+                    local_pf(i,j,nzb+1) = tmrt_grid(j,i)
+                 ENDDO
+              ENDDO
+           ELSE
+              DO  i = nxl, nxr
+                 DO  j = nys, nyn
+                    local_pf(i,j,nzb+1) = tmrt_av_grid(j,i)
+                 ENDDO
+              ENDDO
+           ENDIF
+
 
         CASE ( 'bio_perct*_xy' )        ! 2d-array
@@ -1372,4 +1405,7 @@
           READ ( 13 )  average_trigger_pet
 
+       CASE ( 'average_trigger_mrt' )
+          READ ( 13 )  average_trigger_mrt
+
 
        CASE DEFAULT
@@ -1451,4 +1487,6 @@
     CALL wrd_write_string( 'average_trigger_pet' )
     WRITE ( 14 )  average_trigger_pet
+    CALL wrd_write_string( 'average_trigger_mrt' )
+    WRITE ( 14 )  average_trigger_mrt
 
  END SUBROUTINE bio_wrd_global
@@ -1694,5 +1732,6 @@
     IF ( do_calculate_perct  .OR.  do_calculate_perct_av  .OR.                 &
        do_calculate_utci  .OR.  do_calculate_utci_av  .OR.                     &
-       do_calculate_pet  .OR.  do_calculate_pet_av )  THEN
+       do_calculate_pet  .OR.  do_calculate_pet_av  .OR.                       &
+       do_calculate_mrt2d )  THEN
 
 !
Index: /palm/trunk/SOURCE/data_output_3d.f90
===================================================================
--- /palm/trunk/SOURCE/data_output_3d.f90	(revision 4126)
+++ /palm/trunk/SOURCE/data_output_3d.f90	(revision 4127)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Adjustment for top boundary index for plant-canopy model outputs 
+! (merge from branch resler)
+! 
+! 4048 2019-06-21 21:00:21Z knoop
 ! Moved tcm_data_output_3d to module_interface
 ! 
@@ -276,6 +280,7 @@
         ONLY:  debug_output_timestep,                                          &
                do3d, do3d_no, do3d_time_count, io_blocks, io_group,            &
-               land_surface, message_string, ntdim_3d, nz_do3d, psolver,       &
-               time_since_reference_point, urban_surface, varnamelength
+               land_surface, message_string, ntdim_3d, nz_do3d, plant_canopy,  &
+               psolver, time_since_reference_point, urban_surface,             &
+               varnamelength
 
     USE cpulog,                                                                &
@@ -313,4 +318,7 @@
 
     USE pegrid
+
+    USE plant_canopy_model_mod,                                                &
+        ONLY:  pch_index
 
     USE radiation_model_mod,                                                   &
@@ -410,4 +418,6 @@
        found = .FALSE.
        resorted = .FALSE.
+       trimvar = TRIM( do3d(av,ivar) )
+
 !
 !--    Temporary solution to account for data output within the new urban
@@ -415,7 +425,13 @@
 !--    Store the array chosen on the temporary array.
        nzb_do   = nzb
-       nzt_do   = nz_do3d
-
-       trimvar = TRIM( do3d(av,ivar) )
+!
+!--    Set top index for 3D output. Note in case of plant-canopy model
+!--    these index is determined by pch_index. 
+       IF ( plant_canopy  .AND.  trimvar(1:4) == 'pcm_' )  THEN
+          nzt_do   = pch_index
+       ELSE
+          nzt_do   = nz_do3d
+       ENDIF
+
 !
 !--    Allocate a temporary array with the desired output dimensions.
Index: /palm/trunk/SOURCE/module_interface.f90
===================================================================
--- /palm/trunk/SOURCE/module_interface.f90	(revision 4126)
+++ /palm/trunk/SOURCE/module_interface.f90	(revision 4127)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Add output of 3D plant-canopy outputs (merge from branch resler)
+! 
+! 4048 2019-06-21 21:00:21Z knoop
 ! Moved turbulence_closure_mod calls into this module_interface
 ! 
@@ -373,4 +376,5 @@
                pcm_init,                                                       &
                pcm_header,                                                     &
+               pcm_3d_data_averaging,                                          &
                pcm_data_output_3d
 
@@ -1297,4 +1301,5 @@
     IF ( land_surface        )  CALL lsm_3d_data_averaging( mode, variable )
     IF ( ocean_mode          )  CALL ocean_3d_data_averaging( mode, variable )
+    IF ( plant_canopy        )  CALL pcm_3d_data_averaging( mode, variable )
     IF ( radiation           )  CALL radiation_3d_data_averaging( mode, variable )
     IF ( salsa               )  CALL salsa_3d_data_averaging( mode, variable )
Index: /palm/trunk/SOURCE/netcdf_interface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 4126)
+++ /palm/trunk/SOURCE/netcdf_interface_mod.f90	(revision 4127)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! -Introduce new vertical dimension for plant-canopy output. 
+! -Temporarlily disable masked output for soil (merge from branch resler)
+! 
+! 4069 2019-07-01 14:05:51Z Giersch
 ! Masked output running index mid has been introduced as a local variable to 
 ! avoid runtime error (Loop variable has been modified) in time_integration
@@ -538,5 +542,6 @@
                     id_dim_y_xz, id_dim_yv_xz, id_dim_y_yz, id_dim_yv_yz, &
                     id_dim_y_3d, id_dim_yv_3d, id_dim_zs_xy, id_dim_zs_xz, &
-                    id_dim_zs_yz, id_dim_zs_3d, id_dim_zu_xy, id_dim_zu1_xy, &
+                    id_dim_zs_yz, id_dim_zs_3d, id_dim_zpc_3d, &
+                    id_dim_zu_xy, id_dim_zu1_xy, &
                     id_dim_zu_xz, id_dim_zu_yz, id_dim_zu_3d, id_dim_zw_xy, &
                     id_dim_zw_xz, id_dim_zw_yz, id_dim_zw_3d, id_set_xy, &
@@ -548,6 +553,6 @@
                     id_var_yv_xy, id_var_y_xz, id_var_yv_xz, id_var_y_yz, &
                     id_var_yv_yz, id_var_y_3d, id_var_yv_3d, id_var_zs_xy, &
-                    id_var_zs_xz, id_var_zs_yz, id_var_zs_3d, id_var_zusi_xy, &
-                    id_var_zusi_3d, id_var_zu_xy, id_var_zu1_xy, id_var_zu_xz, &
+                    id_var_zs_xz, id_var_zs_yz, id_var_zs_3d, id_var_zpc_3d, &
+                    id_var_zusi_xy, id_var_zusi_3d, id_var_zu_xy, id_var_zu1_xy, id_var_zu_xz, &
                     id_var_zu_yz, id_var_zu_3d, id_var_zwwi_xy, id_var_zwwi_3d, &
                     id_var_zw_xy, id_var_zw_xz, id_var_zw_yz, id_var_zw_3d
@@ -738,5 +743,5 @@
 
     USE plant_canopy_model_mod,                                                &
-        ONLY:  pcm_define_netcdf_grid
+        ONLY:  pch_index, pcm_define_netcdf_grid
 
     USE profil_parameter,                                                      &
@@ -1552,21 +1557,22 @@
 
           ENDIF
-
-          IF ( land_surface )  THEN
+!
+!--       soil is not in masked output for now - disable temporary this block
+!          IF ( land_surface )  THEN
 !
 !--          Write zs data (vertical axes for soil model), use negative values
 !--          to indicate soil depth
-             ALLOCATE( netcdf_data(mask_size(mid,3)) )
-
-             netcdf_data = zs( mask_k_global(mid,:mask_size(mid,3)) )
-
-             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zs_mask(mid,av), &
-                                     netcdf_data, start = (/ 1 /), &
-                                     count = (/ mask_size(mid,3) /) )
-             CALL netcdf_handle_error( 'netcdf_define_header', 538 )
-
-             DEALLOCATE( netcdf_data )
-
-          ENDIF
+!             ALLOCATE( netcdf_data(mask_size(mid,3)) )
+!
+!             netcdf_data = zs( mask_k_global(mid,:mask_size(mid,3)) )
+!
+!             nc_stat = NF90_PUT_VAR( id_set_mask(mid,av), id_var_zs_mask(mid,av), &
+!                                     netcdf_data, start = (/ 1 /), &
+!                                     count = (/ mask_size(mid,3) /) )
+!             CALL netcdf_handle_error( 'netcdf_define_header', 538 )
+!
+!             DEALLOCATE( netcdf_data )
+!
+!          ENDIF
 
 !
@@ -1847,4 +1853,16 @@
           ENDIF
 
+          IF ( plant_canopy )  THEN
+!
+!--          Define vertical coordinate grid (zpc grid)
+             CALL netcdf_create_dim( id_set_3d(av), 'zpc_3d',                  &
+                                     pch_index+1, id_dim_zpc_3d(av), 70 )
+             !netcdf_create_dim(ncid, dim_name, ncdim_type, ncdim_id, error_no)
+             CALL netcdf_create_var( id_set_3d(av), (/ id_dim_zpc_3d(av) /),    &
+                                     'zpc_3d', NF90_DOUBLE, id_var_zpc_3d(av),   &
+                                     'meters', '', 71, 72, 00 )
+
+          ENDIF
+
 !
 !--       Define the variables
@@ -2006,4 +2024,6 @@
              ELSEIF ( grid_z == 'zs' )  THEN
                 id_z = id_dim_zs_3d(av)
+             ELSEIF ( grid_z == 'zpc' )  THEN
+                id_z = id_dim_zpc_3d(av)
              ENDIF
 
@@ -2248,4 +2268,13 @@
                                         - zs(nzb_soil:nzt_soil), start = (/ 1 /), &
                                         count = (/ nzt_soil-nzb_soil+1 /) )
+                CALL netcdf_handle_error( 'netcdf_define_header', 86 )
+             ENDIF
+
+             IF ( plant_canopy )  THEN
+!
+!--             Write zpc grid
+                nc_stat = NF90_PUT_VAR( id_set_3d(av), id_var_zpc_3d(av),  &
+                                        zu(nzb:nzb+pch_index), start = (/ 1 /), &
+                                        count = (/ pch_index+1 /) )
                 CALL netcdf_handle_error( 'netcdf_define_header', 86 )
              ENDIF
Index: /palm/trunk/SOURCE/plant_canopy_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 4126)
+++ /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
Index: /palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4126)
+++ /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 4127)
@@ -28,4 +28,7 @@
 ! -----------------
 ! $Id$
+! Remove unused pch_index (merge from branch resler)
+! 
+! 4089 2019-07-11 14:30:27Z suehring
 ! - Correct level 2 initialization of spectral albedos in rrtmg branch, long- and
 !   shortwave albedos were mixed-up. 
@@ -6194,10 +6197,7 @@
        USE control_parameters,                                                 &
            ONLY:  dz_stretch_level_start
-           
-       USE netcdf_data_input_mod,                                              &
-           ONLY:  leaf_area_density_f
 
        USE plant_canopy_model_mod,                                             &
-           ONLY:  pch_index, lad_s
+           ONLY:  lad_s
 
        IMPLICIT NONE
@@ -6259,9 +6259,4 @@
            nzutl = MAX( nzutl, MAXVAL( pct ) )
            nzptl = MAXVAL( pct )
-!--        code of plant canopy model uses parameter pch_index
-!--        we need to setup it here to right value
-!--        (pch_index, lad_s and other arrays in PCM are defined flat)
-           pch_index = MERGE( leaf_area_density_f%nz - 1, MAXVAL( pch ),       &
-                              leaf_area_density_f%from_file )
 
            prototype_lad = MAXVAL( lad_s ) * .9_wp  !< better be *1.0 if lad is either 0 or maxval(lad) everywhere
Index: /palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/urban_surface_mod.f90	(revision 4126)
+++ /palm/trunk/SOURCE/urban_surface_mod.f90	(revision 4127)
@@ -28,4 +28,8 @@
 ! -----------------
 ! $Id$
+! Do not add anthopogenic energy during wall/soil spin-up 
+! (merge from branch resler)
+! 
+! 4077 2019-07-09 13:27:11Z gronemeier
 ! Set roughness length z0 and z0h/q at ground-floor level to same value as
 ! those above ground-floor level
@@ -8590,5 +8594,5 @@
 !
 !--     Add-up anthropogenic heat, for now only at upward-facing surfaces
-         IF ( usm_anthropogenic_heat  .AND.  &
+         IF ( usm_anthropogenic_heat  .AND.  .NOT. during_spinup  .AND. &
               intermediate_timestep_count == intermediate_timestep_count_max )  THEN
 !
