Index: palm/trunk/SOURCE/basic_constants_and_equations_mod.f90
===================================================================
--- palm/trunk/SOURCE/basic_constants_and_equations_mod.f90	(revision 3448)
+++ palm/trunk/SOURCE/basic_constants_and_equations_mod.f90	(revision 3449)
@@ -25,4 +25,7 @@
 ! -----------------
 ! $Id$
+! +degc_to_k 
+! 
+! 3361 2018-10-16 20:39:37Z knoop
 ! New module (introduced with modularization of bulk cloud physics model)
 !
@@ -43,7 +46,6 @@
     IMPLICIT NONE
 
-    REAL(wp), PARAMETER ::  pi = 3.141592654_wp  !< PI
-
     REAL(wp), PARAMETER ::  c_p = 1005.0_wp                           !< heat capacity of dry air (J kg-1 K-1)
+    REAL(wp), PARAMETER ::  degc_to_k = 273.15_wp                     !< temperature (in K) of 0 deg C (K)
     REAL(wp), PARAMETER ::  g = 9.81_wp                               !< gravitational acceleration (m s-2)
     REAL(wp), PARAMETER ::  kappa = 0.4_wp                            !< von Karman constant
@@ -55,4 +57,5 @@
     REAL(wp), PARAMETER ::  molecular_weight_of_nh4no3 = 0.08004_wp   !< mol. m. ammonium sulfate (kg mol-1)
     REAL(wp), PARAMETER ::  molecular_weight_of_water = 0.01801528_wp !< mol. m. H2O (kg mol-1)
+    REAL(wp), PARAMETER ::  pi = 3.141592654_wp                       !< PI
     REAL(wp), PARAMETER ::  rho_l = 1.0E3_wp                          !< density of water (kg m-3)
     REAL(wp), PARAMETER ::  rho_nacl = 2165.0_wp                      !< density of NaCl (kg m-3)
@@ -61,5 +64,5 @@
     REAL(wp), PARAMETER ::  r_d = 287.0_wp                            !< sp. gas const. dry air (J kg-1 K-1)
     REAL(wp), PARAMETER ::  r_v = 461.51_wp                           !< sp. gas const. water vapor (J kg-1 K-1)
-    REAL(wp), PARAMETER ::  sigma_sb = 5.67E-08_wp                    !< Stefan-Boltzmann constant
+    REAL(wp), PARAMETER ::  sigma_sb = 5.67037E-08_wp                 !< Stefan-Boltzmann constant
     REAL(wp), PARAMETER ::  solar_constant = 1368.0_wp                !< solar constant at top of atmosphere
     REAL(wp), PARAMETER ::  vanthoff_nacl = 2.0_wp                    !< van't Hoff factor for NaCl
@@ -143,5 +146,5 @@
 !
 !--    Saturation vapor pressure for a specific temperature:
-       magnus_0d =  611.2_wp * EXP( 17.62_wp * ( t - 273.15_wp ) /             &
+       magnus_0d =  611.2_wp * EXP( 17.62_wp * ( t - degc_to_k ) /             &
                                                    ( t - 29.65_wp  ) )
 
@@ -163,5 +166,5 @@
 !
 !--    Saturation vapor pressure for a specific temperature:
-       magnus_1d =  611.2_wp * EXP( 17.62_wp * ( t - 273.15_wp ) /             &
+       magnus_1d =  611.2_wp * EXP( 17.62_wp * ( t - degc_to_k ) /             &
                                                ( t - 29.65_wp  ) )
 
Index: palm/trunk/SOURCE/chemistry_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 3448)
+++ palm/trunk/SOURCE/chemistry_model_mod.f90	(revision 3449)
@@ -27,4 +27,7 @@
 ! -----------------
 ! $Id$
+! additional output - merged from branch resler
+! 
+! 3438 2018-10-28 19:31:42Z pavelkrc
 ! Add terrain-following masked output
 ! 
@@ -262,5 +265,5 @@
                                  !< neumann = zero gradient
 
-    REAL(kind=wp), PUBLIC :: cs_time_step
+    REAL(kind=wp), PUBLIC :: cs_time_step = 0._wp
     CHARACTER(10), PUBLIC :: photolysis_scheme
                                          ! 'constant',
@@ -410,7 +413,7 @@
        IF ( mode == 'allocate' )  THEN
           DO lsp = 1, nspec
-             IF (TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
+             IF ( TRIM(variable(1:3)) == 'kc_' .AND. &
+                  TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
                 chem_species(lsp)%conc_av = 0.0_wp
-                
              ENDIF
           ENDDO
@@ -418,7 +421,7 @@
        ELSEIF ( mode == 'sum' )  THEN
 
-
           DO lsp = 1, nspec
-             IF (TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
+             IF ( TRIM(variable(1:3)) == 'kc_' .AND. &
+                  TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
@@ -462,5 +465,6 @@
        ELSEIF ( mode == 'average' )  THEN
           DO lsp = 1, nspec
-             IF (TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
+             IF ( TRIM(variable(1:3)) == 'kc_' .AND. &
+                  TRIM(variable(4:)) == TRIM(chem_species(lsp)%name)) THEN
                 DO  i = nxlg, nxrg
                    DO  j = nysg, nyng
@@ -811,5 +815,20 @@
        unit = 'illegal'
 
-       spec_name = TRIM(var(4:))             !< var 1:3 is 'kc_'.
+       spec_name = TRIM(var(4:))             !< var 1:3 is 'kc_' or 'em_'.
+
+       IF ( TRIM(var(1:3)) == 'em_' )  THEN
+          DO lsp=1,nspec
+             IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name))   THEN
+                unit = 'mol m-2 s-1'
+             ENDIF
+             !-- It is possible to plant PM10 and PM25 into the gasphase chemistry code
+             !-- as passive species (e.g. 'passive' in GASPHASE_PREPROC/mechanisms):
+             !-- set unit to micrograms per m**3 for PM10 and PM25 (PM2.5)
+             IF (spec_name(1:2) == 'PM')   THEN
+                unit = 'kg m-2 s-1'
+             ENDIF
+          ENDDO
+
+       ELSE
 
           DO lsp=1,nspec
@@ -821,16 +840,17 @@
 !            set unit to micrograms per m**3 for PM10 and PM25 (PM2.5)
              IF (spec_name(1:2) == 'PM')   THEN  
-            unit = 'kg m-3'
+               unit = 'kg m-3'
              ENDIF
           ENDDO
 
-          DO lsp=1,nphot                                                     
+          DO lsp=1,nphot
              IF (TRIM(spec_name) == TRIM(phot_frequen(lsp)%name))   THEN
                 unit = 'sec-1'
              ENDIF
           ENDDO
-
-
-       RETURN
+       ENDIF
+
+
+      RETURN
     END SUBROUTINE chem_check_data_output
 !
@@ -933,4 +953,10 @@
           CALL message( 'chem_check_parameters', 'CM0426', 1, 2, 0, 6, 0 )
        ENDIF
+
+!---------------------
+!--    chem_check_parameters is called before the array chem_species is allocated!
+!--    temporary switch of this part of the check
+       RETURN
+!---------------------
 
 !--    check for initial chem species input
@@ -1076,4 +1102,7 @@
        USE kinds
 
+       USE surface_mod
+
+       !USE chem_modules, ONLY: nvar
 
        IMPLICIT NONE
@@ -1091,50 +1120,83 @@
 
        !-- local variables
-
-       INTEGER              ::  i, j, k, lsp
        CHARACTER(len=16)    ::  spec_name
+       INTEGER(iwp)         ::  i, j, k
+       INTEGER(iwp)         ::  m, l    !< running indices for surfaces
+       INTEGER(iwp)         ::  lsp  !< running index for chem spcs
 
 
        found = .FALSE.
+       IF ( .NOT. (variable(1:3) == 'kc_' .OR. variable(1:3) == 'em_' ) ) THEN
+          RETURN
+       ENDIF
 
        spec_name = TRIM(variable(4:))
 
-       DO lsp=1,nspec
-          IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name))   THEN
+       IF ( variable(1:3) == 'em_' ) THEN
+
+         local_pf = 0.0_wp
+
+         DO lsp = 1, nvar  !!! cssws - nvar species, chem_species - nspec species !!!
+            IF ( TRIM(spec_name) == TRIM(chem_species(lsp)%name) )   THEN
+               ! no average for now
+               DO m = 1, surf_usm_h%ns
+                  local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) = &
+                     local_pf(surf_usm_h%i(m),surf_usm_h%j(m),surf_usm_h%k(m)) + surf_usm_h%cssws(lsp,m)
+               ENDDO
+               DO m = 1, surf_lsm_h%ns
+                  local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) = &
+                     local_pf(surf_lsm_h%i(m),surf_lsm_h%j(m),surf_lsm_h%k(m)) + surf_lsm_h%cssws(lsp,m)
+               ENDDO
+               DO l = 0, 3
+                  DO m = 1, surf_usm_v(l)%ns
+                     local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) = &
+                        local_pf(surf_usm_v(l)%i(m),surf_usm_v(l)%j(m),surf_usm_v(l)%k(m)) + surf_usm_v(l)%cssws(lsp,m)
+                  ENDDO
+                  DO m = 1, surf_lsm_v(l)%ns
+                     local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) = &
+                        local_pf(surf_lsm_v(l)%i(m),surf_lsm_v(l)%j(m),surf_lsm_v(l)%k(m)) + surf_lsm_v(l)%cssws(lsp,m)
+                  ENDDO
+               ENDDO
+               found = .TRUE.
+            ENDIF
+         ENDDO
+       ELSE
+         DO lsp=1,nspec
+            IF (TRIM(spec_name) == TRIM(chem_species(lsp)%name))   THEN
 !
 !--   todo: remove or replace by "CALL message" mechanism (kanani)
 !              IF(myid == 0 .AND. chem_debug0 )  WRITE(6,*) 'Output of species ' // TRIM(variable)  //  &
 !                                                           TRIM(chem_species(lsp)%name)       
-             
-             IF (av == 0) THEN
-                DO  i = nxl, nxr
-                   DO  j = nys, nyn
-                      DO  k = nzb_do, nzt_do
-                          local_pf(i,j,k) = MERGE(                             &
-                                              chem_species(lsp)%conc(k,j,i),   &
-                                              REAL( fill_value, KIND = wp ),   &
-                                              BTEST( wall_flags_0(k,j,i), 0 ) )
-                      ENDDO
-                   ENDDO
-                ENDDO
-           
-             ELSE
-                DO  i = nxl, nxr
-                   DO  j = nys, nyn
-                      DO  k = nzb_do, nzt_do
-                          local_pf(i,j,k) = MERGE(                             &
-                                              chem_species(lsp)%conc_av(k,j,i),&
-                                              REAL( fill_value, KIND = wp ),   &
-                                              BTEST( wall_flags_0(k,j,i), 0 ) )
-                      ENDDO
-                   ENDDO
-                ENDDO
-             ENDIF
-
-             found = .TRUE.
-          ENDIF
-       ENDDO
+               IF (av == 0) THEN
+                  DO  i = nxl, nxr
+                     DO  j = nys, nyn
+                        DO  k = nzb_do, nzt_do
+                            local_pf(i,j,k) = MERGE(                             &
+                                                chem_species(lsp)%conc(k,j,i),   &
+                                                REAL( fill_value, KIND = wp ),   &
+                                                BTEST( wall_flags_0(k,j,i), 0 ) )
+                        ENDDO
+                     ENDDO
+                  ENDDO
+
+               ELSE
+                  DO  i = nxl, nxr
+                     DO  j = nys, nyn
+                        DO  k = nzb_do, nzt_do
+                            local_pf(i,j,k) = MERGE(                             &
+                                                chem_species(lsp)%conc_av(k,j,i),&
+                                                REAL( fill_value, KIND = wp ),   &
+                                                BTEST( wall_flags_0(k,j,i), 0 ) )
+                        ENDDO
+                     ENDDO
+                  ENDDO
+               ENDIF
+               found = .TRUE.
+            ENDIF
+         ENDDO
+       ENDIF
 
        RETURN
