Changeset 3337 for palm

Timestamp:

Oct 12, 2018 3:17:09 PM (6 years ago)

Author:

kanani

Message:

reintegrate branch resler to trunk

Location:

palm/trunk/SOURCE

Files:

• . (modified) (1 prop)
• Makefile (modified) (8 diffs, 1 prop)
• average_3d_data.f90 (modified) (5 diffs)
• biometeorology_mod.f90 (copied) (copied from palm/branches/resler/SOURCE/biometeorology_mod.f90)
• check_parameters.f90 (modified) (5 diffs, 1 prop)
• chem_emissions_mod.f90 (modified) (17 diffs)
• data_output_3d.f90 (modified) (3 diffs)
• modules.f90 (modified) (2 diffs)
• netcdf_data_input_mod.f90 (modified) (7 diffs)
• netcdf_interface_mod.f90 (modified) (7 diffs)
• palm.f90 (modified) (2 diffs)
• plant_canopy_model_mod.f90 (modified) (4 diffs)
• prognostic_equations.f90 (modified) (2 diffs)
• radiation_model_mod.f90 (modified) (125 diffs)
• sum_up_3d_data.f90 (modified) (8 diffs)
• time_integration_spinup.f90 (modified) (3 diffs, 1 prop)
• urban_surface_mod.f90 (modified) (70 diffs)

palm/trunk/SOURCE/Makefile

@@ -25 +25 @@
 # -----------------
 # $Id$
+# (from branch resler)
+# Add biometeorology
+# 
+# 
+# 3322 2018-10-09 10:02:39Z kanani
 # Formatting and cleanup
 # 
@@ -514 +519 @@
         gust_mod.f90 \
         header.f90 \
+        biometeorology_mod.f90 \
         inflow_turbulence.f90 \
         init_3d_model.f90 \
@@ -712 +718 @@
         modules.o
 average_3d_data.o: \
+        biometeorology_mod.o \
         bulk_cloud_model_mod.o \
         chemistry_model_mod.o \
@@ -726 +733 @@
 basic_constants_and_equations_mod.o: \
         mod_kinds.o
+biometeorology_mod.o: \
+        mod_kinds.o \
+        radiation_model_mod.o
 boundary_conds.o: \
         basic_constants_and_equations_mod.o \
@@ -764 +774 @@
 check_parameters.o: \
         basic_constants_and_equations_mod.o \
+        biometeorology_mod.o \
         bulk_cloud_model_mod.o \
         chemistry_model_mod.o \
@@ -904 +915 @@
 data_output_3d.o: \
         basic_constants_and_equations_mod.o \
+        biometeorology_mod.o \
         bulk_cloud_model_mod.o \
         chemistry_model_mod.o \
@@ -1264 +1276 @@
 netcdf_interface_mod.o: \
         basic_constants_and_equations_mod.o \
+        biometeorology_mod.o \
         chemistry_model_mod.o \
         gust_mod.o \
@@ -1513 +1526 @@
 sum_up_3d_data.o: \
         basic_constants_and_equations_mod.o \
+        biometeorology_mod.o \
         bulk_cloud_model_mod.o \
         chemistry_model_mod.o \

palm/trunk/SOURCE/average_3d_data.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! (from resler branch)
+! Add biometeorology output,
+! fix for chemistry output,
+! move of urban surface output
+! 
+! 3294 2018-10-01 02:37:10Z raasch
 ! changes concerning modularization of ocean option
 ! 
@@ -174 +180 @@
         ONLY:  radiation, radiation_3d_data_averaging
 
+    USE biometeorology_mod,                                                    &
+        ONLY:  biometeorology_3d_data_averaging
+
     USE turbulence_closure_mod,                                                &
         ONLY:  tcm_3d_data_averaging
@@ -203 +212 @@
     DO  ii = 1, doav_n
 
-!
-!--    Temporary solution to account for data output within the new urban 
-!--    surface model (urban_surface_mod.f90), see also SELECT CASE ( trimvar )
        trimvar = TRIM( doav(ii) )
-       IF ( urban_surface  .AND.  trimvar(1:4) == 'usm_' )  THEN
-          trimvar = 'usm_output'
-       ENDIF
 
 !
@@ -533 +536 @@
              ENDIF
 
-          CASE ( 'usm_output' )
-!             
-!--          Block of urban surface model outputs
-             CALL usm_average_3d_data( 'average', doav(ii) )
-
           CASE DEFAULT
 !
 !--          Averaging of data from other modules
-             IF ( air_chemistry )  THEN
+             IF ( air_chemistry .AND. &
+                  (trimvar(1:3) == 'kc_' .OR. trimvar(1:3) == 'em_')  )  THEN
                 CALL chem_3d_data_averaging( 'average', doav(ii) )
              ENDIF
@@ -565 +564 @@
              ENDIF
 
+             IF ( radiation  .AND.  trimvar(1:4) == 'bio_' )  THEN
+                CALL biometeorology_3d_data_averaging( 'average', doav(ii) )
+             ENDIF
+
              CALL tcm_3d_data_averaging( 'average', doav(ii) )
+
+             IF ( urban_surface  .AND.  trimvar(1:4) == 'usm_' )  THEN
+                CALL usm_average_3d_data( 'average', doav(ii) )
+             ENDIF
+
 !
 !--          User-defined quantities

palm/trunk/SOURCE/check_parameters.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! (from branch resler)
+! Add biometeorology,
+! fix for chemistry output,
+! increase of uv_heights dimension
+! 
+! 3312 2018-10-06 14:15:46Z knoop
 ! small code rearrangements
 ! 
@@ -756 +762 @@
         ONLY:  radiation, radiation_check_data_output,                         &
                radiation_check_data_output_pr, radiation_check_parameters
+
+    USE biometeorology_mod,                                                    &
+        ONLY:  biometeorology_check_data_output
 
     USE spectra_mod,                                                           &
@@ -1662 +1671 @@
           DO  k = 1, nz+1
 
-             IF ( kk < 100 )  THEN
+             IF ( kk < 200 )  THEN
                 DO  WHILE ( uv_heights(kk+1) <= zu(k) )
                    kk = kk + 1
-                   IF ( kk == 100 )  EXIT
+                   IF ( kk == 200 )  EXIT
                 ENDDO
              ENDIF
 
-             IF ( kk < 100  .AND.  uv_heights(kk+1) /= 9999999.9_wp )  THEN
+             IF ( kk < 200  .AND.  uv_heights(kk+1) /= 9999999.9_wp )  THEN
                 u_init(k) = u_profile(kk) + ( zu(k) - uv_heights(kk) ) /       &
                                        ( uv_heights(kk+1) - uv_heights(kk) ) * &
@@ -3156 +3165 @@
 
              IF ( unit == 'illegal'  .AND.  air_chemistry                      &
-                                     .AND.  var(1:3) == 'kc_' )  THEN
+                  .AND.  (var(1:3) == 'kc_' .OR. var(1:3) == 'em_') )  THEN
                 CALL chem_check_data_output( var, unit, i, ilen, k )
              ENDIF
@@ -3162 +3171 @@
              IF ( unit == 'illegal' )  THEN
                 CALL lsm_check_data_output ( var, unit, i, ilen, k )
+             ENDIF
+
+             IF ( unit == 'illegal' )  THEN
+                CALL biometeorology_check_data_output( var, unit, i, ilen, k )
              ENDIF
 

palm/trunk/SOURCE/chem_emissions_mod.f90

@@ -27 +27 @@
 ! -----------------
 ! $Id$
+! (from branch resler)
+! Formatting
+! 
+! 3312 2018-10-06 14:15:46Z knoop
 ! Initial revision
 !
@@ -719 +723 @@
                 DO  ispec = 1 , len_index
 
-                   IF ( emiss_factor_main(match_spec_input(ispec)) .LT. 0 .AND. emiss_factor_side(match_spec_input(ispec)) .LT. 0 ) THEN
-
-                      message_string = 'PARAMETERIZED emissions mode selected:'            //          &
+                   IF ( emiss_factor_main(match_spec_input(ispec)) .LT. 0 .AND.                         &
+                        emiss_factor_side(match_spec_input(ispec)) .LT. 0 ) THEN
+
+                      message_string = 'PARAMETERIZED emissions mode selected:'            //           &
                                        ' EMISSIONS POSSIBLE ONLY ON STREET SURFACES'        //          &
                                        ' but values of scaling factors for street types'    //          &
@@ -1246 +1251 @@
                 IF (TRIM(spc_names(match_spec_model(ispec)))=="NO") THEN
                 !>             Kg/m2                       kg/m2*s                                                    
-                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%nox_comp(icat,1)*con_factor
+                   delta_emis(nys:nyn,nxl:nxr) =                                                             &
+                         emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%nox_comp(icat,1)*con_factor
                    
-                   emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
+                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                              &
+                         emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
 
                 ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="NO2") THEN
    
-                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%nox_comp(icat,2)*con_factor
-
-                   emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
+                   delta_emis(nys:nyn,nxl:nxr) =                                                             &
+                         emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%nox_comp(icat,2)*con_factor
+
+                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                              &
+                         emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
  
              !> SOX Compositions
@@ -1260 +1269 @@
                 ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="SO2") THEN
                       
-                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,1)*con_factor
-
-                   emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
+                   delta_emis(nys:nyn,nxl:nxr) =                                                             &
+                         emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,1)*con_factor
+
+                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                              &
+                         emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
 
                 ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="SO4") THEN
    
-                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,2)*con_factor
-
-                   emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
+                   delta_emis(nys:nyn,nxl:nxr) =                                                             &
+                         emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,2)*con_factor
+
+                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                              &
+                         emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
  
 
@@ -1278 +1291 @@
                    DO i_pm_comp= 1,SIZE(emt_att%pm_comp(1,:,1))
 
-                      delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,1)*con_factor 
+                      delta_emis(nys:nyn,nxl:nxr) =                                                          &
+                            emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,1)*con_factor
                                                                                          
 
-                      emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
+                      emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                           &
+                            emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
  
                    ENDDO
@@ -1291 +1306 @@
                    DO i_pm_comp= 1,SIZE(emt_att%pm_comp(1,:,2))
 
-                      delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,2)*con_factor 
+                      delta_emis(nys:nyn,nxl:nxr) =                                                          &
+                            emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,2)*con_factor
                                                                                          
 
-                      emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
+                      emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                           &
+                            emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
  
                    ENDDO
@@ -1304 +1321 @@
                    DO i_pm_comp= 1,SIZE(emt_att%pm_comp(1,:,3))  
                        
-                      delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,3)*con_factor 
+                      delta_emis(nys:nyn,nxl:nxr) =                                                          &
+                            emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,3)*con_factor
                                                                                                  
 
-                      emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 
+                      emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                           &
+                            emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
 
                    ENDDO
@@ -1319 +1338 @@
                       IF (TRIM(spc_names(match_spec_model(ispec)))==TRIM(emt_att%voc_name(ivoc))) THEN    
 
-                         delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*                               &
+                         delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*               &
                                                        emt_att%voc_comp(icat,match_spec_voc_input(ivoc))*con_factor
 
-                         emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 
+                         emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                         &
+                               emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
 
                       ENDIF                       
@@ -1333 +1353 @@
                    delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*con_factor
  
-                   emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr) 
+                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                               &
+                         emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
 
                 ENDIF  ! IF (spc_names==)
@@ -1375 +1396 @@
 
                          !> PMs are already in mass units:micrograms:they have to be converted to kilograms 
-                         IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"  &
+                         IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1"         &
+                             .OR.  TRIM(spc_names(match_spec_model(ispec)))=="PM25"  &
                              .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
 
@@ -1389 +1411 @@
                          !> Other Species: inputs are micromoles: they have to be converted in moles
                          !                 ppm/s *m *kg/m^3                
-                            surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_main(match_spec_input(ispec))*   &
+                            surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_main(match_spec_input(ispec))*  &
                          !                                                    micromoles/(m^2*s)
-                                                                          emis_distribution(1,j,i,ispec) *              &
+                                                                          emis_distribution(1,j,i,ispec) *             &
                          !                                                    m**3/Nmole
-                                                                          conv_to_ratio(k,j,i)*                         &       
+                                                                          conv_to_ratio(k,j,i)*                        &
                          !                                                    kg/m^3
                                                                           rho_air(k)   
@@ -1403 +1425 @@
 
 
-                   ELSEIF ( street_type_f%var(j,i) >= side_street_id  .AND. street_type_f%var(j,i) < main_street_id )  THEN
+                   ELSEIF ( street_type_f%var(j,i) >= side_street_id  .AND.                                            &
+                            street_type_f%var(j,i) < main_street_id )  THEN
 
                    !> Cycle over already matched species
@@ -1409 +1432 @@
 
                          !> PMs are already in mass units: micrograms
-                         IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"  &
+                         IF (     TRIM(spc_names(match_spec_model(ispec)))=="PM1"   &
+                             .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"  &
                              .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
 
                             !              kg/(m^2*s) *kg/m^3                               
-                            surf_lsm_h%cssws(match_spec_model(ispec),m)= emiss_factor_side(match_spec_input(ispec)) *   &
+                            surf_lsm_h%cssws(match_spec_model(ispec),m)= emiss_factor_side(match_spec_input(ispec)) *  &
                             !                                                       kg/(m^2*s)
                                                                                 emis_distribution(1,j,i,ispec)*        &
@@ -1423 +1447 @@
                          !>Other Species: inputs are micromoles
                          !                 ppm/s *m *kg/m^3                
-                            surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_side(match_spec_input(ispec)) *   &
+                            surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_side(match_spec_input(ispec)) * &
                          !                                                    micromoles/(m^2*s)
-                                                                          emis_distribution(1,j,i,ispec) *              &
+                                                                          emis_distribution(1,j,i,ispec) *             &
                          !                                                    m**3/Nmole
-                                                                          conv_to_ratio(k,j,i)*                         &       
+                                                                          conv_to_ratio(k,j,i)*                        &
                          !                                                    kg/m^3
                                                                           rho_air(k)   
@@ -1466 +1490 @@
 
                    !> PMs
-                   IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"  &
-                          .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
+                   IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1"                                   &
+                         .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"                           &
+                         .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
                    
                       !            kg/(m^2*s) *kg/m^3                         kg/(m^2*s)                 
@@ -1522 +1547 @@
 
                    !> PMs
-                   IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"  &
-                          .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
+                   IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1"                                            &
+                         .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"                                    &
+                         .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
 
                       !         kg/(m^2*s) * kg/m^3                           kg/(m^2*s)            
@@ -1573 +1599 @@
 
                    !> PMs
-                   IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"  &
-                          .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
+                   IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1"                                     &
+                         .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"                             &
+                         .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
                    
                       !          kg/(m^2*s) *kg/m^3                             kg/(m^2*s)                     

palm/trunk/SOURCE/data_output_3d.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! (from branch resler)
+! Add Biometeorology
+! 
+! 3294 2018-10-01 02:37:10Z raasch
 ! changes concerning modularization of ocean option
 ! 
@@ -269 +273 @@
     USE radiation_model_mod,                                                   &
         ONLY:  nzub, nzut, radiation, radiation_data_output_3d
+
+    USE biometeorology_mod,                                                    &
+        ONLY: biometeorology_data_output_3d
 
     USE turbulence_closure_mod,                                                &
@@ -744 +751 @@
              IF ( .NOT. found  .AND.  radiation )  THEN
                 CALL radiation_data_output_3d( av, do3d(av,if), found,         &
+                                               local_pf, nzb_do, nzt_do )
+                resorted = .TRUE.
+             ENDIF
+
+!
+!--          Biometeorology quantity
+             IF ( .NOT. found  .AND.  radiation )  THEN
+                CALL biometeorology_data_output_3d( av, do3d(av,if), found,    &
                                                local_pf, nzb_do, nzt_do )
                 resorted = .TRUE.

palm/trunk/SOURCE/modules.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! (from branch resler)
+! Increase dimension of uv_heights etc.
+! 
+! 3302 2018-10-03 02:39:40Z raasch
 ! +u/v_stokes_zu, u/v_stokes_zw
 ! 
@@ -1566 +1570 @@
     REAL(wp) ::  tsc(10) = (/ 1.0_wp, 1.0_wp, 0.0_wp, 0.0_wp, &    !< array used for controlling time-integration at different substeps
                  0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp /)
-    REAL(wp) ::  u_profile(100) = 9999999.9_wp                     !< namelist parameter
-    REAL(wp) ::  uv_heights(100) = 9999999.9_wp                    !< namelist parameter
-    REAL(wp) ::  v_profile(100) = 9999999.9_wp                     !< namelist parameter
+    REAL(wp) ::  u_profile(200) = 9999999.9_wp                     !< namelist parameter
+    REAL(wp) ::  uv_heights(200) = 9999999.9_wp                    !< namelist parameter
+    REAL(wp) ::  v_profile(200) = 9999999.9_wp                     !< namelist parameter
     REAL(wp) ::  ug_vertical_gradient(10) = 0.0_wp                 !< namelist parameter
     REAL(wp) ::  ug_vertical_gradient_level(10) = -9999999.9_wp    !< namelist parameter

palm/trunk/SOURCE/netcdf_data_input_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! (from branch resler)
+! Formatting
+! 
+! 3298 2018-10-02 12:21:11Z kanani
 ! Modified EXPERT mode to PRE-PROCESSED mode (Russo)
 ! Introduced Chemistry static netcdf file (Russo)
@@ -1030 +1034 @@
                 ALLOCATE(emt_att%hourly_emis_time_factor(1:emt_att%ncat,1:emt_att%nhoursyear))
                 !Read-in Variable
-                CALL get_variable(id_emis,"emission_time_factors",emt_att%hourly_emis_time_factor,1,emt_att%ncat,1,emt_att%nhoursyear)
+                CALL get_variable(id_emis,"emission_time_factors",emt_att%hourly_emis_time_factor,1,   &
+                                  emt_att%ncat,1,emt_att%nhoursyear)
 
                 !-- MDH
@@ -1038 +1043 @@
                 ALLOCATE(emt_att%mdh_emis_time_factor(1:emt_att%ncat,1:emt_att%nmonthdayhour))
                 !-- Read-in Variable
-                CALL get_variable(id_emis,"emission_time_factors",emt_att%mdh_emis_time_factor,1,emt_att%ncat,1,emt_att%nmonthdayhour)
+                CALL get_variable(id_emis,"emission_time_factors",emt_att%mdh_emis_time_factor,1,       &
+                                  emt_att%ncat,1,emt_att%nmonthdayhour)
 
              ELSE
@@ -1055 +1061 @@
              DO ispec=1,emt_att%nspec
 
-                IF ( .NOT. ALLOCATED( emt(ispec)%default_emission_data ) ) ALLOCATE(emt(ispec)%default_emission_data(1:emt_att%ncat,1:ny+1,1:nx+1))
+                IF ( .NOT. ALLOCATED( emt(ispec)%default_emission_data ) )                              &
+                    ALLOCATE(emt(ispec)%default_emission_data(1:emt_att%ncat,1:ny+1,1:nx+1))
 
                 ALLOCATE(dum_var_3d(1:emt_att%ncat,nys+1:nyn+1,nxl+1:nxr+1))
@@ -1061 +1068 @@
                 CALL get_variable(id_emis,"emission_values",dum_var_3d,ispec,1,emt_att%ncat,nys,nyn,nxl,nxr)          
 
-                emt(ispec)%default_emission_data(:,nys+1:nyn+1,nxl+1:nxr+1)=dum_var_3d(1:emt_att%ncat,nys+1:nyn+1,nxl+1:nxr+1)
+                emt(ispec)%default_emission_data(:,nys+1:nyn+1,nxl+1:nxr+1) =                           &
+                    dum_var_3d(1:emt_att%ncat,nys+1:nyn+1,nxl+1:nxr+1)
 
                 DEALLOCATE (dum_var_3d)
@@ -1105 +1113 @@
 
                 !Allocation for the entire domain has to be done only for the first processor between all the subdomains     
-                IF ( .NOT. ALLOCATED( emt(ispec)%preproc_emission_data ) ) ALLOCATE(emt(ispec)%preproc_emission_data(emt_att%dt_emission,1,1:ny+1,1:nx+1)) 
+                IF ( .NOT. ALLOCATED( emt(ispec)%preproc_emission_data ) )                              &
+                    ALLOCATE(emt(ispec)%preproc_emission_data(emt_att%dt_emission,1,1:ny+1,1:nx+1))
 
                 !> allocate variable where to pass emission values read from netcdf
@@ -1114 +1123 @@
 
       
-                emt(ispec)%preproc_emission_data(:,1,nys+1:nyn+1,nxl+1:nxr+1)=dum_var_4d(1:emt_att%dt_emission,1,nys+1:nyn+1,nxl+1:nxr+1)
+                emt(ispec)%preproc_emission_data(:,1,nys+1:nyn+1,nxl+1:nxr+1) =                         &
+                      dum_var_4d(1:emt_att%dt_emission,1,nys+1:nyn+1,nxl+1:nxr+1)
 
                 DEALLOCATE ( dum_var_4d )

palm/trunk/SOURCE/netcdf_interface_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! (from branch resler)
+! Add biometeorology
+! 
+! 3294 2018-10-01 02:37:10Z raasch
 ! changes concerning modularization of ocean option
 ! 
@@ -610 +614 @@
         ONLY: radiation, radiation_define_netcdf_grid
 
+    USE biometeorology_mod,                                                    &
+        ONLY: biometeorology_define_netcdf_grid
+
     USE spectra_mod,                                                           &
         ONLY:  averaging_interval_sp, comp_spectra_level, data_output_sp, dosp_time_count, spectra_direction
@@ -1016 +1023 @@
                                                          grid_z )
                    ENDIF
+!
+!--                Check for biometeorology quantities
+                   IF ( .NOT. found  .AND.  radiation )  THEN
+                      CALL biometeorology_define_netcdf_grid( domask(mid,av,i),&
+                                                         found, grid_x, grid_y,&
+                                                         grid_z )
+                   ENDIF
 
                    CALL tcm_define_netcdf_grid( domask( mid,av,i), found,      &
@@ -1545 +1559 @@
                                                    grid_y, grid_z )
                    ENDIF
-                   
+
 !
 !--                Check for radiation quantities
                    IF ( .NOT. found  .AND.  radiation )  THEN
                       CALL radiation_define_netcdf_grid( do3d(av,i), found,    &
+                                                         grid_x, grid_y,       &
+                                                         grid_z )
+                   ENDIF
+
+!
+!--                Check for biometeorology quantities
+                   IF ( .NOT. found  .AND.  radiation )  THEN
+                      CALL biometeorology_define_netcdf_grid( do3d(av,i), found,&
                                                          grid_x, grid_y,       &
                                                          grid_z )
@@ -2322 +2344 @@
 
 !
+!--                      Check for biometeorology quantities
+                         IF ( .NOT. found  .AND.  radiation )  THEN
+                            CALL biometeorology_define_netcdf_grid( do2d(av,i),&
+                                                         found, grid_x, grid_y,&
+                                                         grid_z )
+                         ENDIF
+
+!
 !--                      Check for gust module quantities
                          IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
@@ -3034 +3064 @@
 
 !
+!--                   Check for biometeorology quantities
+                      IF ( .NOT. found  .AND.  radiation )  THEN
+                         CALL biometeorology_define_netcdf_grid( do2d(av,i),   &
+                                                            found,             &
+                                                            grid_x, grid_y,    &
+                                                            grid_z )
+                      ENDIF
+
+!
 !--                   Check for gust module quantities
                       IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
@@ -3694 +3733 @@
 
 !
+!--                   Check for biometeorology quantities
+                      IF ( .NOT. found  .AND.  radiation )  THEN
+                         CALL biometeorology_define_netcdf_grid( do2d(av,i),   &
+                                                            found,             &
+                                                            grid_x, grid_y,    &
+                                                            grid_z )
+                      ENDIF
+
+!
 !--                   Check for gust module quantities
                       IF ( .NOT. found  .AND.  gust_module_enabled )  THEN

palm/trunk/SOURCE/palm.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! (from branch resler)
+! Fix chemistry call
+! 
+! 3298 2018-10-02 12:21:11Z kanani
 ! - Minor formatting (kanani)
 ! - Added Call of date_and_time_init (Russo)
@@ -425 +429 @@
 !-- --> Needs to be moved!! What is the dependency about?
 ! IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
-    IF ( air_chemistry  .AND.  do_emis )  THEN
+    IF ( air_chemistry )  THEN
 
        IF ( do_emis ) THEN

palm/trunk/SOURCE/plant_canopy_model_mod.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! Fix reading plant canopy over buildings
+! 
+! 3294 2018-10-01 02:37:10Z raasch
 ! ocean renamed ocean_mode
 ! 
@@ -1131 +1134 @@
        INTEGER(iwp)                        ::  i, j      !< running index
        INTEGER(iwp)                        ::  nzp       !< number of vertical layers of plant canopy
+       INTEGER(iwp)                        ::  nzpltop   !<
+       INTEGER(iwp)                        ::  nzpl      !<
+       INTEGER(iwp)                        ::  kk        !<
        
        REAL(wp), DIMENSION(:), ALLOCATABLE ::  col   !< vertical column of input data
@@ -1144 +1150 @@
                  FORM='FORMATTED', ERR=515)
        READ(152, *, ERR=516, END=517)  nzp   !< read first line = number of vertical layers
-       
+       nzpltop = MIN(nzt+1, nzb+nzp-1)
+       nzpl = nzpltop - nzb + 1    !< no. of layers to assign
        ALLOCATE( col(0:nzp-1) )
 
@@ -1164 +1171 @@
                    CALL message( 'pcm_read_plant_canopy_3d', 'PA0349', 1, 2, 0, 6, 0 )
                 ENDIF
-                lad_s(0:nzp-1,j,i) = col(0:nzp-1) * lad_type_coef(pctype)
-                
+                kk = get_topography_top_index_ji( j, i, 's' )
+                lad_s(nzb:nzpltop-kk, j, i) = col(kk:nzpl-1)*lad_type_coef(pctype)
              CASE DEFAULT
                 WRITE(message_string, '(a,i2,a)')   &

palm/trunk/SOURCE/prognostic_equations.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! (from branch resler)
+! Fix for chemistry call
+! 
+! 3302 2018-10-03 02:39:40Z raasch
 ! Stokes drift + wave breaking term added
 ! 
@@ -482 +486 @@
        lsp_usr = 1
 !
-!--    If required, calculate photolysis frequencies - 
-!--    UNFINISHED: Why not before the intermediate timestep loop?
-       IF ( intermediate_timestep_count ==  1 )  THEN
-          CALL photolysis_control
-       ENDIF
-!
 !--    Chemical reactions
        CALL cpu_log( log_point(82), '(chem react + exch_h)', 'start' )
   
        IF ( chem_gasphase_on ) THEN
+!
+!--       If required, calculate photolysis frequencies - 
+!--       UNFINISHED: Why not before the intermediate timestep loop?
+          IF ( intermediate_timestep_count ==  1 )  THEN
+             CALL photolysis_control
+          ENDIF
           DO  i = nxl, nxr
              DO  j = nys, nyn

palm/trunk/SOURCE/radiation_model_mod.f90

@@ -15 +15 @@
 ! PALM. If not, see <http://www.gnu.org/licenses/>.
 !
-! Copyright 2015-2018 Czech Technical University in Prague
 ! Copyright 2015-2018 Institute of Computer Science of the 
 !                     Czech Academy of Sciences, Prague
+! Copyright 2015-2018 Czech Technical University in Prague
 ! Copyright 1997-2018 Leibniz Universitaet Hannover
 !------------------------------------------------------------------------------!
@@ -28 +28 @@
 ! -----------------
 ! $Id$
