Index: /palm/trunk/SOURCE/chemistry_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 3013)
+++ /palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 3014)
@@ -27,4 +27,7 @@
 ! -----------------
 ! $Id$
+! Bugfix: nzb_do and nzt_do were not used for 3d data output
+! 
+! 3004 2018-04-27 12:33:25Z Giersch
 ! Comment concerning averaged data output added 
 ! 
@@ -1152,5 +1155,5 @@
 
 
-   SUBROUTINE chem_data_output_3d( av, variable, found, local_pf, fill_value )
+   SUBROUTINE chem_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
 
 
@@ -1163,9 +1166,13 @@
 
       CHARACTER (LEN=*) ::  variable !<
+
+      INTEGER(iwp) ::  av    !<
+      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 !<
-      INTEGER(iwp) ::  av    !<
 
       REAL(wp) ::  fill_value !<
-      REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf 
+      REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf 
 
 
@@ -1189,5 +1196,5 @@
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
-                      DO  k = nzb, nzt+1
+                      DO  k = nzb_do, nzt_do
                           local_pf(i,j,k) = MERGE(                             &
                                               chem_species(lsp)%conc(k,j,i),   &
@@ -1201,5 +1208,5 @@
                 DO  i = nxl, nxr
                    DO  j = nys, nyn
-                      DO  k = nzb, nzt+1
+                      DO  k = nzb_do, nzt_do
                           local_pf(i,j,k) = MERGE(                             &
                                               chem_species(lsp)%conc_av(k,j,i),&
Index: /palm/trunk/SOURCE/data_output_2d.f90
===================================================================
--- /palm/trunk/SOURCE/data_output_2d.f90	(revision 3013)
+++ /palm/trunk/SOURCE/data_output_2d.f90	(revision 3014)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Added nzb_do and nzt_do for some modules for 2d output
+! 
+! 3004 2018-04-27 12:33:25Z Giersch
 ! precipitation_rate removed, case prr*_xy removed, to_be_resorted have to point
 ! to ql_vp_av and not to ql_vp, allocation checks implemented (averaged data 
@@ -1395,5 +1398,6 @@
                 IF ( .NOT. found  .AND.  radiation )  THEN
                    CALL radiation_data_output_2d( av, do2d(av,if), found, grid,&
-                                                  mode, local_pf, two_d  )
+                                                  mode, local_pf, two_d,       &
+                                                  nzb_do, nzt_do  )
                 ENDIF
 
Index: /palm/trunk/SOURCE/data_output_3d.f90
===================================================================
--- /palm/trunk/SOURCE/data_output_3d.f90	(revision 3013)
+++ /palm/trunk/SOURCE/data_output_3d.f90	(revision 3014)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Added nzb_do and nzt_do for some modules for 3d output
+! 
+! 3004 2018-04-27 12:33:25Z Giersch
 ! Allocation checks implemented (averaged data will be assigned to fill values 
 ! if no allocation happened so far) 
@@ -779,5 +782,6 @@
 
              IF ( .NOT. found )  THEN
-                CALL tcm_data_output_3d( av, do3d(av,if), found, local_pf )
+                CALL tcm_data_output_3d( av, do3d(av,if), found, local_pf,     &
+                                         nzb_do, nzt_do )
                 resorted = .TRUE.
              ENDIF
@@ -787,5 +791,5 @@
              IF ( .NOT. found  .AND.  radiation )  THEN
                 CALL radiation_data_output_3d( av, do3d(av,if), found,         &
-                                               local_pf )
+                                               local_pf, nzb_do, nzt_do )
                 resorted = .TRUE.
              ENDIF
@@ -803,5 +807,5 @@
              IF ( .NOT. found  .AND.  air_chemistry )  THEN
                 CALL chem_data_output_3d( av, do3d(av,if), found,              &
-                                          local_pf, fill_value )
+                                          local_pf, fill_value, nzb_do, nzt_do )
                 resorted = .TRUE.
              ENDIF
@@ -811,5 +815,5 @@
              IF ( .NOT. found  .AND.  plant_canopy )  THEN
                 CALL pcm_data_output_3d( av, do3d(av,if), found, local_pf,     &
-                                         fill_value )
+                                         fill_value, nzb_do, nzt_do )
                 resorted = .TRUE.
              ENDIF
Index: /palm/trunk/SOURCE/gust_mod.f90
===================================================================
--- /palm/trunk/SOURCE/gust_mod.f90	(revision 3013)
+++ /palm/trunk/SOURCE/gust_mod.f90	(revision 3014)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Bugfix: domain bounds of local_pf corrected
+! 
+! 
 ! Interfaces concerning data output updated
 ! 
@@ -408,5 +411,5 @@
        REAL(wp) ::  fill_value = -999.0_wp    !< value for the _FillValue attribute
 
-       REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf !<
+       REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
 
 
Index: /palm/trunk/SOURCE/init_3d_model.f90
===================================================================
--- /palm/trunk/SOURCE/init_3d_model.f90	(revision 3013)
+++ /palm/trunk/SOURCE/init_3d_model.f90	(revision 3014)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Bugfix: initialization of ts_value missing
+! 
+! 3011 2018-05-07 14:38:42Z schwenkel
 ! removed redundant if statement
 ! 
@@ -1127,4 +1130,8 @@
 
 !
+!-- Initialize time series
+    ts_value = 0.0_wp
+
+!
 !-- Initialize local summation arrays for routine flow_statistics.
 !-- This is necessary because they may not yet have been initialized when they