+
     END SUBROUTINE chem_data_output_3d
 !
@@ -1177,5 +1239,5 @@
 
        spec_name = TRIM( variable(4:) )
-   !av == 0
+       !av == 0
 
        DO lsp=1,nspec
@@ -1292,5 +1354,5 @@
        found  = .TRUE.
 
-       IF (var(1:3) == 'kc_')   THEN                   !< always the same grid for chemistry variables
+       IF ( var(1:3) == 'kc_' .OR. var(1:3) == 'em_' )   THEN                   !< always the same grid for chemistry variables
           grid_x = 'x'
           grid_y = 'y'
@@ -1873,11 +1935,11 @@
  !--   Enable chemistry model
        air_chemistry = .TRUE.                    
-      GOTO 20
-
- 10   BACKSPACE( 11 )
-      READ( 11 , '(A)') line
-      CALL parin_fail_message( 'chemistry_parameters', line )
-
- 20   CONTINUE
+       GOTO 20
+
+ 10    BACKSPACE( 11 )
+       READ( 11 , '(A)') line
+       CALL parin_fail_message( 'chemistry_parameters', line )
+
+ 20    CONTINUE
 
 !
Index: palm/trunk/SOURCE/plant_canopy_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 3448)
+++ palm/trunk/SOURCE/plant_canopy_model_mod.f90	(revision 3449)
@@ -16,4 +16,6 @@
 !
 ! Copyright 1997-2018 Leibniz Universitaet Hannover
+! Copyright 2018 Institute of Computer Science of the
+!                Czech Academy of Sciences, Prague
 !------------------------------------------------------------------------------!
 !
@@ -25,4 +27,10 @@
 ! -----------------
 ! $Id$
+! Add calculation of transpiration for resolved plant canopy (trees, shrubs)
+! (V. Fuka, MFF UK Prague, J.Resler, ICS AS, Prague)
+!
+! Fix reading plant canopy over buildings
+! 
+! 3337 2018-10-12 15:17:09Z kanani
 ! Fix reading plant canopy over buildings
 ! 
@@ -214,5 +222,11 @@
  
     USE arrays_3d,                                                             &
-        ONLY:  dzu, dzw, e, q, s, tend, u, v, w, zu, zw 
+        ONLY:  dzu, dzw, e, exner, hyp, pt, q, s, tend, u, v, w, zu, zw
+
+    USE basic_constants_and_equations_mod,                                     &
+        ONLY:  c_p, degc_to_k, l_v, lv_d_cp, r_d, rd_d_rv
+
+    USE control_parameters,                                                   &
+        ONLY: humidity
 
     USE indices,                                                               &
@@ -222,4 +236,6 @@
     USE kinds
 
+    USE pegrid
+
     USE surface_mod,                                                           &
         ONLY:  get_topography_top_index_ji
@@ -229,12 +245,15 @@
 
 
-    CHARACTER (LEN=30)   ::  canopy_mode = 'block' !< canopy coverage
-
-    INTEGER(iwp) ::  pch_index = 0                               !< plant canopy height/top index
-    INTEGER(iwp) ::  lad_vertical_gradient_level_ind(10) = -9999 !< lad-profile levels (index)
-
-    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  pch_index_ji   !< local plant canopy top
-
-    LOGICAL ::  calc_beta_lad_profile = .FALSE. !< switch for calc. of lad from beta func.
+    CHARACTER (LEN=30)   ::  canopy_mode = 'block'                 !< canopy coverage
+    LOGICAL              ::  plant_canopy_transpiration = .FALSE.  !< flag to switch calculation of transpiration and corresponding latent heat
+                                                                   !< for resolved plant canopy inside radiation model
+                                                                   !< (calls subroutine pcm_calc_transpiration_rate from module plant_canopy_mod)
+
+    INTEGER(iwp) ::  pch_index = 0                                 !< plant canopy height/top index
+    INTEGER(iwp) ::  lad_vertical_gradient_level_ind(10) = -9999   !< lad-profile levels (index)
+
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE ::  pch_index_ji     !< local plant canopy top
+
+    LOGICAL ::  calc_beta_lad_profile = .FALSE.   !< switch for calc. of lad from beta func.
 
     REAL(wp) ::  alpha_lad = 9999999.9_wp         !< coefficient for lad calculation
@@ -261,8 +280,9 @@
     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 ::  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
 
     SAVE
@@ -273,12 +293,17 @@
 !
 !-- Public functions
-    PUBLIC pcm_check_data_output, pcm_check_parameters, pcm_data_output_3d,    &
-           pcm_define_netcdf_grid, pcm_header, pcm_init, pcm_parin, pcm_tendency
+    PUBLIC pcm_calc_transpiration_rate, pcm_check_data_output,                &
+           pcm_check_parameters, pcm_data_output_3d, pcm_define_netcdf_grid,  &
+           pcm_header, pcm_init, pcm_parin, pcm_tendency
 
 !
 !-- Public variables and constants
-    PUBLIC pc_heating_rate, pc_transpiration_rate, canopy_mode, cthf, dt_plant_canopy, lad, lad_s,   &
-           pch_index
-           
+    PUBLIC pc_heating_rate, pc_transpiration_rate, pc_latent_rate,             &
+            canopy_mode, cthf, dt_plant_canopy, lad, lad_s, pch_index,         &
+            plant_canopy_transpiration
+
+    INTERFACE pcm_calc_transpiration_rate
+       MODULE PROCEDURE pcm_calc_transpiration_rate
+    END INTERFACE pcm_calc_transpiration_rate
 
     INTERFACE pcm_check_data_output