+! - New RTM version 2.9: (Pavel Krc, Jaroslav Resler, ICS, Prague)
+!   added calculation of the MRT inside the RTM module
+!   MRT fluxes are consequently used in the new biometeorology module
+!   for calculation of biological indices (MRT, PET)
+!   Fixes of v. 2.5 and SVN trunk:
+!    - proper initialization of rad_net_l
+!    - make arrays nsurfs and surfstart TARGET to prevent some MPI problems
+!    - initialization of arrays used in MPI one-sided operation as 1-D arrays
+!      to prevent problems with some MPI/compiler combinations
+!    - fix indexing of target displacement in subroutine request_itarget to 
+!      consider nzub
+!    - fix LAD dimmension range in PCB calculation
+!    - check ierr in all MPI calls
+!    - use proper per-gridbox sky and diffuse irradiance
+!    - fix shading for reflected irradiance
+!    - clear away the residuals of "atmospheric surfaces" implementation
+!    - fix rounding bug in raytrace_2d introduced in SVN trunk
+! - New RTM version 2.5: (Pavel Krc, Jaroslav Resler, ICS, Prague)
+!   can use angular discretization for all SVF
+!   (i.e. reflected and emitted radiation in addition to direct and diffuse),
+!   allowing for much better scaling wih high resoltion and/or complex terrain
+! - Unite array grow factors
+! - Fix slightly shifted terrain height in raytrace_2d
+! - Use more efficient MPI_Win_allocate for reverse gridsurf index
+! - Fix random MPI RMA bugs on Intel compilers
+! - Fix approx. double plant canopy sink values for reflected radiation
+! - Fix mostly missing plant canopy sinks for direct radiation
+! - Fix discretization errors for plant canopy sink in diffuse radiation
+! - Fix rounding errors in raytrace_2d
+!
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -432 +463 @@
 !> @todo Output of other rrtm arrays (such as volume mixing ratios)
 !> @todo Check for mis-used NINT() calls in raytrace_2d
+!>       RESULT: Original was correct (carefully verified formula), the change
+!>               to INT broke raytracing      -- P. Krc
 !> @todo Optimize radiation_tendency routines
 !>
@@ -440 +473 @@
  
     USE arrays_3d,                                                             &
-        ONLY:  dzw, hyp, nc, pt, q, ql, zu, zw, exner, d_exner
+        ONLY:  dzw, hyp, nc, pt, p, q, ql, u, v, w, zu, zw, exner, d_exner
 
     USE basic_constants_and_equations_mod,                                     &
@@ -735 +768 @@
                                              rrtm_swhr,      & !< RRTM output of shortwave radiation heating rate (K/d)
                                              rrtm_swhrc,     & !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d) 
-                                             rrtm_dirdflux,  & !< RRTM output of incoming direct shortwave (W/m²) 
-                                             rrtm_difdflux     !< RRTM output of incoming diffuse shortwave (W/m²) 
+                                             rrtm_dirdflux,  & !< RRTM output of incoming direct shortwave (W/m) 
+                                             rrtm_difdflux     !< RRTM output of incoming diffuse shortwave (W/m) 
 
     REAL(wp), DIMENSION(1) ::                rrtm_aldif,     & !< surface albedo for longwave diffuse radiation
@@ -771 +804 @@
     INTEGER(iwp), PARAMETER                        ::  ndsvf = 2                          !< number of dimensions of real values in SVF
     INTEGER(iwp), PARAMETER                        ::  idsvf = 2                          !< number of dimensions of integer values in SVF
-    INTEGER(iwp), PARAMETER                        ::  ndcsf = 2                          !< number of dimensions of real values in CSF
+    INTEGER(iwp), PARAMETER                        ::  ndcsf = 1                          !< number of dimensions of real values in CSF
     INTEGER(iwp), PARAMETER                        ::  idcsf = 2                          !< number of dimensions of integer values in CSF
     INTEGER(iwp), PARAMETER                        ::  kdcsf = 4                          !< number of dimensions of integer values in CSF calculation array
@@ -779 +812 @@
     INTEGER(iwp), PARAMETER                        ::  ix = 4                             !< position of i-index in surfl and surf
 
-    INTEGER(iwp), PARAMETER                        ::  nsurf_type = 16                    !< number of surf types incl. phys.(land+urban) & (atm.,sky,boundary) surfaces - 1
+    INTEGER(iwp), PARAMETER                        ::  nsurf_type = 10                    !< number of surf types incl. phys.(land+urban) & (atm.,sky,boundary) surfaces - 1
 
     INTEGER(iwp), PARAMETER                        ::  iup_u    = 0                       !< 0 - index of urban upward surface (ground or roof)
@@ -794 +827 @@
     INTEGER(iwp), PARAMETER                        ::  iwest_l  = 10                      !< 10- index of land westward facing wall
 
-    INTEGER(iwp), PARAMETER                        ::  iup_a    = 11                      !< 11- index of atm. cell ubward virtual surface
-    INTEGER(iwp), PARAMETER                        ::  idown_a  = 12                      !< 12- index of atm. cell downward virtual surface
-    INTEGER(iwp), PARAMETER                        ::  inorth_a = 13                      !< 13- index of atm. cell northward facing virtual surface
-    INTEGER(iwp), PARAMETER                        ::  isouth_a = 14                      !< 14- index of atm. cell southward facing virtual surface
-    INTEGER(iwp), PARAMETER                        ::  ieast_a  = 15                      !< 15- index of atm. cell eastward facing virtual surface
-    INTEGER(iwp), PARAMETER                        ::  iwest_a  = 16                      !< 16- index of atm. cell westward facing virtual surface
-
-    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  idir = (/0, 0,0, 0,1,-1,0,0, 0,1,-1,0, 0,0, 0,1,-1/)   !< surface normal direction x indices
-    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  jdir = (/0, 0,1,-1,0, 0,0,1,-1,0, 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,1,-1,0, 0,0, 0/)   !< surface normal direction z indices
+    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  idir = (/0, 0,0, 0,1,-1,0,0, 0,1,-1/)   !< surface normal direction x indices
+    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
 
@@ -816 +842 @@
 
 !-- indices and sizes of urban and land surface models
-    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  surfl            !< coordinates of i-th local surface in local grid - surfl[:,k] = [d, z, y, x]
-    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  surf             !< coordinates of i-th surface in grid - surf[:,k] = [d, z, y, x]
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfl_l          !< coordinates of i-th local surface in local grid - surfl[:,k] = [d, z, y, x]
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surf_l           !< coordinates of i-th surface in grid - surf[:,k] = [d, z, y, x]
+    INTEGER(iwp), DIMENSION(:,:), POINTER          ::  surfl            !< coordinates of i-th local surface in local grid - surfl[:,k] = [d, z, y, x]
+    INTEGER(iwp), DIMENSION(:,:), POINTER          ::  surf             !< coordinates of i-th surface in grid - surf[:,k] = [d, z, y, x]
     INTEGER(iwp)                                   ::  nsurfl           !< number of all surfaces in local processor
-    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  nsurfs           !< array of number of all surfaces in individual processors
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  nsurfs           !< array of number of all surfaces in individual processors
     INTEGER(iwp)                                   ::  nsurf            !< global number of surfaces in index array of surfaces (nsurf = proc nsurfs)
-    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  surfstart        !< starts of blocks of surfaces for individual processors in array surf
-                                                                        !< respective block for particular processor is surfstart[iproc]+1 : surfstart[iproc+1]
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfstart        !< starts of blocks of surfaces for individual processors in array surf (indexed from 1)
+                                                                        !< respective block for particular processor is surfstart[iproc+1]+1 : surfstart[iproc+1]+nsurfs[iproc+1]
 
 !-- block variables needed for calculation of the plant canopy model inside the urban surface model
@@ -835 +863 @@
 
 !-- configuration parameters (they can be setup in PALM config)
-    LOGICAL                                        ::  rma_lad_raytrace = .FALSE.         !< use MPI RMA to access LAD for raytracing (instead of global array)
-    LOGICAL                                        ::  mrt_factors = .FALSE.              !< whether to generate MRT factor files during init
+    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 
+                                                                                          !< reflected radiation (as opposed to all mutually visible pairs)
+    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
     REAL(wp), PARAMETER                            ::  ext_coef = 0.6_wp                  !< extinction coefficient (a.k.a. alpha)
@@ -874 +908 @@
         INTEGER(iwp)                               :: isurfs            !< 
         REAL(wp)                                   :: rsvf              !< 
-        REAL(wp)                                   :: rtransp           !< 
     END TYPE
 
 !-- arrays storing the values of USM
-    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  svfsurf          !< svfsurf[:,isvf] = index of source and target surface for svf[isvf]
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  svfsurf          !< svfsurf[:,isvf] = index of target and source surface for svf[isvf]
     REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  svf              !< array of shape view factors+direct irradiation factors for local surfaces
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfins          !< array of sw radiation falling to local surface after i-th reflection
@@ -892 +925 @@
     INTEGER(iwp)                                   ::  ndsidir          !< number of apparent solar directions used
     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  dsidir_rev       !< dsidir_rev[ielev,iazim] = i for dsidir or -1 if not present
+
+    INTEGER(iwp)                                   ::  nmrtbl           !< No. of local grid boxes for which MRT is calculated
+    INTEGER(iwp)                                   ::  nmrtf            !< number of MRT factors for local processor
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  mrtbl            !< coordinates of i-th local MRT box - surfl[:,i] = [z, y, x]
+    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  mrtfsurf         !< mrtfsurf[:,imrtf] = index of target MRT box and source surface for mrtf[imrtf]
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtf             !< array of MRT factors for each local MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtft            !< array of MRT factors including transparency for each local MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtsky           !< array of sky view factor for each local MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtskyt          !< array of sky view factor including transparency for each local MRT box
+    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  mrtdsit          !< array of direct solar transparencies for each local MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtinsw          !< mean SW radiant flux for each MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtinlw          !< mean LW radiant flux for each MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrt              !< mean radiant temperature for each MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtinsw_av       !< time average mean SW radiant flux for each MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtinlw_av       !< time average mean LW radiant flux for each MRT box
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrt_av           !< time average mean radiant temperature for each MRT box
 
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinsw         !< array of sw radiation falling to local surface including radiation from reflections
@@ -920 +969 @@
     TYPE(t_svf), DIMENSION(:), POINTER             ::  asvf             !< pointer to growing svc array
     TYPE(t_csf), DIMENSION(:), POINTER             ::  acsf             !< pointer to growing csf array
+    TYPE(t_svf), DIMENSION(:), POINTER             ::  amrtf            !< pointer to growing mrtf array
     TYPE(t_svf), DIMENSION(:), ALLOCATABLE, TARGET ::  asvf1, asvf2     !< realizations of svf array
     TYPE(t_csf), DIMENSION(:), ALLOCATABLE, TARGET ::  acsf1, acsf2     !< realizations of csf array
+    TYPE(t_svf), DIMENSION(:), ALLOCATABLE, TARGET ::  amrtf1, amrtf2   !< realizations of mftf array
     INTEGER(iwp)                                   ::  nsvfla           !< dimmension of array allocated for storage of svf in local processor
     INTEGER(iwp)                                   ::  ncsfla           !< dimmension of array allocated for storage of csf in local processor
-    INTEGER(iwp)                                   ::  msvf, mcsf       !< mod for swapping the growing array
-    INTEGER(iwp), PARAMETER                        ::  gasize = 10000   !< initial size of growing arrays
+    INTEGER(iwp)                                   ::  nmrtfa           !< dimmension of array allocated for storage of mrt
+    INTEGER(iwp)                                   ::  msvf, mcsf, mmrtf!< mod for swapping the growing array
+    INTEGER(iwp), PARAMETER                        ::  gasize = 100000  !< initial size of growing arrays
+    REAL(wp), PARAMETER                            ::  grow_factor = 1.4_wp !< growth factor of growing arrays
     REAL(wp)                                       ::  dist_max_svf = -9999.0 !< maximum distance to calculate the minimum svf to be considered. It is
                                                                         !< used to avoid very small SVFs resulting from too far surfaces with mutual visibility
@@ -932 +985 @@
                                                                         !< needed only during calc_svf but must be here because it is
                                                                         !< shared between subroutines calc_svf and raytrace
-    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE    ::  gridpcbl         !< index of local pcb[k,j,i]
+    INTEGER(iwp), DIMENSION(:,:,:,:), POINTER      ::  gridsurf         !< reverse index of local surfl[d,k,j,i] (for case rad_angular_discretization)
+    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE    ::  gridpcbl         !< reverse index of local pcbl[k,j,i]
+    INTEGER(iwp), PARAMETER                        ::  nsurf_type_u = 6 !< number of urban surf types (used in gridsurf)
 
 !-- temporary arrays for calculation of csf in raytracing
@@ -942 +997 @@
     INTEGER(kind=MPI_ADDRESS_KIND), &
                   DIMENSION(:), ALLOCATABLE        ::  lad_disp         !< array of displaycements of lad in local array of proc lad_ip
+    INTEGER(iwp)                                   ::  win_lad          !< MPI RMA window for leaf area density
+    INTEGER(iwp)                                   ::  win_gridsurf     !< MPI RMA window for reverse grid surface index
 #endif
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  lad_s_ray        !< array of received lad_s for appropriate gridboxes crossed by ray
+    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  target_surfl
     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  rt2_track
     REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  rt2_track_lad
@@ -1083 +1141 @@
 !-- Public variables and constants / NEEDS SORTING
     PUBLIC albedo, albedo_type, decl_1, decl_2, decl_3, dots_rad, dt_radiation,&
-           emissivity, force_radiation_call,                                   &
-           lat, lon, rad_net_av, radiation, radiation_scheme, rad_lw_in,       &
+           emissivity, force_radiation_call, lat, lon,                         &
+           mrt_include_sw, mrt_nlevels, mrtbl, mrtinsw, mrtinlw, nmrtbl,       &
+           rad_net_av, radiation, radiation_scheme, rad_lw_in,                 &
            rad_lw_in_av, rad_lw_out, rad_lw_out_av,                            &
            rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in,  &
@@ -1091 +1150 @@
            skip_time_do_radiation, time_radiation, unscheduled_radiation_calls,&
            zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon,       &
-           nrefsteps, mrt_factors, dist_max_svf, nsvfl, svf,                   &
+           write_svf_on_init, read_svf_on_init,                                &
+           nrefsteps, dist_max_svf, nsvfl, svf,                                &
            svfsurf, surfinsw, surfinlw, surfins, surfinl, surfinswdir,         &
            surfinswdif, surfoutsw, surfoutlw, surfinlwdif, rad_sw_in_dir,      &
            rad_sw_in_diff, rad_lw_in_diff, surfouts, surfoutl, surfoutsl,      &
-           surfoutll, idir, jdir, kdir, id, iz, iy, ix, nsurfs, surfstart,     &
+           surfoutll, idir, jdir, kdir, id, iz, iy, ix,                        &
            surf, surfl, nsurfl, pcbinswdir, pcbinswdif, pcbinsw, pcbinlw,      &
-           iup_u, inorth_u, isouth_u, ieast_u, iwest_u,           &
+           iup_u, inorth_u, isouth_u, ieast_u, iwest_u,                        &
            iup_l, inorth_l, isouth_l, ieast_l, iwest_l,                        &
-           nsurf_type, nzub, nzut, nzu, pch, nsurf,                            &
-           iup_a, idown_a, inorth_a, isouth_a, ieast_a, iwest_a,               &
+           nsurf_type, nzub, nzut, nzpt, nzu, pch, nsurf,                      &
            idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct,                             &
            radiation_interactions, startwall, startland, endland, endwall,     &
-           skyvf, skyvft, radiation_interactions_on, average_radiation
+           skyvf, skyvft, radiation_interactions_on, average_radiation, npcbl, &
+           pcbl
 
 #if defined ( __rrtmg )
@@ -1163 +1223 @@
        SELECT CASE ( TRIM( var ) )
 
-         CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr', 'rad_sw_in' )
+          CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr', 'rad_sw_in' )
              IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
                 message_string = '"output of "' // TRIM( var ) // '" requi' // &
@@ -1205 +1265 @@
              IF ( TRIM( var ) == 'rrtm_asdir*'   ) unit = '' 
 
+          CASE ( 'rad_mrt', 'rad_mrt_sw', 'rad_mrt_lw'  )
+             IF ( .NOT.  radiation ) THEN
+                message_string = 'output of "' // TRIM( var ) // '" require'&
+                                 // 's radiation = .TRUE.'
+                CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
+             ENDIF
+             IF ( mrt_nlevels == 0 ) THEN
+                message_string = 'output of "' // TRIM( var ) // '" require'&
+                                 // 's mrt_nlevels > 0'
+                CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
+             ENDIF
+             IF ( TRIM( var ) == 'rad_mrt_sw'  .AND.  .NOT. mrt_include_sw ) THEN
+                message_string = 'output of "' // TRIM( var ) // '" require'&
+                                 // 's rad_mrt_sw = .TRUE.'
+                CALL message( 'check_parameters', 'PA0511', 1, 2, 0, 6, 0 )
+             ENDIF
+             IF ( TRIM( var ) == 'rad_mrt' ) THEN
+                unit = 'K'
+             ELSE
+                unit = 'W m-2'
+             ENDIF
+
           CASE DEFAULT
              unit = 'illegal'
@@ -1457 +1539 @@
           ENDIF
        ENDIF