Index: /palm/trunk/SOURCE/land_surface_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 3013)
+++ /palm/trunk/SOURCE/land_surface_model_mod.f90	(revision 3014)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Bugfix: set some initial values
+! Bugfix: domain bounds of local_pf corrected
+! 
+! 3004 2018-04-27 12:33:25Z Giersch
 ! Further allocation checks implemented (averaged data will be assigned to fill
 ! values if no allocation happened so far) 
@@ -4639,4 +4643,9 @@
        surf_lsm_h%pavement_surface     = .FALSE.
        surf_lsm_h%vegetation_surface   = .FALSE. 
+
+!
+!--    Set default values
+       surf_lsm_h%r_canopy_min = 0.0_wp
+
 !
 !--    Vertical surfaces
@@ -4671,4 +4680,9 @@
           surf_lsm_v(l)%vegetation_surface   = .FALSE. 
           
+
+!
+!--       Set default values
+          surf_lsm_v(l)%r_canopy_min = 0.0_wp
+        
        ENDDO
 
@@ -5632,5 +5646,5 @@
     REAL(wp) ::  fill_value = -999.0_wp    !< value for the _FillValue attribute
 
-    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf !< 
+    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< 
 
 
Index: /palm/trunk/SOURCE/modules.f90
===================================================================
--- /palm/trunk/SOURCE/modules.f90	(revision 3013)
+++ /palm/trunk/SOURCE/modules.f90	(revision 3014)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Added default values of u_max, v_max, and w_max to avoid floating invalid
+! during spinup
+! 
+! 3004 2018-04-27 12:33:25Z Giersch
 ! precipitation_rate removed
 ! 
@@ -2043,7 +2047,7 @@
                                                   !< (after each timestep)
     
-    REAL(wp) ::  u_max  !< maximum of absolute u-veloctiy in entire domain
-    REAL(wp) ::  v_max  !< maximum of absolute v-veloctiy in entire domain
-    REAL(wp) ::  w_max  !< maximum of absolute w-veloctiy in entire domain
+    REAL(wp) ::  u_max = 0.0_wp  !< maximum of absolute u-veloctiy in entire domain
+    REAL(wp) ::  v_max = 0.0_wp  !< maximum of absolute v-veloctiy in entire domain
+    REAL(wp) ::  w_max = 0.0_wp  !< maximum of absolute w-veloctiy in entire domain
 
     REAL(wp), DIMENSION(2) ::  z_i  !< inversion height
Index: /palm/trunk/SOURCE/plant_canopy_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 3013)
+++ /palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 3014)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Bugfix: nzb_do and nzt_do were not used for 3d data output
+! Added pc_transpiration_rate
+! 
+! 2977 2018-04-17 10:27:57Z kanani
 ! Implement changes from branch radiation (r2948-2971) with minor modifications,
 ! plus some formatting.
@@ -230,4 +234,5 @@
     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
 
     SAVE
@@ -243,5 +248,5 @@
 !
 !-- Public variables and constants
-    PUBLIC pc_heating_rate, canopy_mode, cthf, dt_plant_canopy, lad, lad_s,   &
+    PUBLIC pc_heating_rate, pc_transpiration_rate, canopy_mode, cthf, dt_plant_canopy, lad, lad_s,   &
            pch_index
            
@@ -315,4 +320,7 @@
           unit = 'K s-1'
     
+       CASE ( 'pcm_transpirationrate' )
+          unit = 'kg kg-1 s-1'
+
        CASE ( 'pcm_lad' )
           unit = 'm2 m-3'
@@ -404,9 +412,7 @@
 !> Subroutine defining 3D output variables
 !------------------------------------------------------------------------------!
- SUBROUTINE pcm_data_output_3d( av, variable, found, local_pf, fill_value )
- 
-    USE control_parameters,                                                    &
-        ONLY :  nz_do3d
- 
+ SUBROUTINE pcm_data_output_3d( av, variable, found, local_pf, fill_value,     &
+                                nzb_do, nzt_do )
+  
     USE indices
 
@@ -423,9 +429,11 @@
     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)
 
     LOGICAL      ::  found !< 
 
     REAL(wp)     ::  fill_value
-    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb:nz_do3d) ::  local_pf !< 
+    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< 
 
 
@@ -449,5 +457,18 @@
             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_lad' )
@@ -497,5 +518,5 @@
      SELECT CASE ( TRIM( var ) )
 
-        CASE ( 'pcm_heatrate', 'pcm_lad' )
+        CASE ( 'pcm_heatrate', 'pcm_lad', 'pcm_transpirationrate')
            grid_x = 'x'
            grid_y = 'y'
@@ -873,5 +894,6 @@
 
           ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                 &