@@ -323,4 +348,118 @@
 
 
+
+!------------------------------------------------------------------------------!
+! Description:
+! ------------
+!> Calculation of the plant canopy transpiration rate based on the Jarvis-Stewart
+!> with parametrizations described in Daudet et al. (1999; Agricult. and Forest
+!> Meteorol. 97) and Ngao, Adam and Saudreau (2017;  Agricult. and Forest Meteorol
+!> 237-238). Model functions f1-f4 were adapted from Stewart (1998; Agric.
+!> and Forest. Meteorol. 43) instead, because they are valid for broader intervals
+!> of values. Funcion f4 used in form present in van Wijk et al. (1998;
+!> Tree Physiology 20).
+!>
+!> This subroutine is called from subroutine radiation_interaction
+!> after the calculation of radiation in plant canopy boxes.
+!> (arrays pcbinsw and pcbinlw).
+!>
+!------------------------------------------------------------------------------!
+ SUBROUTINE pcm_calc_transpiration_rate(i, j, k, kk, pcbsw, pcblw, pcbtr, pcblh)
+
+     USE control_parameters,                                                   &
+        ONLY: dz
+
+     USE grid_variables,                                                       &
+        ONLY:  dx, dy
+
+     IMPLICIT NONE
+!--  input parameters
+     INTEGER(iwp), INTENT(IN)          :: i, j, k, kk        !< indices of the pc gridbox
+     REAL(wp), INTENT(IN)              :: pcbsw              !< sw radiation in gridbox (W)
+     REAL(wp), INTENT(IN)              :: pcblw              !< lw radiation in gridbox (W)
+     REAL(wp), INTENT(OUT)             :: pcbtr              !< transpiration rate dq/dt (kg/kg/s)
+     REAL(wp), INTENT(OUT)             :: pcblh              !< latent heat from transpiration dT/dt (K/s)
+
+!--  variables and parameters for calculation of transpiration rate
+     REAL(wp)                          :: sat_press, sat_press_d, temp, v_lad
+     REAL(wp)                          :: d_fact, g_b, g_s, wind_speed, evapor_rate
+     REAL(wp)                          :: f1, f2, f3, f4, vpd, rswc, e_eq, e_imp, rad
+     REAL(wp), PARAMETER               :: gama_psychr = 66   !< psychrometric constant (Pa/K)
+     REAL(wp), PARAMETER               :: g_s_max = 0.01     !< maximum stomatal conductivity (m/s)
+     REAL(wp), PARAMETER               :: m_soil = 0.4_wp    !< soil water content (needs to adjust or take from LSM)
+     REAL(wp), PARAMETER               :: m_wilt = 0.01_wp   !< wilting point soil water content (needs to adjust or take from LSM)
+     REAL(wp), PARAMETER               :: m_sat = 0.51_wp    !< saturation soil water content (needs to adjust or take from LSM)
+     REAL(wp), PARAMETER               :: t2_min = 0.0_wp    !< minimal temperature for calculation of f2
+     REAL(wp), PARAMETER               :: t2_max = 40.0_wp   !< maximal temperature for calculation of f2
+
+
+!--  Temperature (deg C)
+     temp = pt(k,j,i) * exner(k) - degc_to_k
+!--  Coefficient for conversion of radiation to grid to radiation to unit leaves surface
+     v_lad = 1.0_wp / ( MAX(lad_s(kk,j,i), 1.0e-10) * dx * dy * dz(1) )
+!--  Magnus formula for the saturation pressure (see Ngao, Adam and Saudreau (2017) eq. 1)
+!--  There are updated formulas available, kept consistent with the rest of the parametrization
+     sat_press = 610.8_wp * exp(17.27_wp * temp/(temp + 237.3_wp))
+!--  Saturation pressure derivative (derivative of the above)
+     sat_press_d = sat_press * 17.27_wp * 237.3_wp / (temp + 237.3_wp)**2
+!--  Wind speed
+     wind_speed = SQRT( ((u(k,j,i) + u(k,j,i+1))/2.0_wp )**2 +           &
+                        ((v(k,j,i) + v(k,j,i+1))/2.0_wp )**2 +           &
+                        ((w(k,j,i) + w(k,j,i+1))/2.0_wp )**2 )
+!--  Aerodynamic conductivity (Daudet et al. (1999) eq. 14
+     g_b = 0.01_wp * wind_speed + 0.0071_wp
+!--  Radiation flux per leaf surface unit
+     rad = pcbsw * v_lad
+!--  First function for calculation of stomatal conductivity (radiation dependency)
+!--  Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 17
+     f1 = rad * (1000._wp+42.1_wp) / 1000._wp / (rad+42.1_wp)
+!--  Second function for calculation of stomatal conductivity (temperature dependency)
+!--  Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 21
+     f2 = MAX(t2_min, (temp-t2_min) * (t2_max-temp)**((t2_max-16.9_wp)/(16.9_wp-t2_min)) / &
+              ((16.9_wp-t2_min) * (t2_max-16.9_wp)**((t2_max-16.9_wp)/(16.9_wp-t2_min))) )
+!--  Water pressure deficit
+!--  Ngao, Adam and Saudreau (2017) eq. 6 but with water vapour partial pressure
+     vpd = max( sat_press - q(k,j,i) * hyp(k) / rd_d_rv, 0._wp )
+!--  Third function for calculation of stomatal conductivity (water pressure deficit dependency)
+!--  Ngao, Adam and Saudreau (2017) Table 1, limited from below according to Stewart (1988)
+!--  The coefficients of the linear dependence should better correspond to broad-leaved trees
+!--  than the coefficients from Stewart (1988) which correspond to conifer trees.
+     vpd = MIN(MAX(vpd,770.0_wp),3820.0_wp)
+     f3 = -2e-4_wp * vpd + 1.154_wp
+!--  Fourth function for calculation of stomatal conductivity (soil moisture dependency)
+!--  Residual soil water content
+!--  van Wijk et al. (1998; Tree Physiology 20) eq. 7
+!--  TODO - over LSM surface might be calculated from LSM parameters
+     rswc = ( m_sat - m_soil ) / ( m_sat - m_wilt )
+!--  van Wijk et al. (1998; Tree Physiology 20) eq. 5-6 (it is a reformulation of eq. 22-23 of Stewart(1988))
+     f4 = MAX(0._wp, MIN(1.0_wp - 0.041_wp * EXP(3.2_wp * rswc), 1.0_wp - 0.041_wp))
+!--  Stomatal conductivity
+!--  Stewart (1988; Agric. and Forest. Meteorol. 43) eq. 12
+!--  (notation according to Ngao, Adam and Saudreau (2017) and others)
+     g_s = g_s_max * f1 * f2 * f3 * f4 + 1.0e-10_wp
+!--  Decoupling factor
+!--  Daudet et al. (1999) eq. 6
+     d_fact = (sat_press_d / gama_psychr + 2._wp ) /                        &
+              (sat_press_d / gama_psychr + 2._wp + 2 * g_b / g_s )
+!--  Equilibrium evaporation rate
+!--  Daudet et al. (1999) eq. 4
+     e_eq = (pcbsw + pcblw) * v_lad * sat_press_d /                         &
+                 gama_psychr /( sat_press_d / gama_psychr + 2.0_wp ) / l_v
+!--  Imposed evaporation rate
+!--  Daudet et al. (1999) eq. 5
+     e_imp = r_d * pt(k,j,i) * exner(k) / hyp(k) * c_p * g_s * vpd / gama_psychr / l_v
+!--  Evaporation rate
+!--  Daudet et al. (1999) eq. 3
+!--  (evaporation rate is limited to non-negative values)
+     evapor_rate = MAX(d_fact * e_eq + ( 1.0_wp - d_fact ) * e_imp, 0.0_wp)
+!--  Conversion of evaporation rate to q tendency in gridbox
+!--  dq/dt = E * LAD * V_g / (rho_air * V_g)
+     pcbtr = evapor_rate * r_d * pt(k,j,i) * exner(k) * lad_s(kk,j,i) / hyp(k)  !-- = dq/dt
+!--  latent heat from evaporation
+     pcblh = pcbtr * lv_d_cp  !-- = - dT/dt
+
+ END SUBROUTINE pcm_calc_transpiration_rate
+
+
 !------------------------------------------------------------------------------!
 ! Description:
@@ -352,4 +491,10 @@
        CASE ( 'pcm_transpirationrate' )
           unit = 'kg kg-1 s-1'
+
+       CASE ( 'pcm_latentrate' )
+          unit = 'K s-1'
+
+       CASE ( 'pcm_bowenratio' )
+          unit = 'K s-1'
 
        CASE ( 'pcm_lad' )
@@ -503,5 +648,36 @@
             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
+
+       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
@@ -550,5 +726,5 @@
      SELECT CASE ( TRIM( var ) )
 
-        CASE ( 'pcm_heatrate', 'pcm_lad', 'pcm_transpirationrate')
+        CASE ( 'pcm_heatrate', 'pcm_lad', 'pcm_transpirationrate', 'pcm_latentrate', 'pcm_bowenratio')
            grid_x = 'x'
            grid_y = 'y'
@@ -691,5 +867,5 @@
 
        USE control_parameters,                                                 &
-           ONLY: humidity, message_string, ocean_mode, urban_surface
+           ONLY: message_string, ocean_mode, urban_surface
 
        USE netcdf_data_input_mod,                                              &
@@ -705,4 +881,5 @@
        INTEGER(iwp) ::  k   !< running index
        INTEGER(iwp) ::  m   !< running index
+       INTEGER(iwp) ::  pch_index_l
 
        REAL(wp) ::  int_bpdf        !< vertical integral for lad-profile construction
@@ -903,4 +1080,24 @@
        CALL exchange_horiz_2d_int( pch_index_ji, nys, nyn, nxl, nxr, nbgp )
 
+!--    Calculate global pch_index value (index of top of plant canopy from ground)
+       pch_index = MAXVAL(pch_index_ji)
+!--    Exchange pch_index from all processors
+#if defined( __parallel )
+       CALL MPI_Allreduce(MPI_IN_PLACE, pch_index, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr)
+       IF ( ierr /= 0 )  THEN
+          WRITE (9,*) 'Error MPI_Allreduce pch_index:', ierr
+          FLUSH(9)
+       ENDIF
+#endif
+
+!--    Allocation of arrays pc_heating_rate, pc_transpiration_rate and pc_latent_rate
+       ALLOCATE( pc_heating_rate(0:pch_index,nysg:nyng,nxlg:nxrg) )
+       IF ( humidity )  THEN
+          ALLOCATE( pc_transpiration_rate(0:pch_index,nysg:nyng,nxlg:nxrg) )
+          pc_transpiration_rate = 0.0_wp
+          ALLOCATE( pc_latent_rate(0:pch_index,nysg:nyng,nxlg:nxrg) )
+          pc_latent_rate = 0.0_wp
+       ENDIF
+
 !
 !--    Initialization of the canopy heat source distribution due to heating
@@ -915,10 +1112,12 @@
 !--    73â96). This approach has been applied e.g. by Shaw & Schumann (1992;
 !--    Bound.-Layer Meteorol. 61, 47â64).
-!--    When using the urban surface model (USM), canopy heating (pc_heating_rate)
-!--    by radiation is calculated in the USM.
-       IF ( cthf /= 0.0_wp  .AND. .NOT.  urban_surface )  THEN
-
-          ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg),                 &
-                    pc_heating_rate(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
+!--    When using the radiation_interactions, canopy heating (pc_heating_rate)
+!--    and plant canopy transpiration (pc_transpiration_rate, pc_latent_rate)
+!--    are calculated in the RTM after the calculation of radiation.
+!--    We cannot use variable radiation_interactions here to determine the situation
+!--    as it is assigned in init_3d_model after the call of pcm_init.
+       IF ( cthf /= 0.0_wp )  THEN
+
+          ALLOCATE( cum_lai_hf(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
 !
 !--       Piecewise calculation of the cumulative leaf area index by vertical
@@ -996,9 +1195,4 @@
 
        ENDIF
-!
-!--    Allocate transpiration rate
-       IF ( humidity )                                                         &
-          ALLOCATE( pc_transpiration_rate(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
-
 
 
@@ -1028,5 +1222,6 @@
                                   lai_beta,                                    &
                                   leaf_scalar_exch_coeff,                      &
-                                  leaf_surface_conc, pch_index
+                                  leaf_surface_conc, pch_index,                &
+                                  plant_canopy_transpiration
 
        NAMELIST /canopy_par/      alpha_lad, beta_lad, canopy_drag_coeff,      &
@@ -1037,5 +1232,6 @@
                                   lai_beta,                                    &
                                   leaf_scalar_exch_coeff,                      &
-                                  leaf_surface_conc, pch_index
+                                  leaf_surface_conc, pch_index,                &
+                                  plant_canopy_transpiration
 
        line = ' '
@@ -1434,17 +1630,27 @@
 !--       potential temperature
           CASE ( 4 )
-             DO  i = nxl, nxr
-                DO  j = nys, nyn
-!
-!--                Determine topography-top index on scalar-grid
-                   k_wall = get_topography_top_index_ji( j, i, 's' )
-
-                   DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
-
-                      kk = k - k_wall   !- lad arrays are defined flat
-                      tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i)
+             IF ( humidity ) THEN
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+!--                   Determine topography-top index on scalar-grid
+                      k_wall = get_topography_top_index_ji( j, i, 's' )
+                      DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
+                         kk = k - k_wall   !- lad arrays are defined flat
+                         tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i) - pc_latent_rate(kk,j,i)
+                      ENDDO
                    ENDDO
                 ENDDO
-             ENDDO
+             ELSE
+                DO  i = nxl, nxr
+                   DO  j = nys, nyn
+!--                   Determine topography-top index on scalar-grid
+                      k_wall = get_topography_top_index_ji( j, i, 's' )
+                      DO  k = k_wall+1, k_wall + pch_index_ji(j,i)
+                         kk = k - k_wall   !- lad arrays are defined flat
+                         tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i)
+                      ENDDO
+                   ENDDO
+                ENDDO
+             ENDIF
 
 !
@@ -1460,17 +1666,24 @@
 
                       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) )         &
-                                             )**2 +                            &
-                                             ( 0.5_wp * ( v(k,j,i) +           &
-                                                          v(k,j+1,i) )         &
-                                             )**2 +                            &
-                                             ( 0.5_wp * ( w(k-1,j,i) +         & 
-                                                          w(k,j,i) )           &
-                                             )**2                              &
-                                           ) *                                 &
-                                       ( q(k,j,i) - lsc )
+
+                      IF ( .NOT. plant_canopy_transpiration ) THEN
+                         ! pc_transpiration_rate is calculated in radiation model
+                         ! in case of plant_canopy_transpiration = .T.
+                         ! to include also the dependecy to the radiation
+                         ! in the plant canopy box
+                         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) )         &
+                                                )**2 +                            &
+                                                ( 0.5_wp * ( v(k,j,i) +           &
+                                                             v(k,j+1,i) )         &
+                                                )**2 +                            &
+                                                ( 0.5_wp * ( w(k-1,j,i) +         &
+                                                             w(k,j,i) )           &
+                                                )**2                              &
+                                              ) *                                 &
+                                          ( q(k,j,i) - lsc )
+                      ENDIF
 
                       tend(k,j,i) = tend(k,j,i) + pc_transpiration_rate(kk,j,i)
@@ -1778,33 +1991,44 @@
              k_wall = get_topography_top_index_ji( j, i, 's' )
 
+             IF ( humidity ) THEN
+                DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
+                   kk = k - k_wall  !- lad arrays are defined flat
+                   tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i) - pc_latent_rate(kk,j,i)
+                ENDDO
+             ELSE
+                DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
+                   kk = k - k_wall  !- lad arrays are defined flat
+                   tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i)
+                ENDDO
+             ENDIF
+
+!
+!--       humidity
+          CASE ( 5 )
+!
+!--          Determine topography-top index on scalar grid
+             k_wall = get_topography_top_index_ji( j, i, 's' )
+
              DO  k = k_wall + 1, k_wall + pch_index_ji(j,i)
                 kk = k - k_wall  !- lad arrays are defined flat