+!
+!--    Parallel rad_angular_discretization without raytrace_mpi_rma is not implemented
+#if defined( __parallel )      
+       IF ( rad_angular_discretization  .AND.  .NOT. raytrace_mpi_rma )  THEN
+          message_string = 'rad_angular_discretization can only be used ' //  &
+                           'together with raytrace_mpi_rma or when ' //  &
+                           'no parallelization is applied.'
+          CALL message( 'check_parameters', 'PA0486', 1, 2, 0, 6, 0 )
+       ENDIF
+#endif
 
        IF ( cloud_droplets  .AND.   radiation_scheme == 'rrtmg'  .AND.         &
@@ -2354 +2446 @@
 !--    In case averaged radiation is used, calculate mean temperature and 
 !--    liquid water mixing ratio at the urban-layer top.
-       IF ( average_radiation ) THEN   
+       IF ( average_radiation ) THEN
           pt1   = 0.0_wp
           IF ( bulk_cloud_model  .OR.  cloud_droplets )  ql1   = 0.0_wp
@@ -2364 +2456 @@
           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
           CALL MPI_ALLREDUCE( pt1_l, pt1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
-          IF ( bulk_cloud_model  .OR.  cloud_droplets )                            &
-             CALL MPI_ALLREDUCE( ql1_l, ql1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
+          IF ( ierr /= 0 ) THEN
+              WRITE(9,*) 'Error MPI_AllReduce1:', ierr, pt1_l, pt1
+              FLUSH(9)
+          ENDIF
+
+          IF ( bulk_cloud_model  .OR.  cloud_droplets ) THEN
+              CALL MPI_ALLREDUCE( ql1_l, ql1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
+              IF ( ierr /= 0 ) THEN
+                  WRITE(9,*) 'Error MPI_AllReduce2:', ierr, ql1_l, ql1
+                  FLUSH(9)
+              ENDIF
+          ENDIF
 #else
           pt1 = pt1_l 
@@ -2534 +2636 @@
           IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
           CALL MPI_ALLREDUCE( pt1_l, pt1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
-          IF ( bulk_cloud_model  .OR.  cloud_droplets )                             &
+          IF ( ierr /= 0 ) THEN
+              WRITE(9,*) 'Error MPI_AllReduce3:', ierr, pt1_l, pt1
+              FLUSH(9)
+          ENDIF
+          IF ( bulk_cloud_model  .OR.  cloud_droplets )  THEN
              CALL MPI_ALLREDUCE( ql1_l, ql1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
+             IF ( ierr /= 0 ) THEN
+                 WRITE(9,*) 'Error MPI_AllReduce4:', ierr, ql1_l, ql1
+                 FLUSH(9)
+             ENDIF
+          ENDIF
 #else
           pt1 = pt1_l
@@ -2756 +2867 @@
                                   albedo_lw_dif, albedo_sw_dir, albedo_sw_dif, &
                                   constant_albedo, dt_radiation, emissivity,   &
-                                  lw_radiation, net_radiation,                 &
+                                  lw_radiation, mrt_nlevels, mrt_skip_roof,    &
+                                  mrt_include_sw,  net_radiation,              &
                                   radiation_scheme, skip_time_do_radiation,    &
                                   sw_radiation, unscheduled_radiation_calls,   &
+                                  read_svf_on_init, write_svf_on_init,         &
                                   max_raytracing_dist, min_irrf_value,         &
-                                  nrefsteps, mrt_factors, rma_lad_raytrace,    &
+                                  nrefsteps, raytrace_mpi_rma,                 &
                                   dist_max_svf,                                &
                                   surface_reflections, svfnorm_report_thresh,  &
-                                  radiation_interactions_on
+                                  radiation_interactions_on,                   &
+                                  rad_angular_discretization,                  &
+                                  raytrace_discrete_azims,                     &
+                                  raytrace_discrete_elevs
    
        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, net_radiation,                 &
+                                  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, mrt_factors, rma_lad_raytrace,    &
+                                  nrefsteps, raytrace_mpi_rma,                 &
                                   dist_max_svf,                                &
                                   surface_reflections, svfnorm_report_thresh,  &
-                                  radiation_interactions_on
+                                  radiation_interactions_on,                   &
+                                  rad_angular_discretization,                  &
+                                  raytrace_discrete_azims,                     &
+                                  raytrace_discrete_elevs
    
        line = ' '
@@ -4456 +4576 @@
      INTEGER(iwp)                      :: i, j, k, kk, is, js, d, ku, refstep, m, mm, l, ll
      INTEGER(iwp)                      :: isurf, isurfsrc, isvf, icsf, ipcgb
+     INTEGER(iwp)                      :: imrt, imrtf
      INTEGER(iwp)                      :: isd                !< solar direction number
      INTEGER(iwp)                      :: pc_box_dimshift    !< transform for best accuracy
@@ -4547 +4668 @@
      surfoutsl(:)  = 0.0_wp !start-end
      surfoutll(:)  = 0.0_wp !start-end
+     IF ( nmrtbl > 0 )  THEN
+        mrtinsw(:) = 0._wp
+        mrtinlw(:) = 0._wp
+     ENDIF
+
 
 !--  Set up thermal radiation from surfaces
@@ -4623 +4749 @@
      CALL MPI_AllGatherv(surfoutll, nsurfl, MPI_REAL, &
                          surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr) !nsurf global
+     IF ( ierr /= 0 ) THEN
+         WRITE(9,*) 'Error MPI_AllGatherv1:', ierr, SIZE(surfoutll), nsurfl, &
+                     SIZE(surfoutl), nsurfs, surfstart
+         FLUSH(9)
+     ENDIF
 #else
      surfoutl(:) = surfoutll(:) !nsurf global
@@ -4639 +4770 @@
         ENDDO
      ENDIF
-
-     !-- diffuse radiation using sky view factor, TODO: homogeneous rad_*w_in_diff because now it depends on no. of processors
-     surfinswdif(:) = rad_sw_in_diff(nyn,nxl) * skyvft(:)
-     surfinlwdif(:) = rad_lw_in_diff(nyn,nxl) * skyvf(:)
+!
+!--  diffuse radiation using sky view factor
+     DO isurf = 1, nsurfl
+        j = surfl(iy, isurf)
+        i = surfl(ix, isurf)
+        surfinswdif(isurf) = rad_sw_in_diff(j,i) * skyvft(isurf)
+        surfinlwdif(isurf) = rad_lw_in_diff(j,i) * skyvf(isurf)
+     ENDDO
+!
+!--  MRT diffuse irradiance
+     DO  imrt = 1, nmrtbl
+        j = mrtbl(iy, imrt)
+        i = mrtbl(ix, imrt)
+        mrtinsw(imrt) = mrtskyt(imrt) * rad_sw_in_diff(j,i)
+        mrtinlw(imrt) = mrtsky(imrt) * rad_lw_in_diff(j,i)
+     ENDDO
 
      !-- direct radiation
@@ -4654 +4797 @@
         isd = dsidir_rev(j, i)
         DO isurf = 1, nsurfl
-           surfinswdir(isurf) = rad_sw_in_dir(nyn,nxl) * costheta(surfl(id, isurf)) * dsitrans(isurf, isd) / zenith(0)
+           j = surfl(iy, isurf)
+           i = surfl(ix, isurf)
+           surfinswdir(isurf) = rad_sw_in_dir(j,i) * costheta(surfl(id, isurf)) * dsitrans(isurf, isd) / zenith(0)
+        ENDDO
+!
+!--     MRT direct irradiance
+        DO  imrt = 1, nmrtbl
+           j = mrtbl(iy, imrt)
+           i = mrtbl(ix, imrt)
+           mrtinsw(imrt) = mrtinsw(imrt) + mrtdsit(imrt, isd) * rad_sw_in_dir(j,i) &
+                                     / zenith(0) / 4._wp ! normal to sphere
         ENDDO
      ENDIF
+!
+!--  MRT first pass thermal
+     DO  imrtf = 1, nmrtf
+        imrt = mrtfsurf(1, imrtf)
+        isurfsrc = mrtfsurf(2, imrtf)
+        mrtinlw(imrt) = mrtinlw(imrt) + mrtf(imrtf) * surfoutl(isurfsrc)
+     ENDDO
 
      IF ( npcbl > 0 )  THEN
@@ -4674 +4834 @@
              IF ( isurfsrc == -1 )  THEN
 !--                 Diffuse rad from sky.
-                 pcbinswdif(ipcgb) = csf(1,icsf) * csf(2,icsf) * rad_sw_in_diff(j,i)
+                 pcbinswdif(ipcgb) = csf(1,icsf) * rad_sw_in_diff(j,i)
 
                  !--Direct rad
@@ -4685 +4845 @@
                                         * pc_abs_frac * dsitransc(ipcgb, isd)
                  ENDIF
-
-                 EXIT ! only isurfsrc=-1 is processed here
              ENDIF
          ENDDO
@@ -4722 +4880 @@
          CALL MPI_AllGatherv(surfoutsl, nsurfl, MPI_REAL, &
              surfouts, nsurfs, surfstart, MPI_REAL, comm2d, ierr)
+         IF ( ierr /= 0 ) THEN
+             WRITE(9,*) 'Error MPI_AllGatherv2:', ierr, SIZE(surfoutsl), nsurfl, &
+                        SIZE(surfouts), nsurfs, surfstart
+             FLUSH(9)
+         ENDIF
+
          CALL MPI_AllGatherv(surfoutll, nsurfl, MPI_REAL, &
              surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr)
+         IF ( ierr /= 0 ) THEN
+             WRITE(9,*) 'Error MPI_AllGatherv3:', ierr, SIZE(surfoutll), nsurfl, &
+                        SIZE(surfoutl), nsurfs, surfstart
+             FLUSH(9)
+         ENDIF
+
 #else
          surfouts = surfoutsl
@@ -4747 +4917 @@
              IF ( isurfsrc == -1 )  CYCLE ! sky->face only in 1st pass, not here
 
-             pcbinsw(ipcgb) = pcbinsw(ipcgb) + csf(1,icsf) * csf(2,icsf) * surfouts(isurfsrc)
+             pcbinsw(ipcgb) = pcbinsw(ipcgb) + csf(1,icsf) * surfouts(isurfsrc)
+         ENDDO
+!
+!--      MRT reflected
+         DO  imrtf = 1, nmrtf
+            imrt = mrtfsurf(1, imrtf)
+            isurfsrc = mrtfsurf(2, imrtf)
+            mrtinsw(imrt) = mrtinsw(imrt) + mrtft(imrtf) * surfouts(isurfsrc)
+            mrtinlw(imrt) = mrtinlw(imrt) + mrtf(imrtf) * surfoutl(isurfsrc)
          ENDDO
 
@@ -4755 +4933 @@
          surfoutlw = surfoutlw + surfoutll
 
-     ENDDO
+     ENDDO ! refstep
 
 !--  push heat flux absorbed by plant canopy to respective 3D arrays
@@ -4778 +4956 @@
      ENDIF
 !
+!--  Calculate black body MRT (after all reflections)
+     IF ( nmrtbl > 0 )  THEN
+        IF ( mrt_include_sw )  THEN
+           mrt(:) = ((mrtinsw(:) + mrtinlw(:)) / sigma_sb) ** .25_wp
+        ELSE
+           mrt(:) = (mrtinlw(:) / sigma_sb) ** .25_wp
+        ENDIF
+     ENDIF
+!
 !--     Transfer radiation arrays required for energy balance to the respective data types
      DO  i = 1, nsurfl
@@ -4791 +4978 @@
            surf_usm_h%rad_net(m)    = surfinsw(i) - surfoutsw(i) +          &
                                       surfinlw(i) - surfoutlw(i)
+           surf_usm_h%rad_net_l(m)  = surf_usm_h%rad_net(m)
 !
 !--     northward-facding
@@ -4800 +4988 @@
            surf_usm_v(0)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
                                          surfinlw(i) - surfoutlw(i)
+           surf_usm_v(0)%rad_net_l(m)  = surf_usm_v(0)%rad_net(m)
 !
 !--     southward-facding
@@ -4809 +4998 @@
            surf_usm_v(1)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
                                          surfinlw(i) - surfoutlw(i)
+           surf_usm_v(1)%rad_net_l(m)  = surf_usm_v(1)%rad_net(m)
 !
 !--     eastward-facing
@@ -4818 +5008 @@
            surf_usm_v(2)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
                                          surfinlw(i) - surfoutlw(i)
+           surf_usm_v(2)%rad_net_l(m)  = surf_usm_v(2)%rad_net(m)
 !
 !--     westward-facding
@@ -4827 +5018 @@
            surf_usm_v(3)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
                                          surfinlw(i) - surfoutlw(i)
+           surf_usm_v(3)%rad_net_l(m)  = surf_usm_v(3)%rad_net(m)
 !
 !--     (2) land surfaces
@@ -4957 +5149 @@
 #if defined( __parallel )
      CALL MPI_ALLREDUCE( pinswl, pinsw, 1, MPI_REAL, MPI_SUM, comm2d, ierr)
+     IF ( ierr /= 0 ) THEN
+         WRITE(9,*) 'Error MPI_AllReduce5:', ierr, pinswl, pinsw
+         FLUSH(9)
+     ENDIF
      CALL MPI_ALLREDUCE( pinlwl, pinlw, 1, MPI_REAL, MPI_SUM, comm2d, ierr)
+     IF ( ierr /= 0 ) THEN
+         WRITE(9,*) 'Error MPI_AllReduce6:', ierr, pinlwl, pinlw
+         FLUSH(9)
+     ENDIF
      CALL MPI_ALLREDUCE( pabsswl, pabssw, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
+     IF ( ierr /= 0 ) THEN
+         WRITE(9,*) 'Error MPI_AllReduce7:', ierr, pabsswl, pabssw
+         FLUSH(9)
+     ENDIF
      CALL MPI_ALLREDUCE( pabslwl, pabslw, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
+     IF ( ierr /= 0 ) THEN
+         WRITE(9,*) 'Error MPI_AllReduce8:', ierr, pabslwl, pabslw
+         FLUSH(9)
+     ENDIF
      CALL MPI_ALLREDUCE( pemitlwl, pemitlw, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
+     IF ( ierr /= 0 ) THEN
+         WRITE(9,*) 'Error MPI_AllReduce8:', ierr, pemitlwl, pemitlw
+         FLUSH(9)
+     ENDIF
      CALL MPI_ALLREDUCE( emiss_sum_surfl, emiss_sum_surf, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
+     IF ( ierr /= 0 ) THEN
+         WRITE(9,*) 'Error MPI_AllReduce9:', ierr, emiss_sum_surfl, emiss_sum_surf
+         FLUSH(9)
+     ENDIF
      CALL MPI_ALLREDUCE( area_surfl, area_surf, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
+     IF ( ierr /= 0 ) THEN
+         WRITE(9,*) 'Error MPI_AllReduce10:', ierr, area_surfl, area_surf
+         FLUSH(9)
+     ENDIF
 #else
      pinsw = pinswl
@@ -5038 +5258 @@
 !      t_rad_urb = pt_surf_urb / count_surfaces
      
+
     CONTAINS
 
@@ -5201 +5422 @@
        INTEGER(iwp) :: i, j, k, l, m
        INTEGER(iwp) :: k_topo     !< vertical index indicating topography top for given (j,i)
-       INTEGER(iwp) :: nzptl, nzubl, nzutl, isurf, ipcgb
+       INTEGER(iwp) :: nzptl, nzubl, nzutl, isurf, ipcgb, imrt
        REAL(wp)     :: mrl
+#if defined( __parallel )
+       INTEGER(iwp), DIMENSION(:), POINTER       ::  gridsurf_rma   !< fortran pointer, but lower bounds are 1
+       TYPE(c_ptr)                               ::  gridsurf_rma_p !< allocated c pointer
+       INTEGER(iwp)                              ::  minfo          !< MPI RMA window info handle
+#endif
 
 
@@ -5263 +5489 @@
 #if defined( __parallel )
        CALL MPI_AllReduce(nzubl, nzub, 1, MPI_INTEGER, MPI_MIN, comm2d, ierr )
+       IF ( ierr /= 0 ) THEN
+           WRITE(9,*) 'Error MPI_AllReduce11:', ierr, nzubl, nzub
+           FLUSH(9)
+       ENDIF
        CALL MPI_AllReduce(nzutl, nzut, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )
+       IF ( ierr /= 0 ) THEN
+           WRITE(9,*) 'Error MPI_AllReduce12:', ierr, nzutl, nzut
+           FLUSH(9)
+       ENDIF
        CALL MPI_AllReduce(nzptl, nzpt, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )
+       IF ( ierr /= 0 ) THEN
+           WRITE(9,*) 'Error MPI_AllReduce13:', ierr, nzptl, nzpt
+           FLUSH(9)
+       ENDIF
 #else
        nzub = nzubl
@@ -5353 +5591 @@
 !--    cycles must not be altered, certain file input routines may depend
 !--    on it)
-       ALLOCATE(surfl(5,nsurfl))  ! is it mecessary to allocate it with (5,nsurfl)?
+       ALLOCATE(surfl_l(5*nsurfl))  ! is it necessary to allocate it with (5,nsurfl)?
+       surfl(1:5,1:nsurfl) => surfl_l(1:5*nsurfl)
        isurf = 0
+       IF ( rad_angular_discretization )  THEN
+!
+!--       Allocate and fill the reverse indexing array gridsurf
+#if defined( __parallel )
+!
+!--       raytrace_mpi_rma is asserted
+
+          CALL MPI_Info_create(minfo, ierr)
+          IF ( ierr /= 0 ) THEN
+              WRITE(9,*) 'Error MPI_Info_create1:', ierr
+              FLUSH(9)
+          ENDIF
+          CALL MPI_Info_set(minfo, 'accumulate_ordering', '', ierr)
+          IF ( ierr /= 0 ) THEN
+              WRITE(9,*) 'Error MPI_Info_set1:', ierr
+              FLUSH(9)
+          ENDIF
+          CALL MPI_Info_set(minfo, 'accumulate_ops', 'same_op', ierr)
+          IF ( ierr /= 0 ) THEN
+              WRITE(9,*) 'Error MPI_Info_set2:', ierr
+              FLUSH(9)
+          ENDIF
+          CALL MPI_Info_set(minfo, 'same_size', 'true', ierr)
+          IF ( ierr /= 0 ) THEN
+              WRITE(9,*) 'Error MPI_Info_set3:', ierr
+              FLUSH(9)
+          ENDIF
+          CALL MPI_Info_set(minfo, 'same_disp_unit', 'true', ierr)
+          IF ( ierr /= 0 ) THEN
+              WRITE(9,*) 'Error MPI_Info_set4:', ierr
+              FLUSH(9)
+          ENDIF
+
+          CALL MPI_Win_allocate(INT(STORAGE_SIZE(1_iwp)/8*nsurf_type_u*nzu*nny*nnx, &
+                                    kind=MPI_ADDRESS_KIND),                         &
+                                INT(STORAGE_SIZE(1_iwp)/8, kind=MPI_ADDRESS_KIND),  &
+                                minfo, comm2d, gridsurf_rma_p, win_gridsurf, ierr)
+          IF ( ierr /= 0 ) THEN
+              WRITE(9,*) 'Error MPI_Win_allocate1:', ierr, &
+                 INT(STORAGE_SIZE(1_iwp)/8*nsurf_type_u*nzu*nny*nnx,kind=MPI_ADDRESS_KIND), &
+                 INT(STORAGE_SIZE(1_iwp)/8, kind=MPI_ADDRESS_KIND), win_gridsurf
+              FLUSH(9)
+          ENDIF
+
+          CALL MPI_Info_free(minfo, ierr)
+          IF ( ierr /= 0 ) THEN
+              WRITE(9,*) 'Error MPI_Info_free1:', ierr
+              FLUSH(9)
+          ENDIF
+
+!
+!--       On Intel compilers, calling c_f_pointer to transform a C pointer
+!--       directly to a multi-dimensional Fotran pointer leads to strange
+!--       errors on dimension boundaries. However, transforming to a 1D
+!--       pointer and then redirecting a multidimensional pointer to it works
+!--       fine.
+          CALL c_f_pointer(gridsurf_rma_p, gridsurf_rma, (/ nsurf_type_u*nzu*nny*nnx /))
+          gridsurf(0:nsurf_type_u-1, nzub:nzut, nys:nyn, nxl:nxr) =>                &
+                     gridsurf_rma(1:nsurf_type_u*nzu*nny*nnx)
+#else
+          ALLOCATE(gridsurf(0:nsurf_type_u-1,nzub:nzut,nys:nyn,nxl:nxr) )
+#endif
+          gridsurf(:,:,:,:) = -999
+       ENDIF
 
 !--    add horizontal surface elements (land and urban surfaces)
@@ -5362 +5665 @@
               DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
                  k = surf_usm_h%k(m)
-
                  isurf = isurf + 1
                  surfl(:,isurf) = (/iup_u,k,j,i,m/)
+                 IF ( rad_angular_discretization ) THEN
+                    gridsurf(iup_u,k,j,i) = isurf
+                 ENDIF
               ENDDO
 
               DO  m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
                  k = surf_lsm_h%k(m)
-
                  isurf = isurf + 1
                  surfl(:,isurf) = (/iup_l,k,j,i,m/)
+                 IF ( rad_angular_discretization ) THEN
+                    gridsurf(iup_u,k,j,i) = isurf
+                 ENDIF
               ENDDO
 
@@ -5384 +5691 @@
               DO  m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
                  k = surf_usm_v(l)%k(m)
-
-                 isurf          = isurf + 1
+                 isurf = isurf + 1
                  surfl(:,isurf) = (/inorth_u,k,j,i,m/)
+                 IF ( rad_angular_discretization ) THEN
+                    gridsurf(inorth_u,k,j,i) = isurf
+                 ENDIF
               ENDDO
               DO  m = surf_lsm_v(l)%start_index(j,i), surf_lsm_v(l)%end_index(j,i)
                  k = surf_lsm_v(l)%k(m)
-
-                 isurf          = isurf + 1
+                 isurf = isurf + 1
                  surfl(:,isurf) = (/inorth_l,k,j,i,m/)
+                 IF ( rad_angular_discretization ) THEN
+                    gridsurf(inorth_u,k,j,i) = isurf
+                 ENDIF
               ENDDO
 
@@ -5398 +5709 @@
               DO  m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
                  k = surf_usm_v(l)%k(m)
-
-                 isurf          = isurf + 1
+                 isurf = isurf + 1
                  surfl(:,isurf) = (/isouth_u,k,j,i,m/)
+                 IF ( rad_angular_discretization ) THEN
+                    gridsurf(isouth_u,k,j,i) = isurf
+                 ENDIF
               ENDDO
               DO  m = surf_lsm_v(l)%start_index(j,i), surf_lsm_v(l)%end_index(j,i)
                  k = surf_lsm_v(l)%k(m)
-
-                 isurf          = isurf + 1
+                 isurf = isurf + 1
                  surfl(:,isurf) = (/isouth_l,k,j,i,m/)
+                 IF ( rad_angular_discretization ) THEN
+                    gridsurf(isouth_u,k,j,i) = isurf
+                 ENDIF
               ENDDO
 
@@ -5412 +5727 @@
               DO  m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
                  k = surf_usm_v(l)%k(m)
-
-                 isurf          = isurf + 1
+                 isurf = isurf + 1
                  surfl(:,isurf) = (/ieast_u,k,j,i,m/)
+                 IF ( rad_angular_discretization ) THEN
+                    gridsurf(ieast_u,k,j,i) = isurf
+                 ENDIF
               ENDDO
               DO  m = surf_lsm_v(l)%start_index(j,i), surf_lsm_v(l)%end_index(j,i)
                  k = surf_lsm_v(l)%k(m)
-
-                 isurf          = isurf + 1
+                 isurf = isurf + 1
                  surfl(:,isurf) = (/ieast_l,k,j,i,m/)
+                 IF ( rad_angular_discretization ) THEN
+                    gridsurf(ieast_u,k,j,i) = isurf
+                 ENDIF
               ENDDO
 
@@ -5426 +5745 @@
               DO  m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
                  k = surf_usm_v(l)%k(m)
-
-                 isurf          = isurf + 1
+                 isurf = isurf + 1
                  surfl(:,isurf) = (/iwest_u,k,j,i,m/)
+                 IF ( rad_angular_discretization ) THEN
+                    gridsurf(iwest_u,k,j,i) = isurf
+                 ENDIF
               ENDDO
               DO  m = surf_lsm_v(l)%start_index(j,i), surf_lsm_v(l)%end_index(j,i)
                  k = surf_lsm_v(l)%k(m)
-
-                 isurf          = isurf + 1
+                 isurf = isurf + 1
                  surfl(:,isurf) = (/iwest_l,k,j,i,m/)
+                 IF ( rad_angular_discretization ) THEN
+                    gridsurf(iwest_u,k,j,i) = isurf
+                 ENDIF
               ENDDO
            ENDDO
        ENDDO
+!
+!--    Add local MRT boxes for specified number of levels
+       nmrtbl = 0
+       IF ( mrt_nlevels > 0 )  THEN
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
+!
+!--                Skip roof if requested
+                   IF ( mrt_skip_roof  .AND.  surf_usm_h%isroof_surf(m) )  CYCLE
+!
+!--                Cycle over specified no of levels
+                   nmrtbl = nmrtbl + mrt_nlevels
+                ENDDO
+!
+!--             Dtto for LSM
+                DO  m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
+                   nmrtbl = nmrtbl + mrt_nlevels
+                ENDDO
+             ENDDO
+          ENDDO
+
+          ALLOCATE( mrtbl(iz:ix,nmrtbl), mrtsky(nmrtbl), mrtskyt(nmrtbl), &
+                    mrtinsw(nmrtbl), mrtinlw(nmrtbl), mrt(nmrtbl) )
+
+          imrt = 0
+          DO  i = nxl, nxr
+             DO  j = nys, nyn
+                DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
+!
+!--                Skip roof if requested
+                   IF ( mrt_skip_roof  .AND.  surf_usm_h%isroof_surf(m) )  CYCLE
+!
+!--                Cycle over specified no of levels
+                   l = surf_usm_h%k(m)
+                   DO  k = l, l + mrt_nlevels - 1
+                      imrt = imrt + 1
+                      mrtbl(:,imrt) = (/k,j,i/)
+                   ENDDO
+                ENDDO
+!
+!--             Dtto for LSM
+                DO  m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
+                   l = surf_lsm_h%k(m)
+                   DO  k = l, l + mrt_nlevels - 1
+                      imrt = imrt + 1
+                      mrtbl(:,imrt) = (/k,j,i/)
+                   ENDDO
+                ENDDO
+             ENDDO
+          ENDDO
+       ENDIF
 
 !
@@ -5458 +5833 @@
 #if defined( __parallel )
        CALL MPI_Allgather(nsurfl,1,MPI_INTEGER,nsurfs,1,MPI_INTEGER,comm2d,ierr)
+       IF ( ierr /= 0 ) THEN
+         WRITE(9,*) 'Error MPI_AllGather1:', ierr, nsurfl, nsurfs
+         FLUSH(9)
+     ENDIF
+
 #else
        nsurfs(0) = nsurfl
@@ -5469 +5849 @@
        surfstart(numprocs) = k
        nsurf = k
-       ALLOCATE(surf(5,nsurf))
+       ALLOCATE(surf_l(5*nsurf))
+       surf(1:5,1:nsurf) => surf_l(1:5*nsurf)
 
 #if defined( __parallel )
-       CALL MPI_AllGatherv(surfl, nsurfl*5, MPI_INTEGER, surf, nsurfs*5, &
+       CALL MPI_AllGatherv(surfl_l, nsurfl*5, MPI_INTEGER, surf_l, nsurfs*5, &
            surfstart(0:numprocs-1)*5, MPI_INTEGER, comm2d, ierr)
+       IF ( ierr /= 0 ) THEN
+           WRITE(9,*) 'Error MPI_AllGatherv4:', ierr, SIZE(surfl_l), nsurfl*5, &
+                      SIZE(surf_l), nsurfs*5, surfstart(0:numprocs-1)*5
+           FLUSH(9)
+       ENDIF
 #else
        surf = surfl
@@ -5534 +5920 @@
         
         INTEGER(iwp)                                  :: i, j, k, d, ip, jp
-        INTEGER(iwp)                                  :: isvf, ksvf, icsf, kcsf, npcsfl, isvf_surflt, imrtt, imrtf, ipcgb
+        INTEGER(iwp)                                  :: isvf, ksvf, icsf, kcsf, npcsfl, isvf_surflt, imrt, imrtf, ipcgb
         INTEGER(iwp)                                  :: sd, td
         INTEGER(iwp)                                  :: iaz, izn      !< azimuth, zenith counters
@@ -5542 +5928 @@
         REAL(wp)                                      :: az1, az2      !< relative azimuth of section borders
         REAL(wp)                                      :: azmid         !< ray (center) azimuth
-        REAL(wp)                                      :: horizon       !< computed horizon height (tangent of elevation)
-        REAL(wp)                                      :: azen          !< zenith angle
         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           :: zbdry         !< zenith angle boundaries
         REAL(wp), DIMENSION(:), ALLOCATABLE           :: vffrac        !< view factor fractions for individual rays
+        REAL(wp), DIMENSION(:), ALLOCATABLE           :: vffrac0       !< dtto (original values)
         REAL(wp), DIMENSION(:), ALLOCATABLE           :: ztransp       !< array of transparency in z steps
+        INTEGER(iwp)                                  :: lowest_free_ray !< index into zdirs
+        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), DIMENSION(:,:), ALLOCATABLE     :: nzterrl
-        REAL(wp),     DIMENSION(:,:), ALLOCATABLE     :: csflt, pcsflt
-        INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE     :: kcsflt,kpcsflt
+        INTEGER(iwp)                                  :: udim
+        INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET:: nzterrl_l
+        INTEGER(iwp), DIMENSION(:,:), POINTER         :: nzterrl
+        REAL(wp),     DIMENSION(:), ALLOCATABLE,TARGET:: csflt_l, pcsflt_l
+        REAL(wp),     DIMENSION(:,:), POINTER         :: csflt, pcsflt
+        INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET:: kcsflt_l,kpcsflt_l
+        INTEGER(iwp), DIMENSION(:,:), POINTER         :: kcsflt,kpcsflt
         INTEGER(iwp), DIMENSION(:), ALLOCATABLE       :: icsflt,dcsflt,ipcsflt,dpcsflt
         REAL(wp), DIMENSION(3)                        :: uv
@@ -5561 +5958 @@
         INTEGER(idp)                                  :: ray_skip_maxdist, ray_skip_minval !< skipped raytracing counts
         INTEGER(iwp)                                  :: max_track_len !< maximum 2d track length
-        INTEGER(iwp)                                  :: win_lad, minfo
-        REAL(wp), DIMENSION(:,:,:), POINTER           :: lad_s_rma       !< fortran pointer, but lower bounds are 1
+        INTEGER(iwp)                                  :: minfo
+        REAL(wp), DIMENSION(:), POINTER               :: lad_s_rma       !< fortran 1D pointer
         TYPE(c_ptr)                                   :: lad_s_rma_p     !< allocated c pointer
 #if defined( __parallel )
@@ -5569 +5966 @@
 !    
         INTEGER(iwp), DIMENSION(0:svfnorm_report_num) :: svfnorm_counts
-!        CHARACTER(200)                                :: msg
+        CHARACTER(200)                                :: msg
 
 !--     calculation of the SVF
         CALL location_message( '    calculation of SVF and CSF', .TRUE. )
-!         CALL radiation_write_debug_log('Start calculation of SVF and CSF')
+        CALL radiation_write_debug_log('Start calculation of SVF and CSF')
 
 !--     precalculate face areas for different face directions using normal vector
@@ -5596 +5993 @@
             acsf => acsf1
         ENDIF
+        nmrtf = 0
+        IF ( mrt_nlevels > 0 )  THEN
+           nmrtfa = gasize
+           mmrtf = 1
+           ALLOCATE ( amrtf1(nmrtfa) )
+           amrtf => amrtf1
+        ENDIF
         ray_skip_maxdist = 0
         ray_skip_minval = 0
@@ -5602 +6006 @@
         ALLOCATE( nzterr(0:(nx+1)*(ny+1)-1) )
 #if defined( __parallel )
-        ALLOCATE( nzterrl(nys:nyn,nxl:nxr) )
+        !ALLOCATE( nzterrl(nys:nyn,nxl:nxr) )
+        ALLOCATE( nzterrl_l((nyn-nys+1)*(nxr-nxl+1)) )
+        nzterrl(nys:nyn,nxl:nxr) => nzterrl_l(1:(nyn-nys+1)*(nxr-nxl+1))
         nzterrl = get_topography_top_index( 's' )
-        CALL MPI_AllGather( nzterrl, nnx*nny, MPI_INTEGER, &
+        CALL MPI_AllGather( nzterrl_l, nnx*nny, MPI_INTEGER, &
                             nzterr, nnx*nny, MPI_INTEGER, comm2d, ierr )
-        DEALLOCATE(nzterrl)
+        IF ( ierr /= 0 ) THEN
+            WRITE(9,*) 'Error MPI_AllGather1:', ierr, SIZE(nzterrl_l), nnx*nny, &
+                       SIZE(nzterr), nnx*nny
+            FLUSH(9)
+        ENDIF
+        DEALLOCATE(nzterrl_l)
 #else
         nzterr = RESHAPE( get_topography_top_index( 's' ), (/(nx+1)*(ny+1)/) )
@@ -5621 +6032 @@
             CALL MPI_AllGather( pct, nnx*nny, MPI_INTEGER, &
                                 plantt, nnx*nny, MPI_INTEGER, comm2d, ierr )
-            
+            IF ( ierr /= 0 ) THEN
+                WRITE(9,*) 'Error MPI_AllGather2:', ierr, SIZE(pct), nnx*nny, &
+                           SIZE(plantt), nnx*nny
+                FLUSH(9)
+            ENDIF
+
 !--         temporary arrays storing values for csf calculation during raytracing
             ALLOCATE( lad_ip(maxboxesg) )
             ALLOCATE( lad_disp(maxboxesg) )
 
-            IF ( rma_lad_raytrace )  THEN
+            IF ( raytrace_mpi_rma )  THEN
                 ALLOCATE( lad_s_ray(maxboxesg) )
                 
                 ! set conditions for RMA communication
                 CALL MPI_Info_create(minfo, ierr)
+                IF ( ierr /= 0 ) THEN
+                    WRITE(9,*) 'Error MPI_Info_create2:', ierr
+                    FLUSH(9)
+                ENDIF
                 CALL MPI_Info_set(minfo, 'accumulate_ordering', '', ierr)
+                IF ( ierr /= 0 ) THEN
+                    WRITE(9,*) 'Error MPI_Info_set5:', ierr
+                    FLUSH(9)
+                ENDIF
                 CALL MPI_Info_set(minfo, 'accumulate_ops', 'same_op', ierr)
+                IF ( ierr /= 0 ) THEN
+                    WRITE(9,*) 'Error MPI_Info_set6:', ierr
+                    FLUSH(9)
+                ENDIF
                 CALL MPI_Info_set(minfo, 'same_size', 'true', ierr)
+                IF ( ierr /= 0 ) THEN
+                    WRITE(9,*) 'Error MPI_Info_set7:', ierr
+                    FLUSH(9)
+                ENDIF
                 CALL MPI_Info_set(minfo, 'same_disp_unit', 'true', ierr)
+                IF ( ierr /= 0 ) THEN
+                    WRITE(9,*) 'Error MPI_Info_set8:', ierr
+                    FLUSH(9)
+                ENDIF
 
 !--             Allocate and initialize the MPI RMA window
@@ -5643 +6079 @@
                 CALL MPI_Win_allocate(size_lad_rma, STORAGE_SIZE(1.0_wp)/8, minfo, comm2d, &
                                         lad_s_rma_p, win_lad, ierr)
-                CALL c_f_pointer(lad_s_rma_p, lad_s_rma, (/ nzp, nny, nnx /))
-                sub_lad(nzub:, nys:, nxl:) => lad_s_rma(:,:,:)
+                IF ( ierr /= 0 ) THEN
+                    WRITE(9,*) 'Error MPI_Win_allocate2:', ierr, size_lad_rma, &
+                                STORAGE_SIZE(1.0_wp)/8, win_lad
+                    FLUSH(9)
+                ENDIF
+                CALL c_f_pointer(lad_s_rma_p, lad_s_rma, (/ nzp*nny*nnx /))
+                sub_lad(nzub:nzpt, nys:nyn, nxl:nxr) => lad_s_rma(1:nzp*nny*nnx)
             ELSE
                 ALLOCATE(sub_lad(nzub:nzpt, nys:nyn, nxl:nxr))
@@ -5666 +6107 @@
 
 #if defined( __parallel )
-            IF ( rma_lad_raytrace )  THEN
+            IF ( raytrace_mpi_rma )  THEN
                 CALL MPI_Info_free(minfo, ierr)
+                IF ( ierr /= 0 ) THEN
+                    WRITE(9,*) 'Error MPI_Info_free2:', ierr
+                    FLUSH(9)
+                ENDIF
                 CALL MPI_Win_lock_all(0, win_lad, ierr)
+                IF ( ierr /= 0 ) THEN
+                    WRITE(9,*) 'Error MPI_Win_lock_all1:', ierr, win_lad
+                    FLUSH(9)
+                ENDIF
+                
             ELSE
                 ALLOCATE( sub_lad_g(0:(nx+1)*(ny+1)*nzp-1) )
                 CALL MPI_AllGather( sub_lad, nnx*nny*nzp, MPI_REAL, &
                                     sub_lad_g, nnx*nny*nzp, MPI_REAL, comm2d, ierr )
+                IF ( ierr /= 0 ) THEN
+                    WRITE(9,*) 'Error MPI_AllGather3:', ierr, SIZE(sub_lad), &
+                                nnx*nny*nzp, SIZE(sub_lad_g), nnx*nny*nzp
+                    FLUSH(9)
+                ENDIF
             ENDIF
 #endif
         ENDIF
 
-        IF ( mrt_factors )  THEN
-            OPEN(153, file='MRT_TARGETS', access='SEQUENTIAL', &
-                    action='READ', status='OLD', form='FORMATTED', err=524)
-            OPEN(154, file='MRT_FACTORS'//myid_char, access='DIRECT', recl=(5*4+2*8), &
-                    action='WRITE', status='REPLACE', form='UNFORMATTED', err=525)
-            imrtf = 1
-            DO
-                READ(153, *, end=526, err=524) imrtt, i, j, k
-                IF ( i < nxl  .OR.  i > nxr &
-                     .OR.  j < nys  .OR.  j > nyn ) CYCLE
-                ta = (/ REAL(k), REAL(j), REAL(i) /)
-
-                DO isurfs = 1, nsurf
-                    IF ( .NOT.  surface_facing(i, j, k, -1, &
-                        surf(ix, isurfs), surf(iy, isurfs), &
-                        surf(iz, isurfs), surf(id, isurfs)) )  THEN
-                        CYCLE
-                    ENDIF
-                      
-                    sd = surf(id, isurfs)
-                    sa = (/ REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd), &
-                            REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd), &
-                            REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd) /)
-
-!--                 unit vector source -> target
-                    uv = (/ (ta(1)-sa(1))*dz(1), (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)
-                    sqdist = SUM(uv(:)**2)
-                    uv = uv / SQRT(sqdist)
-
-!--                 irradiance factor - see svf. Here we consider that target face is always normal,
-!--                 i.e. the second dot product equals 1
-                    rirrf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) &
-                        / (pi * sqdist) * facearea(sd)
-
-!--                 raytrace while not creating any canopy sink factors
-                    CALL raytrace(sa, ta, isurfs, rirrf, 1._wp, .FALSE., &
-                            visible, transparency, win_lad)
-                    IF ( .NOT.  visible ) CYCLE
-
-                    !rsvf = rirrf * transparency
-                    WRITE(154, rec=imrtf, err=525) INT(imrtt, kind=4), &
-                        INT(surf(id, isurfs), kind=4), &
-                        INT(surf(iz, isurfs), kind=4), &
-                        INT(surf(iy, isurfs), kind=4), &
-                        INT(surf(ix, isurfs), kind=4), &
-                        REAL(rirrf, kind=8), REAL(transparency, kind=8)
-                    imrtf = imrtf + 1
-
-                ENDDO !< isurfs
-            ENDDO !< MRT_TARGETS record
-
-524         message_string = 'error reading file MRT_TARGETS'
-            CALL message( 'radiation_calc_svf', 'PA0524', 1, 2, 0, 6, 0 )
-
-525         message_string = 'error writing file MRT_FACTORS'//myid_char
-            CALL message( 'radiation_calc_svf', 'PA0525', 1, 2, 0, 6, 0 )
-
-526         CLOSE(153)
-            CLOSE(154)
-        ENDIF  !< mrt_factors
-        
+!--     prepare the MPI_Win for collecting the surface indices
+!--     from the reverse index arrays gridsurf from processors of target surfaces
+#if defined( __parallel )
+        IF ( rad_angular_discretization )  THEN
+!
+!--         raytrace_mpi_rma is asserted
+            CALL MPI_Win_lock_all(0, win_gridsurf, ierr)
+            IF ( ierr /= 0 ) THEN
+                WRITE(9,*) 'Error MPI_Win_lock_all2:', ierr, win_gridsurf
+                FLUSH(9)
+            ENDIF
+        ENDIF
+#endif
+
+
         !--Directions opposite to face normals are not even calculated,
         !--they must be preset to 0
@@ -5798 +6207 @@
             END SELECT
 
-            ALLOCATE ( zdirs(1:nzn), zbdry(0:nzn), vffrac(1:nzn), ztransp(1:nzn) )
+            ALLOCATE ( zdirs(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
+
+            itarg0 = 1
+            itarg1 = nzn
             zdirs(:) = (/( TAN(pi/2 - (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
-               !-- For horizontal target, vf fractions are constant per azimuth
-               vffrac(:) = (COS(2 * zbdry(0:nzn-1)) - COS(2 * zbdry(1:nzn))) / 2._wp / REAL(naz, wp)
-               !--sum of vffrac for all iaz equals 1, verified
+               vffrac(1:nzn) = (COS(2 * zbdry(0:nzn-1)) - COS(2 * zbdry(1:nzn))) / 2._wp / REAL(naz, wp)
+!
+!--            For horizontal target, vf fractions are constant per azimuth
+               DO iaz = 1, naz-1
+                  vffrac(iaz*nzn+1:(iaz+1)*nzn) = vffrac(1:nzn)
+               ENDDO
+!--            sum of whole vffrac equals 1, verified
             ENDIF
-
-            !--Calculate sky-view factor and direct solar visibility using 2D raytracing
+!
+!--         Calculate sky-view factor and direct solar visibility using 2D raytracing
             DO iaz = 1, naz
                azmid = az0 + (REAL(iaz, wp) - .5_wp) * azs
@@ -5814 +6232 @@
                   az1 = az2 - azs
                   !TODO precalculate after 1st line
-                  vffrac(:) = (SIN(az2) - SIN(az1))                           &
+                  vffrac(itarg0:itarg1) = (SIN(az2) - SIN(az1))               &
                               * (zbdry(1:nzn) - zbdry(0:nzn-1)                &
                                  + SIN(zbdry(0:nzn-1))*COS(zbdry(0:nzn-1))    &
                                  - SIN(zbdry(1:nzn))*COS(zbdry(1:nzn)))       &
                               / (2._wp * pi)
-                  !--sum of vffrac for all iaz equals 1, verified
+!--               sum of whole vffrac equals 1, verified
                ENDIF
                yxdir = (/ COS(azmid), SIN(azmid) /)
-               CALL raytrace_2d(ta, yxdir, zdirs,                              &
-                               surfstart(myid) + isurflt, facearea(td),        &
-                               vffrac, .TRUE., .FALSE., win_lad, horizon,       &
-                               ztransp)
-               
-
-               azen = pi/2 - ATAN(horizon)
-               IF ( td == iup_u  .OR.  td == iup_l )  THEN
-                  azen = MIN(azen, pi/2) !only above horizontal direction
-                  skyvf(isurflt) = skyvf(isurflt) + (1._wp - COS(2*azen)) /   &
-                     (2._wp * raytrace_discrete_azims)
-               ELSE
-                  skyvf(isurflt) = skyvf(isurflt) + (SIN(az2) - SIN(az1)) *   &
-                              (azen - SIN(azen)*COS(azen)) / (2._wp*pi)
-               ENDIF
-               skyvft(isurflt) = skyvft(isurflt) + SUM(ztransp(:) * vffrac(:))
+               CALL raytrace_2d(ta, yxdir, nzn, zdirs,                        &
+                                    surfstart(myid) + isurflt, facearea(td),  &
+                                    vffrac(itarg0:itarg1), .TRUE., .TRUE.,    &
+                                    .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
+               skyvf(isurflt) = skyvf(isurflt) + &
+                                SUM(vffrac(itarg0:itarg0+lowest_free_ray-1))
+               skyvft(isurflt) = skyvft(isurflt) + &
+                                 SUM(ztransp(itarg0:itarg0+lowest_free_ray-1) &
+                                             * vffrac(itarg0:itarg0+lowest_free_ray-1))
  
-               !--Save direct solar transparency
+!--            Save direct solar transparency
                j = MODULO(NINT(azmid/                                          &
                                (2._wp*pi)*raytrace_discrete_azims-.5_wp, iwp), &
@@ -5846 +6261 @@
                DO k = 1, raytrace_discrete_elevs/2
                   i = dsidir_rev(k-1, j)
-                  IF ( i /= -1 )  dsitrans(isurflt, i) = ztransp(k)
+                  IF ( i /= -1  .AND.  k <= lowest_free_ray )  &
+                     dsitrans(isurflt, i) = ztransp(itarg0+k-1)
                ENDDO
+
+               !
+               !--Advance itarget indices
+               itarg0 = itarg1 + 1
+               itarg1 = itarg1 + nzn
             ENDDO
 
-            DEALLOCATE ( zdirs, zbdry, vffrac, ztransp )
-!
-!--         Following calculations only required for surface_reflections
-            IF ( surface_reflections )  THEN
-
-               DO  isurfs = 1, nsurf
-                  IF ( .NOT.  surface_facing(surfl(ix, isurflt), surfl(iy, isurflt), &
-                     surfl(iz, isurflt), surfl(id, isurflt), &
-                     surf(ix, isurfs), surf(iy, isurfs), &
-                     surf(iz, isurfs), surf(id, isurfs)) )  THEN
-                     CYCLE
+            IF ( rad_angular_discretization )  THEN
+!--            sort itarget by face id
+               CALL quicksort_itarget(itarget,vffrac,ztransp,1,nzn*naz)
+!
+!--            find the first valid position
+               itarg0 = 1
+               DO WHILE ( itarg0 <= nzn*naz )
+                  IF ( itarget(itarg0) /= -1 )  EXIT
+                  itarg0 = itarg0 + 1
+               ENDDO
+
+               DO  i = itarg0, nzn*naz
+!
+!--               For duplicate values, only sum up vf fraction value
+                  IF ( i < nzn*naz )  THEN
+                     IF ( itarget(i+1) == itarget(i) )  THEN
+                        vffrac(i+1) = vffrac(i+1) + vffrac(i)
+                        CYCLE
+                     ENDIF
                   ENDIF
-                  
-                  sd = surf(id, isurfs)
-                  sa = (/ REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd),  &
-                          REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd),  &
-                          REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd)  /)
-
-!--               unit vector source -> target
-                  uv = (/ (ta(1)-sa(1))*dz(1), (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)
-                  sqdist = SUM(uv(:)**2)
-                  uv = uv / SQRT(sqdist)
-
-!--               reject raytracing above max distance
-                  IF ( SQRT(sqdist) > max_raytracing_dist ) THEN
-                     ray_skip_maxdist = ray_skip_maxdist + 1
-                     CYCLE
-                  ENDIF
-                 
-!--               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)
-
-!--               reject raytracing for potentially too small view factor values
-                  IF ( rirrf < min_irrf_value ) THEN
-                      ray_skip_minval = ray_skip_minval + 1
-                      CYCLE
-                  ENDIF
-
-!--               raytrace + process plant canopy sinks within
-                  CALL raytrace(sa, ta, isurfs, rirrf, facearea(td), .TRUE., &
-                                visible, transparency, win_lad)
-
-                  IF ( .NOT.  visible ) CYCLE
-                 ! rsvf = rirrf * transparency
-
+!
 !--               write to the svf array
                   nsvfl = nsvfl + 1
 !--               check dimmension of asvf array and enlarge it if needed
                   IF ( nsvfla < nsvfl )  THEN
-                     k = nsvfla * 2
+                     k = CEILING(REAL(nsvfla, kind=wp) * grow_factor)
                      IF ( msvf == 0 )  THEN
                         msvf = 1
@@ -5917 +6311 @@
                      ENDIF
 
-!                      WRITE(msg,'(A,3I12)') 'Grow asvf:',nsvfl,nsvfla,k
-!                      CALL radiation_write_debug_log( msg )
+                     WRITE (msg,'(A,3I12)') 'Grow asvf:',nsvfl,nsvfla,k
+                     CALL radiation_write_debug_log( msg )
+                     
+                     nsvfla = k
+                  ENDIF
+!--               write svf values into the array
+                  asvf(nsvfl)%isurflt = isurflt
+                  asvf(nsvfl)%isurfs = itarget(i)
+                  asvf(nsvfl)%rsvf = vffrac(i)
+                  asvf(nsvfl)%rtransp = ztransp(i)
+               END DO
+
+            ENDIF ! rad_angular_discretization
+
+            DEALLOCATE ( zdirs, zbdry, vffrac, ztransp, itarget ) !FIXME itarget shall be allocated only
+                                                                  !in case of rad_angular_discretization
+!
+!--         Following calculations only required for surface_reflections
+            IF ( surface_reflections  .AND.  .NOT. rad_angular_discretization )  THEN
+
+               DO  isurfs = 1, nsurf
+                  IF ( .NOT.  surface_facing(surfl(ix, isurflt), surfl(iy, isurflt), &
+                     surfl(iz, isurflt), surfl(id, isurflt), &
+                     surf(ix, isurfs), surf(iy, isurfs), &
+                     surf(iz, isurfs), surf(id, isurfs)) )  THEN
+                     CYCLE
+                  ENDIF
+                  
+                  sd = surf(id, isurfs)
+                  sa = (/ REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd),  &
+                          REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd),  &
+                          REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd)  /)
+
+!--               unit vector source -> target
+                  uv = (/ (ta(1)-sa(1))*dz(1), (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)
+                  sqdist = SUM(uv(:)**2)
+                  uv = uv / SQRT(sqdist)
+
+!--               reject raytracing above max distance
+                  IF ( SQRT(sqdist) > max_raytracing_dist ) THEN
+                     ray_skip_maxdist = ray_skip_maxdist + 1
+                     CYCLE
+                  ENDIF
+                 
+!--               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)
+
+!--               reject raytracing for potentially too small view factor values
+                  IF ( rirrf < min_irrf_value ) THEN
+                      ray_skip_minval = ray_skip_minval + 1
+                      CYCLE
+                  ENDIF
+
+!--               raytrace + process plant canopy sinks within
+                  CALL raytrace(sa, ta, isurfs, rirrf, facearea(td), .TRUE., &
+                                visible, transparency)
+
+                  IF ( .NOT.  visible ) CYCLE
+                 ! rsvf = rirrf * transparency
+
+!--               write to the svf array
+                  nsvfl = nsvfl + 1
+!--               check dimmension of asvf array and enlarge it if needed
+                  IF ( nsvfla < nsvfl )  THEN
+                     k = CEILING(REAL(nsvfla, kind=wp) * grow_factor)
+                     IF ( msvf == 0 )  THEN
+                        msvf = 1
+                        ALLOCATE( asvf1(k) )
+                        asvf => asvf1
+                        asvf1(1:nsvfla) = asvf2
+                        DEALLOCATE( asvf2 )
+                     ELSE
+                        msvf = 0
+                        ALLOCATE( asvf2(k) )
+                        asvf => asvf2
+                        asvf2(1:nsvfla) = asvf1
+                        DEALLOCATE( asvf1 )
+                     ENDIF
+
+                     WRITE(msg,'(A,3I12)') 'Grow asvf:',nsvfl,nsvfla,k
+                     CALL radiation_write_debug_log( msg )
                      
                      nsvfla = k
@@ -5931 +6407 @@
         ENDDO
 
-        !--Raytrace to canopy boxes to fill dsitransc TODO optimize
-        !--
-        dsitransc(:,:) = -999._wp !FIXME
+!--
+!--     Raytrace to canopy boxes to fill dsitransc TODO optimize
+        dsitransc(:,:) = 0._wp
         az0 = 0._wp
         naz = raytrace_discrete_azims
@@ -5940 +6416 @@
         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), vffrac(1:nzn), ztransp(1:nzn), &
+               itarget(1:nzn) )
         zdirs(:) = (/( TAN(pi/2 - (zn0+(REAL(izn,wp)-.5_wp)*zns)), izn=1, nzn )/)
         vffrac(:) = 0._wp
 
-        DO ipcgb = 1, npcbl
+        DO  ipcgb = 1, npcbl
            ta = (/ REAL(pcbl(iz, ipcgb), wp),  &
                    REAL(pcbl(iy, ipcgb), wp),  &
                    REAL(pcbl(ix, ipcgb), wp) /)
-           !--Calculate sky-view factor and direct solar visibility using 2D raytracing
-           DO iaz = 1, naz
+!--        Calculate direct solar visibility using 2D raytracing
+           DO  iaz = 1, naz
               azmid = az0 + (REAL(iaz, wp) - .5_wp) * azs
               yxdir = (/ COS(azmid), SIN(azmid) /)
-              CALL raytrace_2d(ta, yxdir, zdirs,     &
-                                   -999, -999._wp, vffrac, .FALSE., .TRUE., &
-                                   win_lad, horizon, ztransp)
-
-              !--Save direct solar transparency
+              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
+
+!--           Save direct solar transparency
               j = MODULO(NINT(azmid/                                         &
                              (2._wp*pi)*raytrace_discrete_azims-.5_wp, iwp), &
                          raytrace_discrete_azims)
-              DO k = 1, raytrace_discrete_elevs/2
+              DO  k = 1, raytrace_discrete_elevs/2
                  i = dsidir_rev(k-1, j)
-                 IF ( i /= -1 )  dsitransc(ipcgb, i) = ztransp(k)
+                 IF ( i /= -1  .AND.  k <= lowest_free_ray ) &
+                    dsitransc(ipcgb, i) = ztransp(k)
               ENDDO
            ENDDO
         ENDDO
-        DEALLOCATE ( zdirs, vffrac, ztransp )
-
-!         CALL radiation_write_debug_log( 'End of calculation SVF' )
-!         WRITE(msg, *) 'Raytracing skipped for maximum distance of ', &
-!             max_raytracing_dist, ' m on ', ray_skip_maxdist, ' pairs.'
-!         CALL radiation_write_debug_log( msg )
-!         WRITE(msg, *) 'Raytracing skipped for minimum potential value of ', &
-!             min_irrf_value , ' on ', ray_skip_minval, ' pairs.'
-!         CALL radiation_write_debug_log( msg )
+        DEALLOCATE ( zdirs, vffrac, ztransp, itarget )
+!--
+!--     Raytrace to MRT boxes
+        IF ( nmrtbl > 0 )  THEN
+           mrtdsit(:,:) = 0._wp
+           mrtsky(:) = 0._wp
+           mrtskyt(:) = 0._wp
+           az0 = 0._wp
+           naz = raytrace_discrete_azims
+           azs = 2._wp * pi / REAL(naz, wp)
+           zn0 = 0._wp
+           nzn = raytrace_discrete_elevs
+           zns = pi / REAL(nzn, wp)
+           ALLOCATE ( zdirs(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 )/)
+           zbdry(:) = (/( zn0+REAL(izn,wp)*zns, izn=0, nzn )/)
+           vffrac0(:) = (COS(zbdry(0:nzn-1)) - COS(zbdry(1:nzn))) / 2._wp / REAL(naz, wp)
+
+           DO  imrt = 1, nmrtbl
+              ta = (/ REAL(mrtbl(iz, imrt), wp),  &
+                      REAL(mrtbl(iy, imrt), wp),  &
+                      REAL(mrtbl(ix, imrt), wp) /)
+!
+!--           vf fractions are constant per azimuth
+              DO iaz = 0, naz-1
+                 vffrac(iaz*nzn+1:(iaz+1)*nzn) = vffrac0(:)
+              ENDDO
+!--           sum of whole vffrac equals 1, verified
+              itarg0 = 1
+              itarg1 = nzn
+!
+!--           Calculate sky-view factor and direct solar visibility using 2D raytracing
+              DO  iaz = 1, naz
+                 azmid = az0 + (REAL(iaz, wp) - .5_wp) * azs
+                 CALL raytrace_2d(ta, (/ COS(azmid), SIN(azmid) /), 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
+
+!--              Sky view factors for MRT
+                 mrtsky(imrt) = mrtsky(imrt) + &
+                                  SUM(vffrac(itarg0:itarg0+lowest_free_ray-1))
+                 mrtskyt(imrt) = mrtskyt(imrt) + &
+                                   SUM(ztransp(itarg0:itarg0+lowest_free_ray-1) &
+                                               * vffrac(itarg0:itarg0+lowest_free_ray-1))
+!--              Direct solar transparency for MRT
+                 j = MODULO(NINT(azmid/                                         &
+                                (2._wp*pi)*raytrace_discrete_azims-.5_wp, iwp), &
+                            raytrace_discrete_azims)
+                 DO  k = 1, raytrace_discrete_elevs/2
+                    i = dsidir_rev(k-1, j)
+                    IF ( i /= -1  .AND.  k <= lowest_free_ray ) &
+                       mrtdsit(imrt, i) = ztransp(itarg0+k-1)
+                 ENDDO
+!
+!--              Advance itarget indices
+                 itarg0 = itarg1 + 1
+                 itarg1 = itarg1 + nzn
+              ENDDO
+
+!--           sort itarget by face id
+              CALL quicksort_itarget(itarget,vffrac,ztransp,1,nzn*naz)
+!
+!--           find the first valid position
+              itarg0 = 1
+              DO WHILE ( itarg0 <= nzn*naz )
+                 IF ( itarget(itarg0) /= -1 )  EXIT
+                 itarg0 = itarg0 + 1
+              ENDDO
+
+              DO  i = itarg0, nzn*naz
+!
+!--              For duplicate values, only sum up vf fraction value
+                 IF ( i < nzn*naz )  THEN
+                    IF ( itarget(i+1) == itarget(i) )  THEN
+                       vffrac(i+1) = vffrac(i+1) + vffrac(i)
+                       CYCLE
+                    ENDIF
+                 ENDIF
+!
+!--              write to the mrtf array
+                 nmrtf = nmrtf + 1
+!--              check dimmension of mrtf array and enlarge it if needed
+                 IF ( nmrtfa < nmrtf )  THEN
+                    k = CEILING(REAL(nmrtfa, kind=wp) * grow_factor)
+                    IF ( mmrtf == 0 )  THEN
+                       mmrtf = 1
+                       ALLOCATE( amrtf1(k) )
+                       amrtf => amrtf1
+                       amrtf1(1:nmrtfa) = amrtf2
+                       DEALLOCATE( amrtf2 )
+                    ELSE
+                       mmrtf = 0
+                       ALLOCATE( amrtf2(k) )
+                       amrtf => amrtf2
+                       amrtf2(1:nmrtfa) = amrtf1
+                       DEALLOCATE( amrtf1 )
+                    ENDIF
+
+                    WRITE (msg,'(A,3I12)') 'Grow amrtf:', nmrtf, nmrtfa, k
+                    CALL radiation_write_debug_log( msg )
+
+                    nmrtfa = k
+                 ENDIF
+!--              write mrtf values into the array
+                 amrtf(nmrtf)%isurflt = imrt
+                 amrtf(nmrtf)%isurfs = itarget(i)
+                 amrtf(nmrtf)%rsvf = vffrac(i)
+                 amrtf(nmrtf)%rtransp = ztransp(i)
+              ENDDO ! itarg
+
+           ENDDO ! imrt
+           DEALLOCATE ( zdirs, zbdry, vffrac, vffrac0, ztransp, itarget )
+!
+!--        Move MRT factors to final arrays
+           ALLOCATE ( mrtf(nmrtf), mrtft(nmrtf), mrtfsurf(2,nmrtf) )
+           DO  imrtf = 1, nmrtf
+              mrtf(imrtf) = amrtf(imrtf)%rsvf
+              mrtft(imrtf) = amrtf(imrtf)%rsvf * amrtf(imrtf)%rtransp
+              mrtfsurf(:,imrtf) = (/amrtf(imrtf)%isurflt, amrtf(imrtf)%isurfs /)
+           ENDDO
+           IF ( ALLOCATED(amrtf1) )  DEALLOCATE( amrtf1 )
+           IF ( ALLOCATED(amrtf2) )  DEALLOCATE( amrtf2 )
+        ENDIF ! nmrtbl > 0
+
+        IF ( rad_angular_discretization )  THEN
+#if defined( __parallel )
+!--        finalize MPI_RMA communication established to get global index of the surface from grid indices
+!--        flush all MPI window pending requests
+           CALL MPI_Win_flush_all(win_gridsurf, ierr)
+           IF ( ierr /= 0 ) THEN
+               WRITE(9,*) 'Error MPI_Win_flush_all1:', ierr, win_gridsurf
+               FLUSH(9)
+           ENDIF
+!--        unlock MPI window
+           CALL MPI_Win_unlock_all(win_gridsurf, ierr)
+           IF ( ierr /= 0 ) THEN
+               WRITE(9,*) 'Error MPI_Win_unlock_all1:', ierr, win_gridsurf
+               FLUSH(9)
+           ENDIF
+!--        free MPI window
+           CALL MPI_Win_free(win_gridsurf, ierr)
+           IF ( ierr /= 0 ) THEN
+               WRITE(9,*) 'Error MPI_Win_free1:', ierr, win_gridsurf
+               FLUSH(9)
+           ENDIF
+#else
+           DEALLOCATE ( gridsurf )
+#endif
+        ENDIF
+
+        CALL radiation_write_debug_log( 'End of calculation SVF' )
+        WRITE(msg, *) 'Raytracing skipped for maximum distance of ', &
+           max_raytracing_dist, ' m on ', ray_skip_maxdist, ' pairs.'
+        CALL radiation_write_debug_log( msg )
+        WRITE(msg, *) 'Raytracing skipped for minimum potential value of ', &
+           min_irrf_value , ' on ', ray_skip_minval, ' pairs.'
+        CALL radiation_write_debug_log( msg )
 
         CALL location_message( '    waiting for completion of SVF and CSF calculation in all processes', .TRUE. )
@@ -5983 +6616 @@
 !--         finalize mpi_rma communication and deallocate temporary arrays
 #if defined( __parallel )
-            IF ( rma_lad_raytrace )  THEN
+            IF ( raytrace_mpi_rma )  THEN
                 CALL MPI_Win_flush_all(win_lad, ierr)
+                IF ( ierr /= 0 ) THEN
+                    WRITE(9,*) 'Error MPI_Win_flush_all2:', ierr, win_lad
+                    FLUSH(9)
+                ENDIF
 !--             unlock MPI window
                 CALL MPI_Win_unlock_all(win_lad, ierr)
+                IF ( ierr /= 0 ) THEN
+                    WRITE(9,*) 'Error MPI_Win_unlock_all2:', ierr, win_lad
+                    FLUSH(9)
+                ENDIF
 !--             free MPI window
                 CALL MPI_Win_free(win_lad, ierr)
-                
+                IF ( ierr /= 0 ) THEN
+                    WRITE(9,*) 'Error MPI_Win_free2:', ierr, win_lad
+                    FLUSH(9)
+                ENDIF
 !--             deallocate temporary arrays storing values for csf calculation during raytracing
                 DEALLOCATE( lad_s_ray )
-!--             sub_lad is the pointer to lad_s_rma in case of rma_lad_raytrace
+!--             sub_lad is the pointer to lad_s_rma in case of raytrace_mpi_rma
 !--             and must not be deallocated here
             ELSE
@@ -6009 +6653 @@
         CALL location_message( '    calculation of the complete SVF array', .TRUE. )
 
-!         CALL radiation_write_debug_log( 'Start SVF sort' )
-!--     sort svf ( a version of quicksort )
-        CALL quicksort_svf(asvf,1,nsvfl)
-
-        !< load svf from the structure array to plain arrays
-!         CALL radiation_write_debug_log( 'Load svf from the structure array to plain arrays' )
-        ALLOCATE( svf(ndsvf,nsvfl) )
-        ALLOCATE( svfsurf(idsvf,nsvfl) )
-        svfnorm_counts(:) = 0._wp
-        isurflt_prev = -1
-        ksvf = 1
-        svfsum = 0._wp
-        DO isvf = 1, nsvfl
-!--         normalize svf per target face
-            IF ( asvf(ksvf)%isurflt /= isurflt_prev )  THEN
-                IF ( isurflt_prev /= -1  .AND.  svfsum /= 0._wp )  THEN
-                    !< update histogram of logged svf normalization values
-                    i = searchsorted(svfnorm_report_thresh, svfsum / (1._wp-skyvf(isurflt_prev)))
-                    svfnorm_counts(i) = svfnorm_counts(i) + 1
-
-                    svf(1, isvf_surflt:isvf-1) = svf(1, isvf_surflt:isvf-1) / svfsum * (1._wp-skyvf(isurflt_prev))
-                ENDIF
-                isurflt_prev = asvf(ksvf)%isurflt
-                isvf_surflt = isvf
-                svfsum = asvf(ksvf)%rsvf !?? / asvf(ksvf)%rtransp
-            ELSE
-                svfsum = svfsum + asvf(ksvf)%rsvf !?? / asvf(ksvf)%rtransp
-            ENDIF
-
-            svf(:, isvf) = (/ asvf(ksvf)%rsvf, asvf(ksvf)%rtransp /)
-            svfsurf(:, isvf) = (/ asvf(ksvf)%isurflt, asvf(ksvf)%isurfs /)
-
-!--         next element
-            ksvf = ksvf + 1
-        ENDDO
-
-        IF ( isurflt_prev /= -1  .AND.  svfsum /= 0._wp )  THEN
-            i = searchsorted(svfnorm_report_thresh, svfsum / (1._wp-skyvf(isurflt_prev)))
-            svfnorm_counts(i) = svfnorm_counts(i) + 1
-
-            svf(1, isvf_surflt:nsvfl) = svf(1, isvf_surflt:nsvfl) / svfsum * (1._wp-skyvf(isurflt_prev))
-        ENDIF
-        !TODO we should be able to deallocate skyvf, from now on we only need skyvft
+        IF ( rad_angular_discretization )  THEN
+           CALL radiation_write_debug_log( 'Load svf from the structure array to plain arrays' )
+           ALLOCATE( svf(ndsvf,nsvfl) )
+           ALLOCATE( svfsurf(idsvf,nsvfl) )
+
+           DO isvf = 1, nsvfl
+               svf(:, isvf) = (/ asvf(isvf)%rsvf, asvf(isvf)%rtransp /)
+               svfsurf(:, isvf) = (/ asvf(isvf)%isurflt, asvf(isvf)%isurfs /)
+           ENDDO
+        ELSE
+           CALL radiation_write_debug_log( 'Start SVF sort' )
+!--        sort svf ( a version of quicksort )
+           CALL quicksort_svf(asvf,1,nsvfl)
+
+           !< load svf from the structure array to plain arrays
+           CALL radiation_write_debug_log( 'Load svf from the structure array to plain arrays' )
+           ALLOCATE( svf(ndsvf,nsvfl) )
+           ALLOCATE( svfsurf(idsvf,nsvfl) )
+           svfnorm_counts(:) = 0._wp
+           isurflt_prev = -1
+           ksvf = 1
+           svfsum = 0._wp
+           DO isvf = 1, nsvfl
+!--            normalize svf per target face
+               IF ( asvf(ksvf)%isurflt /= isurflt_prev )  THEN
+                   IF ( isurflt_prev /= -1  .AND.  svfsum /= 0._wp )  THEN
+                       !< update histogram of logged svf normalization values
+                       i = searchsorted(svfnorm_report_thresh, svfsum / (1._wp-skyvf(isurflt_prev)))
+                       svfnorm_counts(i) = svfnorm_counts(i) + 1
+
+                       svf(1, isvf_surflt:isvf-1) = svf(1, isvf_surflt:isvf-1) / svfsum * (1._wp-skyvf(isurflt_prev))
+                   ENDIF
+                   isurflt_prev = asvf(ksvf)%isurflt
+                   isvf_surflt = isvf
+                   svfsum = asvf(ksvf)%rsvf !?? / asvf(ksvf)%rtransp
+               ELSE
+                   svfsum = svfsum + asvf(ksvf)%rsvf !?? / asvf(ksvf)%rtransp
+               ENDIF
+
+               svf(:, isvf) = (/ asvf(ksvf)%rsvf, asvf(ksvf)%rtransp /)
+               svfsurf(:, isvf) = (/ asvf(ksvf)%isurflt, asvf(ksvf)%isurfs /)
+
+!--            next element
+               ksvf = ksvf + 1
+           ENDDO
+
+           IF ( isurflt_prev /= -1  .AND.  svfsum /= 0._wp )  THEN
+               i = searchsorted(svfnorm_report_thresh, svfsum / (1._wp-skyvf(isurflt_prev)))
+               svfnorm_counts(i) = svfnorm_counts(i) + 1
+
+               svf(1, isvf_surflt:nsvfl) = svf(1, isvf_surflt:nsvfl) / svfsum * (1._wp-skyvf(isurflt_prev))
+           ENDIF
+           WRITE(9, *) 'SVF normalization histogram:', svfnorm_counts,  &
+                       'on thresholds:', svfnorm_report_thresh(1:svfnorm_report_num), '(val < thresh <= val)'
+           !TODO we should be able to deallocate skyvf, from now on we only need skyvft
+        ENDIF ! rad_angular_discretization
 
 !--     deallocate temporary asvf array
@@ -6067 +6724 @@
 
             CALL location_message( '    calculation of the complete CSF array', .TRUE. )
-!             CALL radiation_write_debug_log( 'Calculation of the complete CSF array' )
+            CALL radiation_write_debug_log( 'Calculation of the complete CSF array' )
 !--         sort and merge csf for the last time, keeping the array size to minimum
             CALL merge_and_grow_csf(-1)
@@ -6073 +6730 @@
 !--         aggregate csb among processors
 !--         allocate necessary arrays
-            ALLOCATE( csflt(ndcsf,max(ncsfl,ndcsf)) )
-            ALLOCATE( kcsflt(kdcsf,max(ncsfl,kdcsf)) )
+            udim = max(ncsfl,1)
+            ALLOCATE( csflt_l(ndcsf*udim) )
+            csflt(1:ndcsf,1:udim) => csflt_l(1:ndcsf*udim)
+            ALLOCATE( kcsflt_l(kdcsf*udim) )
+            kcsflt(1:kdcsf,1:udim) => kcsflt_l(1:kdcsf*udim)
             ALLOCATE( icsflt(0:numprocs-1) )
             ALLOCATE( dcsflt(0:numprocs-1) )
@@ -6108 +6768 @@
 !--             fill out real values of rsvf, rtransp
                 csflt(1,kcsf) = acsf(kcsf)%rsvf
-                csflt(2,kcsf) = acsf(kcsf)%rtransp
 !--             fill out integer values of itz,ity,itx,isurfs
                 kcsflt(1,kcsf) = acsf(kcsf)%itz
@@ -6138 +6797 @@
 !--         scatter and gather the number of elements to and from all processor
 !--         and calculate displacements
-!             CALL radiation_write_debug_log( 'Scatter and gather the number of elements to and from all processor' )
+            CALL radiation_write_debug_log( 'Scatter and gather the number of elements to and from all processor' )
             CALL MPI_AlltoAll(icsflt,1,MPI_INTEGER,ipcsflt,1,MPI_INTEGER,comm2d, ierr)
-            
+            IF ( ierr /= 0 ) THEN
+                WRITE(9,*) 'Error MPI_AlltoAll1:', ierr, SIZE(icsflt), SIZE(ipcsflt)
+                FLUSH(9)
+            ENDIF
+
             npcsfl = SUM(ipcsflt)
             d = 0
@@ -6149 +6812 @@
 
 !--         exchange csf fields between processors
-!             CALL radiation_write_debug_log( 'Exchange csf fields between processors' )
-            ALLOCATE( pcsflt(ndcsf,max(npcsfl,ndcsf)) )
-            ALLOCATE( kpcsflt(kdcsf,max(npcsfl,kdcsf)) )
-            CALL MPI_AlltoAllv(csflt, ndcsf*icsflt, ndcsf*dcsflt, MPI_REAL, &
-                pcsflt, ndcsf*ipcsflt, ndcsf*dpcsflt, MPI_REAL, comm2d, ierr)
-            CALL MPI_AlltoAllv(kcsflt, kdcsf*icsflt, kdcsf*dcsflt, MPI_INTEGER, &
-                kpcsflt, kdcsf*ipcsflt, kdcsf*dpcsflt, MPI_INTEGER, comm2d, ierr)
+            CALL radiation_write_debug_log( 'Exchange csf fields between processors' )
+            udim = max(npcsfl,1)
+            ALLOCATE( pcsflt_l(ndcsf*udim) )
+            pcsflt(1:ndcsf,1:udim) => pcsflt_l(1:ndcsf*udim)
+            ALLOCATE( kpcsflt_l(kdcsf*udim) )
+            kpcsflt(1:kdcsf,1:udim) => kpcsflt_l(1:kdcsf*udim)
+            CALL MPI_AlltoAllv(csflt_l, ndcsf*icsflt, ndcsf*dcsflt, MPI_REAL, &
+                pcsflt_l, ndcsf*ipcsflt, ndcsf*dpcsflt, MPI_REAL, comm2d, ierr)
+            IF ( ierr /= 0 ) THEN
+                WRITE(9,*) 'Error MPI_AlltoAllv1:', ierr, SIZE(ipcsflt), ndcsf*icsflt, &
+                            ndcsf*dcsflt, SIZE(pcsflt_l),ndcsf*ipcsflt, ndcsf*dpcsflt
+                FLUSH(9)
+            ENDIF
+
+            CALL MPI_AlltoAllv(kcsflt_l, kdcsf*icsflt, kdcsf*dcsflt, MPI_INTEGER, &
+                kpcsflt_l, kdcsf*ipcsflt, kdcsf*dpcsflt, MPI_INTEGER, comm2d, ierr)
+            IF ( ierr /= 0 ) THEN
+                WRITE(9,*) 'Error MPI_AlltoAllv2:', ierr, SIZE(kcsflt_l),kdcsf*icsflt, &
+                           kdcsf*dcsflt, SIZE(kpcsflt_l), kdcsf*ipcsflt, kdcsf*dpcsflt
+                FLUSH(9)
+            ENDIF
             
 #else
@@ -6166 +6843 @@
 
 !--         deallocate temporary arrays
-            DEALLOCATE( csflt )
-            DEALLOCATE( kcsflt )
+            DEALLOCATE( csflt_l )
+            DEALLOCATE( kcsflt_l )
             DEALLOCATE( icsflt )
             DEALLOCATE( dcsflt )
@@ -6174 +6851 @@
 
 !--         sort csf ( a version of quicksort )
-!             CALL radiation_write_debug_log( 'Sort csf' )
+            CALL radiation_write_debug_log( 'Sort csf' )
             CALL quicksort_csf2(kpcsflt, pcsflt, 1, npcsfl)
 
 !--         aggregate canopy sink factor records with identical box & source
 !--         againg across all values from all processors
-!             CALL radiation_write_debug_log( 'Aggregate canopy sink factor records with identical box' )
+            CALL radiation_write_debug_log( 'Aggregate canopy sink factor records with identical box' )
 
             IF ( npcsfl > 0 )  THEN
@@ -6190 +6867 @@
                          kpcsflt(1,icsf) == kpcsflt(1,icsf+1)  .AND.  &
                          kpcsflt(4,icsf) == kpcsflt(4,icsf+1) )  THEN
-!--                     We could simply take either first or second rtransp, both are valid. As a very simple heuristic about which ray
-!--                     probably passes nearer the center of the target box, we choose DIF from the entry with greater CSF, since that
-!--                     might mean that the traced beam passes longer through the canopy box.
-                        IF ( pcsflt(1,kcsf) < pcsflt(1,icsf+1) )  THEN
-                            pcsflt(2,kcsf) = pcsflt(2,icsf+1)
-                        ENDIF
+
                         pcsflt(1,kcsf) = pcsflt(1,kcsf) + pcsflt(1,icsf+1)
 
@@ -6224 +6896 @@
             
 !--         deallocation of temporary arrays
-            DEALLOCATE( pcsflt )
-            DEALLOCATE( kpcsflt )
-!             CALL radiation_write_debug_log( 'End of aggregate csf' )
+            IF ( npcbl > 0 )  DEALLOCATE( gridpcbl )
+            DEALLOCATE( pcsflt_l )
+            DEALLOCATE( kpcsflt_l )
+            CALL radiation_write_debug_log( 'End of aggregate csf' )
             
         ENDIF
@@ -6233 +6906 @@
         CALL MPI_BARRIER( comm2d, ierr )
 #endif
-!         CALL radiation_write_debug_log( 'End of radiation_calc_svf (after mpi_barrier)' )
+        CALL radiation_write_debug_log( 'End of radiation_calc_svf (after mpi_barrier)' )
 
         RETURN
@@ -6261 +6934 @@
 !>    doesn't need to be the same in all three dimensions).
 !------------------------------------------------------------------------------!
-    SUBROUTINE raytrace(src, targ, isrc, rirrf, atarg, create_csf, visible, transparency, win_lad)
+    SUBROUTINE raytrace(src, targ, isrc, rirrf, atarg, create_csf, visible, transparency)
         IMPLICIT NONE
 
@@ -6271 +6944 @@
         LOGICAL, INTENT(out)                   :: visible
         REAL(wp), INTENT(out)                  :: transparency !< along whole path
-        INTEGER(iwp), INTENT(in)               :: win_lad
         INTEGER(iwp)                           :: i, k, d
         INTEGER(iwp)                           :: seldim       !< dimension to be incremented
@@ -6294 +6966 @@
         INTEGER(iwp)                           :: ig           !< 1D index of gridbox in global 2D array
         REAL(wp)                               :: lad_s_target !< recieved lad_s of particular grid box
-        REAL(wp), PARAMETER                    :: grow_factor = 1.5_wp !< factor of expansion of grow arrays
 
 !
@@ -6385 +7056 @@
         IF ( plant_canopy )  THEN
 #if defined( __parallel )
-            IF ( rma_lad_raytrace )  THEN
+            IF ( raytrace_mpi_rma )  THEN
 !--             send requests for lad_s to appropriate processor
-                CALL cpu_log( log_point_s(77), 'rad_init_rma', 'start' )
+                CALL cpu_log( log_point_s(77), 'rad_rma_lad', 'start' )
                 DO i = 1, ncsb
                     CALL MPI_Get(lad_s_ray(i), 1, MPI_REAL, lad_ip(i), lad_disp(i), &
                                  1, MPI_REAL, win_lad, ierr)
                     IF ( ierr /= 0 )  THEN
-                        WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Get'
-                        CALL message( 'raytrace', 'PA0519', 1, 2, 0, 6, 0 )
+                        WRITE(9,*) 'Error MPI_Get1:', ierr, lad_s_ray(i), &
+                                   lad_ip(i), lad_disp(i), win_lad
+                        FLUSH(9)
                     ENDIF
                 ENDDO
@@ -6400 +7072 @@
                 CALL MPI_Win_flush_local_all(win_lad, ierr)
                 IF ( ierr /= 0 )  THEN
-                    WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Win_flush_local_all'
-                    CALL message( 'raytrace', 'PA0519', 1, 2, 0, 6, 0 )
+                    WRITE(9,*) 'Error MPI_Win_flush_local_all1:', ierr, win_lad
+                    FLUSH(9)
                 ENDIF
-                CALL cpu_log( log_point_s(77), 'rad_init_rma', 'stop' )
+                CALL cpu_log( log_point_s(77), 'rad_rma_lad', 'stop' )
                 
             ENDIF
@@ -6411 +7083 @@
             DO i = 1, ncsb
 #if defined( __parallel )
-                IF ( rma_lad_raytrace )  THEN
+                IF ( raytrace_mpi_rma )  THEN
                     lad_s_target = lad_s_ray(i)
                 ELSE
@@ -6429 +7101 @@
                     acsf(ncsfl)%itz = boxes(1,i)
                     acsf(ncsfl)%isurfs = isrc
-                    acsf(ncsfl)%rsvf = REAL(cursink*rirrf*atarg, wp) !-- we postpone multiplication by transparency
-                    acsf(ncsfl)%rtransp = REAL(transparency, wp)
+                    acsf(ncsfl)%rsvf = cursink*transparency*rirrf*atarg
                 ENDIF  !< create_csf
 
@@ -6452 +7123 @@
 !> vertical_delta / horizontal_distance
 !------------------------------------------------------------------------------!
-   SUBROUTINE raytrace_2d(origin, yxdir, zdirs, iorig, aorig, vffrac, &
-                              create_csf, skip_1st_pcb, win_lad, horizon, &
-                              transparency)
+   SUBROUTINE raytrace_2d(origin, yxdir, nrays, zdirs, iorig, aorig, vffrac,  &
+                              calc_svf, create_csf, skip_1st_pcb,             &
+                              lowest_free_ray, transparency, itarget)
       IMPLICIT NONE
 
       REAL(wp), DIMENSION(3), INTENT(IN)     ::  origin        !< z,y,x coordinates of ray origin
       REAL(wp), DIMENSION(2), INTENT(IN)     ::  yxdir         !< y,x *unit* vector of ray direction (in grid units)
-      REAL(wp), DIMENSION(:), INTENT(IN)     ::  zdirs         !< list of z directions to raytrace (z/hdist, in grid)
+      INTEGER(iwp)                           ::  nrays         !< number of rays (z directions) to raytrace
+      REAL(wp), DIMENSION(nrays), INTENT(IN) ::  zdirs         !< list of z directions to raytrace (z/hdist, grid, zenith->nadir)
       INTEGER(iwp), INTENT(in)               ::  iorig         !< index of origin face for csf
       REAL(wp), INTENT(in)                   ::  aorig         !< origin face area for csf
-      REAL(wp), DIMENSION(LBOUND(zdirs, 1):UBOUND(zdirs, 1)), INTENT(in) ::  vffrac !< 
-                                                               !< view factor fractions of each ray for csf
-      LOGICAL, INTENT(in)                    ::  create_csf    !< whether to generate new CSFs during raytracing
+      REAL(wp), DIMENSION(nrays), INTENT(in) ::  vffrac        !< view factor fractions of each ray for csf
+      LOGICAL, INTENT(in)                    ::  calc_svf      !< whether to calculate SFV (identify obstacle surfaces)
+      LOGICAL, INTENT(in)                    ::  create_csf    !< whether to create canopy sink factors
       LOGICAL, INTENT(in)                    ::  skip_1st_pcb  !< whether to skip first plant canopy box during raytracing
-      INTEGER(iwp), INTENT(in)               ::  win_lad       !< leaf area density MPI window
-      REAL(wp), INTENT(OUT)                  ::  horizon       !< highest horizon found after raytracing (z/hdist)
-      REAL(wp), DIMENSION(LBOUND(zdirs, 1):UBOUND(zdirs, 1)), INTENT(OUT) ::  transparency !<
-                                                                !< transparencies of zdirs paths
-      !--INTEGER(iwp), DIMENSION(3, LBOUND(zdirs, 1):UBOUND(zdirs, 1)), INTENT(OUT) :: itarget !<
-                                                                !< (z,y,x) coordinates of target faces for zdirs
+      INTEGER(iwp), INTENT(out)              ::  lowest_free_ray !< index into zdirs
+      REAL(wp), DIMENSION(nrays), INTENT(OUT) ::  transparency !< transparencies of zdirs paths
+      INTEGER(iwp), DIMENSION(nrays), INTENT(OUT) ::  itarget  !< global indices of target faces for zdirs
+
+      INTEGER(iwp), DIMENSION(nrays)         ::  target_procs
+      REAL(wp)                               ::  horizon       !< highest horizon found after raytracing (z/hdist)
       INTEGER(iwp)                           ::  i, k, l, d
       INTEGER(iwp)                           ::  seldim       !< dimension to be incremented
@@ -6506 +7178 @@
       INTEGER(iwp)                           ::  iz
       INTEGER(iwp)                           ::  zsgn
-      REAL(wp), PARAMETER                    ::  grow_factor = 1.5_wp !< factor of expansion of grow arrays
+      INTEGER(iwp)                           ::  lowest_lad   !< lowest column cell for which we need LAD
+      INTEGER(iwp)                           ::  lastdir      !< wall direction before hitting this column
+      INTEGER(iwp), DIMENSION(2)             ::  lastcolumn
 
 #if defined( __parallel )
@@ -6515 +7189 @@
       transparency(:) = 1._wp !-- Pre-set the all rays to transparent before reducing
       horizon = -HUGE(1._wp)
-
-      !--Determine distance to boundary (in 2D xy)
+      lowest_free_ray = nrays
+      IF ( rad_angular_discretization  .AND.  calc_svf )  THEN
+         ALLOCATE(target_surfl(nrays))
+         target_surfl(:) = -1
+         lastdir = -999
+         lastcolumn(:) = -999
+      ENDIF
+
+!--   Determine distance to boundary (in 2D xy)
       IF ( yxdir(1) > 0._wp )  THEN
          bdydim = ny + .5_wp !< north global boundary
@@ -6548 +7229 @@
 !--   Since all face coordinates have values *.5 and we'd like to use
 !--   integers, all these have .5 added
-      DO d = 1, 2
+      DO  d = 1, 2
           IF ( yxdir(d) == 0._wp )  THEN
               dimnext(d) = HUGE(1_iwp)
@@ -6569 +7250 @@
          seldim = minloc(dimnextdist, 1)
          nextdist = dimnextdist(seldim)
-         IF ( nextdist > distance ) nextdist = distance
+         IF ( nextdist > distance )  nextdist = distance
 
          IF ( nextdist > lastdist )  THEN
             ntrack = ntrack + 1
             crmid = (lastdist + nextdist) * .5_wp
-            column = INT(yxorigin(:) + yxdir(:) * crmid, iwp)  !NINT(yxorigin(:) + yxdir(:) * crmid, iwp) 
+            column = NINT(yxorigin(:) + yxdir(:) * crmid, iwp)
 
 !--         calculate index of the grid with global indices (column(1),column(2))
@@ -6586 +7267 @@
                horz_entry = -HUGE(1._wp)
             ELSE
-               horz_entry = (nzterr(ig) - origin(1)) / lastdist
+               horz_entry = (REAL(nzterr(ig), wp) + .5_wp - origin(1)) / lastdist
             ENDIF
-            horz_exit = (nzterr(ig) - origin(1)) / nextdist
+            horz_exit = (REAL(nzterr(ig), wp) + .5_wp - origin(1)) / nextdist
+
+            IF ( rad_angular_discretization  .AND.  calc_svf )  THEN
+!
+!--            Identify vertical obstacles hit by rays in current column
+               DO WHILE ( lowest_free_ray > 0 )
+                  IF ( zdirs(lowest_free_ray) > horz_entry )  EXIT
+!
+!--               This may only happen after 1st column, so lastdir and lastcolumn are valid
+                  CALL request_itarget(lastdir,                                         &
+                        CEILING(-0.5_wp + origin(1) + zdirs(lowest_free_ray)*lastdist), &
+                        lastcolumn(1), lastcolumn(2),                                   &
+                        target_surfl(lowest_free_ray), target_procs(lowest_free_ray))
+                  lowest_free_ray = lowest_free_ray - 1
+               ENDDO
+!
+!--            Identify horizontal obstacles hit by rays in current column
+               DO WHILE ( lowest_free_ray > 0 )
+                  IF ( zdirs(lowest_free_ray) > horz_exit )  EXIT
+                  CALL request_itarget(iup_u, nzterr(ig)+1, column(1), column(2), &
+                                       target_surfl(lowest_free_ray),           &
+                                       target_procs(lowest_free_ray))
+                  lowest_free_ray = lowest_free_ray - 1
+               ENDDO
+            ENDIF
+
             horizon = MAX(horizon, horz_entry, horz_exit)
 
@@ -6598 +7304 @@
 
          IF ( nextdist >= distance )  EXIT
+
+         IF ( rad_angular_discretization  .AND.  calc_svf )  THEN
+!
+!--         Save wall direction of coming building column (= this air column)
+            IF ( seldim == 1 )  THEN
+               IF ( dimdelta(seldim) == 1 )  THEN
+                  lastdir = isouth_u
+               ELSE
+                  lastdir = inorth_u
+               ENDIF
+            ELSE
+               IF ( dimdelta(seldim) == 1 )  THEN
+                  lastdir = iwest_u
+               ELSE
+                  lastdir = ieast_u
+               ENDIF
+            ENDIF
+            lastcolumn = column
+         ENDIF
          lastdist = nextdist
          dimnext(seldim) = dimnext(seldim) + dimdelta(seldim)
@@ -6604 +7329 @@
 
       IF ( plant_canopy )  THEN
-         !--Request LAD WHERE applicable
-         !--
+!--      Request LAD WHERE applicable
+!--      
 #if defined( __parallel )
-         IF ( rma_lad_raytrace )  THEN
+         IF ( raytrace_mpi_rma )  THEN
 !--         send requests for lad_s to appropriate processor
             !CALL cpu_log( log_point_s(77), 'usm_init_rma', 'start' )
-            DO i = 1, ntrack
+            DO  i = 1, ntrack
                px = rt2_track(2,i)/nnx
                py = rt2_track(1,i)/nny
                ip = px*pdims(2)+py
                ig = ip*nnx*nny + (rt2_track(2,i)-px*nnx)*nny + rt2_track(1,i)-py*nny
-               IF ( plantt(ig) <= nzterr(ig) )  CYCLE
-               wdisp = (rt2_track(2,i)-px*nnx)*(nny*nzp) + (rt2_track(1,i)-py*nny)*nzp + nzterr(ig)+1-nzub
-               wcount = plantt(ig)-nzterr(ig)
+
+               IF ( rad_angular_discretization  .AND.  calc_svf )  THEN
+!
+!--               For fixed view resolution, we need plant canopy even for rays
+!--               to opposing surfaces
+                  lowest_lad = nzterr(ig) + 1
+               ELSE
+!
+!--               We only need LAD for rays directed above horizon (to sky)
+                  lowest_lad = CEILING( -0.5_wp + origin(1) +            &
+                                    MIN( horizon * rt2_track_dist(i-1),  & ! entry
+                                         horizon * rt2_track_dist(i)   ) ) ! exit
+               ENDIF
+!
+!--            Skip asking for LAD where all plant canopy is under requested level
+               IF ( plantt(ig) < lowest_lad )  CYCLE
+
+               wdisp = (rt2_track(2,i)-px*nnx)*(nny*nzp) + (rt2_track(1,i)-py*nny)*nzp + lowest_lad-nzub
+               wcount = plantt(ig)-lowest_lad+1
                ! TODO send request ASAP - even during raytracing
-               CALL MPI_Get(rt2_track_lad(nzterr(ig)+1:plantt(ig), i), wcount, MPI_REAL, ip,    &
+               CALL MPI_Get(rt2_track_lad(lowest_lad:plantt(ig), i), wcount, MPI_REAL, ip,    &
                             wdisp, wcount, MPI_REAL, win_lad, ierr)
                IF ( ierr /= 0 )  THEN
-                  WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Get'
-                  CALL message( 'raytrace_2d', 'PA0526', 1, 2, 0, 6, 0 )
+                  WRITE(9,*) 'Error MPI_Get2:', ierr, rt2_track_lad(lowest_lad:plantt(ig), i), &
+                             wcount, ip, wdisp, win_lad
+                  FLUSH(9)
                ENDIF
             ENDDO
@@ -6631 +7373 @@
             CALL MPI_Win_flush_local_all(win_lad, ierr)
             IF ( ierr /= 0 )  THEN
-               WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Win_flush_local_all'
-               CALL message( 'raytrace', 'PA0527', 1, 2, 0, 6, 0 )
+               WRITE(9,*) 'Error MPI_Win_flush_local_all2:', ierr, win_lad
+               FLUSH(9)
             ENDIF
             !CALL cpu_log( log_point_s(77), 'usm_init_rma', 'stop' )
-         ELSE ! rma_lad_raytrace
-            DO i = 1, ntrack
+
+         ELSE ! raytrace_mpi_rma = .F.
+            DO  i = 1, ntrack
                px = rt2_track(2,i)/nnx
                py = rt2_track(1,i)/nny
@@ -6645 +7388 @@
          ENDIF
 #else
-         DO i = 1, ntrack
+         DO  i = 1, ntrack
             rt2_track_lad(nzub:plantt_max, i) = sub_lad(rt2_track(1,i), rt2_track(2,i), nzub:plantt_max)
          ENDDO
 #endif
-
-         !--Skip the PCB around origin if requested
-         !--
-         IF ( skip_1st_pcb )  THEN
+      ENDIF ! plant_canopy
+
+      IF ( rad_angular_discretization  .AND.  calc_svf )  THEN
+#if defined( __parallel )
+!--      wait for all gridsurf requests to complete
+         CALL MPI_Win_flush_local_all(win_gridsurf, ierr)
+         IF ( ierr /= 0 )  THEN
+            WRITE(9,*) 'Error MPI_Win_flush_local_all3:', ierr, win_gridsurf
+            FLUSH(9)
+         ENDIF
+#endif
+!
+!--      recalculate local surf indices into global ones
+         DO i = 1, nrays
+            IF ( target_surfl(i) == -1 )  THEN
+               itarget(i) = -1
+            ELSE
+               itarget(i) = target_surfl(i) + surfstart(target_procs(i))
+            ENDIF
+         ENDDO
+         
+         DEALLOCATE( target_surfl )
+         
+      ELSE
+         itarget(:) = -1
+      ENDIF ! rad_angular_discretization
+
+      IF ( plant_canopy )  THEN
+!--      Skip the PCB around origin if requested (for MRT, the PCB might not be there)
+!--      
+         IF ( skip_1st_pcb  .AND.  NINT(origin(1)) <= plantt_max )  THEN
             rt2_track_lad(NINT(origin(1), iwp), 1) = 0._wp
          ENDIF
 
-         !--Assert that we have space allocated for CSFs
-         !--
-         maxboxes = (ntrack + MAX(origin(1) - nzub, nzut - origin(1))) * SIZE(zdirs, 1)
+!--      Assert that we have space allocated for CSFs
+!--      
+         maxboxes = (ntrack + MAX(CEILING(origin(1)-.5_wp) - nzub,          &
+                                  nzpt - CEILING(origin(1)-.5_wp))) * nrays
          IF ( ncsfl + maxboxes > ncsfla )  THEN
 !--         use this code for growing by fixed exponential increments (equivalent to case where ncsfl always increases by 1)
@@ -6668 +7439 @@
          ENDIF
 
-         !--Calculate transparencies and store new CSFs
-         !--
+!--      Calculate transparencies and store new CSFs
+!--      
          zbottom = REAL(nzub, wp) - .5_wp
          ztop = REAL(plantt_max, wp) + .5_wp
 
-         !--Reverse direction of radiation (face->sky), only when create_csf
-         !--
-         IF ( create_csf )  THEN
-            DO i = 1, ntrack ! for each column
+!--      Reverse direction of radiation (face->sky), only when calc_svf
+!--      
+         IF ( calc_svf )  THEN
+            DO  i = 1, ntrack ! for each column
                dxxyy = ((dy*yxdir(1))**2 + (dx*yxdir(2))**2) * (rt2_track_dist(i)-rt2_track_dist(i-1))**2
                px = rt2_track(2,i)/nnx
@@ -6682 +7453 @@
                ip = px*pdims(2)+py
 
-               DO k = LBOUND(zdirs, 1), UBOUND(zdirs, 1) ! for each ray
-                  IF ( zdirs(k) <= horizon )  THEN
-                     CYCLE
+               DO  k = 1, nrays ! for each ray
+!
+!--               NOTE 6778:
+!--               With traditional svf discretization, CSFs under the horizon
+!--               (i.e. for surface to surface radiation)  are created in
+!--               raytrace(). With rad_angular_discretization, we must create
+!--               CSFs under horizon only for one direction, otherwise we would
+!--               have duplicate amount of energy. Although we could choose
+!--               either of the two directions (they differ only by
+!--               discretization error with no bias), we choose the the backward
+!--               direction, because it tends to cumulate high canopy sink
+!--               factors closer to raytrace origin, i.e. it should potentially
+!--               cause less moiree.
+                  IF ( .NOT. rad_angular_discretization )  THEN
+                     IF ( zdirs(k) <= horizon )  CYCLE
                   ENDIF
 
-                  zorig = REAL(origin(1), wp) + zdirs(k) * rt2_track_dist(i-1)
-                  IF ( zorig <= zbottom .OR. zorig >= ztop )  CYCLE
+                  zorig = origin(1) + zdirs(k) * rt2_track_dist(i-1)
+                  IF ( zorig <= zbottom  .OR.  zorig >= ztop )  CYCLE
 
                   zsgn = INT(SIGN(1._wp, zdirs(k)), iwp)
@@ -6695 +7478 @@
                      nz = 2
                      rt2_dist(nz) = SQRT(dxxyy)
-                     iz = NINT(zorig, iwp)
+                     iz = CEILING(-.5_wp + zorig, iwp)
                   ELSE
-                     zexit = MIN(MAX(REAL(origin(1), wp) + zdirs(k) * rt2_track_dist(i), zbottom), ztop)
+                     zexit = MIN(MAX(origin(1) + zdirs(k) * rt2_track_dist(i), zbottom), ztop)
 
                      zb0 = FLOOR(  zorig * zsgn - .5_wp) + 1  ! because it must be greater than orig
@@ -6708 +7491 @@
                   ENDIF
 
-                  DO l = 2, nz
+                  DO  l = 2, nz
                      IF ( rt2_track_lad(iz, i) > 0._wp )  THEN
                         curtrans = exp(-ext_coef * rt2_track_lad(iz, i) * (rt2_dist(l)-rt2_dist(l-1)))
-                        ncsfl = ncsfl + 1
-                        acsf(ncsfl)%ip = ip
-                        acsf(ncsfl)%itx = rt2_track(2,i)
-                        acsf(ncsfl)%ity = rt2_track(1,i)
-                        acsf(ncsfl)%itz = iz
-                        acsf(ncsfl)%isurfs = iorig
-                        acsf(ncsfl)%rsvf = REAL((1._wp - curtrans)*aorig*vffrac(k), wp) ! we postpone multiplication by transparency
-                        acsf(ncsfl)%rtransp = REAL(transparency(k), wp)
+
+                        IF ( create_csf )  THEN
+                           ncsfl = ncsfl + 1
+                           acsf(ncsfl)%ip = ip
+                           acsf(ncsfl)%itx = rt2_track(2,i)
+                           acsf(ncsfl)%ity = rt2_track(1,i)
+                           acsf(ncsfl)%itz = iz
+                           acsf(ncsfl)%isurfs = iorig
+                           acsf(ncsfl)%rsvf = (1._wp - curtrans)*transparency(k)*aorig*vffrac(k)
+                        ENDIF
 
                         transparency(k) = transparency(k) * curtrans
@@ -6724 +7509 @@
                      iz = iz + zsgn
                   ENDDO ! l = 1, nz - 1
-               ENDDO ! k = LBOUND(zdirs, 1), UBOUND(zdirs, 1)
+               ENDDO ! k = 1, nrays
             ENDDO ! i = 1, ntrack
 
-            transparency(:) = 1._wp !-- Reset all rays to transparent
+            transparency(1:lowest_free_ray) = 1._wp !-- Reset rays above horizon to transparent (see NOTE 6778)
          ENDIF
 
-         !-- Forward direction of radiation (sky->face), always
-         !--
-         DO i = ntrack, 1, -1 ! for each column backwards
+!--      Forward direction of radiation (sky->face), always
+!--      
+         DO  i = ntrack, 1, -1 ! for each column backwards
             dxxyy = ((dy*yxdir(1))**2 + (dx*yxdir(2))**2) * (rt2_track_dist(i)-rt2_track_dist(i-1))**2
             px = rt2_track(2,i)/nnx
@@ -6738 +7523 @@
             ip = px*pdims(2)+py
 
-            DO k = LBOUND(zdirs, 1), UBOUND(zdirs, 1) ! for each ray
-               IF ( zdirs(k) <= horizon )  THEN
-                  transparency(k) = 0._wp
-                  CYCLE
-               ENDIF
-
-               zexit = REAL(origin(1), wp) + zdirs(k) * rt2_track_dist(i-1)
-               IF ( zexit <= zbottom .OR. zexit >= ztop )  CYCLE
+            DO  k = 1, nrays ! for each ray
+!
+!--            See NOTE 6778 above
+               IF ( zdirs(k) <= horizon )  CYCLE
+
+               zexit = origin(1) + zdirs(k) * rt2_track_dist(i-1)
+               IF ( zexit <= zbottom  .OR.  zexit >= ztop )  CYCLE
 
                zsgn = -INT(SIGN(1._wp, zdirs(k)), iwp)
@@ -6754 +7538 @@
                   iz = NINT(zexit, iwp)
                ELSE
-                  zorig = MIN(MAX(REAL(origin(1), wp) + zdirs(k) * rt2_track_dist(i), zbottom), ztop)
+                  zorig = MIN(MAX(origin(1) + zdirs(k) * rt2_track_dist(i), zbottom), ztop)
 
                   zb0 = FLOOR(  zorig * zsgn - .5_wp) + 1  ! because it must be greater than orig
@@ -6765 +7549 @@
                ENDIF
 
-               DO l = 2, nz
+               DO  l = 2, nz
                   IF ( rt2_track_lad(iz, i) > 0._wp )  THEN
                      curtrans = exp(-ext_coef * rt2_track_lad(iz, i) * (rt2_dist(l)-rt2_dist(l-1)))
@@ -6775 +7559 @@
                         acsf(ncsfl)%ity = rt2_track(1,i)
                         acsf(ncsfl)%itz = iz
-                        acsf(ncsfl)%isurfs = -1 ! a special ID indicating sky
-                        acsf(ncsfl)%rsvf = REAL((1._wp - curtrans)*aorig*vffrac(k), wp) ! we postpone multiplication by transparency
-                        acsf(ncsfl)%rtransp = REAL(transparency(k), wp)
-                     ENDIF  !< create_csf
+                        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)
+                     ENDIF  ! create_csf
 
                      transparency(k) = transparency(k) * curtrans
@@ -6784 +7568 @@
                   iz = iz + zsgn
                ENDDO ! l = 1, nz - 1
-            ENDDO ! k = LBOUND(zdirs, 1), UBOUND(zdirs, 1)
+            ENDDO ! k = 1, nrays
          ENDDO ! i = 1, ntrack
-
-      ELSE ! not plant_canopy
-         DO k = UBOUND(zdirs, 1), LBOUND(zdirs, 1), -1 ! TODO make more generic
-            IF ( zdirs(k) > horizon )  EXIT
-            transparency(k) = 0._wp
+      ENDIF ! plant_canopy
+
+      IF ( .NOT. (rad_angular_discretization  .AND.  calc_svf) )  THEN
+!
+!--      Just update lowest_free_ray according to horizon
+         DO WHILE ( lowest_free_ray > 0 )
+            IF ( zdirs(lowest_free_ray) > horizon )  EXIT
+            lowest_free_ray = lowest_free_ray - 1
          ENDDO
       ENDIF
+
+   CONTAINS
+      SUBROUTINE request_itarget(d, z, y, x, isurfl, iproc)
+         INTEGER(iwp), INTENT(in)            ::  d, z, y, x
+         INTEGER(iwp), TARGET, INTENT(out)   ::  isurfl
+         INTEGER(iwp), INTENT(out)           ::  iproc
+         INTEGER(kind=MPI_ADDRESS_KIND)      ::  target_displ  !< index of the grid in the local gridsurf array
+         INTEGER(iwp)                        ::  px, py        !< number of processors in x and y direction
+                                                               !< before the processor in the question
+
+!--      calculate target processor and index in the remote local target gridsurf array
+         px = x/nnx
+         py = y/nny
+         iproc = px*pdims(2)+py
+         target_displ = ((x-px*nnx)*nny + y-py*nny)*nzu*nsurf_type_u + (z-nzub)*nsurf_type_u + d
+
+#if defined( __parallel )
+!--      send MPI_Get request to obtain index target_surfl(i)
+         CALL MPI_Get(isurfl, 1, MPI_INTEGER, iproc, target_displ, &
+                        1, MPI_INTEGER, win_gridsurf, ierr)
+         IF ( ierr /= 0 )  THEN
+            WRITE(9,*) 'Error MPI_Get3:', ierr, isurfl, iproc, target_displ, win_gridsurf
+            FLUSH(9)
+         ENDIF
+#else
+!--      set index target_surfl(i)
+         isurfl = gridsurf(d,z,y,x)
+#endif
+      END SUBROUTINE request_itarget
 
    END SUBROUTINE raytrace_2d
@@ -6849 +7665 @@
       ALLOCATE ( dsitrans(nsurfl, ndsidir) )
       ALLOCATE ( dsitransc(npcbl, ndsidir) )
+      IF ( nmrtbl > 0 )  ALLOCATE ( mrtdsit(nmrtbl, ndsidir) )
 
       WRITE ( message_string, * ) 'Precalculated', ndsidir, ' solar positions', &
@@ -6905 +7722 @@
 
 !-- first check: are the two surfaces directed in the same direction
-        IF ( (d==iup_u  .OR.  d==iup_l  .OR.  d==iup_a )                             &
+        IF ( (d==iup_u  .OR.  d==iup_l )                             &
              .AND. (d2==iup_u  .OR. d2==iup_l) ) RETURN
-        IF ( (d==isouth_u  .OR.  d==isouth_l  .OR.  d==isouth_a ) &
+        IF ( (d==isouth_u  .OR.  d==isouth_l ) &
              .AND.  (d2==isouth_u  .OR.  d2==isouth_l) ) RETURN
-        IF ( (d==inorth_u  .OR.  d==inorth_l  .OR.  d==inorth_a ) &
+        IF ( (d==inorth_u  .OR.  d==inorth_l ) &
              .AND.  (d2==inorth_u  .OR.  d2==inorth_l) ) RETURN
-        IF ( (d==iwest_u  .OR.  d==iwest_l  .OR.  d==iwest_a )     &
+        IF ( (d==iwest_u  .OR.  d==iwest_l )     &
              .AND.  (d2==iwest_u  .OR.  d2==iwest_l ) ) RETURN
-        IF ( (d==ieast_u  .OR.  d==ieast_l  .OR.  d==ieast_a )     &
+        IF ( (d==ieast_u  .OR.  d==ieast_l )     &
              .AND.  (d2==ieast_u  .OR.  d2==ieast_l ) ) RETURN
 
 !-- second check: are surfaces facing away from each other
         SELECT CASE (d)
-            CASE (iup_u, iup_l, iup_a)              !< upward facing surfaces
+            CASE (iup_u, iup_l)                     !< upward facing surfaces
                 IF ( z2 < z ) RETURN
-            CASE (idown_a)                          !< downward facing surfaces
-                IF ( z2 > z ) RETURN
-            CASE (isouth_u, isouth_l, isouth_a)     !< southward facing surfaces
+            CASE (isouth_u, isouth_l)               !< southward facing surfaces
                 IF ( y2 > y ) RETURN
-            CASE (inorth_u, inorth_l, inorth_a)     !< northward facing surfaces
+            CASE (inorth_u, inorth_l)               !< northward facing surfaces
                 IF ( y2 < y ) RETURN
-            CASE (iwest_u, iwest_l, iwest_a)        !< westward facing surfaces
+            CASE (iwest_u, iwest_l)                 !< westward facing surfaces
                 IF ( x2 > x ) RETURN
-            CASE (ieast_u, ieast_l, ieast_a)        !< eastward facing surfaces
+            CASE (ieast_u, ieast_l)                 !< eastward facing surfaces
                 IF ( x2 < x ) RETURN
         END SELECT
@@ -7143 +7958 @@
 !>    array for csf
 !------------------------------------------------------------------------------!
+!-- quicksort.f -*-f90-*-
+!-- Author: t-nissie, adaptation J.Resler
+!-- License: GPLv3
+!-- Gist: https://gist.github.com/t-nissie/479f0f16966925fa29ea
+    RECURSIVE SUBROUTINE quicksort_itarget(itarget, vffrac, ztransp, first, last)
+        IMPLICIT NONE
+        INTEGER(iwp), DIMENSION(:), INTENT(INOUT)   :: itarget
+        REAL(wp), DIMENSION(:), INTENT(INOUT)       :: vffrac, ztransp
+        INTEGER(iwp), INTENT(IN)                    :: first, last
+        INTEGER(iwp)                                :: x, t
+        INTEGER(iwp)                                :: i, j
+        REAL(wp)                                    :: tr
+
+        IF ( first>=last ) RETURN
+        x = itarget((first+last)/2)
+        i = first
+        j = last
+        DO
+            DO WHILE ( itarget(i) < x )
+               i=i+1
+            ENDDO
+            DO WHILE ( x < itarget(j) )
+                j=j-1
+            ENDDO
+            IF ( i >= j ) EXIT
+            t = itarget(i);  itarget(i) = itarget(j);  itarget(j) = t
+            tr = vffrac(i);  vffrac(i) = vffrac(j);  vffrac(j) = tr
+            tr = ztransp(i);  ztransp(i) = ztransp(j);  ztransp(j) = tr
+            i=i+1
+            j=j-1
+        ENDDO
+        IF ( first < i-1 ) CALL quicksort_itarget(itarget, vffrac, ztransp, first, i-1)
+        IF ( j+1 < last )  CALL quicksort_itarget(itarget, vffrac, ztransp, j+1, last)
+    END SUBROUTINE quicksort_itarget
+
     PURE FUNCTION svf_lt(svf1,svf2) result (res)
       TYPE (t_svf), INTENT(in) :: svf1,svf2
@@ -7153 +8003 @@
       ENDIF
     END FUNCTION svf_lt
-    
- 
+
+
 !-- quicksort.f -*-f90-*-
 !-- Author: t-nissie, adaptation J.Resler
@@ -7186 +8036 @@
     END SUBROUTINE quicksort_svf
 
-    
     PURE FUNCTION csf_lt(csf1,csf2) result (res)
       TYPE (t_csf), INTENT(in) :: csf1,csf2
@@ -7241 +8090 @@
         INTEGER(iwp)                            :: iread, iwrite
         TYPE(t_csf), DIMENSION(:), POINTER      :: acsfnew
-!        CHARACTER(100)                          :: msg
+        CHARACTER(100)                          :: msg
 
         IF ( newsize == -1 )  THEN
@@ -7270 +8119 @@
                          .AND.  acsfnew(iwrite)%itz == acsf(iread)%itz &
                          .AND.  acsfnew(iwrite)%isurfs == acsf(iread)%isurfs )  THEN
-!--                 We could simply take either first or second rtransp, both are valid. As a very simple heuristic about which ray
-!--                 probably passes nearer the center of the target box, we choose DIF from the entry with greater CSF, since that
-!--                 might mean that the traced beam passes longer through the canopy box.
-                    IF ( acsfnew(iwrite)%rsvf < acsf(iread)%rsvf )  THEN
-                        acsfnew(iwrite)%rtransp = acsf(iread)%rtransp
-                    ENDIF
+
                     acsfnew(iwrite)%rsvf = acsfnew(iwrite)%rsvf + acsf(iread)%rsvf
 !--                 advance reading index, keep writing index
@@ -7310 +8154 @@
         ncsfla = newsize
 
-!         WRITE(msg,'(A,2I12)') 'Grow acsf2:',ncsfl,ncsfla
-!         CALL radiation_write_debug_log( msg )
+        WRITE(msg,'(A,2I12)') 'Grow acsf2:',ncsfl,ncsfla
+        CALL radiation_write_debug_log( msg )
 
     END SUBROUTINE merge_and_grow_csf
@@ -7443 +8287 @@
                SELECT CASE ( d )
 
-               CASE (iup_u,iup_l,iup_a)                    !- gridbox up_facing face
+               CASE (iup_u,iup_l)                          !- gridbox up_facing face
                   sw_gridbox(1) = surfinsw(l)
                   lw_gridbox(1) = surfinlw(l)
                   swd_gridbox(1) = surfinswdif(l)
 
-               CASE (idown_a)                         !- gridbox down_facing face
-                  sw_gridbox(2) = surfinsw(l)
-                  lw_gridbox(2) = surfinlw(l)
-                  swd_gridbox(2) = surfinswdif(l)
-
-               CASE (inorth_u,inorth_l,inorth_a)  !- gridbox north_facing face
+               CASE (inorth_u,inorth_l)                    !- gridbox north_facing face
                   sw_gridbox(3) = surfinsw(l)
                   lw_gridbox(3) = surfinlw(l)
                   swd_gridbox(3) = surfinswdif(l)
 
-               CASE (isouth_u,isouth_l,isouth_a)  !- gridbox south_facing face
+               CASE (isouth_u,isouth_l)                    !- gridbox south_facing face
                   sw_gridbox(4) = surfinsw(l)
                   lw_gridbox(4) = surfinlw(l)
                   swd_gridbox(4) = surfinswdif(l)
 
-               CASE (ieast_u,ieast_l,ieast_a)      !- gridbox east_facing face
+               CASE (ieast_u,ieast_l)                      !- gridbox east_facing face
                   sw_gridbox(5) = surfinsw(l)
                   lw_gridbox(5) = surfinlw(l)
                   swd_gridbox(5) = surfinswdif(l)
 
-               CASE (iwest_u,iwest_l,iwest_a)      !- gridbox west_facing face
+               CASE (iwest_u,iwest_l)                      !- gridbox west_facing face
                   sw_gridbox(6) = surfinsw(l)
                   lw_gridbox(6) = surfinlw(l)
@@ -7520 +8359 @@
     INTEGER(iwp) ::  j !< 
     INTEGER(iwp) ::  k !< 
-    INTEGER(iwp) ::  m !< index of current surface element 
+    INTEGER(iwp) ::  l, m !< index of current surface element
 
     IF ( mode == 'allocate' )  THEN
@@ -7604 +8443 @@
                 rad_sw_hr_av = 0.0_wp
 
+             CASE ( 'rad_mrt_sw' )
+                IF ( .NOT. ALLOCATED( mrtinsw_av ) )  THEN
+                   ALLOCATE( mrtinsw_av(nmrtbl) )
+                ENDIF
+                mrtinsw_av = 0.0_wp
+
+             CASE ( 'rad_mrt_lw' )
+                IF ( .NOT. ALLOCATED( mrtinlw_av ) )  THEN
+                   ALLOCATE( mrtinlw_av(nmrtbl) )
+                ENDIF
+                mrtinlw_av = 0.0_wp
+
+             CASE ( 'rad_mrt' )
+                IF ( .NOT. ALLOCATED( mrt_av ) )  THEN
+                   ALLOCATE( mrt_av(nmrtbl) )
+                ENDIF
+                mrt_av = 0.0_wp
+
           CASE DEFAULT
              CONTINUE
@@ -7819 +8676 @@
              ENDIF
 
+          CASE ( 'rad_mrt_sw' )
+             IF ( ALLOCATED( mrtinsw_av ) )  THEN
+                mrtinsw_av(:) = mrtinsw_av(:) + mrtinsw(:)
+             ENDIF
+
+          CASE ( 'rad_mrt_lw' )
+             IF ( ALLOCATED( mrtinlw_av ) )  THEN
+                mrtinlw_av(:) = mrtinlw_av(:) + mrtinlw(:)
+             ENDIF
+
+          CASE ( 'rad_mrt' )
+             IF ( ALLOCATED( mrt_av ) )  THEN
+                mrt_av(:) = mrt_av(:) + mrt(:)
+             ENDIF
+
           CASE DEFAULT
              CONTINUE
@@ -7828 +8700 @@
        SELECT CASE ( TRIM( variable ) )
 
-         CASE ( 'rad_net*' )
+          CASE ( 'rad_net*' )
              IF ( ALLOCATED( rad_net_av ) ) THEN
                 DO  i = nxlg, nxrg
@@ -7974 +8846 @@
              ENDIF
 
+          CASE ( 'rad_mrt_sw' )
+             IF ( ALLOCATED( mrtinsw_av ) )  THEN
+                mrtinsw_av(:) = mrtinsw_av(:)  / REAL( average_count_3d, KIND=wp )
+             ENDIF
+
+          CASE ( 'rad_mrt_lw' )
+             IF ( ALLOCATED( mrtinlw_av ) )  THEN
+                mrtinlw_av(:) = mrtinlw_av(:) / REAL( average_count_3d, KIND=wp )
+             ENDIF
+
+          CASE ( 'rad_mrt' )
+             IF ( ALLOCATED( mrt_av ) )  THEN
+                mrt_av(:) = mrt_av(:) / REAL( average_count_3d, KIND=wp )
+             ENDIF
+
        END SELECT
 
@@ -8009 +8896 @@
               'rad_sw_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_hr_xz',               &
               'rad_sw_cs_hr_xz', 'rad_sw_hr_xz', 'rad_lw_cs_hr_yz',            &
-              'rad_lw_hr_yz', 'rad_sw_cs_hr_yz', 'rad_sw_hr_yz' )
+              'rad_lw_hr_yz', 'rad_sw_cs_hr_yz', 'rad_sw_hr_yz',               &
+              'rad_mrt', 'rad_mrt_sw', 'rad_mrt_lw' )
           grid_x = 'x'
           grid_y = 'y'
@@ -8037 +8925 @@
 ! Description:
 ! ------------
-!> Subroutine defining 3D output variables
+!> Subroutine defining 2D output variables
 !------------------------------------------------------------------------------!
  SUBROUTINE radiation_data_output_2d( av, variable, found, grid, mode,         &
@@ -8456 +9344 @@
     CHARACTER (LEN=*) ::  variable !< 
 
-    INTEGER(iwp) ::  av    !< 
-    INTEGER(iwp) ::  i     !< 
-    INTEGER(iwp) ::  j     !< 
-    INTEGER(iwp) ::  k     !< 
-    INTEGER(iwp) ::  nzb_do   !< 
-    INTEGER(iwp) ::  nzt_do   !< 
-
-    LOGICAL      ::  found !< 
-
-    REAL(wp) ::  fill_value = -999.0_wp    !< value for the _FillValue attribute
-
-    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< 
-
+    INTEGER(iwp) ::  av          !<
+    INTEGER(iwp) ::  i, j, k, l  !<
+    INTEGER(iwp) ::  nzb_do      !<
+    INTEGER(iwp) ::  nzt_do      !<
+
+    LOGICAL      ::  found       !<
+
+    REAL(wp)     ::  fill_value = -999.0_wp    !< value for the _FillValue attribute
+
+    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
 
     found = .TRUE.
@@ -8657 +9542 @@
                ENDDO
             ENDDO
+         ENDIF
+
+      CASE ( 'rad_mrt_sw' )
+         local_pf = REAL( fill_value, KIND = wp )
+         IF ( av == 0 )  THEN
+            DO  l = 1, nmrtbl
+               i = mrtbl(ix,l)
+               j = mrtbl(iy,l)
+               k = mrtbl(iz,l)
+               local_pf(i,j,k) = mrtinsw(l)
+            ENDDO
+         ELSE
+            IF ( ALLOCATED( mrtinsw_av ) ) THEN
+               DO  l = 1, nmrtbl
+                  i = mrtbl(ix,l)
+                  j = mrtbl(iy,l)
+                  k = mrtbl(iz,l)
+                  local_pf(i,j,k) = mrtinsw_av(l)
+               ENDDO
+            ENDIF
+         ENDIF
+
+      CASE ( 'rad_mrt_lw' )
+         local_pf = REAL( fill_value, KIND = wp )
+         IF ( av == 0 )  THEN
+            DO  l = 1, nmrtbl
+               i = mrtbl(ix,l)
+               j = mrtbl(iy,l)
+               k = mrtbl(iz,l)
+               local_pf(i,j,k) = mrtinlw(l)
+            ENDDO
+         ELSE
+            IF ( ALLOCATED( mrtinlw_av ) ) THEN
+               DO  l = 1, nmrtbl
+                  i = mrtbl(ix,l)
+                  j = mrtbl(iy,l)
+                  k = mrtbl(iz,l)
+                  local_pf(i,j,k) = mrtinlw_av(l)
+               ENDDO
+            ENDIF
+         ENDIF
+
+      CASE ( 'rad_mrt' )
+         local_pf = REAL( fill_value, KIND = wp )
+         IF ( av == 0 )  THEN
+            DO  l = 1, nmrtbl
+               i = mrtbl(ix,l)
+               j = mrtbl(iy,l)
+               k = mrtbl(iz,l)
+               local_pf(i,j,k) = mrt(l)
+            ENDDO
+         ELSE
+            IF ( ALLOCATED( mrt_av ) ) THEN
+               DO  l = 1, nmrtbl
+                  i = mrtbl(ix,l)
+                  j = mrtbl(iy,l)
+                  k = mrtbl(iz,l)
+                  local_pf(i,j,k) = mrt_av(l)
+               ENDDO
+            ENDIF
          ENDIF
 
@@ -9342 +10287 @@
 !> Subroutine writes debug information
 !------------------------------------------------------------------------------!
-! SUBROUTINE radiation_write_debug_log ( message )
+ SUBROUTINE radiation_write_debug_log ( message )
     !> it writes debug log with time stamp
-!    CHARACTER(*)  :: message
-!    CHARACTER(15) :: dtc
-!    CHARACTER(8)  :: date
-!    CHARACTER(10) :: time
-!    CHARACTER(5)  :: zone
-!    CALL date_and_time(date, time, zone)
-!    dtc = date(7:8)//','//time(1:2)//':'//time(3:4)//':'//time(5:10)
-!    WRITE(9,'(2A)') dtc, TRIM(message)
-!    FLUSH(9)
-! END SUBROUTINE radiation_write_debug_log
+    CHARACTER(*)  :: message
+    CHARACTER(15) :: dtc
+    CHARACTER(8)  :: date
+    CHARACTER(10) :: time
+    CHARACTER(5)  :: zone
+    CALL date_and_time(date, time, zone)
+    dtc = date(7:8)//','//time(1:2)//':'//time(3:4)//':'//time(5:10)
+    WRITE(9,'(2A)') dtc, TRIM(message)
+    FLUSH(9)
+ END SUBROUTINE radiation_write_debug_log
 
  END MODULE radiation_model_mod

palm/trunk/SOURCE/sum_up_3d_data.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! (from branch resler)
+! Add biometeorology,
+! fix chemistry output call,
+! move usm calls
+! 
+! 3294 2018-10-01 02:37:10Z raasch
 ! changes concerning modularization of ocean option
 ! 
@@ -260 +266 @@
     USE radiation_model_mod,                                                   &
         ONLY:  radiation, radiation_3d_data_averaging
+
+    USE biometeorology_mod,                                                    &
+        ONLY:  biometeorology_3d_data_averaging
 
     USE surface_mod,                                                           &
@@ -305 +314 @@
 
        DO  ii = 1, doav_n
-!
-!--       Temporary solution to account for data output within the new urban 
-!--       surface model (urban_surface_mod.f90), see also SELECT CASE ( trimvar )
+
           trimvar = TRIM( doav(ii) )
-          IF ( urban_surface  .AND.  trimvar(1:4) == 'usm_' )  THEN
-             trimvar = 'usm_output'
-          ENDIF
-       
+
           SELECT CASE ( trimvar )
 
@@ -495 +499 @@
                 z0q_av = 0.0_wp
 
-             CASE ( 'usm_output' )
-!             
-!--             Block of urban surface model outputs
-                CALL usm_average_3d_data( 'allocate', doav(ii) )
-              
 
              CASE DEFAULT
@@ -505 +504 @@
 !
 !--             Allocating and initializing data arrays for other modules
+
+                IF ( air_chemistry  .AND. &
+                     (trimvar(1:3) == 'kc_' .OR. trimvar(1:3) == 'em_') )  THEN
+                   CALL chem_3d_data_averaging( 'allocate', doav(ii) )
+                ENDIF
+
                 IF ( bulk_cloud_model )  THEN
                    CALL bcm_3d_data_averaging( 'allocate', doav(ii) )
                 ENDIF
 
-                IF ( air_chemistry  .AND.  trimvar(1:3) == 'kc_')  THEN
-                   CALL chem_3d_data_averaging( 'allocate', doav(ii) )
+                IF ( radiation  .AND.  trimvar(1:4) == 'bio_' )  THEN
+                   CALL biometeorology_3d_data_averaging( 'allocate', doav(ii) )
                 ENDIF
 
@@ -529 +534 @@
                 ENDIF
 
+                CALL tcm_3d_data_averaging( 'allocate', doav(ii) )
+
+                IF ( urban_surface  .AND.  trimvar(1:4) == 'usm_' )  THEN
+                   CALL usm_average_3d_data( 'allocate', doav(ii) )
+                ENDIF
+
                 IF ( uv_exposure  .AND.  trimvar(1:5) == 'uvem_')  THEN
                    CALL uvem_3d_data_averaging( 'allocate', doav(ii) )
                 ENDIF
 
-                CALL tcm_3d_data_averaging( 'allocate', doav(ii) )
-
 !
 !--             User-defined quantities
@@ -548 +557 @@
 !-- Loop of all variables to be averaged.
     DO  ii = 1, doav_n
-!
-!--       Temporary solution to account for data output within the new urban 
-!--       surface model (urban_surface_mod.f90), see also SELECT CASE ( trimvar )
-          trimvar = TRIM( doav(ii) )
-          IF ( urban_surface  .AND.  trimvar(1:4) == 'usm_' )  THEN
-             trimvar = 'usm_output'
-          ENDIF
+
+       trimvar = TRIM( doav(ii) )
 !
 !--    Store the array chosen on the temporary array.
@@ -1146 +1150 @@
              ENDIF
 
-          CASE ( 'usm_output' )
-!             
-!--          Block of urban surface model outputs.
+          CASE DEFAULT
+!
+!--          Summing up data from other modules
+             IF ( bulk_cloud_model )  THEN
+                CALL bcm_3d_data_averaging( 'sum', doav(ii) )
+             ENDIF
+
+             IF ( air_chemistry  .AND. &
+                  (trimvar(1:3) == 'kc_' .OR. trimvar(1:3) == 'em_') )  THEN
+                CALL chem_3d_data_averaging( 'sum',doav(ii) )
+             ENDIF
+
+             IF ( radiation   .AND.  trimvar(1:4) == 'bio_')  THEN
+                CALL biometeorology_3d_data_averaging( 'sum', doav(ii) )
+             ENDIF
+
+             IF ( gust_module_enabled )  THEN
+                CALL gust_3d_data_averaging( 'sum', doav(ii) )
+             ENDIF
+
+             IF ( land_surface )  THEN
+                CALL lsm_3d_data_averaging( 'sum', doav(ii) )
+             ENDIF
+
+             IF ( ocean_mode )  THEN
+                CALL ocean_3d_data_averaging( 'sum', doav(ii) )
+             ENDIF
+
+             IF ( radiation )  THEN
+                CALL radiation_3d_data_averaging( 'sum', doav(ii) )
+             ENDIF
+
+             CALL tcm_3d_data_averaging( 'sum', doav(ii) )
+
 !--          In case of urban surface variables it should be always checked
 !--          if respective arrays are allocated, at least in case of a restart
 !--          run, as averaged usm arrays are not read from file at the moment.
-             CALL usm_average_3d_data( 'allocate', doav(ii) )
-             CALL usm_average_3d_data( 'sum', doav(ii) )
-
-          CASE DEFAULT
-!
-!--          Summing up data from other modules
-             IF ( bulk_cloud_model )  THEN
-                CALL bcm_3d_data_averaging( 'sum', doav(ii) )
-             ENDIF
-
-             IF ( air_chemistry  .AND.  trimvar(1:3) == 'kc_')  THEN
-                CALL chem_3d_data_averaging( 'sum',doav(ii) )
-             ENDIF
-
-             IF ( gust_module_enabled )  THEN
-                CALL gust_3d_data_averaging( 'sum', doav(ii) )
-             ENDIF
-
-             IF ( land_surface )  THEN
-                CALL lsm_3d_data_averaging( 'sum', doav(ii) )
-             ENDIF
-
-             IF ( ocean_mode )  THEN
-                CALL ocean_3d_data_averaging( 'sum', doav(ii) )
-             ENDIF
-
-             IF ( radiation )  THEN
-                CALL radiation_3d_data_averaging( 'sum', doav(ii) )
-             ENDIF
-
-             CALL tcm_3d_data_averaging( 'sum', doav(ii) )
+             IF ( urban_surface  .AND.  trimvar(1:4) == 'usm_' )  THEN
+                CALL usm_average_3d_data( 'allocate', doav(ii) )
+                CALL usm_average_3d_data( 'sum', doav(ii) )
+             ENDIF
 
              IF ( uv_exposure )  THEN

palm/trunk/SOURCE/time_integration_spinup.f90

@@ -25 +25 @@
 ! -----------------
 ! $Id$
+! (from branch resler)
+! Add pt1 initialization
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -295 +299 @@
                 k   = surf_usm_h%k(m)
                 pt(k,j,i) = pt_spinup
+                !!!!!!!!!!!!!!!!HACK!!!!!!!!!!!!!
+                surf_usm_h%pt1 = pt_spinup
+                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
              ENDDO
 
@@ -303 +310 @@
                    k   = surf_usm_v(l)%k(m)
                    pt(k,j,i) = pt_spinup
+                   !!!!!!!!!!!!!!!!HACK!!!!!!!!!!!!!
+                   surf_usm_v(l)%pt1 = pt_spinup
+                   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
                 ENDDO
              ENDDO

palm/trunk/SOURCE/urban_surface_mod.f90

@@ -28 +28 @@
 ! -----------------
 ! $Id$
+! Add output variables usm_rad_pc_inlw, usm_rad_pc_insw*
+! 
+! 3274 2018-09-24 15:42:55Z knoop
 ! Modularization of all bulk cloud physics code components
 ! 
@@ -349 +352 @@
                spinup_pt_mean, spinup_time, time_do3d, dt_do3d,                &
                average_count_3d, varnamelength, urban_surface,                 &
-               plant_canopy
+               plant_canopy, dz
 
     USE cpulog,                                                                &
@@ -381 +384 @@
                surfinl, surfinlwdif, rad_sw_in_dir, rad_sw_in_diff,            &
                rad_lw_in_diff, surfouts, surfoutl, surfoutsl, surfoutll, surf, &
-               surfl, nsurfl, nsurfs, surfstart, pcbinsw, pcbinlw,             &
-               iup_u, inorth_u, isouth_u, ieast_u, iwest_u, iup_l,             &
+               surfl, nsurfl, pcbinsw, pcbinlw, pcbinswdir,                    &
+               pcbinswdif, iup_u, inorth_u, isouth_u, ieast_u, iwest_u, iup_l, &
                inorth_l, isouth_l, ieast_l, iwest_l, id,                       &
-               iz, iy, ix, idir, jdir, kdir,  nsurf_type, nsurf, idsvf, ndsvf, &
-               iup_a, idown_a, inorth_a, isouth_a, ieast_a, iwest_a,           &
+               iz, iy, ix,  nsurf, idsvf, ndsvf,                               &
                idcsf, ndcsf, kdcsf, pct,                                       &
-               startland, endland, startwall, endwall, skyvf, skyvft
+               startland, endland, startwall, endwall, skyvf, skyvft, nzub,    &
+               nzut, nzpt, npcbl, pcbl
 
     USE statistics,                                                            &
@@ -414 +417 @@
     INTEGER(iwp) ::  pedestrian_category = 2         !< default category for wall surface in pedestrian zone
     INTEGER(iwp) ::  roof_category = 2               !< default category for root surface
+    REAL(wp)     ::  roughness_concrete = 0.001_wp   !< roughness length of average concrete surface
 !
 !-- Indices of input attributes for (above) ground floor level 
@@ -605 +609 @@
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfins_av       !< average of array of residua of sw radiation absorbed in surface after last reflection
     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinl_av       !< average of array of residua of lw radiation absorbed in surface after last reflection
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinlw_av       !< Average of pcbinlw
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinsw_av       !< Average of pcbinsw
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdir_av    !< Average of pcbinswdir
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdif_av    !< Average of pcbinswdif
+    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswref_av    !< Average of pcbinswref
     
 
@@ -1260 +1269 @@
 !--     find the real name of the variable
         ids = -1
+        l = -1
         var = TRIM(variable)
         DO i = 0, nd-1
             k = len(TRIM(var))
             j = len(TRIM(dirname(i)))
-            IF ( var(k-j+1:k) == dirname(i) )  THEN
+            IF ( TRIM(var(k-j+1:k)) == TRIM(dirname(i)) )  THEN
                 ids = i
                 idsint = dirint(ids)
@@ -1271 +1281 @@
             ENDIF
         ENDDO
+        l = idsint - 2  ! horisontal direction index - terible hack !
+        IF ( l < 0 .OR. l > 3 ) THEN
+           l = -1
+        END IF
         IF ( ids == -1 )  THEN
             var = TRIM(variable)
@@ -1311 +1325 @@
                 CASE ( 'usm_rad_net' )
 !--                 array of complete radiation balance
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%rad_net_av) )  THEN
+                    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
+                    ELSE
+                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%rad_net_av) )  THEN
+                            ALLOCATE( surf_usm_v(l)%rad_net_av(1:surf_usm_v(l)%ns) )
+                            surf_usm_v(l)%rad_net_av = 0.0_wp
+                        ENDIF
                     ENDIF
-                    DO  l = 0, 3
-                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%rad_net_av) )  THEN
-                           ALLOCATE( surf_usm_v(l)%rad_net_av(1:surf_usm_v(l)%ns) )
-                           surf_usm_v(l)%rad_net_av = 0.0_wp
-                       ENDIF
-                    ENDDO
                     
                 CASE ( 'usm_rad_insw' )
@@ -1398 +1411 @@
                     ENDIF
                                     
+                CASE ( 'usm_rad_pc_inlw' )
+!--                 array of of lw radiation absorbed in plant canopy
+                    IF ( .NOT.  ALLOCATED(pcbinlw_av) )  THEN
+                        ALLOCATE( pcbinlw_av(1:npcbl) )
+                        pcbinlw_av = 0.0_wp
+                    ENDIF
+                                    
+                CASE ( 'usm_rad_pc_insw' )
+!--                 array of of sw radiation absorbed in plant canopy
+                    IF ( .NOT.  ALLOCATED(pcbinsw_av) )  THEN
+                        ALLOCATE( pcbinsw_av(1:npcbl) )
+                        pcbinsw_av = 0.0_wp
+                    ENDIF
+                                    
+                CASE ( 'usm_rad_pc_inswdir' )
+!--                 array of of direct sw radiation absorbed in plant canopy
+                    IF ( .NOT.  ALLOCATED(pcbinswdir_av) )  THEN
+                        ALLOCATE( pcbinswdir_av(1:npcbl) )
+                        pcbinswdir_av = 0.0_wp
+                    ENDIF
+                                    
+                CASE ( 'usm_rad_pc_inswdif' )
+!--                 array of of diffuse sw radiation absorbed in plant canopy
+                    IF ( .NOT.  ALLOCATED(pcbinswdif_av) )  THEN
+                        ALLOCATE( pcbinswdif_av(1:npcbl) )
+                        pcbinswdif_av = 0.0_wp
+                    ENDIF
+                                    
+                CASE ( 'usm_rad_pc_inswref' )
+!--                 array of of reflected sw radiation absorbed in plant canopy
+                    IF ( .NOT.  ALLOCATED(pcbinswref_av) )  THEN
+                        ALLOCATE( pcbinswref_av(1:npcbl) )
+                        pcbinswref_av = 0.0_wp
+                    ENDIF
+                                   
                 CASE ( 'usm_rad_hf' )
 !--                 array of heat flux from radiation for surfaces after i-th reflection
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%surfhf_av) )  THEN
+                    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
-                    ENDIF
-                    DO  l = 0, 3
+                    ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%surfhf_av) )  THEN
                            ALLOCATE( surf_usm_v(l)%surfhf_av(1:surf_usm_v(l)%ns) )
                            surf_usm_v(l)%surfhf_av = 0.0_wp
                        ENDIF
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_wshf' )
 !--                 array of sensible heat flux from surfaces
 !--                 land surfaces
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%wshf_eb_av) )  THEN
+                    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
-                    ENDIF
-                    DO  l = 0, 3
+                    ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wshf_eb_av) )  THEN
                            ALLOCATE( surf_usm_v(l)%wshf_eb_av(1:surf_usm_v(l)%ns) )
                            surf_usm_v(l)%wshf_eb_av = 0.0_wp
                        ENDIF