-                    pc_heating_rate(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+                    pc_heating_rate(nzb:nzt+1,nysg:nyng,nxlg:nxrg),            & 
+                    pc_transpiration_rate(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 !
 !--       Piecewise calculation of the cumulative leaf area index by vertical
@@ -1397,7 +1419,6 @@
 
                       kk = k - k_wall   !- lad arrays are defined flat
-                      tend(k,j,i) = tend(k,j,i) -                              &
-                                       lsec *                                  &
-                                       lad_s(kk,j,i) *                         &
+                      pc_transpiration_rate(kk,j,i) =  - lsec                  &
+                                       * lad_s(kk,j,i) *                       &
                                        SQRT( ( 0.5_wp * ( u(k,j,i) +           &
                                                           u(k,j,i+1) )         &
@@ -1411,4 +1432,6 @@
                                            ) *                                 &
                                        ( q(k,j,i) - lsc )
+
+                      tend(k,j,i) = tend(k,j,i) + pc_transpiration_rate(kk,j,i)
                    ENDDO
                 ENDDO
@@ -1728,9 +1751,8 @@
 
              DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
-
-                kk = k - k_wall
-                tend(k,j,i) = tend(k,j,i) -                                    &
-                                 lsec *                                        &
-                                 lad_s(kk,j,i) *                               &
+                kk = k - k_wall  !- lad arrays are defined flat
+
+                pc_transpiration_rate(kk,j,i) = - lsec                         &
+                                 * lad_s(kk,j,i) *                             &
                                  SQRT( ( 0.5_wp * ( u(k,j,i) +                 &
                                                     u(k,j,i+1) )               &
@@ -1744,4 +1766,7 @@
                                      ) *                                       &
                                  ( q(k,j,i) - lsc )
+
+                tend(k,j,i) = tend(k,j,i) + pc_transpiration_rate(kk,j,i)
+
              ENDDO   
 
Index: /palm/trunk/SOURCE/prognostic_equations.f90
===================================================================
--- /palm/trunk/SOURCE/prognostic_equations.f90	(revision 3013)
+++ /palm/trunk/SOURCE/prognostic_equations.f90	(revision 3014)
@@ -25,4 +25,8 @@
 ! -----------------
 ! $Id$
+! Fixed a bug in the IF condition to call pcm_tendency in case of 
+! potential temperature
+! 
+! 2815 2018-02-19 11:29:57Z kanani
 ! Rename chem_tendency to chem_prognostic_equations,
 ! implement vector version for air chemistry
@@ -298,5 +302,5 @@
                timestep_scheme, tsc, use_subsidence_tendencies,                &
                use_upstream_for_tke, wind_turbine, ws_scheme_mom,              & 
-               ws_scheme_sca
+               ws_scheme_sca, urban_surface, land_surface
 
     USE cpulog,                                                                &
@@ -759,5 +763,6 @@
 !
 !--          Consideration of heat sources within the plant canopy
-             IF ( plant_canopy  .AND.  cthf /= 0.0_wp )  THEN
+             IF ( plant_canopy  .AND.                                          &
+                (cthf /= 0.0_wp  .OR. urban_surface  .OR.  land_surface) )  THEN
                 CALL pcm_tendency( i, j, 4 )
              ENDIF
@@ -1325,6 +1330,4 @@
 
 
-
-
     CALL cpu_log( log_point(32), 'all progn.equations', 'stop' )
 
@@ -1671,5 +1674,6 @@
 !
 !--    Consideration of heat sources within the plant canopy
-       IF ( plant_canopy .AND. ( cthf /= 0.0_wp ) ) THEN
+       IF ( plant_canopy  .AND.                                          &
+            (cthf /= 0.0_wp  .OR. urban_surface  .OR.  land_surface) )  THEN
           CALL pcm_tendency( 4 )
        ENDIF
@@ -1718,5 +1722,4 @@
           ENDDO
        ENDDO
-
 !
 !--    Calculate tendencies for the next Runge-Kutta step
Index: /palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 3013)
+++ /palm/trunk/SOURCE/radiation_model_mod.f90	(revision 3014)
@@ -28,4 +28,9 @@
 ! -----------------
 ! $Id$
+! Introduced plant canopy height similar to urban canopy height to limit
+! the memory requirement to allocate lad.
+! Deactivated automatic setting of minimum raytracing distance.
+! 
+! 3004 2018-04-27 12:33:25Z Giersch
 ! Further allocation checks implemented (averaged data will be assigned to fill
 ! values if no allocation happened so far) 
@@ -329,5 +334,5 @@
 
     USE cloud_parameters,                                                      &
-        ONLY:  cp, l_d_cp, r_d, rho_l
+        ONLY:  cp, l_d_cp, l_v, r_d, rho_l
 
     USE constants,                                                             &
@@ -372,5 +377,5 @@
 
     USE plant_canopy_model_mod,                                                &
-        ONLY:  pc_heating_rate, lad_s
+        ONLY:  lad_s, pc_heating_rate, pc_transpiration_rate
 
     USE pegrid
@@ -646,5 +651,7 @@
 !-- Parameters of urban and land surface models
     INTEGER(iwp)                                   ::  nzu                                !< number of layers of urban surface (will be calculated)
+    INTEGER(iwp)                                   ::  nzp                                !< number of layers of plant canopy (will be calculated)
     INTEGER(iwp)                                   ::  nzub,nzut                          !< bottom and top layer of urban surface (will be calculated)
+    INTEGER(iwp)                                   ::  nzpt                               !< top layer of plant canopy (will be calculated)
 !-- parameters of urban and land surface models
     INTEGER(iwp), PARAMETER                        ::  nzut_free = 3                      !< number of free layers above top of of topography
@@ -4324,5 +4331,6 @@
      REAL(wp), DIMENSION(3)            :: sunorig_grid       !< grid squashed solar direction unit vector (zyx)
      REAL(wp), DIMENSION(0:nsurf_type) :: costheta           !< direct irradiance factor of solar angle
-     REAL(wp), DIMENSION(nzub:nzut)    :: pchf_prep          !< precalculated factor for canopy temp tendency
+     REAL(wp), DIMENSION(nzub:nzut)    :: pchf_prep          !< precalculated factor for canopy temperature tendency
+     REAL(wp), DIMENSION(nzub:nzut)    :: pctf_prep          !< precalculated factor for canopy transpiration tendency
      REAL(wp), PARAMETER               :: alpha = 0._wp      !< grid rotation (TODO: add to namelist or remove)
      REAL(wp)                          :: pc_box_area, pc_abs_frac, pc_abs_eff
@@ -4345,8 +4353,12 @@
      REAL(wp)                          :: area_surf          !< total area of surfaces in all processor
 
+
+
 #if ! defined( __nopointer )
      IF ( plant_canopy )  THEN
          pchf_prep(:) = r_d * (hyp(nzub:nzut) / 100000.0_wp)**0.286_wp &
                      / (cp * hyp(nzub:nzut) * dx*dy*dz) !< equals to 1 / (rho * c_p * Vbox * T)
+         pctf_prep(:) = r_d * (hyp(nzub:nzut) / 100000.0_wp)**0.286_wp &
+                     / (l_v * hyp(nzub:nzut) * dx*dy*dz)
      ENDIF
 #endif
@@ -4619,5 +4631,6 @@
 !--  push heat flux absorbed by plant canopy to respective 3D arrays
      IF ( npcbl > 0 )  THEN
-         pc_heating_rate(:,:,:) = 0._wp
+         pc_heating_rate(:,:,:) = 0.0_wp
+         pc_transpiration_rate(:,:,:) = 0.0_wp
          DO ipcgb = 1, npcbl
                  
@@ -4630,4 +4643,8 @@
              pc_heating_rate(kk, j, i) = (pcbinsw(ipcgb) + pcbinlw(ipcgb)) &
                  * pchf_prep(k) * pt(k, j, i) !-- = dT/dt
+
+!             pc_transpiration_rate(kk,j,i) = 0.75_wp* (pcbinsw(ipcgb) + pcbinlw(ipcgb)) &
+!                 * pctf_prep(k) * pt(k, j, i) !-- = dq/dt
+
          ENDDO
      ENDIF
@@ -4995,5 +5012,5 @@
        INTEGER(iwp) :: k_topo     !< vertical index indicating topography top for given (j,i)
        INTEGER(iwp) :: k_topo2    !< vertical index indicating topography top for given (j,i)
-       INTEGER(iwp) :: nzubl, nzutl, isurf, ipcgb
+       INTEGER(iwp) :: nzptl, nzubl, nzutl, isurf, ipcgb
        INTEGER(iwp) :: procid
        REAL(wp)     :: mrl
@@ -5040,4 +5057,5 @@
 
            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
@@ -5058,23 +5076,26 @@
        CALL MPI_AllReduce(nzubl, nzub, 1, MPI_INTEGER, MPI_MIN, comm2d, ierr )
        CALL MPI_AllReduce(nzutl, nzut, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )
+       CALL MPI_AllReduce(nzptl, nzpt, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )
 #else
        nzub = nzubl
        nzut = nzutl
+       nzpt = nzptl
 #endif
 !
-!--    global number of urban layers
+!--    global number of urban and plant layers
        nzu = nzut - nzub + 1
+       nzp = nzpt - nzub + 1
 !
 !--    check max_raytracing_dist relative to urban surface layer height
-       mrl = 2.0_wp * nzu * dz
-       IF ( max_raytracing_dist <= mrl ) THEN
-          IF ( max_raytracing_dist /= -999.0_wp ) THEN
-!--          max_raytracing_dist too low
-             WRITE(message_string, '(a,f6.1)') 'Max_raytracing_dist too low, ' &
-                   // 'override to value ', mrl
-             CALL message('radiation_interaction_init', 'PA0521', 0, 0, -1, 6, 0)
-          ENDIF
-          max_raytracing_dist = mrl
-       ENDIF
+!        mrl = 2.0_wp * nzu * dz
+!        IF ( max_raytracing_dist <= mrl ) THEN
+!           IF ( max_raytracing_dist /= -999.0_wp ) THEN
+! !--          max_raytracing_dist too low
+!              WRITE(message_string, '(a,f6.1)') 'Max_raytracing_dist too low, ' &
+!                    // 'override to value ', mrl
+!              CALL message('radiation_interaction_init', 'PA0521', 0, 0, -1, 6, 0)
+!           ENDIF
+!           max_raytracing_dist = mrl
+!        ENDIF
 !
 !--    allocate urban surfaces grid
@@ -5103,5 +5124,5 @@
        IF ( npcbl > 0 )  THEN
            ALLOCATE( pcbl(iz:ix, 1:npcbl) )
-           ALLOCATE( gridpcbl(nzub:nzut,nys:nyn,nxl:nxr) )
+           ALLOCATE( gridpcbl(nzub:nzpt,nys:nyn,nxl:nxr) )
            pcbl = -1
            gridpcbl(:,:,:) = 0
@@ -5388,5 +5409,5 @@
         IF ( plant_canopy )  THEN
             ALLOCATE( plantt(0:(nx+1)*(ny+1)-1) )
-            maxboxesg = nx + ny + nzu + 1
+            maxboxesg = nx + ny + nzp + 1
             max_track_len = nx + ny + 1
 !--         temporary arrays storing values for csf calculation during raytracing
@@ -5416,15 +5437,15 @@
 !--             optimization of memory should be done
 !--             Argument X of function STORAGE_SIZE(X) needs arbitrary REAL(wp) value, set to 1.0_wp for now
-                size_lad_rma = STORAGE_SIZE(1.0_wp)/8*nnx*nny*nzu
+                size_lad_rma = STORAGE_SIZE(1.0_wp)/8*nnx*nny*nzp
                 CALL MPI_Win_allocate(size_lad_rma, STORAGE_SIZE(1.0_wp)/8, minfo, comm2d, &
                                         lad_s_rma_p, win_lad, ierr)
-                CALL c_f_pointer(lad_s_rma_p, lad_s_rma, (/ nzu, nny, nnx /))
+                CALL c_f_pointer(lad_s_rma_p, lad_s_rma, (/ nzp, nny, nnx /))
                 sub_lad(nzub:, nys:, nxl:) => lad_s_rma(:,:,:)
             ELSE
-                ALLOCATE(sub_lad(nzub:nzut, nys:nyn, nxl:nxr))
+                ALLOCATE(sub_lad(nzub:nzpt, nys:nyn, nxl:nxr))
             ENDIF
 #else
             plantt = RESHAPE( pct(nys:nyn,nxl:nxr), (/(nx+1)*(ny+1)/) )
-            ALLOCATE(sub_lad(nzub:nzut, nys:nyn, nxl:nxr))
+            ALLOCATE(sub_lad(nzub:nzpt, nys:nyn, nxl:nxr))
 #endif
             plantt_max = MAXVAL(plantt)
@@ -5437,5 +5458,5 @@
                     k = get_topography_top_index_ji( j, i, 's' )
 
-                    sub_lad(k:nzut, j, i) = lad_s(0:nzut-k, j, i)
+                    sub_lad(k:nzpt, j, i) = lad_s(0:nzpt-k, j, i)
                 ENDDO
             ENDDO
@@ -5446,7 +5467,7 @@
                 CALL MPI_Win_lock_all(0, win_lad, ierr)
             ELSE
-                ALLOCATE( sub_lad_g(0:(nx+1)*(ny+1)*nzu-1) )
-                CALL MPI_AllGather( sub_lad, nnx*nny*nzu, MPI_REAL, &
-                                    sub_lad_g, nnx*nny*nzu, MPI_REAL, comm2d, ierr )
+                ALLOCATE( sub_lad_g(0:(nx+1)*(ny+1)*nzp-1) )
+                CALL MPI_AllGather( sub_lad, nnx*nny*nzp, MPI_REAL, &
+                                    sub_lad_g, nnx*nny*nzp, MPI_REAL, comm2d, ierr )
             ENDIF
 #endif
@@ -6142,5 +6163,5 @@
 #if defined( __parallel )
                         lad_ip(ncsb) = ip
-                        lad_disp(ncsb) = (box(3)-px*nnx)*(nny*nzu) + (box(2)-py*nny)*nzu + box(1)-nzub
+                        lad_disp(ncsb) = (box(3)-px*nnx)*(nny*nzp) + (box(2)-py*nny)*nzp + box(1)-nzub
 #endif
                     ENDIF
@@ -6185,5 +6206,5 @@
                     lad_s_target = lad_s_ray(i)
                 ELSE
-                    lad_s_target = sub_lad_g(lad_ip(i)*nnx*nny*nzu + lad_disp(i))
+                    lad_s_target = sub_lad_g(lad_ip(i)*nnx*nny*nzp + lad_disp(i))
                 ENDIF
 #else
@@ -6387,5 +6408,5 @@
                ig = ip*nnx*nny + (rt2_track(2,i)-px*nnx)*nny + rt2_track(1,i)-py*nny
                IF ( plantt(ig) <= nzterr(ig) )  CYCLE
-               wdisp = (rt2_track(2,i)-px*nnx)*(nny*nzu) + (rt2_track(1,i)-py*nny)*nzu + nzterr(ig)+1-nzub
+               wdisp = (rt2_track(2,i)-px*nnx)*(nny*nzp) + (rt2_track(1,i)-py*nny)*nzp + nzterr(ig)+1-nzub
                wcount = plantt(ig)-nzterr(ig)
                ! TODO send request ASAP - even during raytracing
@@ -6411,5 +6432,5 @@
                py = rt2_track(1,i)/nny
                ip = px*pdims(2)+py
-               ig = ip*nnx*nny*nzu + (rt2_track(2,i)-px*nnx)*(nny*nzu) + (rt2_track(1,i)-py*nny)*nzu
+               ig = ip*nnx*nny*nzp + (rt2_track(2,i)-px*nnx)*(nny*nzp) + (rt2_track(1,i)-py*nny)*nzp
                rt2_track_lad(nzub:plantt_max, i) = sub_lad_g(ig:ig+nly-1)
             ENDDO
@@ -6429,5 +6450,5 @@
          !--Assert that we have space allocated for CSFs
          !--
-         maxboxes = (ntrack + MAX(origin(1) - nzub, nzut - origin(1))) * SIZE(zdirs, 1)
+         maxboxes = (ntrack + MAX(origin(1) - nzub, nzpt - origin(1))) * SIZE(zdirs, 1)
          IF ( ncsfl + maxboxes > ncsfla )  THEN
 !--         use this code for growing by fixed exponential increments (equivalent to case where ncsfl always increases by 1)
@@ -6814,5 +6835,5 @@
               ENDIF       
 !
-!--           Close binary file                 
+!--           Close binary file               
               CALL close_file( fsvf )
                 
@@ -7617,5 +7638,5 @@
 !------------------------------------------------------------------------------!
  SUBROUTINE radiation_data_output_2d( av, variable, found, grid, mode,         &
-                                      local_pf, two_d )
+                                      local_pf, two_d, nzb_do, nzt_do )
  
     USE indices
@@ -7635,4 +7656,6 @@
     INTEGER(iwp) ::  k  !< 
     INTEGER(iwp) ::  m  !< index of surface element at grid point (j,i)
+    INTEGER(iwp) ::  nzb_do   !< 
+    INTEGER(iwp) ::  nzt_do   !< 
 
     LOGICAL      ::  found !< 
@@ -7641,5 +7664,5 @@
     REAL(wp) ::  fill_value = -999.0_wp    !< value for the _FillValue attribute
 
-    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf !< 
+    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< 
 
     found = .TRUE.
@@ -7685,5 +7708,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = rad_lw_in(k,j,i)
                    ENDDO
@@ -7697,5 +7720,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn 
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = rad_lw_in_av(k,j,i)
                    ENDDO
@@ -7709,5 +7732,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = rad_lw_out(k,j,i)
                    ENDDO
@@ -7721,5 +7744,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn 
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = rad_lw_out_av(k,j,i)
                    ENDDO
@@ -7733,5 +7756,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = rad_lw_cs_hr(k,j,i)
                    ENDDO
@@ -7745,5 +7768,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn 
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = rad_lw_cs_hr_av(k,j,i)
                    ENDDO
@@ -7757,5 +7780,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = rad_lw_hr(k,j,i)
                    ENDDO
@@ -7769,5 +7792,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn 
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = rad_lw_hr_av(k,j,i)
                    ENDDO
@@ -7781,5 +7804,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = rad_sw_in(k,j,i)
                    ENDDO
@@ -7793,5 +7816,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn 
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = rad_sw_in_av(k,j,i)
                    ENDDO
@@ -7805,5 +7828,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = rad_sw_out(k,j,i)
                    ENDDO
@@ -7829,5 +7852,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = rad_sw_cs_hr(k,j,i)
                    ENDDO
@@ -7841,5 +7864,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn 
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = rad_sw_cs_hr_av(k,j,i)
                    ENDDO
@@ -7853,5 +7876,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = rad_sw_hr(k,j,i)
                    ENDDO
@@ -7865,5 +7888,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn 
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = rad_sw_hr_av(k,j,i)
                    ENDDO
@@ -7888,5 +7911,5 @@
 !> Subroutine defining 3D output variables
 !------------------------------------------------------------------------------!
- SUBROUTINE radiation_data_output_3d( av, variable, found, local_pf )
+ SUBROUTINE radiation_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
  
 
@@ -7904,4 +7927,6 @@
     INTEGER(iwp) ::  j     !< 
     INTEGER(iwp) ::  k     !< 
+    INTEGER(iwp) ::  nzb_do   !< 
+    INTEGER(iwp) ::  nzt_do   !< 
 
     LOGICAL      ::  found !< 
@@ -7909,5 +7934,5 @@
     REAL(wp) ::  fill_value = -999.0_wp    !< value for the _FillValue attribute
 
-    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf !< 
+    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< 
 
 
@@ -7921,5 +7946,5 @@
             DO  i = nxl, nxr
                DO  j = nys, nyn
-                  DO  k = nzb, nzt+1
+                  DO  k = nzb_do, nzt_do
                      local_pf(i,j,k) = rad_sw_in(k,j,i)
                   ENDDO
@@ -7933,5 +7958,5 @@
             DO  i = nxl, nxr
                DO  j = nys, nyn
-                  DO  k = nzb, nzt+1
+                  DO  k = nzb_do, nzt_do
                      local_pf(i,j,k) = rad_sw_in_av(k,j,i)
                   ENDDO
@@ -7944,5 +7969,5 @@
             DO  i = nxl, nxr
                DO  j = nys, nyn
-                  DO  k = nzb, nzt+1
+                  DO  k = nzb_do, nzt_do
                      local_pf(i,j,k) = rad_sw_out(k,j,i)
                   ENDDO
@@ -7956,5 +7981,5 @@
             DO  i = nxl, nxr
                DO  j = nys, nyn
-                  DO  k = nzb, nzt+1
+                  DO  k = nzb_do, nzt_do
                      local_pf(i,j,k) = rad_sw_out_av(k,j,i)
                   ENDDO
@@ -7967,5 +7992,5 @@
             DO  i = nxl, nxr
                DO  j = nys, nyn
-                  DO  k = nzb, nzt+1
+                  DO  k = nzb_do, nzt_do
                      local_pf(i,j,k) = rad_sw_cs_hr(k,j,i)
                   ENDDO
@@ -7979,5 +8004,5 @@
             DO  i = nxl, nxr
                DO  j = nys, nyn
-                  DO  k = nzb, nzt+1
+                  DO  k = nzb_do, nzt_do
                      local_pf(i,j,k) = rad_sw_cs_hr_av(k,j,i)
                   ENDDO
@@ -7990,5 +8015,5 @@
             DO  i = nxl, nxr
                DO  j = nys, nyn
-                  DO  k = nzb, nzt+1
+                  DO  k = nzb_do, nzt_do
                      local_pf(i,j,k) = rad_sw_hr(k,j,i)
                   ENDDO
@@ -8002,5 +8027,5 @@
             DO  i = nxl, nxr
                DO  j = nys, nyn
-                  DO  k = nzb, nzt+1
+                  DO  k = nzb_do, nzt_do
                      local_pf(i,j,k) = rad_sw_hr_av(k,j,i)
                   ENDDO
@@ -8013,5 +8038,5 @@
             DO  i = nxl, nxr
                DO  j = nys, nyn
-                  DO  k = nzb, nzt+1
+                  DO  k = nzb_do, nzt_do
                      local_pf(i,j,k) = rad_lw_in(k,j,i)
                   ENDDO
@@ -8025,5 +8050,5 @@
             DO  i = nxl, nxr
                DO  j = nys, nyn
-                  DO  k = nzb, nzt+1
+                  DO  k = nzb_do, nzt_do
                      local_pf(i,j,k) = rad_lw_in_av(k,j,i)
                   ENDDO
@@ -8036,5 +8061,5 @@
             DO  i = nxl, nxr
                DO  j = nys, nyn
-                  DO  k = nzb, nzt+1
+                  DO  k = nzb_do, nzt_do
                      local_pf(i,j,k) = rad_lw_out(k,j,i)
                   ENDDO
@@ -8048,5 +8073,5 @@
             DO  i = nxl, nxr
                DO  j = nys, nyn
-                  DO  k = nzb, nzt+1
+                  DO  k = nzb_do, nzt_do
                      local_pf(i,j,k) = rad_lw_out_av(k,j,i)
                   ENDDO
@@ -8059,5 +8084,5 @@
             DO  i = nxl, nxr
                DO  j = nys, nyn
-                  DO  k = nzb, nzt+1
+                  DO  k = nzb_do, nzt_do
                      local_pf(i,j,k) = rad_lw_cs_hr(k,j,i)
                   ENDDO
@@ -8071,5 +8096,5 @@
             DO  i = nxl, nxr
                DO  j = nys, nyn
-                  DO  k = nzb, nzt+1
+                  DO  k = nzb_do, nzt_do
                      local_pf(i,j,k) = rad_lw_cs_hr_av(k,j,i)
                   ENDDO
@@ -8082,5 +8107,5 @@
             DO  i = nxl, nxr
                DO  j = nys, nyn
-                  DO  k = nzb, nzt+1
+                  DO  k = nzb_do, nzt_do
                      local_pf(i,j,k) = rad_lw_hr(k,j,i)
                   ENDDO
@@ -8094,5 +8119,5 @@
             DO  i = nxl, nxr
                DO  j = nys, nyn
-                  DO  k = nzb, nzt+1
+                  DO  k = nzb_do, nzt_do
                      local_pf(i,j,k) = rad_lw_hr_av(k,j,i)
                   ENDDO
Index: /palm/trunk/SOURCE/time_integration.f90
===================================================================
--- /palm/trunk/SOURCE/time_integration.f90	(revision 3013)
+++ /palm/trunk/SOURCE/time_integration.f90	(revision 3014)
@@ -25,4 +25,9 @@
 ! -----------------
 ! $Id$
+! Fixed bug in IF statement
+! Ensure that the time when calling the radiation to be the time step of the
+! pre-calculated time when first calculate the positions of the sun
+! 
+! 3004 2018-04-27 12:33:25Z Giersch
 ! First call of flow_statistics has been removed. It is already called in 
 ! run_control itself
@@ -491,9 +496,12 @@
 
     CHARACTER (LEN=9) ::  time_to_string          !< 
-    INTEGER           ::  n
-    INTEGER           ::  lsp
+    INTEGER(iwp)      ::  it
+    INTEGER(iwp)      ::  lsp
+    INTEGER(iwp)      ::  n
+
 
     REAL(wp) ::  dt_3d_old  !< temporary storage of timestep to be used for
                             !< steering of run control output interval
+    REAL(wp) ::  tsrp_org   !< original value of time_since_reference_point
 !
 !-- At beginning determine the first time step
@@ -1040,4 +1048,12 @@
                 ENDIF
 
+!
+!--             Adjust the current_ time to the time step of the radiation model.
+!--             Needed since radiation is pre-calculated and stored only on apparent 
+!--             solar positions
+                it = FLOOR(time_since_reference_point/dt_radiation)
+                tsrp_org = time_since_reference_point
+                time_since_reference_point = REAL(it,wp) * dt_radiation
+
                 CALL radiation_control
 
@@ -1050,4 +1066,8 @@
                    CALL cpu_log( log_point(75), 'radiation_interaction', 'stop' )
                 ENDIF
+    
+!
+!--             Return the current time to its original value
+                time_since_reference_point = tsrp_org
 
              ENDIF
Index: /palm/trunk/SOURCE/turbulence_closure_mod.f90
===================================================================
--- /palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3013)
+++ /palm/trunk/SOURCE/turbulence_closure_mod.f90	(revision 3014)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Bugfix: nzb_do and nzt_do were not used for 3d data output
+! 
+! 3004 2018-04-27 12:33:25Z Giersch
 ! Further allocation checks implemented
 ! 
@@ -626,5 +629,5 @@
     REAL(wp) ::  fill_value = -999.0_wp  !< value for the _FillValue attribute
 
-    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf !< local
+    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< local
        !< array to which output data is resorted to
 
@@ -723,5 +726,5 @@
 !> Define 3D output variables.
 !------------------------------------------------------------------------------!
- SUBROUTINE tcm_data_output_3d( av, variable, found, local_pf )
+ SUBROUTINE tcm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
  
 
@@ -733,8 +736,10 @@
     CHARACTER (LEN=*) ::  variable   !< 
 
-    INTEGER(iwp) ::  av   !< 
-    INTEGER(iwp) ::  i    !< 
-    INTEGER(iwp) ::  j    !< 
-    INTEGER(iwp) ::  k    !< 
+    INTEGER(iwp) ::  av     !< 
+    INTEGER(iwp) ::  i      !< 
+    INTEGER(iwp) ::  j      !< 
+    INTEGER(iwp) ::  k      !< 
+    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   !< 
@@ -742,5 +747,5 @@
     REAL(wp) ::  fill_value = -999.0_wp  !< value for the _FillValue attribute
 
-    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf   !< local
+    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf   !< local
        !< array to which output data is resorted to
 
@@ -756,5 +761,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = diss(k,j,i)
                    ENDDO
@@ -768,5 +773,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = diss_av(k,j,i)
                    ENDDO
@@ -779,5 +784,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = kh(k,j,i)
                    ENDDO
@@ -791,5 +796,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = kh_av(k,j,i)
                    ENDDO
@@ -802,5 +807,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = km(k,j,i)
                    ENDDO
@@ -814,5 +819,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = km_av(k,j,i)
                    ENDDO
@@ -825,5 +830,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = dummy1(k,j,i)
                    ENDDO
@@ -836,5 +841,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = dummy2(k,j,i)
                    ENDDO
@@ -847,5 +852,5 @@
              DO  i = nxl, nxr
                 DO  j = nys, nyn
-                   DO  k = nzb, nzt+1
+                   DO  k = nzb_do, nzt_do
                       local_pf(i,j,k) = dummy3(k,j,i)
                    ENDDO
Index: /palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- /palm/trunk/SOURCE/urban_surface_mod.f90	(revision 3013)
+++ /palm/trunk/SOURCE/urban_surface_mod.f90	(revision 3014)
@@ -28,4 +28,7 @@
 ! -----------------
 ! $Id$
+! Added pc_transpiration_rate
+! 
+! 2977 2018-04-17 10:27:57Z kanani
 ! Implement changes from branch radiation (r2948-2971) with minor modifications.
 ! (moh.hefny):
@@ -264,4 +267,5 @@
 !>       by zero, e.g. in case fraq(0,m) + fraq(1,m) = 0?!
 !> @todo Use unit 90 for OPEN/CLOSE of input files (FK)
+!> @todo Move plant canopy stuff into plant canopy code
 !------------------------------------------------------------------------------!
  MODULE urban_surface_mod
@@ -309,5 +313,5 @@
     
     USE plant_canopy_model_mod,                                                &
-        ONLY:  pc_heating_rate
+        ONLY:  pc_heating_rate, pc_transpiration_rate
     
     USE radiation_model_mod,                                                   &
@@ -2361,5 +2365,5 @@
         INTEGER(iwp)                                           :: ids,idsint,idsidx,isurf,isvf,isurfs,isurflt
         INTEGER(iwp)                                           :: is,js,ks,i,j,k,iwl,istat, l, m
-        INTEGER(iwp)                                           ::  k_topo    !< topography top index
+        INTEGER(iwp)                                           :: k_topo    !< topography top index
 
         dirstart = (/ startland, startwall, startwall, startwall, startwall /)
@@ -4539,4 +4543,17 @@
 !--             in case of cthf /= 0 => we need to allocate it for our use here
                 ALLOCATE( pc_heating_rate(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+
+                pc_heating_rate = 0.0_wp
+
+            ENDIF
+
+            IF ( .NOT.  ALLOCATED( pc_transpiration_rate) )  THEN
+!--             then pc_heating_rate is allocated in init_plant_canopy
+!--             in case of cthf /= 0 => we need to allocate it for our use here
+                ALLOCATE( pc_transpiration_rate(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+
+                pc_transpiration_rate = 0.0_wp
+
+
             ENDIF
         ENDIF
Index: /palm/trunk/SOURCE/user_data_output_2d.f90
===================================================================
--- /palm/trunk/SOURCE/user_data_output_2d.f90	(revision 3013)
+++ /palm/trunk/SOURCE/user_data_output_2d.f90	(revision 3014)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Bugfix: domain bounds of local_pf corrected
+! 
+! 3004 2018-04-27 12:33:25Z Giersch
 ! Further allocation checks implemented (averaged data will be assigned to fill
 ! values if no allocation happened so far)
@@ -100,5 +103,5 @@
     REAL(wp) ::  fill_value = -999.0_wp    !< value for the _FillValue attribute
 
-    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf !< 
+    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< 
 
 
Index: /palm/trunk/SOURCE/uv_exposure_model_mod.f90
===================================================================
--- /palm/trunk/SOURCE/uv_exposure_model_mod.f90	(revision 3013)
+++ /palm/trunk/SOURCE/uv_exposure_model_mod.f90	(revision 3014)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! Bugfix: domain bounds of local_pf corrected
+! 
+! 3004 2018-04-27 12:33:25Z Giersch
 ! Further allocation checks implemented (averaged data will be assigned to fill
 ! values if no allocation happened so far) 
@@ -337,5 +340,5 @@
     REAL(wp) ::  fill_value = -999.0_wp    !< value for the _FillValue attribute
 
-    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb:nzt+1) ::  local_pf !< 
+    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< 
 
 