-                tend(k,j,i) = tend(k,j,i) + pc_heating_rate(kk,j,i)
-             ENDDO
-
-
-!
-!--       humidity
-          CASE ( 5 )
-!
-!--          Determine topography-top index on scalar grid
-             k_wall = get_topography_top_index_ji( j, i, 's' )
-
-             DO  k = k_wall + 1, k_wall + pch_index_ji(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) )               &
-                                       )**2  +                                 &
-                                       ( 0.5_wp * ( v(k,j,i) +                 &
-                                                    v(k,j+1,i) )               &
-                                       )**2 +                                  &
-                                       ( 0.5_wp * ( w(k-1,j,i) +               &
-                                                    w(k,j,i) )                 &
-                                       )**2                                    &
-                                     ) *                                       &
-                                 ( q(k,j,i) - lsc )
+                IF ( .NOT. plant_canopy_transpiration ) THEN
+                   ! pc_transpiration_rate is calculated in radiation model
+                   ! in case of plant_canopy_transpiration = .T.
+                   ! to include also the dependecy to the radiation
+                   ! in the plant canopy box
+                   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) )               &
+                                          )**2  +                                 &
+                                          ( 0.5_wp * ( v(k,j,i) +                 &
+                                                       v(k,j+1,i) )               &
+                                          )**2 +                                  &
+                                          ( 0.5_wp * ( w(k-1,j,i) +               &
+                                                       w(k,j,i) )                 &
+                                          )**2                                    &
+                                        ) *                                       &
+                                    ( q(k,j,i) - lsc )
+                ENDIF
 
                 tend(k,j,i) = tend(k,j,i) + pc_transpiration_rate(kk,j,i)
Index: palm/trunk/SOURCE/radiation_model_mod.f90
===================================================================
--- palm/trunk/SOURCE/radiation_model_mod.f90	(revision 3448)
+++ palm/trunk/SOURCE/radiation_model_mod.f90	(revision 3449)
@@ -28,4 +28,11 @@
 ! -----------------
 ! $Id$
+! New RTM version 3.0: (Pavel Krc, Jaroslav Resler, ICS, Prague)
+!   - Interaction of plant canopy with LW radiation
+!   - Transpiration from resolved plant canopy dependent on radiation
+!     called from RTM
+! 
+! 
+! 3435 2018-10-26 18:25:44Z gronemeier
 ! - workaround: return unit=illegal in check_data_output for certain variables
 !   when check called from init_masks
@@ -499,5 +506,5 @@
 
     USE basic_constants_and_equations_mod,                                     &
-        ONLY:  c_p, g, lv_d_cp, l_v, pi, r_d, rho_l, solar_constant,           &
+        ONLY:  c_p, g, lv_d_cp, l_v, pi, r_d, rho_l, solar_constant,      &
                barometric_formula
 
@@ -544,5 +551,6 @@
 
     USE plant_canopy_model_mod,                                                &
-        ONLY:  lad_s, pc_heating_rate, pc_transpiration_rate
+        ONLY:  lad_s, pc_heating_rate, pc_transpiration_rate, pc_latent_rate,  &
+               plant_canopy_transpiration, pcm_calc_transpiration_rate
 
     USE pegrid
@@ -853,5 +861,6 @@
     INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  jdir = (/0, 0,1,-1,0, 0,0,1,-1,0, 0/)   !< surface normal direction y indices
     INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  kdir = (/1,-1,0, 0,0, 0,1,0, 0,0, 0/)   !< surface normal direction z indices
-                                                                                          !< parameter but set in the code
+    REAL(wp),     DIMENSION(0:nsurf_type)          :: facearea                            !< area of single face in respective
+                                                                                          !< direction (will be calc'd)
 
 
@@ -887,15 +896,14 @@
 !-- configuration parameters (they can be setup in PALM config)
     LOGICAL                                        ::  raytrace_mpi_rma = .TRUE.          !< use MPI RMA to access LAD and gridsurf from remote processes during raytracing
-    LOGICAL                                        ::  read_svf_on_init = .FALSE.         !< flag parameter indicating wheather SVFs will be read from a file at initialization
-    LOGICAL                                        ::  write_svf_on_init = .FALSE.        !< flag parameter indicating wheather SVFs will be written out to a file
-    LOGICAL                                        ::  rad_angular_discretization = .TRUE.!< whether to use fixed resolution discretization of view factors for 
+    LOGICAL                                        ::  rad_angular_discretization = .TRUE.!< whether to use fixed resolution discretization of view factors for
                                                                                           !< reflected radiation (as opposed to all mutually visible pairs)
+    LOGICAL                                        ::  plant_lw_interact = .TRUE.         !< whether plant canopy interacts with LW radiation (in addition to SW)
     INTEGER(iwp)                                   ::  mrt_nlevels = 0                    !< number of vertical boxes above surface for which to calculate MRT
     LOGICAL                                        ::  mrt_skip_roof = .TRUE.             !< do not calculate MRT above roof surfaces
     LOGICAL                                        ::  mrt_include_sw = .TRUE.            !< should MRT calculation include SW radiation as well?
-    INTEGER(iwp)                                   ::  nrefsteps = 0                      !< number of reflection steps to perform
+    INTEGER(iwp)                                   ::  nrefsteps = 3                      !< number of reflection steps to perform
     REAL(wp), PARAMETER                            ::  ext_coef = 0.6_wp                  !< extinction coefficient (a.k.a. alpha)
     INTEGER(iwp), PARAMETER                        ::  rad_version_len = 10               !< length of identification string of rad version
-    CHARACTER(rad_version_len), PARAMETER          ::  rad_version = 'RAD v. 1.1'         !< identification of version of binary svf and restart files
+    CHARACTER(rad_version_len), PARAMETER          ::  rad_version = 'RAD v. 3.0'         !< identification of version of binary svf and restart files
     INTEGER(iwp)                                   ::  raytrace_discrete_elevs = 40       !< number of discretization steps for elevation (nadir to zenith)
     INTEGER(iwp)                                   ::  raytrace_discrete_azims = 80       !< number of discretization steps for azimuth (out of 360 degrees)
@@ -929,6 +937,7 @@
         INTEGER(iwp)                               :: ity               !< 
         INTEGER(iwp)                               :: itz               !< 
-        INTEGER(iwp)                               :: isurfs            !< 
-        REAL(wp)                                   :: rsvf              !< 
+        INTEGER(iwp)                               :: isurfs            !< Idx of source face / -1 for sky
+        REAL(wp)                                   :: rcvf              !< Canopy view factor for faces /
+                                                                        !< canopy sink factor for sky (-1)
     END TYPE
 
@@ -976,4 +985,5 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfouts         !< array of reflected sw radiation for all surfaces in i-th reflection
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutl         !< array of reflected + emitted lw radiation for all surfaces in i-th reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlg         !< global array of incoming lw radiation from plant canopy
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutsw        !< array of total sw radiation outgoing from nonvirtual surfaces surfaces after all reflection
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutlw        !< array of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
@@ -2928,33 +2938,32 @@
        CHARACTER (LEN=80) ::  line  !< dummy string that contains the current line of the parameter file 
        
-       NAMELIST /radiation_par/   albedo, albedo_type, albedo_lw_dir,          &
-                                  albedo_lw_dif, albedo_sw_dir, albedo_sw_dif, &
-                                  constant_albedo, dt_radiation, emissivity,   &
-                                  lw_radiation, mrt_nlevels, mrt_skip_roof,    &
-                                  mrt_include_sw,  net_radiation,              &
-                                  radiation_scheme, skip_time_do_radiation,    &
-                                  sw_radiation, unscheduled_radiation_calls,   &
-                                  max_raytracing_dist, min_irrf_value,         &
-                                  nrefsteps, raytrace_mpi_rma,                 &
-                                  surface_reflections, svfnorm_report_thresh,  &
-                                  radiation_interactions_on,                   &
-                                  rad_angular_discretization,                  &
-                                  raytrace_discrete_azims,                     &
-                                  raytrace_discrete_elevs
+       NAMELIST /radiation_par/   albedo, albedo_lw_dif, albedo_lw_dir,         &
+                                  albedo_sw_dif, albedo_sw_dir, albedo_type,    &
+                                  constant_albedo, dt_radiation, emissivity,    &
+                                  lw_radiation, max_raytracing_dist,            &
+                                  min_irrf_value, mrt_include_sw, mrt_nlevels,  &
+                                  mrt_skip_roof, net_radiation, nrefsteps,      &
+                                  plant_lw_interact, rad_angular_discretization,&
+                                  radiation_interactions_on, radiation_scheme,  &
+                                  raytrace_discrete_azims,                      &
+                                  raytrace_discrete_elevs, raytrace_mpi_rma,    &
+                                  skip_time_do_radiation, surface_reflections,  &
+                                  svfnorm_report_thresh, sw_radiation,          &
+                                  unscheduled_radiation_calls
+
    
-       NAMELIST /radiation_parameters/   albedo, albedo_type, albedo_lw_dir,   &
-                                  albedo_lw_dif, albedo_sw_dir, albedo_sw_dif, &
-                                  constant_albedo, dt_radiation, emissivity,   &
-                                  lw_radiation, mrt_nlevels, mrt_skip_roof,    &
-                                  mrt_include_sw,  net_radiation,              &
-                                  radiation_scheme, skip_time_do_radiation,    &
-                                  sw_radiation, unscheduled_radiation_calls,   &
-                                  max_raytracing_dist, min_irrf_value,         &
-                                  nrefsteps, raytrace_mpi_rma,                 &
-                                  surface_reflections, svfnorm_report_thresh,  &
-                                  radiation_interactions_on,                   &
-                                  rad_angular_discretization,                  &
-                                  raytrace_discrete_azims,                     &
-                                  raytrace_discrete_elevs
+       NAMELIST /radiation_parameters/ albedo, albedo_lw_dif, albedo_lw_dir,    &
+                                  albedo_sw_dif, albedo_sw_dir, albedo_type,    &
+                                  constant_albedo, dt_radiation, emissivity,    &
+                                  lw_radiation, max_raytracing_dist,            &
+                                  min_irrf_value, mrt_include_sw, mrt_nlevels,  &
+                                  mrt_skip_roof, net_radiation, nrefsteps,      &
+                                  plant_lw_interact, rad_angular_discretization,&
+                                  radiation_interactions_on, radiation_scheme,  &
+                                  raytrace_discrete_azims,                      &
+                                  raytrace_discrete_elevs, raytrace_mpi_rma,    &
+                                  skip_time_do_radiation, surface_reflections,  &
+                                  svfnorm_report_thresh, sw_radiation,          &
+                                  unscheduled_radiation_calls
    
        line = ' '
@@ -4648,8 +4657,8 @@
      REAL(wp), DIMENSION(0:nsurf_type) :: costheta           !< direct irradiance factor of solar angle
      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
-     REAL(wp), DIMENSION(0:nsurf_type) :: facearea
+     REAL(wp)                          :: asrc               !< area of source face
+     REAL(wp)                          :: pcrad              !< irradiance from plant canopy
      REAL(wp)                          :: pabsswl  = 0.0_wp  !< total absorbed SW radiation energy in local processor (W)
      REAL(wp)                          :: pabssw   = 0.0_wp  !< total absorbed SW radiation energy in all processors (W)
@@ -4672,6 +4681,4 @@
          pchf_prep(:) = r_d * exner(nzub:nzut)                                 &
                      / (c_p * hyp(nzub:nzut) * dx*dy*dz(1)) !< equals to 1 / (rho * c_p * Vbox * T)
-         pctf_prep(:) = r_d * exner(nzub:nzut)                                 &
-                     / (l_v * hyp(nzub:nzut) * dx*dy*dz(1))
      ENDIF
 #endif
@@ -4729,4 +4736,5 @@
         mrtinlw(:) = 0._wp
      ENDIF
+     surfinlg(:)  = 0._wp !global
 
 
@@ -4824,5 +4832,9 @@
 !
 !--        For surface-to-surface factors we calculate thermal radiation in 1st pass
-           surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc)
+           IF ( plant_lw_interact )  THEN
+              surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * svf(2,isvf) * surfoutl(isurfsrc)
+           ELSE
+              surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc)
+           ENDIF
         ENDDO
      ENDIF