-                    ENDDO
+                    ENDIF
 !
 !--             Please note, the following output quantities belongs to the 
@@ -1431 +1477 @@
                 CASE ( 'usm_wghf' )
 !--                 array of heat flux from ground (wall, roof, land)
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%wghf_eb_av) )  THEN
+                    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
-                    ENDIF
-                    DO  l = 0, 3
+                    ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wghf_eb_av) )  THEN
                            ALLOCATE( surf_usm_v(l)%wghf_eb_av(1:surf_usm_v(l)%ns) )
                            surf_usm_v(l)%wghf_eb_av = 0.0_wp
                        ENDIF
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_wghf_window' )
 !--                 array of heat flux from window ground (wall, roof, land)
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%wghf_eb_window_av) )  THEN
+                    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
-                    ENDIF
-                    DO  l = 0, 3
+                    ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wghf_eb_window_av) )  THEN
                            ALLOCATE( surf_usm_v(l)%wghf_eb_window_av(1:surf_usm_v(l)%ns) )
                            surf_usm_v(l)%wghf_eb_window_av = 0.0_wp
                        ENDIF
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_wghf_green' )
 !--                 array of heat flux from green ground (wall, roof, land)
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%wghf_eb_green_av) )  THEN
+                    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
-                    ENDIF
-                    DO  l = 0, 3
+                    ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wghf_eb_green_av) )  THEN
                            ALLOCATE( surf_usm_v(l)%wghf_eb_green_av(1:surf_usm_v(l)%ns) )
                            surf_usm_v(l)%wghf_eb_green_av = 0.0_wp
                        ENDIF
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_iwghf' )
 !--                 array of heat flux from indoor ground (wall, roof, land)
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%iwghf_eb_av) )  THEN
+                    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
-                    ENDIF
-                    DO  l = 0, 3
+                    ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%iwghf_eb_av) )  THEN
                            ALLOCATE( surf_usm_v(l)%iwghf_eb_av(1:surf_usm_v(l)%ns) )
                            surf_usm_v(l)%iwghf_eb_av = 0.0_wp
                        ENDIF
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_iwghf_window' )
 !--                 array of heat flux from indoor window ground (wall, roof, land)
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%iwghf_eb_window_av) )  THEN
+                    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
-                    ENDIF
-                    DO  l = 0, 3
+                    ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%iwghf_eb_window_av) )  THEN
                            ALLOCATE( surf_usm_v(l)%iwghf_eb_window_av(1:surf_usm_v(l)%ns) )
                            surf_usm_v(l)%iwghf_eb_window_av = 0.0_wp
                        ENDIF
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_t_surf' )
 !--                 surface temperature for surfaces
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_av) )  THEN
+                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%t_surf_av) )  THEN
                         ALLOCATE( surf_usm_h%t_surf_av(1:surf_usm_h%ns) )
                         surf_usm_h%t_surf_av = 0.0_wp
-                    ENDIF
-                    DO  l = 0, 3
+                    ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_av) )  THEN
                            ALLOCATE( surf_usm_v(l)%t_surf_av(1:surf_usm_v(l)%ns) )
                            surf_usm_v(l)%t_surf_av = 0.0_wp
                        ENDIF
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_t_surf_window' )
 !--                 surface temperature for window surfaces
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_window_av) )  THEN
+                    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
-                    ENDIF
-                    DO  l = 0, 3
+                    ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_window_av) )  THEN
                            ALLOCATE( surf_usm_v(l)%t_surf_window_av(1:surf_usm_v(l)%ns) )
                            surf_usm_v(l)%t_surf_window_av = 0.0_wp
                        ENDIF
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_t_surf_green' )
 !--                 surface temperature for green surfaces
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_green_av) )  THEN
+                    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
-                    ENDIF
-                    DO  l = 0, 3
+                    ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_green_av) )  THEN
                            ALLOCATE( surf_usm_v(l)%t_surf_green_av(1:surf_usm_v(l)%ns) )
                            surf_usm_v(l)%t_surf_green_av = 0.0_wp
                        ENDIF
-                    ENDDO
+                    ENDIF
                 
                 CASE ( 'usm_t_surf_10cm' )
 !--                 near surface temperature for whole surfaces
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_10cm_av) )  THEN
+                    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
-                    ENDIF
-                    DO  l = 0, 3
+                    ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_10cm_av) )  THEN
                            ALLOCATE( surf_usm_v(l)%t_surf_10cm_av(1:surf_usm_v(l)%ns) )
                            surf_usm_v(l)%t_surf_10cm_av = 0.0_wp
                        ENDIF
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_t_wall' )
 !--                 wall temperature for iwl layer of walls and land
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_wall_av) )  THEN
+                    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
-                    ENDIF
-                    DO  l = 0, 3
+                    ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_wall_av) )  THEN
                            ALLOCATE( surf_usm_v(l)%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
                            surf_usm_v(l)%t_wall_av = 0.0_wp
                        ENDIF
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_t_window' )
 !--                 window temperature for iwl layer of walls and land
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_window_av) )  THEN
+                    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
-                    ENDIF
-                    DO  l = 0, 3
+                    ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_window_av) )  THEN
                            ALLOCATE( surf_usm_v(l)%t_window_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
                            surf_usm_v(l)%t_window_av = 0.0_wp
                        ENDIF
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_t_green' )
 !--                 green temperature for iwl layer of walls and land
-                    IF ( .NOT.  ALLOCATED(surf_usm_h%t_green_av) )  THEN
+                    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
-                    ENDIF
-                    DO  l = 0, 3
+                    ELSE
                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_green_av) )  THEN
                            ALLOCATE( surf_usm_v(l)%t_green_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
                            surf_usm_v(l)%t_green_av = 0.0_wp
                        ENDIF
-                    ENDDO
+                    ENDIF
 
                CASE DEFAULT
@@ -1596 +1630 @@
                 CASE ( 'usm_rad_net' )
 !--                 array of complete radiation balance
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%rad_net_av(m) =                              &
-                                          surf_usm_h%rad_net_av(m) +           &
-                                          surf_usm_h%rad_net_l(m)
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%rad_net_av(m) =                              &
+                                             surf_usm_h%rad_net_av(m) +           &
+                                             surf_usm_h%rad_net_l(m)
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%rad_net_av(m) =                        &
@@ -1607 +1642 @@
                                           surf_usm_v(l)%rad_net_l(m)
                        ENDDO
-                    ENDDO
-                    
+                    ENDIF
+
                 CASE ( 'usm_rad_insw' )
 !--                 array of sw radiation falling to surface after i-th reflection
@@ -1700 +1735 @@
                     ENDDO
                     
+                CASE ( 'usm_rad_pc_inlw' )
+                    pcbinlw_av(:) = pcbinlw_av(:) + pcbinlw(:)
+                    
+                CASE ( 'usm_rad_pc_insw' )
+                    pcbinsw_av(:) = pcbinsw_av(:) + pcbinsw(:)
+                    
+                CASE ( 'usm_rad_pc_inswdir' )
+                    pcbinswdir_av(:) = pcbinswdir_av(:) + pcbinswdir(:)
+                    
+                CASE ( 'usm_rad_pc_inswdif' )
+                    pcbinswdif_av(:) = pcbinswdif_av(:) + pcbinswdif(:)
+                    
+                CASE ( 'usm_rad_pc_inswref' )
+                    pcbinswref_av(:) = pcbinswref_av(:) + pcbinsw(:)     &
+                                                        - pcbinswdir(:)  &
+                                                        - pcbinswdif(:)
+                    
                 CASE ( 'usm_rad_hf' )
 !--                 array of heat flux from radiation for surfaces after i-th reflection
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%surfhf_av(m) =                               &
-                                          surf_usm_h%surfhf_av(m) +            &
-                                          surf_usm_h%surfhf(m)
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%surfhf_av(m) =                               &
+                                             surf_usm_h%surfhf_av(m) +            &
+                                             surf_usm_h%surfhf(m)
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%surfhf_av(m) =                         &
@@ -1713 +1766 @@
                                           surf_usm_v(l)%surfhf(m)
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_wshf' )
 !--                 array of sensible heat flux from surfaces (land, roof, wall)
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%wshf_eb_av(m) =                              &
-                                          surf_usm_h%wshf_eb_av(m) +           &
-                                          surf_usm_h%wshf_eb(m)
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%wshf_eb_av(m) =                              &
+                                             surf_usm_h%wshf_eb_av(m) +           &
+                                             surf_usm_h%wshf_eb(m)
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%wshf_eb_av(m) =                        &
@@ -1728 +1782 @@
                                           surf_usm_v(l)%wshf_eb(m)
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_wghf' )
 !--                 array of heat flux from ground (wall, roof, land)
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%wghf_eb_av(m) =                              &
-                                          surf_usm_h%wghf_eb_av(m) +           &
-                                          surf_usm_h%wghf_eb(m)
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%wghf_eb_av(m) =                              &
+                                             surf_usm_h%wghf_eb_av(m) +           &
+                                             surf_usm_h%wghf_eb(m)
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%wghf_eb_av(m) =                        &
@@ -1743 +1798 @@
                                           surf_usm_v(l)%wghf_eb(m)
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_wghf_window' )
 !--                 array of heat flux from window ground (wall, roof, land)
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%wghf_eb_window_av(m) =                              &
-                                          surf_usm_h%wghf_eb_window_av(m) +           &
-                                          surf_usm_h%wghf_eb_window(m)
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%wghf_eb_window_av(m) =                              &
+                                             surf_usm_h%wghf_eb_window_av(m) +           &
+                                             surf_usm_h%wghf_eb_window(m)
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%wghf_eb_window_av(m) =                        &
@@ -1758 +1814 @@
                                           surf_usm_v(l)%wghf_eb_window(m)
                        ENDDO
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_wghf_green' )
 !--                 array of heat flux from green ground (wall, roof, land)
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%wghf_eb_green_av(m) =                              &
-                                          surf_usm_h%wghf_eb_green_av(m) +           &
-                                          surf_usm_h%wghf_eb_green(m)
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%wghf_eb_green_av(m) =                              &
+                                             surf_usm_h%wghf_eb_green_av(m) +           &
+                                             surf_usm_h%wghf_eb_green(m)
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%wghf_eb_green_av(m) =                        &
@@ -1773 +1830 @@
                                           surf_usm_v(l)%wghf_eb_green(m)
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_iwghf' )
 !--                 array of heat flux from indoor ground (wall, roof, land)
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%iwghf_eb_av(m) =                              &
-                                          surf_usm_h%iwghf_eb_av(m) +           &
-                                          surf_usm_h%iwghf_eb(m)
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%iwghf_eb_av(m) =                              &
+                                             surf_usm_h%iwghf_eb_av(m) +           &
+                                             surf_usm_h%iwghf_eb(m)
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%iwghf_eb_av(m) =                        &
@@ -1788 +1846 @@
                                           surf_usm_v(l)%iwghf_eb(m)
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_iwghf_window' )
 !--                 array of heat flux from indoor window ground (wall, roof, land)
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%iwghf_eb_window_av(m) =                              &
-                                          surf_usm_h%iwghf_eb_window_av(m) +           &
-                                          surf_usm_h%iwghf_eb_window(m)
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%iwghf_eb_window_av(m) =                              &
+                                             surf_usm_h%iwghf_eb_window_av(m) +           &
+                                             surf_usm_h%iwghf_eb_window(m)
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%iwghf_eb_window_av(m) =                        &
@@ -1803 +1862 @@
                                           surf_usm_v(l)%iwghf_eb_window(m)
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_t_surf' )
 !--                 surface temperature for surfaces
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%t_surf_av(m) =                               & 
-                                          surf_usm_h%t_surf_av(m) +            &
-                                          t_surf_h(m)
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%t_surf_av(m) =                               &
+                                             surf_usm_h%t_surf_av(m) +            &
+                                             t_surf_h(m)
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_surf_av(m) =                         &
@@ -1818 +1878 @@
                                           t_surf_v(l)%t(m)
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_t_surf_window' )
 !--                 surface temperature for window surfaces
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%t_surf_window_av(m) =                               & 
-                                          surf_usm_h%t_surf_window_av(m) +            &
-                                          t_surf_window_h(m)
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%t_surf_window_av(m) =                               &
+                                             surf_usm_h%t_surf_window_av(m) +            &
+                                             t_surf_window_h(m)
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_surf_window_av(m) =                         &
@@ -1833 +1894 @@
                                           t_surf_window_v(l)%t(m)
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_t_surf_green' )
 !--                 surface temperature for green surfaces
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%t_surf_green_av(m) =                               & 
-                                          surf_usm_h%t_surf_green_av(m) +            &
-                                          t_surf_green_h(m)
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%t_surf_green_av(m) =                               &
+                                             surf_usm_h%t_surf_green_av(m) +            &
+                                             t_surf_green_h(m)
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_surf_green_av(m) =                         &
@@ -1848 +1910 @@
                                           t_surf_green_v(l)%t(m)
                        ENDDO
-                    ENDDO
+                    ENDIF
                 
                 CASE ( 'usm_t_surf_10cm' )
 !--                 near surface temperature for whole surfaces
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%t_surf_10cm_av(m) =                               & 
-                                          surf_usm_h%t_surf_10cm_av(m) +            &
-                                          t_surf_10cm_h(m)
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%t_surf_10cm_av(m) =                               &
+                                             surf_usm_h%t_surf_10cm_av(m) +            &
+                                             t_surf_10cm_h(m)
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_surf_10cm_av(m) =                         &
@@ -1863 +1926 @@
                                           t_surf_10cm_v(l)%t(m)
                        ENDDO
-                    ENDDO
+                    ENDIF
 
                     
                 CASE ( 'usm_t_wall' )
 !--                 wall temperature for  iwl layer of walls and land
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%t_wall_av(iwl,m) =                           &
-                                          surf_usm_h%t_wall_av(iwl,m) +        &
-                                          t_wall_h(iwl,m)
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%t_wall_av(iwl,m) =                           &
+                                             surf_usm_h%t_wall_av(iwl,m) +        &
+                                             t_wall_h(iwl,m)
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_wall_av(iwl,m) =                     &
@@ -1879 +1943 @@
                                           t_wall_v(l)%t(iwl,m)
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_t_window' )
 !--                 window temperature for  iwl layer of walls and land
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%t_window_av(iwl,m) =                           &
-                                          surf_usm_h%t_window_av(iwl,m) +        &
-                                          t_window_h(iwl,m)
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%t_window_av(iwl,m) =                           &
+                                             surf_usm_h%t_window_av(iwl,m) +        &
+                                             t_window_h(iwl,m)
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_window_av(iwl,m) =                     &
@@ -1894 +1959 @@
                                           t_window_v(l)%t(iwl,m)
                        ENDDO
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_t_green' )
 !--                 green temperature for  iwl layer of walls and land
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%t_green_av(iwl,m) =                           &
-                                          surf_usm_h%t_green_av(iwl,m) +        &
-                                          t_green_h(iwl,m)
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%t_green_av(iwl,m) =                           &
+                                             surf_usm_h%t_green_av(iwl,m) +        &
+                                             t_green_h(iwl,m)
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_green_av(iwl,m) =                     &
@@ -1909 +1975 @@
                                           t_green_v(l)%t(iwl,m)
                        ENDDO
-                    ENDDO
+                    ENDIF
 
                 CASE DEFAULT
@@ -1922 +1988 @@
                 CASE ( 'usm_rad_net' )
 !--                 array of complete radiation balance
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%rad_net_av(m) =                              &
-                                          surf_usm_h%rad_net_av(m) /           &
-                                          REAL( average_count_3d, kind=wp )
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%rad_net_av(m) =                              &
+                                             surf_usm_h%rad_net_av(m) /           &
+                                             REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%rad_net_av(m) =                        &
@@ -1933 +2000 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_rad_insw' )
@@ -2023 +2090 @@
                     ENDDO
                     
+                CASE ( 'usm_rad_pc_inlw' )
+                    pcbinlw_av(:) = pcbinlw_av(:) / REAL( average_count_3d, kind=wp )
+                    
+                CASE ( 'usm_rad_pc_insw' )
+                    pcbinsw_av(:) = pcbinsw_av(:) / REAL( average_count_3d, kind=wp )
+                    
+                CASE ( 'usm_rad_pc_inswdir' )
+                    pcbinswdir_av(:) = pcbinswdir_av(:) / REAL( average_count_3d, kind=wp )
+                    
+                CASE ( 'usm_rad_pc_inswdif' )
+                    pcbinswdif_av(:) = pcbinswdif_av(:) / REAL( average_count_3d, kind=wp )
+                    
+                CASE ( 'usm_rad_pc_inswref' )
+                    pcbinswref_av(:) = pcbinswref_av(:) / REAL( average_count_3d, kind=wp )
+                    
                 CASE ( 'usm_rad_hf' )
 !--                 array of heat flux from radiation for surfaces after i-th reflection
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%surfhf_av(m) =                               &
-                                          surf_usm_h%surfhf_av(m) /            &
-                                          REAL( average_count_3d, kind=wp )
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%surfhf_av(m) =                               &
+                                             surf_usm_h%surfhf_av(m) /            &
+                                             REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%surfhf_av(m) =                         &
@@ -2036 +2119 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_wshf' )
 !--                 array of sensible heat flux from surfaces (land, roof, wall)
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%wshf_eb_av(m) =                              &
-                                          surf_usm_h%wshf_eb_av(m) /           &
-                                          REAL( average_count_3d, kind=wp )
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%wshf_eb_av(m) =                              &
+                                             surf_usm_h%wshf_eb_av(m) /           &
+                                             REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%wshf_eb_av(m) =                        &
@@ -2051 +2135 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_wghf' )
 !--                 array of heat flux from ground (wall, roof, land)
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%wghf_eb_av(m) =                              &
-                                          surf_usm_h%wghf_eb_av(m) /           &
-                                          REAL( average_count_3d, kind=wp )
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%wghf_eb_av(m) =                              &
+                                             surf_usm_h%wghf_eb_av(m) /           &
+                                             REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%wghf_eb_av(m) =                        &
@@ -2066 +2151 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_wghf_window' )
 !--                 array of heat flux from window ground (wall, roof, land)
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%wghf_eb_window_av(m) =                              &
-                                          surf_usm_h%wghf_eb_window_av(m) /           &
-                                          REAL( average_count_3d, kind=wp )
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%wghf_eb_window_av(m) =                              &
+                                             surf_usm_h%wghf_eb_window_av(m) /           &
+                                             REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%wghf_eb_window_av(m) =                        &
@@ -2081 +2167 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_wghf_green' )
 !--                 array of heat flux from green ground (wall, roof, land)
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%wghf_eb_green_av(m) =                              &
-                                          surf_usm_h%wghf_eb_green_av(m) /           &
-                                          REAL( average_count_3d, kind=wp )
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%wghf_eb_green_av(m) =                              &
+                                             surf_usm_h%wghf_eb_green_av(m) /           &
+                                             REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%wghf_eb_green_av(m) =                        &
@@ -2096 +2183 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_iwghf' )
 !--                 array of heat flux from indoor ground (wall, roof, land)
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%iwghf_eb_av(m) =                              &
-                                          surf_usm_h%iwghf_eb_av(m) /           &
-                                          REAL( average_count_3d, kind=wp )
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%iwghf_eb_av(m) =                              &
+                                             surf_usm_h%iwghf_eb_av(m) /           &
+                                             REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%iwghf_eb_av(m) =                        &
@@ -2111 +2199 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_iwghf_window' )
 !--                 array of heat flux from indoor window ground (wall, roof, land)
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%iwghf_eb_window_av(m) =                              &
-                                          surf_usm_h%iwghf_eb_window_av(m) /           &
-                                          REAL( average_count_3d, kind=wp )
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%iwghf_eb_window_av(m) =                              &
+                                             surf_usm_h%iwghf_eb_window_av(m) /           &
+                                             REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%iwghf_eb_window_av(m) =                        &
@@ -2126 +2215 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_t_surf' )
 !--                 surface temperature for surfaces
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%t_surf_av(m) =                               & 
-                                          surf_usm_h%t_surf_av(m) /            &
-                                          REAL( average_count_3d, kind=wp )
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%t_surf_av(m) =                               &
+                                             surf_usm_h%t_surf_av(m) /            &
+                                             REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_surf_av(m) =                         &
@@ -2141 +2231 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_t_surf_window' )
 !--                 surface temperature for window surfaces
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%t_surf_window_av(m) =                               & 
-                                          surf_usm_h%t_surf_window_av(m) /            &
-                                          REAL( average_count_3d, kind=wp )
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%t_surf_window_av(m) =                               &
+                                             surf_usm_h%t_surf_window_av(m) /            &
+                                             REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_surf_window_av(m) =                         &
@@ -2156 +2247 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_t_surf_green' )
 !--                 surface temperature for green surfaces
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%t_surf_green_av(m) =                               & 
-                                          surf_usm_h%t_surf_green_av(m) /            &
-                                          REAL( average_count_3d, kind=wp )
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%t_surf_green_av(m) =                               &
+                                             surf_usm_h%t_surf_green_av(m) /            &
+                                             REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_surf_green_av(m) =                         &
@@ -2171 +2263 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_t_surf_10cm' )
 !--                 near surface temperature for whole surfaces
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%t_surf_10cm_av(m) =                               & 
-                                          surf_usm_h%t_surf_10cm_av(m) /            &
-                                          REAL( average_count_3d, kind=wp )
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%t_surf_10cm_av(m) =                               &
+                                             surf_usm_h%t_surf_10cm_av(m) /            &
+                                             REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_surf_10cm_av(m) =                         &
@@ -2186 +2279 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
                     
                 CASE ( 'usm_t_wall' )
 !--                 wall temperature for  iwl layer of walls and land
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%t_wall_av(iwl,m) =                           &
-                                          surf_usm_h%t_wall_av(iwl,m) /        &
-                                          REAL( average_count_3d, kind=wp )
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%t_wall_av(iwl,m) =                           &
+                                             surf_usm_h%t_wall_av(iwl,m) /        &
+                                             REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_wall_av(iwl,m) =                     &
@@ -2201 +2295 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_t_window' )
 !--                 window temperature for  iwl layer of walls and land
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%t_window_av(iwl,m) =                           &
-                                          surf_usm_h%t_window_av(iwl,m) /        &
-                                          REAL( average_count_3d, kind=wp )
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%t_window_av(iwl,m) =                           &
+                                             surf_usm_h%t_window_av(iwl,m) /        &
+                                             REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_window_av(iwl,m) =                     &
@@ -2216 +2311 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
 
                 CASE ( 'usm_t_green' )
 !--                 green temperature for  iwl layer of walls and land
-                    DO  m = 1, surf_usm_h%ns
-                       surf_usm_h%t_green_av(iwl,m) =                           &
-                                          surf_usm_h%t_green_av(iwl,m) /        &
-                                          REAL( average_count_3d, kind=wp )
-                    ENDDO
-                    DO  l = 0, 3
+                    IF ( l == -1 ) THEN
+                       DO  m = 1, surf_usm_h%ns
+                          surf_usm_h%t_green_av(iwl,m) =                           &
+                                             surf_usm_h%t_green_av(iwl,m) /        &
+                                             REAL( average_count_3d, kind=wp )
+                       ENDDO
+                    ELSE
                        DO  m = 1, surf_usm_v(l)%ns
                           surf_usm_v(l)%t_green_av(iwl,m) =                     &
@@ -2231 +2327 @@
                                           REAL( average_count_3d, kind=wp )
                        ENDDO
-                    ENDDO
+                    ENDIF
 
 
@@ -2309 +2405 @@
              var(1:9)  == 'usm_wshf_'  .OR.  var(1:9) == 'usm_wghf_' .OR.                 &
              var(1:16) == 'usm_wghf_window_' .OR. var(1:15) == 'usm_wghf_green_' .OR.     &
-             var(1:10) == 'usm_iwghf_' .OR. var(1:17) == 'usm_iwghf_window_' )  THEN
+             var(1:10) == 'usm_iwghf_' .OR. var(1:17) == 'usm_iwghf_window_' .OR.         &
+             var(1:17) == 'usm_surfwintrans_' )  THEN
             unit = 'W/m2'