@@ -4833,5 +4845,9 @@
         i = surfl(ix, isurf)
         surfinswdif(isurf) = rad_sw_in_diff(j,i) * skyvft(isurf)
-        surfinlwdif(isurf) = rad_lw_in_diff(j,i) * skyvf(isurf)
+        IF ( plant_lw_interact )  THEN
+           surfinlwdif(isurf) = rad_lw_in_diff(j,i) * skyvft(isurf)
+        ELSE
+           surfinlwdif(isurf) = rad_lw_in_diff(j,i) * skyvf(isurf)
+        ENDIF
      ENDDO
 !
@@ -4880,7 +4896,7 @@
          pcbinswdir(:) = 0._wp
          pcbinswdif(:) = 0._wp
-         pcbinlw(:) = 0._wp  !< will stay always 0 since we don't absorb lw anymore
-!
-!--         pcsf first pass
+         pcbinlw(:) = 0._wp
+!
+!--      pcsf first pass
          DO icsf = 1, ncsfl
              ipcgb = csfsurf(1, icsf)
@@ -4891,16 +4907,38 @@
 
              IF ( isurfsrc == -1 )  THEN
-!--                 Diffuse rad from sky.
-                 pcbinswdif(ipcgb) = csf(1,icsf) * rad_sw_in_diff(j,i)
-
-                 !--Direct rad
-                 IF ( zenith(0) > 0 )  THEN
-                    !--Estimate directed box absorption
-                    pc_abs_frac = 1._wp - exp(pc_abs_eff * lad_s(k,j,i))
-
-                    !--isd has already been established, see 1)
-                    pcbinswdir(ipcgb) = rad_sw_in_dir(j, i) * pc_box_area &
-                                        * pc_abs_frac * dsitransc(ipcgb, isd)
-                 ENDIF
+!
+!--             Diffuse rad from sky.
+                pcbinswdif(ipcgb) = csf(1,icsf) * rad_sw_in_diff(j,i)
+!
+!--             Absorbed diffuse LW from sky minus emitted to sky
+                IF ( plant_lw_interact )  THEN
+                   pcbinlw(ipcgb) = csf(1,icsf)                                  &
+                                       * (rad_lw_in_diff(j, i)                   &
+                                          - sigma_sb * (pt(k,j,i)*exner(k))**4)
+                ENDIF
+!
+!--             Direct rad
+                IF ( zenith(0) > 0 )  THEN
+!--                Estimate directed box absorption
+                   pc_abs_frac = 1._wp - exp(pc_abs_eff * lad_s(k,j,i))
+!
+!--                isd has already been established, see 1)
+                   pcbinswdir(ipcgb) = rad_sw_in_dir(j, i) * pc_box_area &
+                                       * pc_abs_frac * dsitransc(ipcgb, isd)
+                ENDIF
+             ELSE
+                IF ( plant_lw_interact )  THEN
+!
+!--                Thermal emission from plan canopy towards respective face
+                   pcrad = sigma_sb * (pt(k,j,i) * exner(k))**4 * csf(1,icsf)
+                   surfinlg(isurfsrc) = surfinlg(isurfsrc) + pcrad
+!
+!--                Remove the flux above + absorb LW from first pass from surfaces
+                   asrc = facearea(surf(id, isurfsrc))
+                   pcbinlw(ipcgb) = pcbinlw(ipcgb)                      &
+                                    + (csf(1,icsf) * surfoutl(isurfsrc) & ! Absorb from first pass surf emit
+                                       - pcrad)                         & ! Remove emitted heatflux
+                                    * asrc
+                ENDIF
              ENDIF
          ENDDO
@@ -4908,4 +4946,20 @@
          pcbinsw(:) = pcbinswdir(:) + pcbinswdif(:)
      ENDIF
+
+     IF ( plant_lw_interact )  THEN
+!
+!--     Exchange incoming lw radiation from plant canopy
+#if defined( __parallel )
+        CALL MPI_Allreduce(MPI_IN_PLACE, surfinlg, nsurf, MPI_REAL, MPI_SUM, comm2d, ierr)
+        IF ( ierr /= 0 )  THEN
+           WRITE (9,*) 'Error MPI_Allreduce:', ierr
+           FLUSH(9)
+        ENDIF
+        surfinl(:) = surfinl(:) + surfinlg(surfstart(myid)+1:surfstart(myid+1))
+#else
+        surfinl(:) = surfinl(:) + surfinlg(:)
+#endif
+     ENDIF
+
      surfins = surfinswdir + surfinswdif
      surfinl = surfinl + surfinlwdif
@@ -4966,14 +5020,30 @@
              isurfsrc = svfsurf(2, isvf)
              surfins(isurf) = surfins(isurf) + svf(1,isvf) * svf(2,isvf) * surfouts(isurfsrc)
-             surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc)
+             IF ( plant_lw_interact )  THEN
+                surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * svf(2,isvf) * surfoutl(isurfsrc)
+             ELSE
+                surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc)
+             ENDIF
          ENDDO
-
-!--         radiation absorbed by plant canopy
-         DO icsf = 1, ncsfl
+!
+!--      NOTE: PC absorbtion and MRT from reflected can both be done at once
+!--      after all reflections if we do one more MPI_ALLGATHERV on surfout.
+!--      Advantage: less local computation. Disadvantage: one more collective
+!--      MPI call.
+!
+!--      Radiation absorbed by plant canopy
+         DO  icsf = 1, ncsfl
              ipcgb = csfsurf(1, icsf)
              isurfsrc = csfsurf(2, icsf)
              IF ( isurfsrc == -1 )  CYCLE ! sky->face only in 1st pass, not here
-
-             pcbinsw(ipcgb) = pcbinsw(ipcgb) + csf(1,icsf) * surfouts(isurfsrc)
+!
+!--          Calculate source surface area. If the `surf' array is removed
+!--          before timestepping starts (future version), then asrc must be
+!--          stored within `csf'
+             asrc = facearea(surf(id, isurfsrc))
+             pcbinsw(ipcgb) = pcbinsw(ipcgb) + csf(1,icsf) * surfouts(isurfsrc) * asrc
+             IF ( plant_lw_interact )  THEN
+                pcbinlw(ipcgb) = pcbinlw(ipcgb) + csf(1,icsf) * surfoutl(isurfsrc) * asrc
+             ENDIF
          ENDDO
 !
@@ -4996,20 +5066,28 @@
      IF ( npcbl > 0 )  THEN
          pc_heating_rate(:,:,:) = 0.0_wp
-         pc_transpiration_rate(:,:,:) = 0.0_wp
          DO ipcgb = 1, npcbl
-                 
              j = pcbl(iy, ipcgb)
              i = pcbl(ix, ipcgb)
              k = pcbl(iz, ipcgb)
 !
-!--             Following expression equals former kk = k - nzb_s_inner(j,i)
+!--          Following expression equals former kk = k - nzb_s_inner(j,i)
              kk = k - get_topography_top_index_ji( j, i, 's' )  !- lad arrays are defined flat
              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
+
+         IF ( humidity .AND. plant_canopy_transpiration ) THEN
+!--          Calculation of plant canopy transpiration rate and correspondidng latent heat rate
+             pc_transpiration_rate(:,:,:) = 0.0_wp
+             pc_latent_rate(:,:,:) = 0.0_wp
+             DO ipcgb = 1, npcbl
+                 i = pcbl(ix, ipcgb)
+                 j = pcbl(iy, ipcgb)
+                 k = pcbl(iz, ipcgb)
+                 kk = k - get_topography_top_index_ji( j, i, 's' )  !- lad arrays are defined flat
+                 CALL pcm_calc_transpiration_rate( i, j, k, kk, pcbinsw(ipcgb), pcbinlw(ipcgb), &
+                                                   pc_transpiration_rate(kk,j,i), pc_latent_rate(kk,j,i) )
+              ENDDO
+         ENDIF
      ENDIF
 !
@@ -5161,11 +5239,4 @@
 !--  domain when using average_radiation in the respective radiation model
 
-!--  Precalculate face areas for all face directions using normal vector
-     DO d = 0, nsurf_type
-        facearea(d) = 1._wp
-        IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx
-        IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy
-        IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz(1)
-     ENDDO
 !--  calculate horizontal area
 ! !!! ATTENTION!!! uniform grid is assumed here
@@ -5429,5 +5500,5 @@
        IMPLICIT NONE
 
-       INTEGER(iwp) :: i, j, k, l, m
+       INTEGER(iwp) :: i, j, k, l, m, d
        INTEGER(iwp) :: k_topo     !< vertical index indicating topography top for given (j,i)
        INTEGER(iwp) :: nzptl, nzubl, nzutl, isurf, ipcgb, imrt
@@ -5439,4 +5510,12 @@
 #endif
 
+!
+!--     precalculate face areas for different face directions using normal vector
+        DO d = 0, nsurf_type
+            facearea(d) = 1._wp
+            IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx
+            IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy
+            IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz(1)
+        ENDDO
 !
 !--    Find nzub, nzut, nzu via wall_flag_0 array (nzb_s_inner will be
@@ -5904,4 +5983,5 @@
        ALLOCATE( surfouts(nsurf) )
        ALLOCATE( surfoutl(nsurf) )
+       ALLOCATE( surfinlg(nsurf) )
        ALLOCATE( skyvf(nsurfl) )
        ALLOCATE( skyvft(nsurfl) )
@@ -5940,6 +6020,8 @@
         REAL(wp)                                      :: az1, az2      !< relative azimuth of section borders
         REAL(wp)                                      :: azmid         !< ray (center) azimuth
+        REAL(wp)                                      :: yxlen         !< |yxdir|
         REAL(wp), DIMENSION(2)                        :: yxdir         !< y,x *unit* vector of ray direction (in grid units)
         REAL(wp), DIMENSION(:), ALLOCATABLE           :: zdirs         !< directions in z (tangent of elevation)
+        REAL(wp), DIMENSION(:), ALLOCATABLE           :: zcent         !< zenith angle centers
         REAL(wp), DIMENSION(:), ALLOCATABLE           :: zbdry         !< zenith angle boundaries
         REAL(wp), DIMENSION(:), ALLOCATABLE           :: vffrac        !< view factor fractions for individual rays
@@ -5949,10 +6031,5 @@
         INTEGER(iwp), DIMENSION(:), ALLOCATABLE       :: itarget       !< face indices of detected obstacles
         INTEGER(iwp)                                  :: itarg0, itarg1
-#if defined( __parallel )
-#endif
-
-
-
-        REAL(wp),     DIMENSION(0:nsurf_type)         :: facearea
+
         INTEGER(iwp)                                  :: udim
         INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET:: nzterrl_l
@@ -5966,4 +6043,5 @@
         LOGICAL                                       :: visible
         REAL(wp), DIMENSION(3)                        :: sa, ta          !< real coordinates z,y,x of source and target
+        REAL(wp)                                      :: difvf           !< differential view factor
         REAL(wp)                                      :: transparency, rirrf, sqdist, svfsum
         INTEGER(iwp)                                  :: isurflt, isurfs, isurflt_prev
@@ -5983,12 +6061,4 @@
         CALL location_message( '    calculation of SVF and CSF', .TRUE. )
         CALL radiation_write_debug_log('Start calculation of SVF and CSF')
-
-!--     precalculate face areas for different face directions using normal vector
-        DO d = 0, nsurf_type
-            facearea(d) = 1._wp
-            IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx
-            IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy
-            IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz(1)
-        ENDDO
 
 !--     initialize variables and temporary arrays for calculation of svf and csf
@@ -6219,5 +6289,5 @@
             END SELECT
 
-            ALLOCATE ( zdirs(1:nzn), zbdry(0:nzn), vffrac(1:nzn*naz), &
+            ALLOCATE ( zdirs(1:nzn), zcent(1:nzn), zbdry(0:nzn), vffrac(1:nzn*naz), &
                        ztransp(1:nzn*naz), itarget(1:nzn*naz) )   !FIXME allocate itarget only
                                                                   !in case of rad_angular_discretization
@@ -6225,5 +6295,5 @@
             itarg0 = 1
             itarg1 = nzn
-            zdirs(:) = (/( TAN(pi/2 - (zn0+(REAL(izn,wp)-.5_wp)*zns)), izn=1, nzn )/)
+            zcent(:) = (/( zn0+(REAL(izn,wp)-.5_wp)*zns, izn=1, nzn )/)
             zbdry(:) = (/( zn0+REAL(izn,wp)*zns, izn=0, nzn )/)
             IF ( td == iup_u  .OR.  td == iup_l )  THEN
@@ -6251,5 +6321,9 @@
 !--               sum of whole vffrac equals 1, verified
                ENDIF
-               yxdir = (/ COS(azmid), SIN(azmid) /)
+               yxdir(:) = (/ COS(azmid) / dy, SIN(azmid) / dx /)
+               yxlen = SQRT(SUM(yxdir(:)**2))
+               zdirs(:) = COS(zcent(:)) / (dz(1) * yxlen * SIN(zcent(:)))
+               yxdir(:) = yxdir(:) / yxlen
+
                CALL raytrace_2d(ta, yxdir, nzn, zdirs,                        &
                                     surfstart(myid) + isurflt, facearea(td),  &
@@ -6257,7 +6331,6 @@
                                     .FALSE., lowest_free_ray,                 &
                                     ztransp(itarg0:itarg1),                   &
-                                    itarget(itarg0:itarg1))   !FIXME unit vect in grid units + zdirs
-                                                              !FIXME itarget available only in
-                                                              !case of rad_angular_discretization
+                                    itarget(itarg0:itarg1))
+
                skyvf(isurflt) = skyvf(isurflt) + &
                                 SUM(vffrac(itarg0:itarg0+lowest_free_ray-1))
@@ -6337,5 +6410,5 @@
             ENDIF ! rad_angular_discretization
 
-            DEALLOCATE ( zdirs, zbdry, vffrac, ztransp, itarget ) !FIXME itarget shall be allocated only
+            DEALLOCATE ( zdirs, zcent, zbdry, vffrac, ztransp, itarget ) !FIXME itarget shall be allocated only
                                                                   !in case of rad_angular_discretization
 !
@@ -6367,9 +6440,10 @@
                   ENDIF
                  
+                  difvf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & ! cosine of source normal and direction
+                      * dot_product((/ kdir(td), jdir(td), idir(td) /), -uv) &  ! cosine of target normal and reverse direction
+                      / (pi * sqdist) ! square of distance between centers
+!
 !--               irradiance factor (our unshaded shape view factor) = view factor per differential target area * source area
-                  rirrf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & ! cosine of source normal and direction
-                      * dot_product((/ kdir(td), jdir(td), idir(td) /), -uv) &  ! cosine of target normal and reverse direction
-                      / (pi * sqdist) & ! square of distance between centers
-                      * facearea(sd)
+                  rirrf = difvf * facearea(sd)
 
 !--               reject raytracing for potentially too small view factor values
@@ -6380,5 +6454,5 @@
 
 !--               raytrace + process plant canopy sinks within
-                  CALL raytrace(sa, ta, isurfs, rirrf, facearea(td), .TRUE., &
+                  CALL raytrace(sa, ta, isurfs, difvf, facearea(td), .TRUE., &
                                 visible, transparency)
 
@@ -6428,7 +6502,7 @@
         nzn = raytrace_discrete_elevs / 2
         zns = pi / 2._wp / REAL(nzn, wp)
-        ALLOCATE ( zdirs(1:nzn), vffrac(1:nzn), ztransp(1:nzn), &
+        ALLOCATE ( zdirs(1:nzn), zcent(1:nzn), vffrac(1:nzn), ztransp(1:nzn), &
                itarget(1:nzn) )
-        zdirs(:) = (/( TAN(pi/2 - (zn0+(REAL(izn,wp)-.5_wp)*zns)), izn=1, nzn )/)
+        zcent(:) = (/( zn0+(REAL(izn,wp)-.5_wp)*zns, izn=1, nzn )/)
         vffrac(:) = 0._wp
 
@@ -6440,8 +6514,11 @@
            DO  iaz = 1, naz
               azmid = az0 + (REAL(iaz, wp) - .5_wp) * azs
-              yxdir = (/ COS(azmid), SIN(azmid) /)
+              yxdir(:) = (/ COS(azmid) / dy, SIN(azmid) / dx /)
+              yxlen = SQRT(SUM(yxdir(:)**2))
+              zdirs(:) = COS(zcent(:)) / (dz(1) * yxlen * SIN(zcent(:)))
+              yxdir(:) = yxdir(:) / yxlen
               CALL raytrace_2d(ta, yxdir, nzn, zdirs,                                &
                                    -999, -999._wp, vffrac, .FALSE., .FALSE., .TRUE., &
-                                   lowest_free_ray, ztransp, itarget) !FIXME unit vect in grid units + zdirs
+                                   lowest_free_ray, ztransp, itarget)
 
 !--           Save direct solar transparency
@@ -6456,5 +6533,5 @@
            ENDDO
         ENDDO
-        DEALLOCATE ( zdirs, vffrac, ztransp, itarget )
+        DEALLOCATE ( zdirs, zcent, vffrac, ztransp, itarget )
 !--
 !--     Raytrace to MRT boxes
@@ -6469,9 +6546,9 @@
            nzn = raytrace_discrete_elevs
            zns = pi / REAL(nzn, wp)
-           ALLOCATE ( zdirs(1:nzn), zbdry(0:nzn), vffrac(1:nzn*naz), vffrac0(1:nzn), &
+           ALLOCATE ( zdirs(1:nzn), zcent(1:nzn), zbdry(0:nzn), vffrac(1:nzn*naz), vffrac0(1:nzn), &
                       ztransp(1:nzn*naz), itarget(1:nzn*naz) )   !FIXME allocate itarget only
                                                                  !in case of rad_angular_discretization
 
-           zdirs(:) = (/( TAN(pi/2 - (zn0+(REAL(izn,wp)-.5_wp)*zns)), izn=1, nzn )/)
+           zcent(:) = (/( zn0+(REAL(izn,wp)-.5_wp)*zns, izn=1, nzn )/)
            zbdry(:) = (/( zn0+REAL(izn,wp)*zns, izn=0, nzn )/)
            vffrac0(:) = (COS(zbdry(0:nzn-1)) - COS(zbdry(1:nzn))) / 2._wp / REAL(naz, wp)
@@ -6493,11 +6570,14 @@
               DO  iaz = 1, naz
                  azmid = az0 + (REAL(iaz, wp) - .5_wp) * azs
-                 CALL raytrace_2d(ta, (/ COS(azmid), SIN(azmid) /), nzn, zdirs,  &
+                 yxdir(:) = (/ COS(azmid) / dy, SIN(azmid) / dx /)
+                 yxlen = SQRT(SUM(yxdir(:)**2))
+                 zdirs(:) = COS(zcent(:)) / (dz(1) * yxlen * SIN(zcent(:)))
+                 yxdir(:) = yxdir(:) / yxlen
+
+                 CALL raytrace_2d(ta, yxdir, nzn, zdirs,                         &
                                   -999, -999._wp, vffrac(itarg0:itarg1), .TRUE., &
                                   .FALSE., .TRUE., lowest_free_ray,              &
                                   ztransp(itarg0:itarg1),                        &
-                                  itarget(itarg0:itarg1))   !FIXME unit vect in grid units + zdirs
-                                                            !FIXME itarget available only in
-                                                            !case of rad_angular_discretization
+                                  itarget(itarg0:itarg1))
 
 !--              Sky view factors for MRT
@@ -6574,5 +6654,5 @@
 
            ENDDO ! imrt
-           DEALLOCATE ( zdirs, zbdry, vffrac, vffrac0, ztransp, itarget )
+           DEALLOCATE ( zdirs, zcent, zbdry, vffrac, vffrac0, ztransp, itarget )
 !
 !--        Move MRT factors to final arrays
@@ -6778,6 +6858,5 @@
                     j = 0
                 ENDIF
-!--             fill out real values of rsvf, rtransp
-                csflt(1,kcsf) = acsf(kcsf)%rsvf
+                csflt(1,kcsf) = acsf(kcsf)%rcvf
 !--             fill out integer values of itz,ity,itx,isurfs
                 kcsflt(1,kcsf) = acsf(kcsf)%itz
@@ -6946,10 +7025,10 @@
 !>    doesn't need to be the same in all three dimensions).
 !------------------------------------------------------------------------------!
-    SUBROUTINE raytrace(src, targ, isrc, rirrf, atarg, create_csf, visible, transparency)
+    SUBROUTINE raytrace(src, targ, isrc, difvf, atarg, create_csf, visible, transparency)
         IMPLICIT NONE
 
         REAL(wp), DIMENSION(3), INTENT(in)     :: src, targ    !< real coordinates z,y,x
         INTEGER(iwp), INTENT(in)               :: isrc         !< index of source face for csf
-        REAL(wp), INTENT(in)                   :: rirrf        !< irradiance factor for csf
+        REAL(wp), INTENT(in)                   :: difvf        !< differential view factor for csf
         REAL(wp), INTENT(in)                   :: atarg        !< target surface area for csf
         LOGICAL, INTENT(in)                    :: create_csf   !< whether to generate new CSFs during raytracing
@@ -7113,5 +7192,5 @@
                     acsf(ncsfl)%itz = boxes(1,i)
                     acsf(ncsfl)%isurfs = isrc
-                    acsf(ncsfl)%rsvf = cursink*transparency*rirrf*atarg
+                    acsf(ncsfl)%rcvf = cursink*transparency*difvf*atarg
                 ENDIF  !< create_csf
 
@@ -7514,5 +7593,5 @@
                            acsf(ncsfl)%itz = iz
                            acsf(ncsfl)%isurfs = iorig
-                           acsf(ncsfl)%rsvf = (1._wp - curtrans)*transparency(k)*aorig*vffrac(k)
+                           acsf(ncsfl)%rcvf = (1._wp - curtrans)*transparency(k)*vffrac(k)
                         ENDIF
 
@@ -7571,7 +7650,7 @@
                         acsf(ncsfl)%ity = rt2_track(1,i)
                         acsf(ncsfl)%itz = iz
-                        acsf(ncsfl)%isurfs = itarget(k) !if above horizon, then itarget(k)==-1, which
-                                                        !is also a special ID indicating sky
-                        acsf(ncsfl)%rsvf = (1._wp - curtrans)*transparency(k)*aorig*vffrac(k)
+                        IF ( itarget(k) /= -1 )  ERROR STOP !FIXME remove after test
+                        acsf(ncsfl)%isurfs = -1
+                        acsf(ncsfl)%rcvf = (1._wp - curtrans)*transparency(k)*aorig*vffrac(k)
                      ENDIF  ! create_csf
 
@@ -7600,4 +7679,8 @@
          INTEGER(iwp), TARGET, INTENT(out)   ::  isurfl
          INTEGER(iwp), INTENT(out)           ::  iproc
+#if defined( __parallel )
+#else
+         INTEGER(iwp)                        ::  target_displ  !< index of the grid in the local gridsurf array
+#endif
          INTEGER(iwp)                        ::  px, py        !< number of processors in x and y direction
                                                                !< before the processor in the question
@@ -8138,5 +8221,5 @@
                          .AND.  acsfnew(iwrite)%isurfs == acsf(iread)%isurfs )  THEN
 
-                    acsfnew(iwrite)%rsvf = acsfnew(iwrite)%rsvf + acsf(iread)%rsvf
+                    acsfnew(iwrite)%rcvf = acsfnew(iwrite)%rcvf + acsf(iread)%rcvf
 !--                 advance reading index, keep writing index
                 ELSE
Index: palm/trunk/SOURCE/urban_surface_mod.f90
===================================================================
--- palm/trunk/SOURCE/urban_surface_mod.f90	(revision 3448)
+++ palm/trunk/SOURCE/urban_surface_mod.f90	(revision 3449)
@@ -28,4 +28,10 @@
 ! -----------------
 ! $Id$
+! Bugfix: Fix average arrays allocations in usm_average_3d_data (J.Resler)
+! Bugfix: Fix reading wall temperatures (J.Resler)
+! Bugfix: Fix treating of outputs for wall temperature and sky view factors (J.Resler)
+!
+! 
+! 3435 2018-10-26 18:25:44Z gronemeier
 ! Bugfix: allocate gamma_w_green_sat until nzt_wall+1
 ! 
@@ -411,5 +417,5 @@
     
     USE plant_canopy_model_mod,                                                &
-        ONLY:  pc_heating_rate, pc_transpiration_rate
+        ONLY:  pc_heating_rate, pc_transpiration_rate, pc_latent_rate
     
     USE radiation_model_mod,                                                   &
@@ -1845,7 +1851,9 @@
                 CASE ( 'usm_rad_net' )
 !--                 array of complete radiation balance
-                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%rad_net_av) )  THEN
-                        ALLOCATE( surf_usm_h%rad_net_av(1:surf_usm_h%ns) )
-                        surf_usm_h%rad_net_av = 0.0_wp
+                    IF ( l == -1 ) THEN
+                        IF ( .NOT.  ALLOCATED(surf_usm_h%rad_net_av) )  THEN
+                            ALLOCATE( surf_usm_h%rad_net_av(1:surf_usm_h%ns) )
+                            surf_usm_h%rad_net_av = 0.0_wp
+                        ENDIF
                     ELSE
                         IF ( .NOT.  ALLOCATED(surf_usm_v(l)%rad_net_av) )  THEN