-        ELSE IF ( var(1:10) == 'usm_t_surf'   .OR.  var(1:10) == 'usm_t_wall' .OR.         &
+        ELSE IF ( var(1:10) == 'usm_t_surf'   .OR.  var(1:10) == 'usm_t_wall' .OR.        &
                   var(1:12) == 'usm_t_window' .OR. var(1:17) == 'usm_t_surf_window' .OR.  &
                   var(1:16) == 'usm_t_surf_green'  .OR.                                   &
                   var(1:11) == 'usm_t_green' .OR.                                         &
-                  var(1:15) == 'usm_t_surf_10cm')  THEN
+                  var(1:15) == 'usm_t_surf_10cm' )  THEN
             unit = 'K'
+        ELSE IF ( var == 'usm_rad_pc_inlw'  .OR.  var == 'usm_rad_pc_insw'  .OR.          &
+                  var == 'usm_rad_pc_inswdir'  .OR.  var == 'usm_rad_pc_inswdif'  .OR.    &
+                  var == 'usm_rad_pc_inswref' )  THEN
+            unit = 'W'
         ELSE IF ( var(1:9) == 'usm_surfz'  .OR.  var(1:7) == 'usm_svf'  .OR.              & 
                   var(1:7) == 'usm_dif'  .OR.  var(1:11) == 'usm_surfcat'  .OR.           &
@@ -2417 +2518 @@
         INTEGER(iwp), DIMENSION(0:nd-1)                        :: dirstart
         INTEGER(iwp), DIMENSION(0:nd-1)                        :: dirend
-        INTEGER(iwp)                                           :: ids,idsint,idsidx,isurf,isvf,isurfs,isurflt
+        INTEGER(iwp)                                           :: ids,idsint,idsidx,isurf,isvf,isurfs,isurflt,ipcgb
         INTEGER(iwp)                                           :: is,js,ks,i,j,k,iwl,istat, l, m
 
@@ -2431 +2532 @@
             k = len(TRIM(var))
             j = len(TRIM(dirname(i)))
-            IF ( var(k-j+1:k) == dirname(i) )  THEN
+            IF ( TRIM(var(k-j+1:k)) == TRIM(dirname(i)) )  THEN
                 ids = i
                 idsint = dirint(ids)
@@ -2802 +2903 @@
                    ENDIF
                  ENDIF
+              ENDDO
+
+          CASE ( 'usm_rad_pc_inlw' )
+!--           array of lw radiation absorbed by plant canopy
+              DO ipcgb = 1, npcbl
+                  IF ( av == 0 )  THEN
+                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinlw(ipcgb)
+                  ELSE
+                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinlw_av(ipcgb)
+                  ENDIF
+              ENDDO
+
+          CASE ( 'usm_rad_pc_insw' )
+!--           array of sw radiation absorbed by plant canopy
+              DO ipcgb = 1, npcbl
+                  IF ( av == 0 )  THEN
+                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinsw(ipcgb)
+                  ELSE
+                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinsw_av(ipcgb)
+                  ENDIF
+              ENDDO
+
+          CASE ( 'usm_rad_pc_inswdir' )
+!--           array of direct sw radiation absorbed by plant canopy
+              DO ipcgb = 1, npcbl
+                  IF ( av == 0 )  THEN
+                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinswdir(ipcgb)
+                  ELSE
+                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinswdir_av(ipcgb)
+                  ENDIF
+              ENDDO
+
+          CASE ( 'usm_rad_pc_inswdif' )
+!--           array of diffuse sw radiation absorbed by plant canopy
+              DO ipcgb = 1, npcbl
+                  IF ( av == 0 )  THEN
+                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinswdif(ipcgb)
+                  ELSE
+                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinswdif_av(ipcgb)
+                  ENDIF
+              ENDDO
+
+          CASE ( 'usm_rad_pc_inswref' )
+!--           array of reflected sw radiation absorbed by plant canopy
+              DO ipcgb = 1, npcbl
+                  IF ( av == 0 )  THEN
+                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinsw(ipcgb)      &
+                                                                              - pcbinswdir(ipcgb) &
+                                                                              - pcbinswdif(ipcgb)
+                  ELSE
+                      temp_pf(pcbl(iz,ipcgb),pcbl(iy,ipcgb),pcbl(ix,ipcgb)) = pcbinswref_av(ipcgb)
+                  ENDIF
               ENDDO
  
@@ -3356 +3509 @@
           CASE DEFAULT
               found = .FALSE.
-              
+              RETURN
         END SELECT
 
@@ -3399 +3552 @@
              var(1:14) == 'usm_rad_outsw_'  .OR.  var(1:14) == 'usm_rad_outlw_'  .OR.       &
              var(1:14) == 'usm_rad_ressw_'  .OR.  var(1:14) == 'usm_rad_reslw_'  .OR.       &
-             var(1:11) == 'usm_rad_hf_'  .OR.                                               &
+             var(1:11) == 'usm_rad_hf_'  .OR.  var == 'usm_rad_pc_inlw'  .OR.               &
+             var == 'usm_rad_pc_insw'  .OR.  var == 'usm_rad_pc_inswdir'  .OR.              &
+             var == 'usm_rad_pc_inswdif'  .OR.  var == 'usm_rad_pc_inswref'  .OR.           &
              var(1:9) == 'usm_wshf_'  .OR.  var(1:9) == 'usm_wghf_'  .OR.                   &
              var(1:16) == 'usm_wghf_window_'  .OR. var(1:15) == 'usm_wghf_green_' .OR.      &
@@ -4965 +5120 @@
                                  * wintend(nzb_wall:nzt_wall) + tsc(3)            &
                                  * surf_usm_h%tt_window_m(nzb_wall:nzt_wall,m) )   
-           
+
 !
 !--        calculate t_wall tendencies for the next Runge-Kutta step
@@ -5168 +5323 @@
               ENDIF
            ENDDO
+!!!!!!!!!!!!!HACK!!!!!!!!!!!!!!!!!!!
+!           t_window_v_p(l)%t = t_wall_v_p(l)%t
+!           surf_usm_v(l)%tt_window_m  = surf_usm_v(l)%tt_wall_m
+!           t_green_v_p(l)%t = t_wall_v_p(l)%t
+!           surf_usm_v(l)%tt_green_m  = surf_usm_v(l)%tt_wall_m
+!!!!!!!!!!!!!HACK!!!!!!!!!!!!!!!!!!!
         ENDDO
 
@@ -5320 +5481 @@
                            naheatlayers,                                       &
                            pedestrian_category,                                &
+                           roughness_concrete,                                 &
                            read_wall_temp_3d,                                  &
                            roof_category,                                      &
@@ -5337 +5499 @@
                            naheatlayers,                                       &
                            pedestrian_category,                                &
+                           roughness_concrete,                                 &
                            read_wall_temp_3d,                                  &
                            roof_category,                                      &
@@ -6937 +7100 @@
            ENDDO
            IF ( ip == -99999 )  THEN
-!--           wall category not found
-              WRITE (message_string, "(A,I5,A,3I5)") 'wall category ', it,     &
-                                     ' not found  for i,j,k=', iw,jw,kw
-              CALL message( 'usm_read_urban_surface', 'PA0506', 1, 2, 0, 6, 0 )
+!--           land/roof category not found
+              WRITE (9,"(A,I5,A,3I5)") 'land/roof category ', it,     &
+                                       ' not found  for i,j,k=', iw,jw,kw
+              FLUSH(9)
+              IF ( surf_usm_h%isroof_surf(m) ) THEN
+                 category = roof_category
+              ELSE
+                 category = land_category
+              ENDIF
+              DO k = 1, n_surface_types
+                 IF ( surface_type_codes(k) == roof_category ) THEN
+                    ip = k
+                    EXIT
+                 ENDIF
+              ENDDO
+              IF ( ip == -99999 )  THEN
+!--              default land/roof category not found
+                 WRITE (9,"(A,I5,A,3I5)") 'Default land/roof category', category, ' not found!'
+                 FLUSH(9)
+                 ip = 1
+              ENDIF
            ENDIF
 !