@@ -1968,7 +1976,9 @@
                 CASE ( 'usm_rad_hf' )
 !--                 array of heat flux from radiation for surfaces after i-th reflection
-                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%surfhf_av) )  THEN
-                        ALLOCATE( surf_usm_h%surfhf_av(1:surf_usm_h%ns) )
-                        surf_usm_h%surfhf_av = 0.0_wp
+                    IF ( l == -1 ) THEN
+                        IF ( .NOT.  ALLOCATED(surf_usm_h%surfhf_av) )  THEN
+                           ALLOCATE( surf_usm_h%surfhf_av(1:surf_usm_h%ns) )
+                           surf_usm_h%surfhf_av = 0.0_wp
+                        ENDIF
                     ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%surfhf_av) )  THEN
@@ -1981,7 +1991,9 @@
 !--                 array of sensible heat flux from surfaces
 !--                 land surfaces
-                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%wshf_eb_av) )  THEN
-                        ALLOCATE( surf_usm_h%wshf_eb_av(1:surf_usm_h%ns) )
-                        surf_usm_h%wshf_eb_av = 0.0_wp
+                    IF ( l == -1 ) THEN
+                       IF ( .NOT.  ALLOCATED(surf_usm_h%wshf_eb_av) )  THEN
+                          ALLOCATE( surf_usm_h%wshf_eb_av(1:surf_usm_h%ns) )
+                          surf_usm_h%wshf_eb_av = 0.0_wp
+                       ENDIF
                     ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wshf_eb_av) )  THEN
@@ -2036,7 +2048,9 @@
                 CASE ( 'usm_wghf' )
 !--                 array of heat flux from ground (wall, roof, land)
-                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%wghf_eb_av) )  THEN
-                        ALLOCATE( surf_usm_h%wghf_eb_av(1:surf_usm_h%ns) )
-                        surf_usm_h%wghf_eb_av = 0.0_wp
+                    IF ( l == -1 ) THEN
+                       IF ( .NOT.  ALLOCATED(surf_usm_h%wghf_eb_av) )  THEN
+                           ALLOCATE( surf_usm_h%wghf_eb_av(1:surf_usm_h%ns) )
+                           surf_usm_h%wghf_eb_av = 0.0_wp
+                       ENDIF
                     ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wghf_eb_av) )  THEN
@@ -2048,7 +2062,9 @@
                 CASE ( 'usm_wghf_window' )
 !--                 array of heat flux from window ground (wall, roof, land)
-                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%wghf_eb_window_av) )  THEN
-                        ALLOCATE( surf_usm_h%wghf_eb_window_av(1:surf_usm_h%ns) )
-                        surf_usm_h%wghf_eb_window_av = 0.0_wp
+                    IF ( l == -1 ) THEN
+                       IF ( .NOT.  ALLOCATED(surf_usm_h%wghf_eb_window_av) )  THEN
+                           ALLOCATE( surf_usm_h%wghf_eb_window_av(1:surf_usm_h%ns) )
+                           surf_usm_h%wghf_eb_window_av = 0.0_wp
+                       ENDIF
                     ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wghf_eb_window_av) )  THEN
@@ -2060,7 +2076,9 @@
                 CASE ( 'usm_wghf_green' )
 !--                 array of heat flux from green ground (wall, roof, land)
-                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%wghf_eb_green_av) )  THEN
-                        ALLOCATE( surf_usm_h%wghf_eb_green_av(1:surf_usm_h%ns) )
-                        surf_usm_h%wghf_eb_green_av = 0.0_wp
+                    IF ( l == -1 ) THEN
+                       IF ( .NOT.  ALLOCATED(surf_usm_h%wghf_eb_green_av) )  THEN
+                           ALLOCATE( surf_usm_h%wghf_eb_green_av(1:surf_usm_h%ns) )
+                           surf_usm_h%wghf_eb_green_av = 0.0_wp
+                       ENDIF
                     ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wghf_eb_green_av) )  THEN
@@ -2072,7 +2090,9 @@
                 CASE ( 'usm_iwghf' )
 !--                 array of heat flux from indoor ground (wall, roof, land)
-                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%iwghf_eb_av) )  THEN
-                        ALLOCATE( surf_usm_h%iwghf_eb_av(1:surf_usm_h%ns) )
-                        surf_usm_h%iwghf_eb_av = 0.0_wp
+                    IF ( l == -1 ) THEN
+                       IF ( .NOT.  ALLOCATED(surf_usm_h%iwghf_eb_av) )  THEN
+                           ALLOCATE( surf_usm_h%iwghf_eb_av(1:surf_usm_h%ns) )
+                           surf_usm_h%iwghf_eb_av = 0.0_wp
+                       ENDIF
                     ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%iwghf_eb_av) )  THEN
@@ -2084,7 +2104,9 @@
                 CASE ( 'usm_iwghf_window' )
 !--                 array of heat flux from indoor window ground (wall, roof, land)
-                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%iwghf_eb_window_av) )  THEN
-                        ALLOCATE( surf_usm_h%iwghf_eb_window_av(1:surf_usm_h%ns) )
-                        surf_usm_h%iwghf_eb_window_av = 0.0_wp
+                    IF ( l == -1 ) THEN
+                       IF ( .NOT.  ALLOCATED(surf_usm_h%iwghf_eb_window_av) )  THEN
+                           ALLOCATE( surf_usm_h%iwghf_eb_window_av(1:surf_usm_h%ns) )
+                           surf_usm_h%iwghf_eb_window_av = 0.0_wp
+                       ENDIF
                     ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%iwghf_eb_window_av) )  THEN
@@ -2093,10 +2115,12 @@
                        ENDIF
                     ENDIF
-                    
+
                 CASE ( 'usm_t_surf_wall' )
 !--                 surface temperature for surfaces
-                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%t_surf_wall_av) )  THEN
-                        ALLOCATE( surf_usm_h%t_surf_wall_av(1:surf_usm_h%ns) )
-                        surf_usm_h%t_surf_wall_av = 0.0_wp
+                    IF ( l == -1 ) THEN
+                       IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_wall_av) )  THEN
+                           ALLOCATE( surf_usm_h%t_surf_wall_av(1:surf_usm_h%ns) )
+                           surf_usm_h%t_surf_wall_av = 0.0_wp
+                       ENDIF
                     ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_wall_av) )  THEN
@@ -2108,7 +2132,9 @@
                 CASE ( 'usm_t_surf_window' )
 !--                 surface temperature for window surfaces
-                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%t_surf_window_av) )  THEN
-                        ALLOCATE( surf_usm_h%t_surf_window_av(1:surf_usm_h%ns) )
-                        surf_usm_h%t_surf_window_av = 0.0_wp
+                    IF ( l == -1 ) THEN
+                       IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_window_av) )  THEN
+                           ALLOCATE( surf_usm_h%t_surf_window_av(1:surf_usm_h%ns) )
+                           surf_usm_h%t_surf_window_av = 0.0_wp
+                       ENDIF
                     ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_window_av) )  THEN
@@ -2120,7 +2146,9 @@
                 CASE ( 'usm_t_surf_green' )
 !--                 surface temperature for green surfaces
-                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%t_surf_green_av) )  THEN
-                        ALLOCATE( surf_usm_h%t_surf_green_av(1:surf_usm_h%ns) )
-                        surf_usm_h%t_surf_green_av = 0.0_wp
+                    IF ( l == -1 ) THEN
+                       IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_green_av) )  THEN
+                           ALLOCATE( surf_usm_h%t_surf_green_av(1:surf_usm_h%ns) )
+                           surf_usm_h%t_surf_green_av = 0.0_wp
+                       ENDIF
                     ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_green_av) )  THEN
@@ -2132,7 +2160,9 @@
                 CASE ( 'usm_t_surf_10cm' )
 !--                 near surface temperature for whole surfaces
-                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%t_surf_10cm_av) )  THEN
-                        ALLOCATE( surf_usm_h%t_surf_10cm_av(1:surf_usm_h%ns) )
-                        surf_usm_h%t_surf_10cm_av = 0.0_wp
+                    IF ( l == -1 ) THEN
+                       IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_10cm_av) )  THEN
+                           ALLOCATE( surf_usm_h%t_surf_10cm_av(1:surf_usm_h%ns) )
+                           surf_usm_h%t_surf_10cm_av = 0.0_wp
+                       ENDIF
                     ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_10cm_av) )  THEN