@@ -7036 +7216 @@
                  surf_usm_v(l)%albedo(:,m)    = -1.0_wp
                  surf_usm_v(l)%thickness_wall(m) = -1.0_wp
+                 surf_usm_v(l)%thickness_window(m)   = -1.0_wp
+                 surf_usm_v(l)%thickness_green(m)    = -1.0_wp
+                 surf_usm_v(l)%transmissivity(m)  = -1.0_wp
               ELSE IF ( kw <= usm_par(ii,jw,iw) )  THEN
 !--                 pedestrian zone
@@ -7089 +7272 @@
               ELSE
 !
+                 WRITE(9,*) 'Problem reading USM data:'
+                 WRITE(9,*) l,i,j,kw,get_topography_top_index_ji( j, i, 's' )
+                 WRITE(9,*) ii,iw,jw,kw,get_topography_top_index_ji( jw, iw, 's' )
+                 WRITE(9,*) usm_par(ii,jw,iw),usm_par(ii+1,jw,iw)
+                 WRITE(9,*) usm_par(ii+2,jw,iw),usm_par(ii+3,jw,iw)
+                 WRITE(9,*) usm_par(ii+4,jw,iw),usm_par(ii+5,jw,iw)
+                 WRITE(9,*) kw,roof_height_limit,wall_category,roof_category
+                 FLUSH(9)
 !--              supply the default category
                  IF ( kw <= roof_height_limit ) THEN
@@ -7097 +7288 @@
                  surf_usm_v(l)%albedo(:,m)    = -1.0_wp
                  surf_usm_v(l)%thickness_wall(m) = -1.0_wp
+                 surf_usm_v(l)%thickness_window(m) = -1.0_wp
+                 surf_usm_v(l)%thickness_green(m) = -1.0_wp
+                 surf_usm_v(l)%transmissivity(m)  = -1.0_wp
               ENDIF
 !
@@ -7110 +7304 @@
               IF ( ip == -99999 )  THEN
 !--              wall category not found
-                 WRITE (message_string, "(A,I7,A,3I5)") 'wall category ', it,  &
-                                        ' not found  for i,j,k=', iw,jw,kw
-                 WRITE(9,*) message_string
+                 WRITE (9, "(A,I7,A,3I5)") 'wall category ', it,  &
+                                           ' not found  for i,j,k=', iw,jw,kw
+                 FLUSH(9)
+                 category = wall_category 
+                 DO k = 1, n_surface_types
+                    IF ( surface_type_codes(k) == category ) THEN
+                       ip = k
+                       EXIT
+                    ENDIF
+                 ENDDO
+                 IF ( ip == -99999 )  THEN
+!--                 default wall category not found
+                    WRITE (9, "(A,I5,A,3I5)") 'Default wall category', category, ' not found!'
+                    FLUSH(9)
+                    ip = 1
+                 ENDIF
               ENDIF
+
 !
 !--           Albedo
@@ -7192 +7400 @@
         ENDDO
 
+        
+        WRITE(9,*) 'Urban surfaces read'
+        FLUSH(9)
+        
         CALL location_message( '    types and parameters of urban surfaces read', .TRUE. )
    
@@ -7473 +7685 @@
            IF ( humidity )  surf_usm_h%vpt_surface(m) =                        &
                                                    surf_usm_h%pt_surface(m)
-                       
+
 !--        calculate true tendency
            stend = ( t_surf_h_p(m) - t_surf_h(m) - dt_3d * tsc(3) *           &
@@ -7508 +7720 @@
 !
 !--        calculate fluxes
-!--        rad_net_l is never used!           
+!--        rad_net_l is never used!
            surf_usm_h%rad_net_l(m) = surf_usm_h%rad_net_l(m) +                           &
                                      surf_usm_h%frac(ind_veg_wall,m) *                   &
@@ -7522 +7734 @@
            surf_usm_h%wghf_eb(m)   = lambda_surface *                         &
                                       ( t_surf_h_p(m) - t_wall_h(nzb_wall,m) )
-           surf_usm_h%wghf_eb_green(m)  = lambda_surface_green *                        &
+           surf_usm_h%wghf_eb_green(m)  = lambda_surface_green *                         &
                                           ( t_surf_green_h_p(m) - t_green_h(nzb_wall,m) )
-           surf_usm_h%wghf_eb_window(m) = lambda_surface_window *                       &
+           surf_usm_h%wghf_eb_window(m) = lambda_surface_window *                        &
                                            ( t_surf_window_h_p(m) - t_window_h(nzb_wall,m) )
 
@@ -7539 +7751 @@
 !--        diffusion_s, surface_layer_fluxes,...
            surf_usm_h%shf(m) = surf_usm_h%wshf_eb(m) / c_p
-
+      
        ENDDO
 !
@@ -7562 +7774 @@
 #if ! defined( __nopointer )          
 !
-!--          calculate rho * c_p coefficient at surface layer
+!--          calculate rho * c_p coefficient at wall layer
              rho_cp  = c_p * hyp(k) / ( r_d * surf_usm_v(l)%pt1(m) * exner(k) )
 #endif
@@ -7589 +7801 @@
 !--          obtained by simple linear interpolation. ( An alternative would
 !--          be an logarithmic interpolation. )
-!--          A roughness lenght of 0.001 is assumed for concrete (the inverse,
-!--          1000 is used in the nominator for scaling)
-             surf_usm_v(l)%r_a(m) = rho_cp / ( surf_usm_v(l)%z0(m) * 1000.0_wp &
-                        * ( 11.8_wp + 4.2_wp *                                 &
+!--          Parameter roughness_concrete (default value = 0.001) is used
+!--          to calculation of roughness relative to concrete
+             surf_usm_v(l)%r_a(m) = rho_cp / ( surf_usm_v(l)%z0(m) /           &
+                        roughness_concrete * ( 11.8_wp + 4.2_wp *              &
                         SQRT( MAX( ( ( u(k,j,i) + u(k,j,i+1) ) * 0.5_wp )**2 + &
                                    ( ( v(k,j,i) + v(k,j+1,i) ) * 0.5_wp )**2 + &
@@ -7606 +7818 @@
              f_shf_window  = rho_cp / surf_usm_v(l)%r_a(m)
              f_shf_green   = rho_cp / surf_usm_v(l)%r_a(m)
-
-
 
 !--          add LW up so that it can be removed in prognostic equation
@@ -7657 +7867 @@
                            ( surf_usm_v(l)%c_surface_green(m) + coef_green_2 * dt_3d * tsc(2) ) 
 
-
-
 !--           add RK3 term
               t_surf_v_p(l)%t(m) = t_surf_v_p(l)%t(m) + dt_3d * tsc(3) *       &
@@ -7677 +7885 @@
               IF ( humidity )  surf_usm_v(l)%vpt_surface(m) =                  &
                                                     surf_usm_v(l)%pt_surface(m)
-            
+
 !--           calculate true tendency
               stend = ( t_surf_v_p(l)%t(m) - t_surf_v(l)%t(m) - dt_3d * tsc(3) *&
@@ -7766 +7974 @@
            dtime = mod(simulated_time + time_utc_init, 24.0_wp*3600.0_wp)
            dhour = INT(dtime/3600.0_wp)
-           DO m = 1, naheatlayers
-!--           Get indices of respective grid point
-              i = surf_usm_h%i(m)
-              j = surf_usm_h%j(m)
-              k = surf_usm_h%k(m)
-              IF ( k > get_topography_top_index_ji( j, i, 's' )  .AND. &
-                   k <= naheatlayers )  THEN
-!--              increase of pt in box i,j,k in time dt_3d 
-!--              given to anthropogenic heat aheat*acoef (W*m-2)
-!--              linear interpolation of coeficient
-                 acoef = (REAL(dhour+1,wp)-dtime/3600.0_wp)*aheatprof(k,dhour) + &
-                         (dtime/3600.0_wp-REAL(dhour,wp))*aheatprof(k,dhour+1)
-                 IF ( aheat(k,j,i) > 0.0_wp )  THEN
-                    pt(k,j,i) = pt(k,j,i) + aheat(k,j,i)*acoef*dt_3d/(exner(k)*rho_cp*dzu(k))
-                 ENDIF
-              ENDIF
+           DO i = nxl, nxr
+              DO j = nys, nyn
+                 DO k = nzub, min(nzut,naheatlayers)
+                    IF ( k > get_topography_top_index_ji( j, i, 's' ) ) THEN
+!--                    increase of pt in box i,j,k in time dt_3d 
+!--                    given to anthropogenic heat aheat*acoef (W*m-2)
+!--                    linear interpolation of coeficient
+                       acoef = (REAL(dhour+1,wp)-dtime/3600.0_wp)*aheatprof(k,dhour) + &
+                               (dtime/3600.0_wp-REAL(dhour,wp))*aheatprof(k,dhour+1)
+                       IF ( aheat(k,j,i) > 0.0_wp )  THEN
+!--                       calculate rho * c_p coefficient at layer k
+                          rho_cp  = c_p * hyp(k) / ( r_d * pt(k+1,j,i) * exner(k) )
+                          pt(k,j,i) = pt(k,j,i) + aheat(k,j,i)*acoef*dt_3d/(exner(k)*rho_cp*dz(1))
+                       ENDIF
+                    ENDIF
+                 ENDDO
+              ENDDO
            ENDDO