@@ -2144,7 +2174,9 @@
                 CASE ( 'usm_t_wall' )
 !--                 wall temperature for iwl layer of walls and land
-                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%t_wall_av) )  THEN
-                        ALLOCATE( surf_usm_h%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
-                        surf_usm_h%t_wall_av = 0.0_wp
+                    IF ( l == -1 ) THEN
+                       IF ( .NOT.  ALLOCATED(surf_usm_h%t_wall_av) )  THEN
+                           ALLOCATE( surf_usm_h%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
+                           surf_usm_h%t_wall_av = 0.0_wp
+                       ENDIF
                     ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_wall_av) )  THEN
@@ -2156,7 +2188,9 @@
                 CASE ( 'usm_t_window' )
 !--                 window temperature for iwl layer of walls and land
-                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%t_window_av) )  THEN
-                        ALLOCATE( surf_usm_h%t_window_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
-                        surf_usm_h%t_window_av = 0.0_wp
+                    IF ( l == -1 ) THEN
+                       IF ( .NOT.  ALLOCATED(surf_usm_h%t_window_av) )  THEN
+                           ALLOCATE( surf_usm_h%t_window_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
+                           surf_usm_h%t_window_av = 0.0_wp
+                       ENDIF
                     ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_window_av) )  THEN
@@ -2168,7 +2202,9 @@
                 CASE ( 'usm_t_green' )
 !--                 green temperature for iwl layer of walls and land
-                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%t_green_av) )  THEN
-                        ALLOCATE( surf_usm_h%t_green_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
-                        surf_usm_h%t_green_av = 0.0_wp
+                    IF ( l == -1 ) THEN
+                       IF ( .NOT.  ALLOCATED(surf_usm_h%t_green_av) )  THEN
+                           ALLOCATE( surf_usm_h%t_green_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
+                           surf_usm_h%t_green_av = 0.0_wp
+                       ENDIF
                     ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_green_av) )  THEN
@@ -3091,7 +3127,8 @@
         CHARACTER(LEN=*),INTENT(OUT)   ::  unit       !<
 
-        INTEGER(iwp)                                  :: i,j          !< index
+        INTEGER(iwp)                                  :: i,j,l        !< index
+        CHARACTER(LEN=2)                              :: ls
         CHARACTER(LEN=varnamelength)                  :: var          !< TRIM(variable)
-        INTEGER(iwp), PARAMETER                       :: nl1 = 32     !< number of directional usm variables
+        INTEGER(iwp), PARAMETER                       :: nl1 = 31     !< number of directional usm variables
         CHARACTER(LEN=varnamelength), DIMENSION(nl1)  :: varlist1 = & !< list of directional usm variables
                   (/'usm_rad_net                   ', &
@@ -3119,18 +3156,15 @@
                     'usm_surfalb                   ', &
                     'usm_surfemis                  ', &
-                    'usm_t_surf                    ', &
+                    'usm_t_surf_wall               ', &
                     'usm_t_surf_window             ', &
                     'usm_t_surf_green              ', &
                     'usm_t_green                   ', &
                     'usm_t_surf_10cm               ', &
-                    'usm_t_wall                    ', &
-                    'usm_t_window                  ', &
-                    'usm_t_green                   '/)
-
-        INTEGER(iwp), PARAMETER                       :: nl2 = 9      !< number of other variables
+                    'usm_skyvf                     ', &
+                    'usm_skyvft                    '/)
+
+        INTEGER(iwp), PARAMETER                       :: nl2 = 7      !< number of other variables
         CHARACTER(LEN=varnamelength), DIMENSION(nl2)  :: varlist2 = & !< list of other usm variables
-                  (/'usm_skyvf                     ', &
-                    'usm_skyvft                    ', &
-                    'usm_svf                       ', &
+                  (/'usm_svf                       ', &
                     'usm_dif                       ', &
                     'usm_rad_pc_inlw               ', &
@@ -3139,4 +3173,10 @@
                     'usm_rad_pc_inswdif            ', &
                     'usm_rad_pc_inswref            '/)
+
+        INTEGER(iwp), PARAMETER                       :: nl3 = 3      !< number of directional layer usm variables
+        CHARACTER(LEN=varnamelength), DIMENSION(nl3)  :: varlist3 = & !< list of directional layer usm variables
+                  (/'usm_t_wall                    ', &
+                    'usm_t_window                  ', &
+                    'usm_t_green                   '/)
 
         INTEGER(iwp), PARAMETER                       :: nd = 5     !< number of directions
@@ -3158,4 +3198,18 @@
                  EXIT
               ENDIF
+              IF ( lfound ) EXIT
+           ENDDO
+        ENDDO
+        IF ( lfound ) GOTO 10
+!       directional layer variables
+        DO i = 1, nl3
+           DO j = 1, nd
+              DO l = nzb_wall, nzt_wall
+                 WRITE(ls,'(A1,I1)') '_',l
+                 IF ( TRIM(var) == TRIM(varlist3(i))//TRIM(ls)//TRIM(dirname(j)) ) THEN
+                    lfound = .TRUE.
+                    EXIT
+                 ENDIF
+              ENDDO
               IF ( lfound ) EXIT
            ENDDO
@@ -4503,5 +4557,5 @@
              var(1:9) == 'usm_qsws_'  .OR.  var(1:13) == 'usm_qsws_veg_'  .OR.              &
              var(1:13) == 'usm_qsws_liq_' .OR.                                              &
-             var(1:10) == 'usm_t_surf_wall'  .OR.  var(1:10) == 'usm_t_wall'  .OR.          &
+             var(1:15) == 'usm_t_surf_wall'  .OR.  var(1:10) == 'usm_t_wall'  .OR.          &
              var(1:17) == 'usm_t_surf_window'  .OR.  var(1:12) == 'usm_t_window'  .OR.      &
              var(1:16) == 'usm_t_surf_green'  .OR. var(1:11) == 'usm_t_green' .OR.          &
@@ -5910,25 +5964,4 @@
         ENDIF
 
-        IF ( plant_canopy )  THEN
-            
-            IF ( .NOT.  ALLOCATED( pc_heating_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_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
 !
 !--    Check for consistent initialization.
@@ -6092,11 +6125,10 @@
                ENDDO
             ENDDO
-        ELSE
-!--         If specified, replace constant wall temperatures with fully 3D values from file
-            IF ( read_wall_temp_3d )  CALL usm_read_wall_temperature()
-!
         ENDIF
-        
-!--    
+
+!--     If specified, replace constant wall temperatures with fully 3D values from file
+        IF ( read_wall_temp_3d )  CALL usm_read_wall_temperature()
+
+!--
 !--     Possibly DO user-defined actions (e.g. define heterogeneous wall surface)
         CALL user_init_urban_surface
@@ -8870,19 +8902,25 @@
 !> x, y, d, z coordinates
 !------------------------------------------------------------------------------!
-    PURE FUNCTION advance_surface(isurfl_start, isurfl_stop, x, y, z, d) &
-            result(isurfl)
-
-        INTEGER(iwp), INTENT(in)                :: isurfl_start, isurfl_stop
+    PURE FUNCTION find_surface( x, y, z, d ) result(isurfl)
+
         INTEGER(iwp), INTENT(in)                :: x, y, z, d
-        INTEGER(iwp)                            :: isx, isy, isz, isd
         INTEGER(iwp)                            :: isurfl
-
-        DO isurfl = isurfl_start, isurfl_stop
-            isx = surfl(ix, isurfl)
-            isy = surfl(iy, isurfl)
-            isz = surfl(iz, isurfl)
-            isd = surfl(id, isurfl)
-            IF ( isx==x .and. isy==y .and. isz==z .and. isd==d )  RETURN
-        ENDDO
+        INTEGER(iwp)                            :: isx, isy, isz
+
+        IF ( d == 0 ) THEN
+           DO  isurfl = 1, surf_usm_h%ns
+              isx = surf_usm_h%i(isurfl)
+              isy = surf_usm_h%j(isurfl)
+              isz = surf_usm_h%k(isurfl)
+              IF ( isx==x .and. isy==y .and. isz==z )  RETURN
+           ENDDO
+        ELSE
+           DO  isurfl = 1, surf_usm_v(d-1)%ns
+              isx = surf_usm_v(d-1)%i(isurfl)
+              isy = surf_usm_v(d-1)%j(isurfl)
+              isz = surf_usm_v(d-1)%k(isurfl)
+              IF ( isx==x .and. isy==y .and. isz==z )  RETURN
+           ENDDO
+        ENDIF
 
 !--     coordinate not found
@@ -8927,5 +8965,5 @@
                         nys <= j .and. j <= nyn)  THEN  !< local processor
 !--                     identify surface id
-                        isurfl = advance_surface(isurfl+1, nsurfl, i, j, k, d)
+                        isurfl = find_surface( i, j, k, d )
                         IF ( isurfl == -1 )  THEN
                             WRITE(message_string, '(a,4i5,a,i5,a)') 'Coordinates (xyzd) ', i, j, k, d, &
@@ -8939,7 +8977,11 @@
                            t_surf_wall_h(isurfl) = rtsurf
                            t_wall_h(:,isurfl) = rtwall(:)
+                           t_window_h(:,isurfl) = rtwall(:)
+                           t_green_h(:,isurfl) = rtwall(:)
                         ELSE
                            t_surf_wall_v(d-1)%t(isurfl) = rtsurf
                            t_wall_v(d-1)%t(:,isurfl) = rtwall(:)
+                           t_window_v(d-1)%t(:,isurfl) = rtwall(:)
+                           t_green_v(d-1)%t(:,isurfl) = rtwall(:)
                         ENDIF
                     ENDIF
@@ -9663,14 +9705,4 @@
                                           surf_usm_v(l)%rad_lw_out(m)
 
-if ((abs(t_surf_wall_v(l)%t(m)-276.).gt.100.).or.(abs(t_surf_window_v(l)%t(m)-276.).gt.300.) &
-.or.(abs(t_surf_green_v(l)%t(m)-276.).gt.100)) then
-print*, "tsurfvvvv",m,t_surf_wall_v(l)%t(m),t_surf_window_v(l)%t(m),t_surf_green_v(l)%t(m)
-print*, "params", surf_usm_v(l)%emissivity(ind_veg_wall,m),lambda_surface, t_wall_v(l)%t(nzb_wall:nzt_wall,m), &
-                  surf_usm_v(l)%emissivity(ind_wat_win,m),lambda_surface_window,t_window_v(l)%t(nzb_wall:nzt_wall,m),t_green_v(l)%t(nzb_wall:nzt_wall,m)
-print*, "dicken",surf_usm_v(l)%zw(nzb_wall:nzt_wall,m),surf_usm_v(l)%zw_window(nzb_wall:nzt_wall,m),surf_usm_v(l)%zw_green(nzb_wall:nzt_wall,m)
-print*, "c",surf_usm_v(l)%c_surface_window(m),surf_usm_v(l)%c_surface(m)
-!if ((abs(t_surf_v(l)%t(m)-276.).gt.10.).or.(abs(t_surf_window_v(l)%t(m)-276.).gt.10.)) then
-stop
-endif
 !--           numerator of the prognostic equation
               coef_1 = surf_usm_v(l)%rad_net_l(m) +                            & ! coef +1 corresponds to -lwout included in calculation of radnet_l
