Changeset 3337


Ignore:
Timestamp:
Oct 12, 2018 3:17:09 PM (3 years ago)
Author:
kanani
Message:

reintegrate branch resler to trunk

Location:
palm/trunk/SOURCE
Files:
16 edited
1 copied

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE

  • palm/trunk/SOURCE/Makefile

    r3322 r3337  
    2525# -----------------
    2626# $Id$
     27# (from branch resler)
     28# Add biometeorology
     29#
     30#
     31# 3322 2018-10-09 10:02:39Z kanani
    2732# Formatting and cleanup
    2833#
     
    514519        gust_mod.f90 \
    515520        header.f90 \
     521        biometeorology_mod.f90 \
    516522        inflow_turbulence.f90 \
    517523        init_3d_model.f90 \
     
    712718        modules.o
    713719average_3d_data.o: \
     720        biometeorology_mod.o \
    714721        bulk_cloud_model_mod.o \
    715722        chemistry_model_mod.o \
     
    726733basic_constants_and_equations_mod.o: \
    727734        mod_kinds.o
     735biometeorology_mod.o: \
     736        mod_kinds.o \
     737        radiation_model_mod.o
    728738boundary_conds.o: \
    729739        basic_constants_and_equations_mod.o \
     
    764774check_parameters.o: \
    765775        basic_constants_and_equations_mod.o \
     776        biometeorology_mod.o \
    766777        bulk_cloud_model_mod.o \
    767778        chemistry_model_mod.o \
     
    904915data_output_3d.o: \
    905916        basic_constants_and_equations_mod.o \
     917        biometeorology_mod.o \
    906918        bulk_cloud_model_mod.o \
    907919        chemistry_model_mod.o \
     
    12641276netcdf_interface_mod.o: \
    12651277        basic_constants_and_equations_mod.o \
     1278        biometeorology_mod.o \
    12661279        chemistry_model_mod.o \
    12671280        gust_mod.o \
     
    15131526sum_up_3d_data.o: \
    15141527        basic_constants_and_equations_mod.o \
     1528        biometeorology_mod.o \
    15151529        bulk_cloud_model_mod.o \
    15161530        chemistry_model_mod.o \
  • palm/trunk/SOURCE/average_3d_data.f90

    r3294 r3337  
    2525! -----------------
    2626! $Id$
     27! (from resler branch)
     28! Add biometeorology output,
     29! fix for chemistry output,
     30! move of urban surface output
     31!
     32! 3294 2018-10-01 02:37:10Z raasch
    2733! changes concerning modularization of ocean option
    2834!
     
    174180        ONLY:  radiation, radiation_3d_data_averaging
    175181
     182    USE biometeorology_mod,                                                    &
     183        ONLY:  biometeorology_3d_data_averaging
     184
    176185    USE turbulence_closure_mod,                                                &
    177186        ONLY:  tcm_3d_data_averaging
     
    203212    DO  ii = 1, doav_n
    204213
    205 !
    206 !--    Temporary solution to account for data output within the new urban
    207 !--    surface model (urban_surface_mod.f90), see also SELECT CASE ( trimvar )
    208214       trimvar = TRIM( doav(ii) )
    209        IF ( urban_surface  .AND.  trimvar(1:4) == 'usm_' )  THEN
    210           trimvar = 'usm_output'
    211        ENDIF
    212215
    213216!
     
    533536             ENDIF
    534537
    535           CASE ( 'usm_output' )
    536 !             
    537 !--          Block of urban surface model outputs
    538              CALL usm_average_3d_data( 'average', doav(ii) )
    539 
    540538          CASE DEFAULT
    541539!
    542540!--          Averaging of data from other modules
    543              IF ( air_chemistry )  THEN
     541             IF ( air_chemistry .AND. &
     542                  (trimvar(1:3) == 'kc_' .OR. trimvar(1:3) == 'em_')  )  THEN
    544543                CALL chem_3d_data_averaging( 'average', doav(ii) )
    545544             ENDIF
     
    565564             ENDIF
    566565
     566             IF ( radiation  .AND.  trimvar(1:4) == 'bio_' )  THEN
     567                CALL biometeorology_3d_data_averaging( 'average', doav(ii) )
     568             ENDIF
     569
    567570             CALL tcm_3d_data_averaging( 'average', doav(ii) )
     571
     572             IF ( urban_surface  .AND.  trimvar(1:4) == 'usm_' )  THEN
     573                CALL usm_average_3d_data( 'average', doav(ii) )
     574             ENDIF
     575
    568576!
    569577!--          User-defined quantities
  • palm/trunk/SOURCE/check_parameters.f90

    r3312 r3337  
    2525! -----------------
    2626! $Id$
     27! (from branch resler)
     28! Add biometeorology,
     29! fix for chemistry output,
     30! increase of uv_heights dimension
     31!
     32! 3312 2018-10-06 14:15:46Z knoop
    2733! small code rearrangements
    2834!
     
    756762        ONLY:  radiation, radiation_check_data_output,                         &
    757763               radiation_check_data_output_pr, radiation_check_parameters
     764
     765    USE biometeorology_mod,                                                    &
     766        ONLY:  biometeorology_check_data_output
    758767
    759768    USE spectra_mod,                                                           &
     
    16621671          DO  k = 1, nz+1
    16631672
    1664              IF ( kk < 100 )  THEN
     1673             IF ( kk < 200 )  THEN
    16651674                DO  WHILE ( uv_heights(kk+1) <= zu(k) )
    16661675                   kk = kk + 1
    1667                    IF ( kk == 100 )  EXIT
     1676                   IF ( kk == 200 )  EXIT
    16681677                ENDDO
    16691678             ENDIF
    16701679
    1671              IF ( kk < 100  .AND.  uv_heights(kk+1) /= 9999999.9_wp )  THEN
     1680             IF ( kk < 200  .AND.  uv_heights(kk+1) /= 9999999.9_wp )  THEN
    16721681                u_init(k) = u_profile(kk) + ( zu(k) - uv_heights(kk) ) /       &
    16731682                                       ( uv_heights(kk+1) - uv_heights(kk) ) * &
     
    31563165
    31573166             IF ( unit == 'illegal'  .AND.  air_chemistry                      &
    3158                                      .AND.  var(1:3) == 'kc_' )  THEN
     3167                  .AND.  (var(1:3) == 'kc_' .OR. var(1:3) == 'em_') )  THEN
    31593168                CALL chem_check_data_output( var, unit, i, ilen, k )
    31603169             ENDIF
     
    31623171             IF ( unit == 'illegal' )  THEN
    31633172                CALL lsm_check_data_output ( var, unit, i, ilen, k )
     3173             ENDIF
     3174
     3175             IF ( unit == 'illegal' )  THEN
     3176                CALL biometeorology_check_data_output( var, unit, i, ilen, k )
    31643177             ENDIF
    31653178
  • palm/trunk/SOURCE/chem_emissions_mod.f90

    r3312 r3337  
    2727! -----------------
    2828! $Id$
     29! (from branch resler)
     30! Formatting
     31!
     32! 3312 2018-10-06 14:15:46Z knoop
    2933! Initial revision
    3034!
     
    719723                DO  ispec = 1 , len_index
    720724
    721                    IF ( emiss_factor_main(match_spec_input(ispec)) .LT. 0 .AND. emiss_factor_side(match_spec_input(ispec)) .LT. 0 ) THEN
    722 
    723                       message_string = 'PARAMETERIZED emissions mode selected:'            //          &
     725                   IF ( emiss_factor_main(match_spec_input(ispec)) .LT. 0 .AND.                         &
     726                        emiss_factor_side(match_spec_input(ispec)) .LT. 0 ) THEN
     727
     728                      message_string = 'PARAMETERIZED emissions mode selected:'            //           &
    724729                                       ' EMISSIONS POSSIBLE ONLY ON STREET SURFACES'        //          &
    725730                                       ' but values of scaling factors for street types'    //          &
     
    12461251                IF (TRIM(spc_names(match_spec_model(ispec)))=="NO") THEN
    12471252                !>             Kg/m2                       kg/m2*s                                                   
    1248                    delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%nox_comp(icat,1)*con_factor
     1253                   delta_emis(nys:nyn,nxl:nxr) =                                                             &
     1254                         emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%nox_comp(icat,1)*con_factor
    12491255                   
    1250                    emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
     1256                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                              &
     1257                         emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
    12511258
    12521259                ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="NO2") THEN
    12531260   
    1254                    delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%nox_comp(icat,2)*con_factor
    1255 
    1256                    emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
     1261                   delta_emis(nys:nyn,nxl:nxr) =                                                             &
     1262                         emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%nox_comp(icat,2)*con_factor
     1263
     1264                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                              &
     1265                         emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
    12571266 
    12581267             !> SOX Compositions
     
    12601269                ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="SO2") THEN
    12611270                     
    1262                    delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,1)*con_factor
    1263 
    1264                    emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
     1271                   delta_emis(nys:nyn,nxl:nxr) =                                                             &
     1272                         emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,1)*con_factor
     1273
     1274                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                              &
     1275                         emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
    12651276
    12661277                ELSE IF (TRIM(spc_names(match_spec_model(ispec)))=="SO4") THEN
    12671278   
    1268                    delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,2)*con_factor
    1269 
    1270                    emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
     1279                   delta_emis(nys:nyn,nxl:nxr) =                                                             &
     1280                         emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%sox_comp(icat,2)*con_factor
     1281
     1282                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                              &
     1283                         emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
    12711284 
    12721285
     
    12781291                   DO i_pm_comp= 1,SIZE(emt_att%pm_comp(1,:,1))
    12791292
    1280                       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
     1293                      delta_emis(nys:nyn,nxl:nxr) =                                                          &
     1294                            emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,1)*con_factor
    12811295                                                                                         
    12821296
    1283                       emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
     1297                      emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                           &
     1298                            emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
    12841299 
    12851300                   ENDDO
     
    12911306                   DO i_pm_comp= 1,SIZE(emt_att%pm_comp(1,:,2))
    12921307
    1293                       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
     1308                      delta_emis(nys:nyn,nxl:nxr) =                                                          &
     1309                            emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,2)*con_factor
    12941310                                                                                         
    12951311
    1296                       emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
     1312                      emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                           &
     1313                            emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
    12971314 
    12981315                   ENDDO
     
    13041321                   DO i_pm_comp= 1,SIZE(emt_att%pm_comp(1,:,3)) 
    13051322                       
    1306                       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
     1323                      delta_emis(nys:nyn,nxl:nxr) =                                                          &
     1324                            emis(nys:nyn,nxl:nxr)*time_factor(icat)*emt_att%pm_comp(icat,i_pm_comp,3)*con_factor
    13071325                                                                                                 
    13081326
    1309                       emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
     1327                      emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                           &
     1328                            emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
    13101329
    13111330                   ENDDO
     
    13191338                      IF (TRIM(spc_names(match_spec_model(ispec)))==TRIM(emt_att%voc_name(ivoc))) THEN   
    13201339
    1321                          delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*                               &
     1340                         delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*               &
    13221341                                                       emt_att%voc_comp(icat,match_spec_voc_input(ivoc))*con_factor
    13231342
    1324                          emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
     1343                         emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                         &
     1344                               emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
    13251345
    13261346                      ENDIF                       
     
    13331353                   delta_emis(nys:nyn,nxl:nxr) = emis(nys:nyn,nxl:nxr)*time_factor(icat)*con_factor
    13341354 
    1335                    emis_distribution(1,nys:nyn,nxl:nxr,ispec)=emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
     1355                   emis_distribution(1,nys:nyn,nxl:nxr,ispec) =                                               &
     1356                         emis_distribution(1,nys:nyn,nxl:nxr,ispec)+delta_emis(nys:nyn,nxl:nxr)
    13361357
    13371358                ENDIF  ! IF (spc_names==)
     
    13751396
    13761397                         !> PMs are already in mass units:micrograms:they have to be converted to kilograms
    1377                          IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"  &
     1398                         IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1"         &
     1399                             .OR.  TRIM(spc_names(match_spec_model(ispec)))=="PM25"  &
    13781400                             .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
    13791401
     
    13891411                         !> Other Species: inputs are micromoles: they have to be converted in moles
    13901412                         !                 ppm/s *m *kg/m^3               
    1391                             surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_main(match_spec_input(ispec))*   &
     1413                            surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_main(match_spec_input(ispec))*  &
    13921414                         !                                                    micromoles/(m^2*s)
    1393                                                                           emis_distribution(1,j,i,ispec) *              &
     1415                                                                          emis_distribution(1,j,i,ispec) *             &
    13941416                         !                                                    m**3/Nmole
    1395                                                                           conv_to_ratio(k,j,i)*                         &       
     1417                                                                          conv_to_ratio(k,j,i)*                        &
    13961418                         !                                                    kg/m^3
    13971419                                                                          rho_air(k)   
     
    14031425
    14041426
    1405                    ELSEIF ( street_type_f%var(j,i) >= side_street_id  .AND. street_type_f%var(j,i) < main_street_id )  THEN
     1427                   ELSEIF ( street_type_f%var(j,i) >= side_street_id  .AND.                                            &
     1428                            street_type_f%var(j,i) < main_street_id )  THEN
    14061429
    14071430                   !> Cycle over already matched species
     
    14091432
    14101433                         !> PMs are already in mass units: micrograms
    1411                          IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"  &
     1434                         IF (     TRIM(spc_names(match_spec_model(ispec)))=="PM1"   &
     1435                             .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"  &
    14121436                             .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
    14131437
    14141438                            !              kg/(m^2*s) *kg/m^3                               
    1415                             surf_lsm_h%cssws(match_spec_model(ispec),m)= emiss_factor_side(match_spec_input(ispec)) *   &
     1439                            surf_lsm_h%cssws(match_spec_model(ispec),m)= emiss_factor_side(match_spec_input(ispec)) *  &
    14161440                            !                                                       kg/(m^2*s)
    14171441                                                                                emis_distribution(1,j,i,ispec)*        &
     
    14231447                         !>Other Species: inputs are micromoles
    14241448                         !                 ppm/s *m *kg/m^3               
    1425                             surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_side(match_spec_input(ispec)) *   &
     1449                            surf_lsm_h%cssws(match_spec_model(ispec),m) = emiss_factor_side(match_spec_input(ispec)) * &
    14261450                         !                                                    micromoles/(m^2*s)
    1427                                                                           emis_distribution(1,j,i,ispec) *              &
     1451                                                                          emis_distribution(1,j,i,ispec) *             &
    14281452                         !                                                    m**3/Nmole
    1429                                                                           conv_to_ratio(k,j,i)*                         &       
     1453                                                                          conv_to_ratio(k,j,i)*                        &
    14301454                         !                                                    kg/m^3
    14311455                                                                          rho_air(k)   
     
    14661490
    14671491                   !> PMs
    1468                    IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"  &
    1469                           .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
     1492                   IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1"                                   &
     1493                         .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"                           &
     1494                         .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
    14701495                   
    14711496                      !            kg/(m^2*s) *kg/m^3                         kg/(m^2*s)                 
     
    15221547
    15231548                   !> PMs
    1524                    IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"  &
    1525                           .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
     1549                   IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1"                                            &
     1550                         .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"                                    &
     1551                         .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
    15261552
    15271553                      !         kg/(m^2*s) * kg/m^3                           kg/(m^2*s)           
     
    15731599
    15741600                   !> PMs
    1575                    IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1" .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"  &
    1576                           .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
     1601                   IF (TRIM(spc_names(match_spec_model(ispec)))=="PM1"                                     &
     1602                         .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM25"                             &
     1603                         .OR. TRIM(spc_names(match_spec_model(ispec)))=="PM10") THEN
    15771604                   
    15781605                      !          kg/(m^2*s) *kg/m^3                             kg/(m^2*s)                     
  • palm/trunk/SOURCE/data_output_3d.f90

    r3294 r3337  
    2525! -----------------
    2626! $Id$
     27! (from branch resler)
     28! Add Biometeorology
     29!
     30! 3294 2018-10-01 02:37:10Z raasch
    2731! changes concerning modularization of ocean option
    2832!
     
    269273    USE radiation_model_mod,                                                   &
    270274        ONLY:  nzub, nzut, radiation, radiation_data_output_3d
     275
     276    USE biometeorology_mod,                                                    &
     277        ONLY: biometeorology_data_output_3d
    271278
    272279    USE turbulence_closure_mod,                                                &
     
    744751             IF ( .NOT. found  .AND.  radiation )  THEN
    745752                CALL radiation_data_output_3d( av, do3d(av,if), found,         &
     753                                               local_pf, nzb_do, nzt_do )
     754                resorted = .TRUE.
     755             ENDIF
     756
     757!
     758!--          Biometeorology quantity
     759             IF ( .NOT. found  .AND.  radiation )  THEN
     760                CALL biometeorology_data_output_3d( av, do3d(av,if), found,    &
    746761                                               local_pf, nzb_do, nzt_do )
    747762                resorted = .TRUE.
  • palm/trunk/SOURCE/modules.f90

    r3302 r3337  
    2525! -----------------
    2626! $Id$
     27! (from branch resler)
     28! Increase dimension of uv_heights etc.
     29!
     30! 3302 2018-10-03 02:39:40Z raasch
    2731! +u/v_stokes_zu, u/v_stokes_zw
    2832!
     
    15661570    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
    15671571                 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp, 0.0_wp /)
    1568     REAL(wp) ::  u_profile(100) = 9999999.9_wp                     !< namelist parameter
    1569     REAL(wp) ::  uv_heights(100) = 9999999.9_wp                    !< namelist parameter
    1570     REAL(wp) ::  v_profile(100) = 9999999.9_wp                     !< namelist parameter
     1572    REAL(wp) ::  u_profile(200) = 9999999.9_wp                     !< namelist parameter
     1573    REAL(wp) ::  uv_heights(200) = 9999999.9_wp                    !< namelist parameter
     1574    REAL(wp) ::  v_profile(200) = 9999999.9_wp                     !< namelist parameter
    15711575    REAL(wp) ::  ug_vertical_gradient(10) = 0.0_wp                 !< namelist parameter
    15721576    REAL(wp) ::  ug_vertical_gradient_level(10) = -9999999.9_wp    !< namelist parameter
  • palm/trunk/SOURCE/netcdf_data_input_mod.f90

    r3298 r3337  
    2525! -----------------
    2626! $Id$
     27! (from branch resler)
     28! Formatting
     29!
     30! 3298 2018-10-02 12:21:11Z kanani
    2731! Modified EXPERT mode to PRE-PROCESSED mode (Russo)
    2832! Introduced Chemistry static netcdf file (Russo)
     
    10301034                ALLOCATE(emt_att%hourly_emis_time_factor(1:emt_att%ncat,1:emt_att%nhoursyear))
    10311035                !Read-in Variable
    1032                 CALL get_variable(id_emis,"emission_time_factors",emt_att%hourly_emis_time_factor,1,emt_att%ncat,1,emt_att%nhoursyear)
     1036                CALL get_variable(id_emis,"emission_time_factors",emt_att%hourly_emis_time_factor,1,   &
     1037                                  emt_att%ncat,1,emt_att%nhoursyear)
    10331038
    10341039                !-- MDH
     
    10381043                ALLOCATE(emt_att%mdh_emis_time_factor(1:emt_att%ncat,1:emt_att%nmonthdayhour))
    10391044                !-- Read-in Variable
    1040                 CALL get_variable(id_emis,"emission_time_factors",emt_att%mdh_emis_time_factor,1,emt_att%ncat,1,emt_att%nmonthdayhour)
     1045                CALL get_variable(id_emis,"emission_time_factors",emt_att%mdh_emis_time_factor,1,       &
     1046                                  emt_att%ncat,1,emt_att%nmonthdayhour)
    10411047
    10421048             ELSE
     
    10551061             DO ispec=1,emt_att%nspec
    10561062
    1057                 IF ( .NOT. ALLOCATED( emt(ispec)%default_emission_data ) ) ALLOCATE(emt(ispec)%default_emission_data(1:emt_att%ncat,1:ny+1,1:nx+1))
     1063                IF ( .NOT. ALLOCATED( emt(ispec)%default_emission_data ) )                              &
     1064                    ALLOCATE(emt(ispec)%default_emission_data(1:emt_att%ncat,1:ny+1,1:nx+1))
    10581065
    10591066                ALLOCATE(dum_var_3d(1:emt_att%ncat,nys+1:nyn+1,nxl+1:nxr+1))
     
    10611068                CALL get_variable(id_emis,"emission_values",dum_var_3d,ispec,1,emt_att%ncat,nys,nyn,nxl,nxr)         
    10621069
    1063                 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)
     1070                emt(ispec)%default_emission_data(:,nys+1:nyn+1,nxl+1:nxr+1) =                           &
     1071                    dum_var_3d(1:emt_att%ncat,nys+1:nyn+1,nxl+1:nxr+1)
    10641072
    10651073                DEALLOCATE (dum_var_3d)
     
    11051113
    11061114                !Allocation for the entire domain has to be done only for the first processor between all the subdomains     
    1107                 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))
     1115                IF ( .NOT. ALLOCATED( emt(ispec)%preproc_emission_data ) )                              &
     1116                    ALLOCATE(emt(ispec)%preproc_emission_data(emt_att%dt_emission,1,1:ny+1,1:nx+1))
    11081117
    11091118                !> allocate variable where to pass emission values read from netcdf
     
    11141123
    11151124     
    1116                 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)
     1125                emt(ispec)%preproc_emission_data(:,1,nys+1:nyn+1,nxl+1:nxr+1) =                         &
     1126                      dum_var_4d(1:emt_att%dt_emission,1,nys+1:nyn+1,nxl+1:nxr+1)
    11171127
    11181128                DEALLOCATE ( dum_var_4d )
  • palm/trunk/SOURCE/netcdf_interface_mod.f90

    r3294 r3337  
    2525! -----------------
    2626! $Id$
     27! (from branch resler)
     28! Add biometeorology
     29!
     30! 3294 2018-10-01 02:37:10Z raasch
    2731! changes concerning modularization of ocean option
    2832!
     
    610614        ONLY: radiation, radiation_define_netcdf_grid
    611615
     616    USE biometeorology_mod,                                                    &
     617        ONLY: biometeorology_define_netcdf_grid
     618
    612619    USE spectra_mod,                                                           &
    613620        ONLY:  averaging_interval_sp, comp_spectra_level, data_output_sp, dosp_time_count, spectra_direction
     
    10161023                                                         grid_z )
    10171024                   ENDIF
     1025!
     1026!--                Check for biometeorology quantities
     1027                   IF ( .NOT. found  .AND.  radiation )  THEN
     1028                      CALL biometeorology_define_netcdf_grid( domask(mid,av,i),&
     1029                                                         found, grid_x, grid_y,&
     1030                                                         grid_z )
     1031                   ENDIF
    10181032
    10191033                   CALL tcm_define_netcdf_grid( domask( mid,av,i), found,      &
     
    15451559                                                   grid_y, grid_z )
    15461560                   ENDIF
    1547                    
     1561
    15481562!
    15491563!--                Check for radiation quantities
    15501564                   IF ( .NOT. found  .AND.  radiation )  THEN
    15511565                      CALL radiation_define_netcdf_grid( do3d(av,i), found,    &
     1566                                                         grid_x, grid_y,       &
     1567                                                         grid_z )
     1568                   ENDIF
     1569
     1570!
     1571!--                Check for biometeorology quantities
     1572                   IF ( .NOT. found  .AND.  radiation )  THEN
     1573                      CALL biometeorology_define_netcdf_grid( do3d(av,i), found,&
    15521574                                                         grid_x, grid_y,       &
    15531575                                                         grid_z )
     
    23222344
    23232345!
     2346!--                      Check for biometeorology quantities
     2347                         IF ( .NOT. found  .AND.  radiation )  THEN
     2348                            CALL biometeorology_define_netcdf_grid( do2d(av,i),&
     2349                                                         found, grid_x, grid_y,&
     2350                                                         grid_z )
     2351                         ENDIF
     2352
     2353!
    23242354!--                      Check for gust module quantities
    23252355                         IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
     
    30343064
    30353065!
     3066!--                   Check for biometeorology quantities
     3067                      IF ( .NOT. found  .AND.  radiation )  THEN
     3068                         CALL biometeorology_define_netcdf_grid( do2d(av,i),   &
     3069                                                            found,             &
     3070                                                            grid_x, grid_y,    &
     3071                                                            grid_z )
     3072                      ENDIF
     3073
     3074!
    30363075!--                   Check for gust module quantities
    30373076                      IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
     
    36943733
    36953734!
     3735!--                   Check for biometeorology quantities
     3736                      IF ( .NOT. found  .AND.  radiation )  THEN
     3737                         CALL biometeorology_define_netcdf_grid( do2d(av,i),   &
     3738                                                            found,             &
     3739                                                            grid_x, grid_y,    &
     3740                                                            grid_z )
     3741                      ENDIF
     3742
     3743!
    36963744!--                   Check for gust module quantities
    36973745                      IF ( .NOT. found  .AND.  gust_module_enabled )  THEN
  • palm/trunk/SOURCE/palm.f90

    r3298 r3337  
    2525! -----------------
    2626! $Id$
     27! (from branch resler)
     28! Fix chemistry call
     29!
     30! 3298 2018-10-02 12:21:11Z kanani
    2731! - Minor formatting (kanani)
    2832! - Added Call of date_and_time_init (Russo)
     
    425429!-- --> Needs to be moved!! What is the dependency about?
    426430! IF (  TRIM( initializing_actions ) /= 'read_restart_data' )  THEN
    427     IF ( air_chemistry  .AND.  do_emis )  THEN
     431    IF ( air_chemistry )  THEN
    428432
    429433       IF ( do_emis ) THEN
  • palm/trunk/SOURCE/plant_canopy_model_mod.f90

    r3294 r3337  
    2525! -----------------
    2626! $Id$
     27! Fix reading plant canopy over buildings
     28!
     29! 3294 2018-10-01 02:37:10Z raasch
    2730! ocean renamed ocean_mode
    2831!
     
    11311134       INTEGER(iwp)                        ::  i, j      !< running index
    11321135       INTEGER(iwp)                        ::  nzp       !< number of vertical layers of plant canopy
     1136       INTEGER(iwp)                        ::  nzpltop   !<
     1137       INTEGER(iwp)                        ::  nzpl      !<
     1138       INTEGER(iwp)                        ::  kk        !<
    11331139       
    11341140       REAL(wp), DIMENSION(:), ALLOCATABLE ::  col   !< vertical column of input data
     
    11441150                 FORM='FORMATTED', ERR=515)
    11451151       READ(152, *, ERR=516, END=517)  nzp   !< read first line = number of vertical layers
    1146        
     1152       nzpltop = MIN(nzt+1, nzb+nzp-1)
     1153       nzpl = nzpltop - nzb + 1    !< no. of layers to assign
    11471154       ALLOCATE( col(0:nzp-1) )
    11481155
     
    11641171                   CALL message( 'pcm_read_plant_canopy_3d', 'PA0349', 1, 2, 0, 6, 0 )
    11651172                ENDIF
    1166                 lad_s(0:nzp-1,j,i) = col(0:nzp-1) * lad_type_coef(pctype)
    1167                
     1173                kk = get_topography_top_index_ji( j, i, 's' )
     1174                lad_s(nzb:nzpltop-kk, j, i) = col(kk:nzpl-1)*lad_type_coef(pctype)
    11681175             CASE DEFAULT
    11691176                WRITE(message_string, '(a,i2,a)')   &
  • palm/trunk/SOURCE/prognostic_equations.f90

    r3302 r3337  
    2525! -----------------
    2626! $Id$
     27! (from branch resler)
     28! Fix for chemistry call
     29!
     30! 3302 2018-10-03 02:39:40Z raasch
    2731! Stokes drift + wave breaking term added
    2832!
     
    482486       lsp_usr = 1
    483487!
    484 !--    If required, calculate photolysis frequencies -
    485 !--    UNFINISHED: Why not before the intermediate timestep loop?
    486        IF ( intermediate_timestep_count ==  1 )  THEN
    487           CALL photolysis_control
    488        ENDIF
    489 !
    490488!--    Chemical reactions
    491489       CALL cpu_log( log_point(82), '(chem react + exch_h)', 'start' )
    492490 
    493491       IF ( chem_gasphase_on ) THEN
     492!
     493!--       If required, calculate photolysis frequencies -
     494!--       UNFINISHED: Why not before the intermediate timestep loop?
     495          IF ( intermediate_timestep_count ==  1 )  THEN
     496             CALL photolysis_control
     497          ENDIF
    494498          DO  i = nxl, nxr
    495499             DO  j = nys, nyn
  • palm/trunk/SOURCE/radiation_model_mod.f90

    r3274 r3337  
    1515! PALM. If not, see <http://www.gnu.org/licenses/>.
    1616!
    17 ! Copyright 2015-2018 Czech Technical University in Prague
    1817! Copyright 2015-2018 Institute of Computer Science of the
    1918!                     Czech Academy of Sciences, Prague
     19! Copyright 2015-2018 Czech Technical University in Prague
    2020! Copyright 1997-2018 Leibniz Universitaet Hannover
    2121!------------------------------------------------------------------------------!
     
    2828! -----------------
    2929! $Id$
     30! - New RTM version 2.9: (Pavel Krc, Jaroslav Resler, ICS, Prague)
     31!   added calculation of the MRT inside the RTM module
     32!   MRT fluxes are consequently used in the new biometeorology module
     33!   for calculation of biological indices (MRT, PET)
     34!   Fixes of v. 2.5 and SVN trunk:
     35!    - proper initialization of rad_net_l
     36!    - make arrays nsurfs and surfstart TARGET to prevent some MPI problems
     37!    - initialization of arrays used in MPI one-sided operation as 1-D arrays
     38!      to prevent problems with some MPI/compiler combinations
     39!    - fix indexing of target displacement in subroutine request_itarget to
     40!      consider nzub
     41!    - fix LAD dimmension range in PCB calculation
     42!    - check ierr in all MPI calls
     43!    - use proper per-gridbox sky and diffuse irradiance
     44!    - fix shading for reflected irradiance
     45!    - clear away the residuals of "atmospheric surfaces" implementation
     46!    - fix rounding bug in raytrace_2d introduced in SVN trunk
     47! - New RTM version 2.5: (Pavel Krc, Jaroslav Resler, ICS, Prague)
     48!   can use angular discretization for all SVF
     49!   (i.e. reflected and emitted radiation in addition to direct and diffuse),
     50!   allowing for much better scaling wih high resoltion and/or complex terrain
     51! - Unite array grow factors
     52! - Fix slightly shifted terrain height in raytrace_2d
     53! - Use more efficient MPI_Win_allocate for reverse gridsurf index
     54! - Fix random MPI RMA bugs on Intel compilers
     55! - Fix approx. double plant canopy sink values for reflected radiation
     56! - Fix mostly missing plant canopy sinks for direct radiation
     57! - Fix discretization errors for plant canopy sink in diffuse radiation
     58! - Fix rounding errors in raytrace_2d
     59!
     60! 3274 2018-09-24 15:42:55Z knoop
    3061! Modularization of all bulk cloud physics code components
    3162!
     
    432463!> @todo Output of other rrtm arrays (such as volume mixing ratios)
    433464!> @todo Check for mis-used NINT() calls in raytrace_2d
     465!>       RESULT: Original was correct (carefully verified formula), the change
     466!>               to INT broke raytracing      -- P. Krc
    434467!> @todo Optimize radiation_tendency routines
    435468!>
     
    440473 
    441474    USE arrays_3d,                                                             &
    442         ONLY:  dzw, hyp, nc, pt, q, ql, zu, zw, exner, d_exner
     475        ONLY:  dzw, hyp, nc, pt, p, q, ql, u, v, w, zu, zw, exner, d_exner
    443476
    444477    USE basic_constants_and_equations_mod,                                     &
     
    735768                                             rrtm_swhr,      & !< RRTM output of shortwave radiation heating rate (K/d)
    736769                                             rrtm_swhrc,     & !< RRTM output of incoming shortwave clear sky radiation heating rate (K/d)
    737                                              rrtm_dirdflux,  & !< RRTM output of incoming direct shortwave (W/m²)
    738                                              rrtm_difdflux     !< RRTM output of incoming diffuse shortwave (W/m²)
     770                                             rrtm_dirdflux,  & !< RRTM output of incoming direct shortwave (W/m)
     771                                             rrtm_difdflux     !< RRTM output of incoming diffuse shortwave (W/m)
    739772
    740773    REAL(wp), DIMENSION(1) ::                rrtm_aldif,     & !< surface albedo for longwave diffuse radiation
     
    771804    INTEGER(iwp), PARAMETER                        ::  ndsvf = 2                          !< number of dimensions of real values in SVF
    772805    INTEGER(iwp), PARAMETER                        ::  idsvf = 2                          !< number of dimensions of integer values in SVF
    773     INTEGER(iwp), PARAMETER                        ::  ndcsf = 2                          !< number of dimensions of real values in CSF
     806    INTEGER(iwp), PARAMETER                        ::  ndcsf = 1                          !< number of dimensions of real values in CSF
    774807    INTEGER(iwp), PARAMETER                        ::  idcsf = 2                          !< number of dimensions of integer values in CSF
    775808    INTEGER(iwp), PARAMETER                        ::  kdcsf = 4                          !< number of dimensions of integer values in CSF calculation array
     
    779812    INTEGER(iwp), PARAMETER                        ::  ix = 4                             !< position of i-index in surfl and surf
    780813
    781     INTEGER(iwp), PARAMETER                        ::  nsurf_type = 16                    !< number of surf types incl. phys.(land+urban) & (atm.,sky,boundary) surfaces - 1
     814    INTEGER(iwp), PARAMETER                        ::  nsurf_type = 10                    !< number of surf types incl. phys.(land+urban) & (atm.,sky,boundary) surfaces - 1
    782815
    783816    INTEGER(iwp), PARAMETER                        ::  iup_u    = 0                       !< 0 - index of urban upward surface (ground or roof)
     
    794827    INTEGER(iwp), PARAMETER                        ::  iwest_l  = 10                      !< 10- index of land westward facing wall
    795828
    796     INTEGER(iwp), PARAMETER                        ::  iup_a    = 11                      !< 11- index of atm. cell ubward virtual surface
    797     INTEGER(iwp), PARAMETER                        ::  idown_a  = 12                      !< 12- index of atm. cell downward virtual surface
    798     INTEGER(iwp), PARAMETER                        ::  inorth_a = 13                      !< 13- index of atm. cell northward facing virtual surface
    799     INTEGER(iwp), PARAMETER                        ::  isouth_a = 14                      !< 14- index of atm. cell southward facing virtual surface
    800     INTEGER(iwp), PARAMETER                        ::  ieast_a  = 15                      !< 15- index of atm. cell eastward facing virtual surface
    801     INTEGER(iwp), PARAMETER                        ::  iwest_a  = 16                      !< 16- index of atm. cell westward facing virtual surface
    802 
    803     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
    804     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
    805     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
     829    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  idir = (/0, 0,0, 0,1,-1,0,0, 0,1,-1/)   !< surface normal direction x indices
     830    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  jdir = (/0, 0,1,-1,0, 0,0,1,-1,0, 0/)   !< surface normal direction y indices
     831    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  kdir = (/1,-1,0, 0,0, 0,1,0, 0,0, 0/)   !< surface normal direction z indices
    806832                                                                                          !< parameter but set in the code
    807833
     
    816842
    817843!-- indices and sizes of urban and land surface models
    818     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  surfl            !< coordinates of i-th local surface in local grid - surfl[:,k] = [d, z, y, x]
    819     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  surf             !< coordinates of i-th surface in grid - surf[:,k] = [d, z, y, x]
     844    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfl_l          !< coordinates of i-th local surface in local grid - surfl[:,k] = [d, z, y, x]
     845    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surf_l           !< coordinates of i-th surface in grid - surf[:,k] = [d, z, y, x]
     846    INTEGER(iwp), DIMENSION(:,:), POINTER          ::  surfl            !< coordinates of i-th local surface in local grid - surfl[:,k] = [d, z, y, x]
     847    INTEGER(iwp), DIMENSION(:,:), POINTER          ::  surf             !< coordinates of i-th surface in grid - surf[:,k] = [d, z, y, x]
    820848    INTEGER(iwp)                                   ::  nsurfl           !< number of all surfaces in local processor
    821     INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  nsurfs           !< array of number of all surfaces in individual processors
     849    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  nsurfs           !< array of number of all surfaces in individual processors
    822850    INTEGER(iwp)                                   ::  nsurf            !< global number of surfaces in index array of surfaces (nsurf = proc nsurfs)
    823     INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  surfstart        !< starts of blocks of surfaces for individual processors in array surf
    824                                                                         !< respective block for particular processor is surfstart[iproc]+1 : surfstart[iproc+1]
     851    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfstart        !< starts of blocks of surfaces for individual processors in array surf (indexed from 1)
     852                                                                        !< respective block for particular processor is surfstart[iproc+1]+1 : surfstart[iproc+1]+nsurfs[iproc+1]
    825853
    826854!-- block variables needed for calculation of the plant canopy model inside the urban surface model
     
    835863
    836864!-- configuration parameters (they can be setup in PALM config)
    837     LOGICAL                                        ::  rma_lad_raytrace = .FALSE.         !< use MPI RMA to access LAD for raytracing (instead of global array)
    838     LOGICAL                                        ::  mrt_factors = .FALSE.              !< whether to generate MRT factor files during init
     865    LOGICAL                                        ::  raytrace_mpi_rma = .TRUE.          !< use MPI RMA to access LAD and gridsurf from remote processes during raytracing
     866    LOGICAL                                        ::  read_svf_on_init = .FALSE.         !< flag parameter indicating wheather SVFs will be read from a file at initialization
     867    LOGICAL                                        ::  write_svf_on_init = .FALSE.        !< flag parameter indicating wheather SVFs will be written out to a file
     868    LOGICAL                                        ::  rad_angular_discretization = .TRUE.!< whether to use fixed resolution discretization of view factors for
     869                                                                                          !< reflected radiation (as opposed to all mutually visible pairs)
     870    INTEGER(iwp)                                   ::  mrt_nlevels = 0                    !< number of vertical boxes above surface for which to calculate MRT
     871    LOGICAL                                        ::  mrt_skip_roof = .TRUE.             !< do not calculate MRT above roof surfaces
     872    LOGICAL                                        ::  mrt_include_sw = .TRUE.            !< should MRT calculation include SW radiation as well?
    839873    INTEGER(iwp)                                   ::  nrefsteps = 0                      !< number of reflection steps to perform
    840874    REAL(wp), PARAMETER                            ::  ext_coef = 0.6_wp                  !< extinction coefficient (a.k.a. alpha)
     
    874908        INTEGER(iwp)                               :: isurfs            !<
    875909        REAL(wp)                                   :: rsvf              !<
    876         REAL(wp)                                   :: rtransp           !<
    877910    END TYPE
    878911
    879912!-- arrays storing the values of USM
    880     INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  svfsurf          !< svfsurf[:,isvf] = index of source and target surface for svf[isvf]
     913    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  svfsurf          !< svfsurf[:,isvf] = index of target and source surface for svf[isvf]
    881914    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  svf              !< array of shape view factors+direct irradiation factors for local surfaces
    882915    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfins          !< array of sw radiation falling to local surface after i-th reflection
     
    892925    INTEGER(iwp)                                   ::  ndsidir          !< number of apparent solar directions used
    893926    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  dsidir_rev       !< dsidir_rev[ielev,iazim] = i for dsidir or -1 if not present
     927
     928    INTEGER(iwp)                                   ::  nmrtbl           !< No. of local grid boxes for which MRT is calculated
     929    INTEGER(iwp)                                   ::  nmrtf            !< number of MRT factors for local processor
     930    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  mrtbl            !< coordinates of i-th local MRT box - surfl[:,i] = [z, y, x]
     931    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  mrtfsurf         !< mrtfsurf[:,imrtf] = index of target MRT box and source surface for mrtf[imrtf]
     932    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtf             !< array of MRT factors for each local MRT box
     933    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtft            !< array of MRT factors including transparency for each local MRT box
     934    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtsky           !< array of sky view factor for each local MRT box
     935    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtskyt          !< array of sky view factor including transparency for each local MRT box
     936    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  mrtdsit          !< array of direct solar transparencies for each local MRT box
     937    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtinsw          !< mean SW radiant flux for each MRT box
     938    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtinlw          !< mean LW radiant flux for each MRT box
     939    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrt              !< mean radiant temperature for each MRT box
     940    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtinsw_av       !< time average mean SW radiant flux for each MRT box
     941    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrtinlw_av       !< time average mean LW radiant flux for each MRT box
     942    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  mrt_av           !< time average mean radiant temperature for each MRT box
    894943
    895944    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinsw         !< array of sw radiation falling to local surface including radiation from reflections
     
    920969    TYPE(t_svf), DIMENSION(:), POINTER             ::  asvf             !< pointer to growing svc array
    921970    TYPE(t_csf), DIMENSION(:), POINTER             ::  acsf             !< pointer to growing csf array
     971    TYPE(t_svf), DIMENSION(:), POINTER             ::  amrtf            !< pointer to growing mrtf array
    922972    TYPE(t_svf), DIMENSION(:), ALLOCATABLE, TARGET ::  asvf1, asvf2     !< realizations of svf array
    923973    TYPE(t_csf), DIMENSION(:), ALLOCATABLE, TARGET ::  acsf1, acsf2     !< realizations of csf array
     974    TYPE(t_svf), DIMENSION(:), ALLOCATABLE, TARGET ::  amrtf1, amrtf2   !< realizations of mftf array
    924975    INTEGER(iwp)                                   ::  nsvfla           !< dimmension of array allocated for storage of svf in local processor
    925976    INTEGER(iwp)                                   ::  ncsfla           !< dimmension of array allocated for storage of csf in local processor
    926     INTEGER(iwp)                                   ::  msvf, mcsf       !< mod for swapping the growing array
    927     INTEGER(iwp), PARAMETER                        ::  gasize = 10000   !< initial size of growing arrays
     977    INTEGER(iwp)                                   ::  nmrtfa           !< dimmension of array allocated for storage of mrt
     978    INTEGER(iwp)                                   ::  msvf, mcsf, mmrtf!< mod for swapping the growing array
     979    INTEGER(iwp), PARAMETER                        ::  gasize = 100000  !< initial size of growing arrays
     980    REAL(wp), PARAMETER                            ::  grow_factor = 1.4_wp !< growth factor of growing arrays
    928981    REAL(wp)                                       ::  dist_max_svf = -9999.0 !< maximum distance to calculate the minimum svf to be considered. It is
    929982                                                                        !< used to avoid very small SVFs resulting from too far surfaces with mutual visibility
     
    932985                                                                        !< needed only during calc_svf but must be here because it is
    933986                                                                        !< shared between subroutines calc_svf and raytrace
    934     INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE    ::  gridpcbl         !< index of local pcb[k,j,i]
     987    INTEGER(iwp), DIMENSION(:,:,:,:), POINTER      ::  gridsurf         !< reverse index of local surfl[d,k,j,i] (for case rad_angular_discretization)
     988    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE    ::  gridpcbl         !< reverse index of local pcbl[k,j,i]
     989    INTEGER(iwp), PARAMETER                        ::  nsurf_type_u = 6 !< number of urban surf types (used in gridsurf)
    935990
    936991!-- temporary arrays for calculation of csf in raytracing
     
    942997    INTEGER(kind=MPI_ADDRESS_KIND), &
    943998                  DIMENSION(:), ALLOCATABLE        ::  lad_disp         !< array of displaycements of lad in local array of proc lad_ip
     999    INTEGER(iwp)                                   ::  win_lad          !< MPI RMA window for leaf area density
     1000    INTEGER(iwp)                                   ::  win_gridsurf     !< MPI RMA window for reverse grid surface index
    9441001#endif
    9451002    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  lad_s_ray        !< array of received lad_s for appropriate gridboxes crossed by ray
     1003    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  target_surfl
    9461004    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  rt2_track
    9471005    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  rt2_track_lad
     
    10831141!-- Public variables and constants / NEEDS SORTING
    10841142    PUBLIC albedo, albedo_type, decl_1, decl_2, decl_3, dots_rad, dt_radiation,&
    1085            emissivity, force_radiation_call,                                   &
    1086            lat, lon, rad_net_av, radiation, radiation_scheme, rad_lw_in,       &
     1143           emissivity, force_radiation_call, lat, lon,                         &
     1144           mrt_include_sw, mrt_nlevels, mrtbl, mrtinsw, mrtinlw, nmrtbl,       &
     1145           rad_net_av, radiation, radiation_scheme, rad_lw_in,                 &
    10871146           rad_lw_in_av, rad_lw_out, rad_lw_out_av,                            &
    10881147           rad_lw_cs_hr, rad_lw_cs_hr_av, rad_lw_hr, rad_lw_hr_av, rad_sw_in,  &
     
    10911150           skip_time_do_radiation, time_radiation, unscheduled_radiation_calls,&
    10921151           zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon,       &
    1093            nrefsteps, mrt_factors, dist_max_svf, nsvfl, svf,                   &
     1152           write_svf_on_init, read_svf_on_init,                                &
     1153           nrefsteps, dist_max_svf, nsvfl, svf,                                &
    10941154           svfsurf, surfinsw, surfinlw, surfins, surfinl, surfinswdir,         &
    10951155           surfinswdif, surfoutsw, surfoutlw, surfinlwdif, rad_sw_in_dir,      &
    10961156           rad_sw_in_diff, rad_lw_in_diff, surfouts, surfoutl, surfoutsl,      &
    1097            surfoutll, idir, jdir, kdir, id, iz, iy, ix, nsurfs, surfstart,     &
     1157           surfoutll, idir, jdir, kdir, id, iz, iy, ix,                        &
    10981158           surf, surfl, nsurfl, pcbinswdir, pcbinswdif, pcbinsw, pcbinlw,      &
    1099            iup_u, inorth_u, isouth_u, ieast_u, iwest_u,           &
     1159           iup_u, inorth_u, isouth_u, ieast_u, iwest_u,                        &
    11001160           iup_l, inorth_l, isouth_l, ieast_l, iwest_l,                        &
    1101            nsurf_type, nzub, nzut, nzu, pch, nsurf,                            &
    1102            iup_a, idown_a, inorth_a, isouth_a, ieast_a, iwest_a,               &
     1161           nsurf_type, nzub, nzut, nzpt, nzu, pch, nsurf,                      &
    11031162           idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct,                             &
    11041163           radiation_interactions, startwall, startland, endland, endwall,     &
    1105            skyvf, skyvft, radiation_interactions_on, average_radiation
     1164           skyvf, skyvft, radiation_interactions_on, average_radiation, npcbl, &
     1165           pcbl
    11061166
    11071167#if defined ( __rrtmg )
     
    11631223       SELECT CASE ( TRIM( var ) )
    11641224
    1165          CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr', 'rad_sw_in' )
     1225          CASE ( 'rad_lw_cs_hr', 'rad_lw_hr', 'rad_sw_cs_hr', 'rad_sw_hr', 'rad_sw_in' )
    11661226             IF (  .NOT.  radiation  .OR.  radiation_scheme /= 'rrtmg' )  THEN
    11671227                message_string = '"output of "' // TRIM( var ) // '" requi' // &
     
    12051265             IF ( TRIM( var ) == 'rrtm_asdir*'   ) unit = ''
    12061266
     1267          CASE ( 'rad_mrt', 'rad_mrt_sw', 'rad_mrt_lw'  )
     1268             IF ( .NOT.  radiation ) THEN
     1269                message_string = 'output of "' // TRIM( var ) // '" require'&
     1270                                 // 's radiation = .TRUE.'
     1271                CALL message( 'check_parameters', 'PA0509', 1, 2, 0, 6, 0 )
     1272             ENDIF
     1273             IF ( mrt_nlevels == 0 ) THEN
     1274                message_string = 'output of "' // TRIM( var ) // '" require'&
     1275                                 // 's mrt_nlevels > 0'
     1276                CALL message( 'check_parameters', 'PA0510', 1, 2, 0, 6, 0 )
     1277             ENDIF
     1278             IF ( TRIM( var ) == 'rad_mrt_sw'  .AND.  .NOT. mrt_include_sw ) THEN
     1279                message_string = 'output of "' // TRIM( var ) // '" require'&
     1280                                 // 's rad_mrt_sw = .TRUE.'
     1281                CALL message( 'check_parameters', 'PA0511', 1, 2, 0, 6, 0 )
     1282             ENDIF
     1283             IF ( TRIM( var ) == 'rad_mrt' ) THEN
     1284                unit = 'K'
     1285             ELSE
     1286                unit = 'W m-2'
     1287             ENDIF
     1288
    12071289          CASE DEFAULT
    12081290             unit = 'illegal'
     
    14571539          ENDIF
    14581540       ENDIF
     1541!
     1542!--    Parallel rad_angular_discretization without raytrace_mpi_rma is not implemented
     1543#if defined( __parallel )     
     1544       IF ( rad_angular_discretization  .AND.  .NOT. raytrace_mpi_rma )  THEN
     1545          message_string = 'rad_angular_discretization can only be used ' //  &
     1546                           'together with raytrace_mpi_rma or when ' //  &
     1547                           'no parallelization is applied.'
     1548          CALL message( 'check_parameters', 'PA0486', 1, 2, 0, 6, 0 )
     1549       ENDIF
     1550#endif
    14591551
    14601552       IF ( cloud_droplets  .AND.   radiation_scheme == 'rrtmg'  .AND.         &
     
    23542446!--    In case averaged radiation is used, calculate mean temperature and
    23552447!--    liquid water mixing ratio at the urban-layer top.
    2356        IF ( average_radiation ) THEN   
     2448       IF ( average_radiation ) THEN
    23572449          pt1   = 0.0_wp
    23582450          IF ( bulk_cloud_model  .OR.  cloud_droplets )  ql1   = 0.0_wp
     
    23642456          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    23652457          CALL MPI_ALLREDUCE( pt1_l, pt1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
    2366           IF ( bulk_cloud_model  .OR.  cloud_droplets )                            &
    2367              CALL MPI_ALLREDUCE( ql1_l, ql1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
     2458          IF ( ierr /= 0 ) THEN
     2459              WRITE(9,*) 'Error MPI_AllReduce1:', ierr, pt1_l, pt1
     2460              FLUSH(9)
     2461          ENDIF
     2462
     2463          IF ( bulk_cloud_model  .OR.  cloud_droplets ) THEN
     2464              CALL MPI_ALLREDUCE( ql1_l, ql1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
     2465              IF ( ierr /= 0 ) THEN
     2466                  WRITE(9,*) 'Error MPI_AllReduce2:', ierr, ql1_l, ql1
     2467                  FLUSH(9)
     2468              ENDIF
     2469          ENDIF
    23682470#else
    23692471          pt1 = pt1_l
     
    25342636          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    25352637          CALL MPI_ALLREDUCE( pt1_l, pt1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
    2536           IF ( bulk_cloud_model  .OR.  cloud_droplets )                             &
     2638          IF ( ierr /= 0 ) THEN
     2639              WRITE(9,*) 'Error MPI_AllReduce3:', ierr, pt1_l, pt1
     2640              FLUSH(9)
     2641          ENDIF
     2642          IF ( bulk_cloud_model  .OR.  cloud_droplets )  THEN
    25372643             CALL MPI_ALLREDUCE( ql1_l, ql1, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
     2644             IF ( ierr /= 0 ) THEN
     2645                 WRITE(9,*) 'Error MPI_AllReduce4:', ierr, ql1_l, ql1
     2646                 FLUSH(9)
     2647             ENDIF
     2648          ENDIF
    25382649#else
    25392650          pt1 = pt1_l
     
    27562867                                  albedo_lw_dif, albedo_sw_dir, albedo_sw_dif, &
    27572868                                  constant_albedo, dt_radiation, emissivity,   &
    2758                                   lw_radiation, net_radiation,                 &
     2869                                  lw_radiation, mrt_nlevels, mrt_skip_roof,    &
     2870                                  mrt_include_sw,  net_radiation,              &
    27592871                                  radiation_scheme, skip_time_do_radiation,    &
    27602872                                  sw_radiation, unscheduled_radiation_calls,   &
     2873                                  read_svf_on_init, write_svf_on_init,         &
    27612874                                  max_raytracing_dist, min_irrf_value,         &
    2762                                   nrefsteps, mrt_factors, rma_lad_raytrace,    &
     2875                                  nrefsteps, raytrace_mpi_rma,                 &
    27632876                                  dist_max_svf,                                &
    27642877                                  surface_reflections, svfnorm_report_thresh,  &
    2765                                   radiation_interactions_on
     2878                                  radiation_interactions_on,                   &
     2879                                  rad_angular_discretization,                  &
     2880                                  raytrace_discrete_azims,                     &
     2881                                  raytrace_discrete_elevs
    27662882   
    27672883       NAMELIST /radiation_parameters/   albedo, albedo_type, albedo_lw_dir,   &
    27682884                                  albedo_lw_dif, albedo_sw_dir, albedo_sw_dif, &
    27692885                                  constant_albedo, dt_radiation, emissivity,   &
    2770                                   lw_radiation, net_radiation,                 &
     2886                                  lw_radiation, mrt_nlevels, mrt_skip_roof,    &
     2887                                  mrt_include_sw,  net_radiation,              &
    27712888                                  radiation_scheme, skip_time_do_radiation,    &
    27722889                                  sw_radiation, unscheduled_radiation_calls,   &
    27732890                                  max_raytracing_dist, min_irrf_value,         &
    2774                                   nrefsteps, mrt_factors, rma_lad_raytrace,    &
     2891                                  nrefsteps, raytrace_mpi_rma,                 &
    27752892                                  dist_max_svf,                                &
    27762893                                  surface_reflections, svfnorm_report_thresh,  &
    2777                                   radiation_interactions_on
     2894                                  radiation_interactions_on,                   &
     2895                                  rad_angular_discretization,                  &
     2896                                  raytrace_discrete_azims,                     &
     2897                                  raytrace_discrete_elevs
    27782898   
    27792899       line = ' '
     
    44564576     INTEGER(iwp)                      :: i, j, k, kk, is, js, d, ku, refstep, m, mm, l, ll
    44574577     INTEGER(iwp)                      :: isurf, isurfsrc, isvf, icsf, ipcgb
     4578     INTEGER(iwp)                      :: imrt, imrtf
    44584579     INTEGER(iwp)                      :: isd                !< solar direction number
    44594580     INTEGER(iwp)                      :: pc_box_dimshift    !< transform for best accuracy
     
    45474668     surfoutsl(:)  = 0.0_wp !start-end
    45484669     surfoutll(:)  = 0.0_wp !start-end
     4670     IF ( nmrtbl > 0 )  THEN
     4671        mrtinsw(:) = 0._wp
     4672        mrtinlw(:) = 0._wp
     4673     ENDIF
     4674
    45494675
    45504676!--  Set up thermal radiation from surfaces
     
    46234749     CALL MPI_AllGatherv(surfoutll, nsurfl, MPI_REAL, &
    46244750                         surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr) !nsurf global
     4751     IF ( ierr /= 0 ) THEN
     4752         WRITE(9,*) 'Error MPI_AllGatherv1:', ierr, SIZE(surfoutll), nsurfl, &
     4753                     SIZE(surfoutl), nsurfs, surfstart
     4754         FLUSH(9)
     4755     ENDIF
    46254756#else
    46264757     surfoutl(:) = surfoutll(:) !nsurf global
     
    46394770        ENDDO
    46404771     ENDIF
    4641 
    4642      !-- diffuse radiation using sky view factor, TODO: homogeneous rad_*w_in_diff because now it depends on no. of processors
    4643      surfinswdif(:) = rad_sw_in_diff(nyn,nxl) * skyvft(:)
    4644      surfinlwdif(:) = rad_lw_in_diff(nyn,nxl) * skyvf(:)
     4772!
     4773!--  diffuse radiation using sky view factor
     4774     DO isurf = 1, nsurfl
     4775        j = surfl(iy, isurf)
     4776        i = surfl(ix, isurf)
     4777        surfinswdif(isurf) = rad_sw_in_diff(j,i) * skyvft(isurf)
     4778        surfinlwdif(isurf) = rad_lw_in_diff(j,i) * skyvf(isurf)
     4779     ENDDO
     4780!
     4781!--  MRT diffuse irradiance
     4782     DO  imrt = 1, nmrtbl
     4783        j = mrtbl(iy, imrt)
     4784        i = mrtbl(ix, imrt)
     4785        mrtinsw(imrt) = mrtskyt(imrt) * rad_sw_in_diff(j,i)
     4786        mrtinlw(imrt) = mrtsky(imrt) * rad_lw_in_diff(j,i)
     4787     ENDDO
    46454788
    46464789     !-- direct radiation
     
    46544797        isd = dsidir_rev(j, i)
    46554798        DO isurf = 1, nsurfl
    4656            surfinswdir(isurf) = rad_sw_in_dir(nyn,nxl) * costheta(surfl(id, isurf)) * dsitrans(isurf, isd) / zenith(0)
     4799           j = surfl(iy, isurf)
     4800           i = surfl(ix, isurf)
     4801           surfinswdir(isurf) = rad_sw_in_dir(j,i) * costheta(surfl(id, isurf)) * dsitrans(isurf, isd) / zenith(0)
     4802        ENDDO
     4803!
     4804!--     MRT direct irradiance
     4805        DO  imrt = 1, nmrtbl
     4806           j = mrtbl(iy, imrt)
     4807           i = mrtbl(ix, imrt)
     4808           mrtinsw(imrt) = mrtinsw(imrt) + mrtdsit(imrt, isd) * rad_sw_in_dir(j,i) &
     4809                                     / zenith(0) / 4._wp ! normal to sphere
    46574810        ENDDO
    46584811     ENDIF
     4812!
     4813!--  MRT first pass thermal
     4814     DO  imrtf = 1, nmrtf
     4815        imrt = mrtfsurf(1, imrtf)
     4816        isurfsrc = mrtfsurf(2, imrtf)
     4817        mrtinlw(imrt) = mrtinlw(imrt) + mrtf(imrtf) * surfoutl(isurfsrc)
     4818     ENDDO
    46594819
    46604820     IF ( npcbl > 0 )  THEN
     
    46744834             IF ( isurfsrc == -1 )  THEN
    46754835!--                 Diffuse rad from sky.
    4676                  pcbinswdif(ipcgb) = csf(1,icsf) * csf(2,icsf) * rad_sw_in_diff(j,i)
     4836                 pcbinswdif(ipcgb) = csf(1,icsf) * rad_sw_in_diff(j,i)
    46774837
    46784838                 !--Direct rad
     
    46854845                                        * pc_abs_frac * dsitransc(ipcgb, isd)
    46864846                 ENDIF
    4687 
    4688                  EXIT ! only isurfsrc=-1 is processed here
    46894847             ENDIF
    46904848         ENDDO
     
    47224880         CALL MPI_AllGatherv(surfoutsl, nsurfl, MPI_REAL, &
    47234881             surfouts, nsurfs, surfstart, MPI_REAL, comm2d, ierr)
     4882         IF ( ierr /= 0 ) THEN
     4883             WRITE(9,*) 'Error MPI_AllGatherv2:', ierr, SIZE(surfoutsl), nsurfl, &
     4884                        SIZE(surfouts), nsurfs, surfstart
     4885             FLUSH(9)
     4886         ENDIF
     4887
    47244888         CALL MPI_AllGatherv(surfoutll, nsurfl, MPI_REAL, &
    47254889             surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr)
     4890         IF ( ierr /= 0 ) THEN
     4891             WRITE(9,*) 'Error MPI_AllGatherv3:', ierr, SIZE(surfoutll), nsurfl, &
     4892                        SIZE(surfoutl), nsurfs, surfstart
     4893             FLUSH(9)
     4894         ENDIF
     4895
    47264896#else
    47274897         surfouts = surfoutsl
     
    47474917             IF ( isurfsrc == -1 )  CYCLE ! sky->face only in 1st pass, not here
    47484918
    4749              pcbinsw(ipcgb) = pcbinsw(ipcgb) + csf(1,icsf) * csf(2,icsf) * surfouts(isurfsrc)
     4919             pcbinsw(ipcgb) = pcbinsw(ipcgb) + csf(1,icsf) * surfouts(isurfsrc)
     4920         ENDDO
     4921!
     4922!--      MRT reflected
     4923         DO  imrtf = 1, nmrtf
     4924            imrt = mrtfsurf(1, imrtf)
     4925            isurfsrc = mrtfsurf(2, imrtf)
     4926            mrtinsw(imrt) = mrtinsw(imrt) + mrtft(imrtf) * surfouts(isurfsrc)
     4927            mrtinlw(imrt) = mrtinlw(imrt) + mrtf(imrtf) * surfoutl(isurfsrc)
    47504928         ENDDO
    47514929
     
    47554933         surfoutlw = surfoutlw + surfoutll
    47564934
    4757      ENDDO
     4935     ENDDO ! refstep
    47584936
    47594937!--  push heat flux absorbed by plant canopy to respective 3D arrays
     
    47784956     ENDIF
    47794957!
     4958!--  Calculate black body MRT (after all reflections)
     4959     IF ( nmrtbl > 0 )  THEN
     4960        IF ( mrt_include_sw )  THEN
     4961           mrt(:) = ((mrtinsw(:) + mrtinlw(:)) / sigma_sb) ** .25_wp
     4962        ELSE
     4963           mrt(:) = (mrtinlw(:) / sigma_sb) ** .25_wp
     4964        ENDIF
     4965     ENDIF
     4966!
    47804967!--     Transfer radiation arrays required for energy balance to the respective data types
    47814968     DO  i = 1, nsurfl
     
    47914978           surf_usm_h%rad_net(m)    = surfinsw(i) - surfoutsw(i) +          &
    47924979                                      surfinlw(i) - surfoutlw(i)
     4980           surf_usm_h%rad_net_l(m)  = surf_usm_h%rad_net(m)
    47934981!
    47944982!--     northward-facding
     
    48004988           surf_usm_v(0)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
    48014989                                         surfinlw(i) - surfoutlw(i)
     4990           surf_usm_v(0)%rad_net_l(m)  = surf_usm_v(0)%rad_net(m)
    48024991!
    48034992!--     southward-facding
     
    48094998           surf_usm_v(1)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
    48104999                                         surfinlw(i) - surfoutlw(i)
     5000           surf_usm_v(1)%rad_net_l(m)  = surf_usm_v(1)%rad_net(m)
    48115001!
    48125002!--     eastward-facing
     
    48185008           surf_usm_v(2)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
    48195009                                         surfinlw(i) - surfoutlw(i)
     5010           surf_usm_v(2)%rad_net_l(m)  = surf_usm_v(2)%rad_net(m)
    48205011!
    48215012!--     westward-facding
     
    48275018           surf_usm_v(3)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
    48285019                                         surfinlw(i) - surfoutlw(i)
     5020           surf_usm_v(3)%rad_net_l(m)  = surf_usm_v(3)%rad_net(m)
    48295021!
    48305022!--     (2) land surfaces
     
    49575149#if defined( __parallel )
    49585150     CALL MPI_ALLREDUCE( pinswl, pinsw, 1, MPI_REAL, MPI_SUM, comm2d, ierr)
     5151     IF ( ierr /= 0 ) THEN
     5152         WRITE(9,*) 'Error MPI_AllReduce5:', ierr, pinswl, pinsw
     5153         FLUSH(9)
     5154     ENDIF
    49595155     CALL MPI_ALLREDUCE( pinlwl, pinlw, 1, MPI_REAL, MPI_SUM, comm2d, ierr)
     5156     IF ( ierr /= 0 ) THEN
     5157         WRITE(9,*) 'Error MPI_AllReduce6:', ierr, pinlwl, pinlw
     5158         FLUSH(9)
     5159     ENDIF
    49605160     CALL MPI_ALLREDUCE( pabsswl, pabssw, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
     5161     IF ( ierr /= 0 ) THEN
     5162         WRITE(9,*) 'Error MPI_AllReduce7:', ierr, pabsswl, pabssw
     5163         FLUSH(9)
     5164     ENDIF
    49615165     CALL MPI_ALLREDUCE( pabslwl, pabslw, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
     5166     IF ( ierr /= 0 ) THEN
     5167         WRITE(9,*) 'Error MPI_AllReduce8:', ierr, pabslwl, pabslw
     5168         FLUSH(9)
     5169     ENDIF
    49625170     CALL MPI_ALLREDUCE( pemitlwl, pemitlw, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
     5171     IF ( ierr /= 0 ) THEN
     5172         WRITE(9,*) 'Error MPI_AllReduce8:', ierr, pemitlwl, pemitlw
     5173         FLUSH(9)
     5174     ENDIF
    49635175     CALL MPI_ALLREDUCE( emiss_sum_surfl, emiss_sum_surf, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
     5176     IF ( ierr /= 0 ) THEN
     5177         WRITE(9,*) 'Error MPI_AllReduce9:', ierr, emiss_sum_surfl, emiss_sum_surf
     5178         FLUSH(9)
     5179     ENDIF
    49645180     CALL MPI_ALLREDUCE( area_surfl, area_surf, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
     5181     IF ( ierr /= 0 ) THEN
     5182         WRITE(9,*) 'Error MPI_AllReduce10:', ierr, area_surfl, area_surf
     5183         FLUSH(9)
     5184     ENDIF
    49655185#else
    49665186     pinsw = pinswl
     
    50385258!      t_rad_urb = pt_surf_urb / count_surfaces
    50395259     
     5260
    50405261    CONTAINS
    50415262
     
    52015422       INTEGER(iwp) :: i, j, k, l, m
    52025423       INTEGER(iwp) :: k_topo     !< vertical index indicating topography top for given (j,i)
    5203        INTEGER(iwp) :: nzptl, nzubl, nzutl, isurf, ipcgb
     5424       INTEGER(iwp) :: nzptl, nzubl, nzutl, isurf, ipcgb, imrt
    52045425       REAL(wp)     :: mrl
     5426#if defined( __parallel )
     5427       INTEGER(iwp), DIMENSION(:), POINTER       ::  gridsurf_rma   !< fortran pointer, but lower bounds are 1
     5428       TYPE(c_ptr)                               ::  gridsurf_rma_p !< allocated c pointer
     5429       INTEGER(iwp)                              ::  minfo          !< MPI RMA window info handle
     5430#endif
    52055431
    52065432
     
    52635489#if defined( __parallel )
    52645490       CALL MPI_AllReduce(nzubl, nzub, 1, MPI_INTEGER, MPI_MIN, comm2d, ierr )
     5491       IF ( ierr /= 0 ) THEN
     5492           WRITE(9,*) 'Error MPI_AllReduce11:', ierr, nzubl, nzub
     5493           FLUSH(9)
     5494       ENDIF
    52655495       CALL MPI_AllReduce(nzutl, nzut, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )
     5496       IF ( ierr /= 0 ) THEN
     5497           WRITE(9,*) 'Error MPI_AllReduce12:', ierr, nzutl, nzut
     5498           FLUSH(9)
     5499       ENDIF
    52665500       CALL MPI_AllReduce(nzptl, nzpt, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr )
     5501       IF ( ierr /= 0 ) THEN
     5502           WRITE(9,*) 'Error MPI_AllReduce13:', ierr, nzptl, nzpt
     5503           FLUSH(9)
     5504       ENDIF
    52675505#else
    52685506       nzub = nzubl
     
    53535591!--    cycles must not be altered, certain file input routines may depend
    53545592!--    on it)
    5355        ALLOCATE(surfl(5,nsurfl))  ! is it mecessary to allocate it with (5,nsurfl)?
     5593       ALLOCATE(surfl_l(5*nsurfl))  ! is it necessary to allocate it with (5,nsurfl)?
     5594       surfl(1:5,1:nsurfl) => surfl_l(1:5*nsurfl)
    53565595       isurf = 0
     5596       IF ( rad_angular_discretization )  THEN
     5597!
     5598!--       Allocate and fill the reverse indexing array gridsurf
     5599#if defined( __parallel )
     5600!
     5601!--       raytrace_mpi_rma is asserted
     5602
     5603          CALL MPI_Info_create(minfo, ierr)
     5604          IF ( ierr /= 0 ) THEN
     5605              WRITE(9,*) 'Error MPI_Info_create1:', ierr
     5606              FLUSH(9)
     5607          ENDIF
     5608          CALL MPI_Info_set(minfo, 'accumulate_ordering', '', ierr)
     5609          IF ( ierr /= 0 ) THEN
     5610              WRITE(9,*) 'Error MPI_Info_set1:', ierr
     5611              FLUSH(9)
     5612          ENDIF
     5613          CALL MPI_Info_set(minfo, 'accumulate_ops', 'same_op', ierr)
     5614          IF ( ierr /= 0 ) THEN
     5615              WRITE(9,*) 'Error MPI_Info_set2:', ierr
     5616              FLUSH(9)
     5617          ENDIF
     5618          CALL MPI_Info_set(minfo, 'same_size', 'true', ierr)
     5619          IF ( ierr /= 0 ) THEN
     5620              WRITE(9,*) 'Error MPI_Info_set3:', ierr
     5621              FLUSH(9)
     5622          ENDIF
     5623          CALL MPI_Info_set(minfo, 'same_disp_unit', 'true', ierr)
     5624          IF ( ierr /= 0 ) THEN
     5625              WRITE(9,*) 'Error MPI_Info_set4:', ierr
     5626              FLUSH(9)
     5627          ENDIF
     5628
     5629          CALL MPI_Win_allocate(INT(STORAGE_SIZE(1_iwp)/8*nsurf_type_u*nzu*nny*nnx, &
     5630                                    kind=MPI_ADDRESS_KIND),                         &
     5631                                INT(STORAGE_SIZE(1_iwp)/8, kind=MPI_ADDRESS_KIND),  &
     5632                                minfo, comm2d, gridsurf_rma_p, win_gridsurf, ierr)
     5633          IF ( ierr /= 0 ) THEN
     5634              WRITE(9,*) 'Error MPI_Win_allocate1:', ierr, &
     5635                 INT(STORAGE_SIZE(1_iwp)/8*nsurf_type_u*nzu*nny*nnx,kind=MPI_ADDRESS_KIND), &
     5636                 INT(STORAGE_SIZE(1_iwp)/8, kind=MPI_ADDRESS_KIND), win_gridsurf
     5637              FLUSH(9)
     5638          ENDIF
     5639
     5640          CALL MPI_Info_free(minfo, ierr)
     5641          IF ( ierr /= 0 ) THEN
     5642              WRITE(9,*) 'Error MPI_Info_free1:', ierr
     5643              FLUSH(9)
     5644          ENDIF
     5645
     5646!
     5647!--       On Intel compilers, calling c_f_pointer to transform a C pointer
     5648!--       directly to a multi-dimensional Fotran pointer leads to strange
     5649!--       errors on dimension boundaries. However, transforming to a 1D
     5650!--       pointer and then redirecting a multidimensional pointer to it works
     5651!--       fine.
     5652          CALL c_f_pointer(gridsurf_rma_p, gridsurf_rma, (/ nsurf_type_u*nzu*nny*nnx /))
     5653          gridsurf(0:nsurf_type_u-1, nzub:nzut, nys:nyn, nxl:nxr) =>                &
     5654                     gridsurf_rma(1:nsurf_type_u*nzu*nny*nnx)
     5655#else
     5656          ALLOCATE(gridsurf(0:nsurf_type_u-1,nzub:nzut,nys:nyn,nxl:nxr) )
     5657#endif
     5658          gridsurf(:,:,:,:) = -999
     5659       ENDIF
    53575660
    53585661!--    add horizontal surface elements (land and urban surfaces)
     
    53625665              DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
    53635666                 k = surf_usm_h%k(m)
    5364 
    53655667                 isurf = isurf + 1
    53665668                 surfl(:,isurf) = (/iup_u,k,j,i,m/)
     5669                 IF ( rad_angular_discretization ) THEN
     5670                    gridsurf(iup_u,k,j,i) = isurf
     5671                 ENDIF
    53675672              ENDDO
    53685673
    53695674              DO  m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
    53705675                 k = surf_lsm_h%k(m)
    5371 
    53725676                 isurf = isurf + 1
    53735677                 surfl(:,isurf) = (/iup_l,k,j,i,m/)
     5678                 IF ( rad_angular_discretization ) THEN
     5679                    gridsurf(iup_u,k,j,i) = isurf
     5680                 ENDIF
    53745681              ENDDO
    53755682
     
    53845691              DO  m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
    53855692                 k = surf_usm_v(l)%k(m)
    5386 
    5387                  isurf          = isurf + 1
     5693                 isurf = isurf + 1
    53885694                 surfl(:,isurf) = (/inorth_u,k,j,i,m/)
     5695                 IF ( rad_angular_discretization ) THEN
     5696                    gridsurf(inorth_u,k,j,i) = isurf
     5697                 ENDIF
    53895698              ENDDO
    53905699              DO  m = surf_lsm_v(l)%start_index(j,i), surf_lsm_v(l)%end_index(j,i)
    53915700                 k = surf_lsm_v(l)%k(m)
    5392 
    5393                  isurf          = isurf + 1
     5701                 isurf = isurf + 1
    53945702                 surfl(:,isurf) = (/inorth_l,k,j,i,m/)
     5703                 IF ( rad_angular_discretization ) THEN
     5704                    gridsurf(inorth_u,k,j,i) = isurf
     5705                 ENDIF
    53955706              ENDDO
    53965707
     
    53985709              DO  m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
    53995710                 k = surf_usm_v(l)%k(m)
    5400 
    5401                  isurf          = isurf + 1
     5711                 isurf = isurf + 1
    54025712                 surfl(:,isurf) = (/isouth_u,k,j,i,m/)
     5713                 IF ( rad_angular_discretization ) THEN
     5714                    gridsurf(isouth_u,k,j,i) = isurf
     5715                 ENDIF
    54035716              ENDDO
    54045717              DO  m = surf_lsm_v(l)%start_index(j,i), surf_lsm_v(l)%end_index(j,i)
    54055718                 k = surf_lsm_v(l)%k(m)
    5406 
    5407                  isurf          = isurf + 1
     5719                 isurf = isurf + 1
    54085720                 surfl(:,isurf) = (/isouth_l,k,j,i,m/)
     5721                 IF ( rad_angular_discretization ) THEN
     5722                    gridsurf(isouth_u,k,j,i) = isurf
     5723                 ENDIF
    54095724              ENDDO
    54105725
     
    54125727              DO  m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
    54135728                 k = surf_usm_v(l)%k(m)
    5414 
    5415                  isurf          = isurf + 1
     5729                 isurf = isurf + 1
    54165730                 surfl(:,isurf) = (/ieast_u,k,j,i,m/)
     5731                 IF ( rad_angular_discretization ) THEN
     5732                    gridsurf(ieast_u,k,j,i) = isurf
     5733                 ENDIF
    54175734              ENDDO
    54185735              DO  m = surf_lsm_v(l)%start_index(j,i), surf_lsm_v(l)%end_index(j,i)
    54195736                 k = surf_lsm_v(l)%k(m)
    5420 
    5421                  isurf          = isurf + 1
     5737                 isurf = isurf + 1
    54225738                 surfl(:,isurf) = (/ieast_l,k,j,i,m/)
     5739                 IF ( rad_angular_discretization ) THEN
     5740                    gridsurf(ieast_u,k,j,i) = isurf
     5741                 ENDIF
    54235742              ENDDO
    54245743
     
    54265745              DO  m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
    54275746                 k = surf_usm_v(l)%k(m)
    5428 
    5429                  isurf          = isurf + 1
     5747                 isurf = isurf + 1
    54305748                 surfl(:,isurf) = (/iwest_u,k,j,i,m/)
     5749                 IF ( rad_angular_discretization ) THEN
     5750                    gridsurf(iwest_u,k,j,i) = isurf
     5751                 ENDIF
    54315752              ENDDO
    54325753              DO  m = surf_lsm_v(l)%start_index(j,i), surf_lsm_v(l)%end_index(j,i)
    54335754                 k = surf_lsm_v(l)%k(m)
    5434 
    5435                  isurf          = isurf + 1
     5755                 isurf = isurf + 1
    54365756                 surfl(:,isurf) = (/iwest_l,k,j,i,m/)
     5757                 IF ( rad_angular_discretization ) THEN
     5758                    gridsurf(iwest_u,k,j,i) = isurf
     5759                 ENDIF
    54375760              ENDDO
    54385761           ENDDO
    54395762       ENDDO
     5763!
     5764!--    Add local MRT boxes for specified number of levels
     5765       nmrtbl = 0
     5766       IF ( mrt_nlevels > 0 )  THEN
     5767          DO  i = nxl, nxr
     5768             DO  j = nys, nyn
     5769                DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
     5770!
     5771!--                Skip roof if requested
     5772                   IF ( mrt_skip_roof  .AND.  surf_usm_h%isroof_surf(m) )  CYCLE
     5773!
     5774!--                Cycle over specified no of levels
     5775                   nmrtbl = nmrtbl + mrt_nlevels
     5776                ENDDO
     5777!
     5778!--             Dtto for LSM
     5779                DO  m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
     5780                   nmrtbl = nmrtbl + mrt_nlevels
     5781                ENDDO
     5782             ENDDO
     5783          ENDDO
     5784
     5785          ALLOCATE( mrtbl(iz:ix,nmrtbl), mrtsky(nmrtbl), mrtskyt(nmrtbl), &
     5786                    mrtinsw(nmrtbl), mrtinlw(nmrtbl), mrt(nmrtbl) )
     5787
     5788          imrt = 0
     5789          DO  i = nxl, nxr
     5790             DO  j = nys, nyn
     5791                DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
     5792!
     5793!--                Skip roof if requested
     5794                   IF ( mrt_skip_roof  .AND.  surf_usm_h%isroof_surf(m) )  CYCLE
     5795!
     5796!--                Cycle over specified no of levels
     5797                   l = surf_usm_h%k(m)
     5798                   DO  k = l, l + mrt_nlevels - 1
     5799                      imrt = imrt + 1
     5800                      mrtbl(:,imrt) = (/k,j,i/)
     5801                   ENDDO
     5802                ENDDO
     5803!
     5804!--             Dtto for LSM
     5805                DO  m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
     5806                   l = surf_lsm_h%k(m)
     5807                   DO  k = l, l + mrt_nlevels - 1
     5808                      imrt = imrt + 1
     5809                      mrtbl(:,imrt) = (/k,j,i/)
     5810                   ENDDO
     5811                ENDDO
     5812             ENDDO
     5813          ENDDO
     5814       ENDIF
    54405815
    54415816!
     
    54585833#if defined( __parallel )
    54595834       CALL MPI_Allgather(nsurfl,1,MPI_INTEGER,nsurfs,1,MPI_INTEGER,comm2d,ierr)
     5835       IF ( ierr /= 0 ) THEN
     5836         WRITE(9,*) 'Error MPI_AllGather1:', ierr, nsurfl, nsurfs
     5837         FLUSH(9)
     5838     ENDIF
     5839
    54605840#else
    54615841       nsurfs(0) = nsurfl
     
    54695849       surfstart(numprocs) = k
    54705850       nsurf = k
    5471        ALLOCATE(surf(5,nsurf))
     5851       ALLOCATE(surf_l(5*nsurf))
     5852       surf(1:5,1:nsurf) => surf_l(1:5*nsurf)
    54725853
    54735854#if defined( __parallel )
    5474        CALL MPI_AllGatherv(surfl, nsurfl*5, MPI_INTEGER, surf, nsurfs*5, &
     5855       CALL MPI_AllGatherv(surfl_l, nsurfl*5, MPI_INTEGER, surf_l, nsurfs*5, &
    54755856           surfstart(0:numprocs-1)*5, MPI_INTEGER, comm2d, ierr)
     5857       IF ( ierr /= 0 ) THEN
     5858           WRITE(9,*) 'Error MPI_AllGatherv4:', ierr, SIZE(surfl_l), nsurfl*5, &
     5859                      SIZE(surf_l), nsurfs*5, surfstart(0:numprocs-1)*5
     5860           FLUSH(9)
     5861       ENDIF
    54765862#else
    54775863       surf = surfl
     
    55345920       
    55355921        INTEGER(iwp)                                  :: i, j, k, d, ip, jp
    5536         INTEGER(iwp)                                  :: isvf, ksvf, icsf, kcsf, npcsfl, isvf_surflt, imrtt, imrtf, ipcgb
     5922        INTEGER(iwp)                                  :: isvf, ksvf, icsf, kcsf, npcsfl, isvf_surflt, imrt, imrtf, ipcgb
    55375923        INTEGER(iwp)                                  :: sd, td
    55385924        INTEGER(iwp)                                  :: iaz, izn      !< azimuth, zenith counters
     
    55425928        REAL(wp)                                      :: az1, az2      !< relative azimuth of section borders
    55435929        REAL(wp)                                      :: azmid         !< ray (center) azimuth
    5544         REAL(wp)                                      :: horizon       !< computed horizon height (tangent of elevation)
    5545         REAL(wp)                                      :: azen          !< zenith angle
    55465930        REAL(wp), DIMENSION(2)                        :: yxdir         !< y,x *unit* vector of ray direction (in grid units)
    55475931        REAL(wp), DIMENSION(:), ALLOCATABLE           :: zdirs         !< directions in z (tangent of elevation)
    55485932        REAL(wp), DIMENSION(:), ALLOCATABLE           :: zbdry         !< zenith angle boundaries
    55495933        REAL(wp), DIMENSION(:), ALLOCATABLE           :: vffrac        !< view factor fractions for individual rays
     5934        REAL(wp), DIMENSION(:), ALLOCATABLE           :: vffrac0       !< dtto (original values)
    55505935        REAL(wp), DIMENSION(:), ALLOCATABLE           :: ztransp       !< array of transparency in z steps
     5936        INTEGER(iwp)                                  :: lowest_free_ray !< index into zdirs
     5937        INTEGER(iwp), DIMENSION(:), ALLOCATABLE       :: itarget       !< face indices of detected obstacles
     5938        INTEGER(iwp)                                  :: itarg0, itarg1
     5939#if defined( __parallel )
     5940#endif
     5941
     5942
     5943
    55515944        REAL(wp),     DIMENSION(0:nsurf_type)         :: facearea
    5552         INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE     :: nzterrl
    5553         REAL(wp),     DIMENSION(:,:), ALLOCATABLE     :: csflt, pcsflt
    5554         INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE     :: kcsflt,kpcsflt
     5945        INTEGER(iwp)                                  :: udim
     5946        INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET:: nzterrl_l
     5947        INTEGER(iwp), DIMENSION(:,:), POINTER         :: nzterrl
     5948        REAL(wp),     DIMENSION(:), ALLOCATABLE,TARGET:: csflt_l, pcsflt_l
     5949        REAL(wp),     DIMENSION(:,:), POINTER         :: csflt, pcsflt
     5950        INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET:: kcsflt_l,kpcsflt_l
     5951        INTEGER(iwp), DIMENSION(:,:), POINTER         :: kcsflt,kpcsflt
    55555952        INTEGER(iwp), DIMENSION(:), ALLOCATABLE       :: icsflt,dcsflt,ipcsflt,dpcsflt
    55565953        REAL(wp), DIMENSION(3)                        :: uv
     
    55615958        INTEGER(idp)                                  :: ray_skip_maxdist, ray_skip_minval !< skipped raytracing counts
    55625959        INTEGER(iwp)                                  :: max_track_len !< maximum 2d track length
    5563         INTEGER(iwp)                                  :: win_lad, minfo
    5564         REAL(wp), DIMENSION(:,:,:), POINTER           :: lad_s_rma       !< fortran pointer, but lower bounds are 1
     5960        INTEGER(iwp)                                  :: minfo
     5961        REAL(wp), DIMENSION(:), POINTER               :: lad_s_rma       !< fortran 1D pointer
    55655962        TYPE(c_ptr)                                   :: lad_s_rma_p     !< allocated c pointer
    55665963#if defined( __parallel )
     
    55695966!   
    55705967        INTEGER(iwp), DIMENSION(0:svfnorm_report_num) :: svfnorm_counts
    5571 !        CHARACTER(200)                                :: msg
     5968        CHARACTER(200)                                :: msg
    55725969
    55735970!--     calculation of the SVF
    55745971        CALL location_message( '    calculation of SVF and CSF', .TRUE. )
    5575 !         CALL radiation_write_debug_log('Start calculation of SVF and CSF')
     5972        CALL radiation_write_debug_log('Start calculation of SVF and CSF')
    55765973
    55775974!--     precalculate face areas for different face directions using normal vector
     
    55965993            acsf => acsf1
    55975994        ENDIF
     5995        nmrtf = 0
     5996        IF ( mrt_nlevels > 0 )  THEN
     5997           nmrtfa = gasize
     5998           mmrtf = 1
     5999           ALLOCATE ( amrtf1(nmrtfa) )
     6000           amrtf => amrtf1
     6001        ENDIF
    55986002        ray_skip_maxdist = 0
    55996003        ray_skip_minval = 0
     
    56026006        ALLOCATE( nzterr(0:(nx+1)*(ny+1)-1) )
    56036007#if defined( __parallel )
    5604         ALLOCATE( nzterrl(nys:nyn,nxl:nxr) )
     6008        !ALLOCATE( nzterrl(nys:nyn,nxl:nxr) )
     6009        ALLOCATE( nzterrl_l((nyn-nys+1)*(nxr-nxl+1)) )
     6010        nzterrl(nys:nyn,nxl:nxr) => nzterrl_l(1:(nyn-nys+1)*(nxr-nxl+1))
    56056011        nzterrl = get_topography_top_index( 's' )
    5606         CALL MPI_AllGather( nzterrl, nnx*nny, MPI_INTEGER, &
     6012        CALL MPI_AllGather( nzterrl_l, nnx*nny, MPI_INTEGER, &
    56076013                            nzterr, nnx*nny, MPI_INTEGER, comm2d, ierr )
    5608         DEALLOCATE(nzterrl)
     6014        IF ( ierr /= 0 ) THEN
     6015            WRITE(9,*) 'Error MPI_AllGather1:', ierr, SIZE(nzterrl_l), nnx*nny, &
     6016                       SIZE(nzterr), nnx*nny
     6017            FLUSH(9)
     6018        ENDIF
     6019        DEALLOCATE(nzterrl_l)
    56096020#else
    56106021        nzterr = RESHAPE( get_topography_top_index( 's' ), (/(nx+1)*(ny+1)/) )
     
    56216032            CALL MPI_AllGather( pct, nnx*nny, MPI_INTEGER, &
    56226033                                plantt, nnx*nny, MPI_INTEGER, comm2d, ierr )
    5623            
     6034            IF ( ierr /= 0 ) THEN
     6035                WRITE(9,*) 'Error MPI_AllGather2:', ierr, SIZE(pct), nnx*nny, &
     6036                           SIZE(plantt), nnx*nny
     6037                FLUSH(9)
     6038            ENDIF
     6039
    56246040!--         temporary arrays storing values for csf calculation during raytracing
    56256041            ALLOCATE( lad_ip(maxboxesg) )
    56266042            ALLOCATE( lad_disp(maxboxesg) )
    56276043
    5628             IF ( rma_lad_raytrace )  THEN
     6044            IF ( raytrace_mpi_rma )  THEN
    56296045                ALLOCATE( lad_s_ray(maxboxesg) )
    56306046               
    56316047                ! set conditions for RMA communication
    56326048                CALL MPI_Info_create(minfo, ierr)
     6049                IF ( ierr /= 0 ) THEN
     6050                    WRITE(9,*) 'Error MPI_Info_create2:', ierr
     6051                    FLUSH(9)
     6052                ENDIF
    56336053                CALL MPI_Info_set(minfo, 'accumulate_ordering', '', ierr)
     6054                IF ( ierr /= 0 ) THEN
     6055                    WRITE(9,*) 'Error MPI_Info_set5:', ierr
     6056                    FLUSH(9)
     6057                ENDIF
    56346058                CALL MPI_Info_set(minfo, 'accumulate_ops', 'same_op', ierr)
     6059                IF ( ierr /= 0 ) THEN
     6060                    WRITE(9,*) 'Error MPI_Info_set6:', ierr
     6061                    FLUSH(9)
     6062                ENDIF
    56356063                CALL MPI_Info_set(minfo, 'same_size', 'true', ierr)
     6064                IF ( ierr /= 0 ) THEN
     6065                    WRITE(9,*) 'Error MPI_Info_set7:', ierr
     6066                    FLUSH(9)
     6067                ENDIF
    56366068                CALL MPI_Info_set(minfo, 'same_disp_unit', 'true', ierr)
     6069                IF ( ierr /= 0 ) THEN
     6070                    WRITE(9,*) 'Error MPI_Info_set8:', ierr
     6071                    FLUSH(9)
     6072                ENDIF
    56376073
    56386074!--             Allocate and initialize the MPI RMA window
     
    56436079                CALL MPI_Win_allocate(size_lad_rma, STORAGE_SIZE(1.0_wp)/8, minfo, comm2d, &
    56446080                                        lad_s_rma_p, win_lad, ierr)
    5645                 CALL c_f_pointer(lad_s_rma_p, lad_s_rma, (/ nzp, nny, nnx /))
    5646                 sub_lad(nzub:, nys:, nxl:) => lad_s_rma(:,:,:)
     6081                IF ( ierr /= 0 ) THEN
     6082                    WRITE(9,*) 'Error MPI_Win_allocate2:', ierr, size_lad_rma, &
     6083                                STORAGE_SIZE(1.0_wp)/8, win_lad
     6084                    FLUSH(9)
     6085                ENDIF
     6086                CALL c_f_pointer(lad_s_rma_p, lad_s_rma, (/ nzp*nny*nnx /))
     6087                sub_lad(nzub:nzpt, nys:nyn, nxl:nxr) => lad_s_rma(1:nzp*nny*nnx)
    56476088            ELSE
    56486089                ALLOCATE(sub_lad(nzub:nzpt, nys:nyn, nxl:nxr))
     
    56666107
    56676108#if defined( __parallel )
    5668             IF ( rma_lad_raytrace )  THEN
     6109            IF ( raytrace_mpi_rma )  THEN
    56696110                CALL MPI_Info_free(minfo, ierr)
     6111                IF ( ierr /= 0 ) THEN
     6112                    WRITE(9,*) 'Error MPI_Info_free2:', ierr
     6113                    FLUSH(9)
     6114                ENDIF
    56706115                CALL MPI_Win_lock_all(0, win_lad, ierr)
     6116                IF ( ierr /= 0 ) THEN
     6117                    WRITE(9,*) 'Error MPI_Win_lock_all1:', ierr, win_lad
     6118                    FLUSH(9)
     6119                ENDIF
     6120               
    56716121            ELSE
    56726122                ALLOCATE( sub_lad_g(0:(nx+1)*(ny+1)*nzp-1) )
    56736123                CALL MPI_AllGather( sub_lad, nnx*nny*nzp, MPI_REAL, &
    56746124                                    sub_lad_g, nnx*nny*nzp, MPI_REAL, comm2d, ierr )
     6125                IF ( ierr /= 0 ) THEN
     6126                    WRITE(9,*) 'Error MPI_AllGather3:', ierr, SIZE(sub_lad), &
     6127                                nnx*nny*nzp, SIZE(sub_lad_g), nnx*nny*nzp
     6128                    FLUSH(9)
     6129                ENDIF
    56756130            ENDIF
    56766131#endif
    56776132        ENDIF
    56786133
    5679         IF ( mrt_factors )  THEN
    5680             OPEN(153, file='MRT_TARGETS', access='SEQUENTIAL', &
    5681                     action='READ', status='OLD', form='FORMATTED', err=524)
    5682             OPEN(154, file='MRT_FACTORS'//myid_char, access='DIRECT', recl=(5*4+2*8), &
    5683                     action='WRITE', status='REPLACE', form='UNFORMATTED', err=525)
    5684             imrtf = 1
    5685             DO
    5686                 READ(153, *, end=526, err=524) imrtt, i, j, k
    5687                 IF ( i < nxl  .OR.  i > nxr &
    5688                      .OR.  j < nys  .OR.  j > nyn ) CYCLE
    5689                 ta = (/ REAL(k), REAL(j), REAL(i) /)
    5690 
    5691                 DO isurfs = 1, nsurf
    5692                     IF ( .NOT.  surface_facing(i, j, k, -1, &
    5693                         surf(ix, isurfs), surf(iy, isurfs), &
    5694                         surf(iz, isurfs), surf(id, isurfs)) )  THEN
    5695                         CYCLE
    5696                     ENDIF
    5697                      
    5698                     sd = surf(id, isurfs)
    5699                     sa = (/ REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd), &
    5700                             REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd), &
    5701                             REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd) /)
    5702 
    5703 !--                 unit vector source -> target
    5704                     uv = (/ (ta(1)-sa(1))*dz(1), (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)
    5705                     sqdist = SUM(uv(:)**2)
    5706                     uv = uv / SQRT(sqdist)
    5707 
    5708 !--                 irradiance factor - see svf. Here we consider that target face is always normal,
    5709 !--                 i.e. the second dot product equals 1
    5710                     rirrf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) &
    5711                         / (pi * sqdist) * facearea(sd)
    5712 
    5713 !--                 raytrace while not creating any canopy sink factors
    5714                     CALL raytrace(sa, ta, isurfs, rirrf, 1._wp, .FALSE., &
    5715                             visible, transparency, win_lad)
    5716                     IF ( .NOT.  visible ) CYCLE
    5717 
    5718                     !rsvf = rirrf * transparency
    5719                     WRITE(154, rec=imrtf, err=525) INT(imrtt, kind=4), &
    5720                         INT(surf(id, isurfs), kind=4), &
    5721                         INT(surf(iz, isurfs), kind=4), &
    5722                         INT(surf(iy, isurfs), kind=4), &
    5723                         INT(surf(ix, isurfs), kind=4), &
    5724                         REAL(rirrf, kind=8), REAL(transparency, kind=8)
    5725                     imrtf = imrtf + 1
    5726 
    5727                 ENDDO !< isurfs
    5728             ENDDO !< MRT_TARGETS record
    5729 
    5730 524         message_string = 'error reading file MRT_TARGETS'
    5731             CALL message( 'radiation_calc_svf', 'PA0524', 1, 2, 0, 6, 0 )
    5732 
    5733 525         message_string = 'error writing file MRT_FACTORS'//myid_char
    5734             CALL message( 'radiation_calc_svf', 'PA0525', 1, 2, 0, 6, 0 )
    5735 
    5736 526         CLOSE(153)
    5737             CLOSE(154)
    5738         ENDIF  !< mrt_factors
    5739        
     6134!--     prepare the MPI_Win for collecting the surface indices
     6135!--     from the reverse index arrays gridsurf from processors of target surfaces
     6136#if defined( __parallel )
     6137        IF ( rad_angular_discretization )  THEN
     6138!
     6139!--         raytrace_mpi_rma is asserted
     6140            CALL MPI_Win_lock_all(0, win_gridsurf, ierr)
     6141            IF ( ierr /= 0 ) THEN
     6142                WRITE(9,*) 'Error MPI_Win_lock_all2:', ierr, win_gridsurf
     6143                FLUSH(9)
     6144            ENDIF
     6145        ENDIF
     6146#endif
     6147
     6148
    57406149        !--Directions opposite to face normals are not even calculated,
    57416150        !--they must be preset to 0
     
    57986207            END SELECT
    57996208
    5800             ALLOCATE ( zdirs(1:nzn), zbdry(0:nzn), vffrac(1:nzn), ztransp(1:nzn) )
     6209            ALLOCATE ( zdirs(1:nzn), zbdry(0:nzn), vffrac(1:nzn*naz), &
     6210                       ztransp(1:nzn*naz), itarget(1:nzn*naz) )   !FIXME allocate itarget only
     6211                                                                  !in case of rad_angular_discretization
     6212
     6213            itarg0 = 1
     6214            itarg1 = nzn
    58016215            zdirs(:) = (/( TAN(pi/2 - (zn0+(REAL(izn,wp)-.5_wp)*zns)), izn=1, nzn )/)
    58026216            zbdry(:) = (/( zn0+REAL(izn,wp)*zns, izn=0, nzn )/)
    58036217            IF ( td == iup_u  .OR.  td == iup_l )  THEN
    5804                !-- For horizontal target, vf fractions are constant per azimuth
    5805                vffrac(:) = (COS(2 * zbdry(0:nzn-1)) - COS(2 * zbdry(1:nzn))) / 2._wp / REAL(naz, wp)
    5806                !--sum of vffrac for all iaz equals 1, verified
     6218               vffrac(1:nzn) = (COS(2 * zbdry(0:nzn-1)) - COS(2 * zbdry(1:nzn))) / 2._wp / REAL(naz, wp)
     6219!
     6220!--            For horizontal target, vf fractions are constant per azimuth
     6221               DO iaz = 1, naz-1
     6222                  vffrac(iaz*nzn+1:(iaz+1)*nzn) = vffrac(1:nzn)
     6223               ENDDO
     6224!--            sum of whole vffrac equals 1, verified
    58076225            ENDIF
    5808 
    5809             !--Calculate sky-view factor and direct solar visibility using 2D raytracing
     6226!
     6227!--         Calculate sky-view factor and direct solar visibility using 2D raytracing
    58106228            DO iaz = 1, naz
    58116229               azmid = az0 + (REAL(iaz, wp) - .5_wp) * azs
     
    58146232                  az1 = az2 - azs
    58156233                  !TODO precalculate after 1st line
    5816                   vffrac(:) = (SIN(az2) - SIN(az1))                           &
     6234                  vffrac(itarg0:itarg1) = (SIN(az2) - SIN(az1))               &
    58176235                              * (zbdry(1:nzn) - zbdry(0:nzn-1)                &
    58186236                                 + SIN(zbdry(0:nzn-1))*COS(zbdry(0:nzn-1))    &
    58196237                                 - SIN(zbdry(1:nzn))*COS(zbdry(1:nzn)))       &
    58206238                              / (2._wp * pi)
    5821                   !--sum of vffrac for all iaz equals 1, verified
     6239!--               sum of whole vffrac equals 1, verified
    58226240               ENDIF
    58236241               yxdir = (/ COS(azmid), SIN(azmid) /)
    5824                CALL raytrace_2d(ta, yxdir, zdirs,                              &
    5825                                surfstart(myid) + isurflt, facearea(td),        &
    5826                                vffrac, .TRUE., .FALSE., win_lad, horizon,       &
    5827                                ztransp)
    5828                
    5829 
    5830                azen = pi/2 - ATAN(horizon)
    5831                IF ( td == iup_u  .OR.  td == iup_l )  THEN
    5832                   azen = MIN(azen, pi/2) !only above horizontal direction
    5833                   skyvf(isurflt) = skyvf(isurflt) + (1._wp - COS(2*azen)) /   &
    5834                      (2._wp * raytrace_discrete_azims)
    5835                ELSE
    5836                   skyvf(isurflt) = skyvf(isurflt) + (SIN(az2) - SIN(az1)) *   &
    5837                               (azen - SIN(azen)*COS(azen)) / (2._wp*pi)
    5838                ENDIF
    5839                skyvft(isurflt) = skyvft(isurflt) + SUM(ztransp(:) * vffrac(:))
     6242               CALL raytrace_2d(ta, yxdir, nzn, zdirs,                        &
     6243                                    surfstart(myid) + isurflt, facearea(td),  &
     6244                                    vffrac(itarg0:itarg1), .TRUE., .TRUE.,    &
     6245                                    .FALSE., lowest_free_ray,                 &
     6246                                    ztransp(itarg0:itarg1),                   &
     6247                                    itarget(itarg0:itarg1))   !FIXME unit vect in grid units + zdirs
     6248                                                              !FIXME itarget available only in
     6249                                                              !case of rad_angular_discretization
     6250               skyvf(isurflt) = skyvf(isurflt) + &
     6251                                SUM(vffrac(itarg0:itarg0+lowest_free_ray-1))
     6252               skyvft(isurflt) = skyvft(isurflt) + &
     6253                                 SUM(ztransp(itarg0:itarg0+lowest_free_ray-1) &
     6254                                             * vffrac(itarg0:itarg0+lowest_free_ray-1))
    58406255 
    5841                !--Save direct solar transparency
     6256!--            Save direct solar transparency
    58426257               j = MODULO(NINT(azmid/                                          &
    58436258                               (2._wp*pi)*raytrace_discrete_azims-.5_wp, iwp), &
     
    58466261               DO k = 1, raytrace_discrete_elevs/2
    58476262                  i = dsidir_rev(k-1, j)
    5848                   IF ( i /= -1 )  dsitrans(isurflt, i) = ztransp(k)
     6263                  IF ( i /= -1  .AND.  k <= lowest_free_ray )  &
     6264                     dsitrans(isurflt, i) = ztransp(itarg0+k-1)
    58496265               ENDDO
     6266
     6267               !
     6268               !--Advance itarget indices
     6269               itarg0 = itarg1 + 1
     6270               itarg1 = itarg1 + nzn
    58506271            ENDDO
    58516272
    5852             DEALLOCATE ( zdirs, zbdry, vffrac, ztransp )
    5853 !
    5854 !--         Following calculations only required for surface_reflections
    5855             IF ( surface_reflections )  THEN
    5856 
    5857                DO  isurfs = 1, nsurf
    5858                   IF ( .NOT.  surface_facing(surfl(ix, isurflt), surfl(iy, isurflt), &
    5859                      surfl(iz, isurflt), surfl(id, isurflt), &
    5860                      surf(ix, isurfs), surf(iy, isurfs), &
    5861                      surf(iz, isurfs), surf(id, isurfs)) )  THEN
    5862                      CYCLE
     6273            IF ( rad_angular_discretization )  THEN
     6274!--            sort itarget by face id
     6275               CALL quicksort_itarget(itarget,vffrac,ztransp,1,nzn*naz)
     6276!
     6277!--            find the first valid position
     6278               itarg0 = 1
     6279               DO WHILE ( itarg0 <= nzn*naz )
     6280                  IF ( itarget(itarg0) /= -1 )  EXIT
     6281                  itarg0 = itarg0 + 1
     6282               ENDDO
     6283
     6284               DO  i = itarg0, nzn*naz
     6285!
     6286!--               For duplicate values, only sum up vf fraction value
     6287                  IF ( i < nzn*naz )  THEN
     6288                     IF ( itarget(i+1) == itarget(i) )  THEN
     6289                        vffrac(i+1) = vffrac(i+1) + vffrac(i)
     6290                        CYCLE
     6291                     ENDIF
    58636292                  ENDIF
    5864                  
    5865                   sd = surf(id, isurfs)
    5866                   sa = (/ REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd),  &
    5867                           REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd),  &
    5868                           REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd)  /)
    5869 
    5870 !--               unit vector source -> target
    5871                   uv = (/ (ta(1)-sa(1))*dz(1), (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)
    5872                   sqdist = SUM(uv(:)**2)
    5873                   uv = uv / SQRT(sqdist)
    5874 
    5875 !--               reject raytracing above max distance
    5876                   IF ( SQRT(sqdist) > max_raytracing_dist ) THEN
    5877                      ray_skip_maxdist = ray_skip_maxdist + 1
    5878                      CYCLE
    5879                   ENDIF
    5880                  
    5881 !--               irradiance factor (our unshaded shape view factor) = view factor per differential target area * source area
    5882                   rirrf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & ! cosine of source normal and direction
    5883                       * dot_product((/ kdir(td), jdir(td), idir(td) /), -uv) &  ! cosine of target normal and reverse direction
    5884                       / (pi * sqdist) & ! square of distance between centers
    5885                       * facearea(sd)
    5886 
    5887 !--               reject raytracing for potentially too small view factor values
    5888                   IF ( rirrf < min_irrf_value ) THEN
    5889                       ray_skip_minval = ray_skip_minval + 1
    5890                       CYCLE
    5891                   ENDIF
    5892 
    5893 !--               raytrace + process plant canopy sinks within
    5894                   CALL raytrace(sa, ta, isurfs, rirrf, facearea(td), .TRUE., &
    5895                                 visible, transparency, win_lad)
    5896 
    5897                   IF ( .NOT.  visible ) CYCLE
    5898                  ! rsvf = rirrf * transparency
    5899 
     6293!
    59006294!--               write to the svf array
    59016295                  nsvfl = nsvfl + 1
    59026296!--               check dimmension of asvf array and enlarge it if needed
    59036297                  IF ( nsvfla < nsvfl )  THEN
    5904                      k = nsvfla * 2
     6298                     k = CEILING(REAL(nsvfla, kind=wp) * grow_factor)
    59056299                     IF ( msvf == 0 )  THEN
    59066300                        msvf = 1
     
    59176311                     ENDIF
    59186312
    5919 !                      WRITE(msg,'(A,3I12)') 'Grow asvf:',nsvfl,nsvfla,k
    5920 !                      CALL radiation_write_debug_log( msg )
     6313                     WRITE (msg,'(A,3I12)') 'Grow asvf:',nsvfl,nsvfla,k
     6314                     CALL radiation_write_debug_log( msg )
     6315                     
     6316                     nsvfla = k
     6317                  ENDIF
     6318!--               write svf values into the array
     6319                  asvf(nsvfl)%isurflt = isurflt
     6320                  asvf(nsvfl)%isurfs = itarget(i)
     6321                  asvf(nsvfl)%rsvf = vffrac(i)
     6322                  asvf(nsvfl)%rtransp = ztransp(i)
     6323               END DO
     6324
     6325            ENDIF ! rad_angular_discretization
     6326
     6327            DEALLOCATE ( zdirs, zbdry, vffrac, ztransp, itarget ) !FIXME itarget shall be allocated only
     6328                                                                  !in case of rad_angular_discretization
     6329!
     6330!--         Following calculations only required for surface_reflections
     6331            IF ( surface_reflections  .AND.  .NOT. rad_angular_discretization )  THEN
     6332
     6333               DO  isurfs = 1, nsurf
     6334                  IF ( .NOT.  surface_facing(surfl(ix, isurflt), surfl(iy, isurflt), &
     6335                     surfl(iz, isurflt), surfl(id, isurflt), &
     6336                     surf(ix, isurfs), surf(iy, isurfs), &
     6337                     surf(iz, isurfs), surf(id, isurfs)) )  THEN
     6338                     CYCLE
     6339                  ENDIF
     6340                 
     6341                  sd = surf(id, isurfs)
     6342                  sa = (/ REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd),  &
     6343                          REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd),  &
     6344                          REAL(surf(ix, isurfs), wp) - 0.5_wp * idir(sd)  /)
     6345
     6346!--               unit vector source -> target
     6347                  uv = (/ (ta(1)-sa(1))*dz(1), (ta(2)-sa(2))*dy, (ta(3)-sa(3))*dx /)
     6348                  sqdist = SUM(uv(:)**2)
     6349                  uv = uv / SQRT(sqdist)
     6350
     6351!--               reject raytracing above max distance
     6352                  IF ( SQRT(sqdist) > max_raytracing_dist ) THEN
     6353                     ray_skip_maxdist = ray_skip_maxdist + 1
     6354                     CYCLE
     6355                  ENDIF
     6356                 
     6357!--               irradiance factor (our unshaded shape view factor) = view factor per differential target area * source area
     6358                  rirrf = dot_product((/ kdir(sd), jdir(sd), idir(sd) /), uv) & ! cosine of source normal and direction
     6359                      * dot_product((/ kdir(td), jdir(td), idir(td) /), -uv) &  ! cosine of target normal and reverse direction
     6360                      / (pi * sqdist) & ! square of distance between centers
     6361                      * facearea(sd)
     6362
     6363!--               reject raytracing for potentially too small view factor values
     6364                  IF ( rirrf < min_irrf_value ) THEN
     6365                      ray_skip_minval = ray_skip_minval + 1
     6366                      CYCLE
     6367                  ENDIF
     6368
     6369!--               raytrace + process plant canopy sinks within
     6370                  CALL raytrace(sa, ta, isurfs, rirrf, facearea(td), .TRUE., &
     6371                                visible, transparency)
     6372
     6373                  IF ( .NOT.  visible ) CYCLE
     6374                 ! rsvf = rirrf * transparency
     6375
     6376!--               write to the svf array
     6377                  nsvfl = nsvfl + 1
     6378!--               check dimmension of asvf array and enlarge it if needed
     6379                  IF ( nsvfla < nsvfl )  THEN
     6380                     k = CEILING(REAL(nsvfla, kind=wp) * grow_factor)
     6381                     IF ( msvf == 0 )  THEN
     6382                        msvf = 1
     6383                        ALLOCATE( asvf1(k) )
     6384                        asvf => asvf1
     6385                        asvf1(1:nsvfla) = asvf2
     6386                        DEALLOCATE( asvf2 )
     6387                     ELSE
     6388                        msvf = 0
     6389                        ALLOCATE( asvf2(k) )
     6390                        asvf => asvf2
     6391                        asvf2(1:nsvfla) = asvf1
     6392                        DEALLOCATE( asvf1 )
     6393                     ENDIF
     6394
     6395                     WRITE(msg,'(A,3I12)') 'Grow asvf:',nsvfl,nsvfla,k
     6396                     CALL radiation_write_debug_log( msg )
    59216397                     
    59226398                     nsvfla = k
     
    59316407        ENDDO
    59326408
    5933         !--Raytrace to canopy boxes to fill dsitransc TODO optimize
    5934         !--
    5935         dsitransc(:,:) = -999._wp !FIXME
     6409!--
     6410!--     Raytrace to canopy boxes to fill dsitransc TODO optimize
     6411        dsitransc(:,:) = 0._wp
    59366412        az0 = 0._wp
    59376413        naz = raytrace_discrete_azims
     
    59406416        nzn = raytrace_discrete_elevs / 2
    59416417        zns = pi / 2._wp / REAL(nzn, wp)
    5942         ALLOCATE ( zdirs(1:nzn), vffrac(1:nzn), ztransp(1:nzn) )
     6418        ALLOCATE ( zdirs(1:nzn), vffrac(1:nzn), ztransp(1:nzn), &
     6419               itarget(1:nzn) )
    59436420        zdirs(:) = (/( TAN(pi/2 - (zn0+(REAL(izn,wp)-.5_wp)*zns)), izn=1, nzn )/)
    59446421        vffrac(:) = 0._wp
    59456422
    5946         DO ipcgb = 1, npcbl
     6423        DO  ipcgb = 1, npcbl
    59476424           ta = (/ REAL(pcbl(iz, ipcgb), wp),  &
    59486425                   REAL(pcbl(iy, ipcgb), wp),  &
    59496426                   REAL(pcbl(ix, ipcgb), wp) /)
    5950            !--Calculate sky-view factor and direct solar visibility using 2D raytracing
    5951            DO iaz = 1, naz
     6427!--        Calculate direct solar visibility using 2D raytracing
     6428           DO  iaz = 1, naz
    59526429              azmid = az0 + (REAL(iaz, wp) - .5_wp) * azs
    59536430              yxdir = (/ COS(azmid), SIN(azmid) /)
    5954               CALL raytrace_2d(ta, yxdir, zdirs,     &
    5955                                    -999, -999._wp, vffrac, .FALSE., .TRUE., &
    5956                                    win_lad, horizon, ztransp)
    5957 
    5958               !--Save direct solar transparency
     6431              CALL raytrace_2d(ta, yxdir, nzn, zdirs,                                &
     6432                                   -999, -999._wp, vffrac, .FALSE., .FALSE., .TRUE., &
     6433                                   lowest_free_ray, ztransp, itarget) !FIXME unit vect in grid units + zdirs
     6434
     6435!--           Save direct solar transparency
    59596436              j = MODULO(NINT(azmid/                                         &
    59606437                             (2._wp*pi)*raytrace_discrete_azims-.5_wp, iwp), &
    59616438                         raytrace_discrete_azims)
    5962               DO k = 1, raytrace_discrete_elevs/2
     6439              DO  k = 1, raytrace_discrete_elevs/2
    59636440                 i = dsidir_rev(k-1, j)
    5964                  IF ( i /= -1 )  dsitransc(ipcgb, i) = ztransp(k)
     6441                 IF ( i /= -1  .AND.  k <= lowest_free_ray ) &
     6442                    dsitransc(ipcgb, i) = ztransp(k)
    59656443              ENDDO
    59666444           ENDDO
    59676445        ENDDO
    5968         DEALLOCATE ( zdirs, vffrac, ztransp )
    5969 
    5970 !         CALL radiation_write_debug_log( 'End of calculation SVF' )
    5971 !         WRITE(msg, *) 'Raytracing skipped for maximum distance of ', &
    5972 !             max_raytracing_dist, ' m on ', ray_skip_maxdist, ' pairs.'
    5973 !         CALL radiation_write_debug_log( msg )
    5974 !         WRITE(msg, *) 'Raytracing skipped for minimum potential value of ', &
    5975 !             min_irrf_value , ' on ', ray_skip_minval, ' pairs.'
    5976 !         CALL radiation_write_debug_log( msg )
     6446        DEALLOCATE ( zdirs, vffrac, ztransp, itarget )
     6447!--
     6448!--     Raytrace to MRT boxes
     6449        IF ( nmrtbl > 0 )  THEN
     6450           mrtdsit(:,:) = 0._wp
     6451           mrtsky(:) = 0._wp
     6452           mrtskyt(:) = 0._wp
     6453           az0 = 0._wp
     6454           naz = raytrace_discrete_azims
     6455           azs = 2._wp * pi / REAL(naz, wp)
     6456           zn0 = 0._wp
     6457           nzn = raytrace_discrete_elevs
     6458           zns = pi / REAL(nzn, wp)
     6459           ALLOCATE ( zdirs(1:nzn), zbdry(0:nzn), vffrac(1:nzn*naz), vffrac0(1:nzn), &
     6460                      ztransp(1:nzn*naz), itarget(1:nzn*naz) )   !FIXME allocate itarget only
     6461                                                                 !in case of rad_angular_discretization
     6462
     6463           zdirs(:) = (/( TAN(pi/2 - (zn0+(REAL(izn,wp)-.5_wp)*zns)), izn=1, nzn )/)
     6464           zbdry(:) = (/( zn0+REAL(izn,wp)*zns, izn=0, nzn )/)
     6465           vffrac0(:) = (COS(zbdry(0:nzn-1)) - COS(zbdry(1:nzn))) / 2._wp / REAL(naz, wp)
     6466
     6467           DO  imrt = 1, nmrtbl
     6468              ta = (/ REAL(mrtbl(iz, imrt), wp),  &
     6469                      REAL(mrtbl(iy, imrt), wp),  &
     6470                      REAL(mrtbl(ix, imrt), wp) /)
     6471!
     6472!--           vf fractions are constant per azimuth
     6473              DO iaz = 0, naz-1
     6474                 vffrac(iaz*nzn+1:(iaz+1)*nzn) = vffrac0(:)
     6475              ENDDO
     6476!--           sum of whole vffrac equals 1, verified
     6477              itarg0 = 1
     6478              itarg1 = nzn
     6479!
     6480!--           Calculate sky-view factor and direct solar visibility using 2D raytracing
     6481              DO  iaz = 1, naz
     6482                 azmid = az0 + (REAL(iaz, wp) - .5_wp) * azs
     6483                 CALL raytrace_2d(ta, (/ COS(azmid), SIN(azmid) /), nzn, zdirs,  &
     6484                                  -999, -999._wp, vffrac(itarg0:itarg1), .TRUE., &
     6485                                  .FALSE., .TRUE., lowest_free_ray,              &
     6486                                  ztransp(itarg0:itarg1),                        &
     6487                                  itarget(itarg0:itarg1))   !FIXME unit vect in grid units + zdirs
     6488                                                            !FIXME itarget available only in
     6489                                                            !case of rad_angular_discretization
     6490
     6491!--              Sky view factors for MRT
     6492                 mrtsky(imrt) = mrtsky(imrt) + &
     6493                                  SUM(vffrac(itarg0:itarg0+lowest_free_ray-1))
     6494                 mrtskyt(imrt) = mrtskyt(imrt) + &
     6495                                   SUM(ztransp(itarg0:itarg0+lowest_free_ray-1) &
     6496                                               * vffrac(itarg0:itarg0+lowest_free_ray-1))
     6497!--              Direct solar transparency for MRT
     6498                 j = MODULO(NINT(azmid/                                         &
     6499                                (2._wp*pi)*raytrace_discrete_azims-.5_wp, iwp), &
     6500                            raytrace_discrete_azims)
     6501                 DO  k = 1, raytrace_discrete_elevs/2
     6502                    i = dsidir_rev(k-1, j)
     6503                    IF ( i /= -1  .AND.  k <= lowest_free_ray ) &
     6504                       mrtdsit(imrt, i) = ztransp(itarg0+k-1)
     6505                 ENDDO
     6506!
     6507!--              Advance itarget indices
     6508                 itarg0 = itarg1 + 1
     6509                 itarg1 = itarg1 + nzn
     6510              ENDDO
     6511
     6512!--           sort itarget by face id
     6513              CALL quicksort_itarget(itarget,vffrac,ztransp,1,nzn*naz)
     6514!
     6515!--           find the first valid position
     6516              itarg0 = 1
     6517              DO WHILE ( itarg0 <= nzn*naz )
     6518                 IF ( itarget(itarg0) /= -1 )  EXIT
     6519                 itarg0 = itarg0 + 1
     6520              ENDDO
     6521
     6522              DO  i = itarg0, nzn*naz
     6523!
     6524!--              For duplicate values, only sum up vf fraction value
     6525                 IF ( i < nzn*naz )  THEN
     6526                    IF ( itarget(i+1) == itarget(i) )  THEN
     6527                       vffrac(i+1) = vffrac(i+1) + vffrac(i)
     6528                       CYCLE
     6529                    ENDIF
     6530                 ENDIF
     6531!
     6532!--              write to the mrtf array
     6533                 nmrtf = nmrtf + 1
     6534!--              check dimmension of mrtf array and enlarge it if needed
     6535                 IF ( nmrtfa < nmrtf )  THEN
     6536                    k = CEILING(REAL(nmrtfa, kind=wp) * grow_factor)
     6537                    IF ( mmrtf == 0 )  THEN
     6538                       mmrtf = 1
     6539                       ALLOCATE( amrtf1(k) )
     6540                       amrtf => amrtf1
     6541                       amrtf1(1:nmrtfa) = amrtf2
     6542                       DEALLOCATE( amrtf2 )
     6543                    ELSE
     6544                       mmrtf = 0
     6545                       ALLOCATE( amrtf2(k) )
     6546                       amrtf => amrtf2
     6547                       amrtf2(1:nmrtfa) = amrtf1
     6548                       DEALLOCATE( amrtf1 )
     6549                    ENDIF
     6550
     6551                    WRITE (msg,'(A,3I12)') 'Grow amrtf:', nmrtf, nmrtfa, k
     6552                    CALL radiation_write_debug_log( msg )
     6553
     6554                    nmrtfa = k
     6555                 ENDIF
     6556!--              write mrtf values into the array
     6557                 amrtf(nmrtf)%isurflt = imrt
     6558                 amrtf(nmrtf)%isurfs = itarget(i)
     6559                 amrtf(nmrtf)%rsvf = vffrac(i)
     6560                 amrtf(nmrtf)%rtransp = ztransp(i)
     6561              ENDDO ! itarg
     6562
     6563           ENDDO ! imrt
     6564           DEALLOCATE ( zdirs, zbdry, vffrac, vffrac0, ztransp, itarget )
     6565!
     6566!--        Move MRT factors to final arrays
     6567           ALLOCATE ( mrtf(nmrtf), mrtft(nmrtf), mrtfsurf(2,nmrtf) )
     6568           DO  imrtf = 1, nmrtf
     6569              mrtf(imrtf) = amrtf(imrtf)%rsvf
     6570              mrtft(imrtf) = amrtf(imrtf)%rsvf * amrtf(imrtf)%rtransp
     6571              mrtfsurf(:,imrtf) = (/amrtf(imrtf)%isurflt, amrtf(imrtf)%isurfs /)
     6572           ENDDO
     6573           IF ( ALLOCATED(amrtf1) )  DEALLOCATE( amrtf1 )
     6574           IF ( ALLOCATED(amrtf2) )  DEALLOCATE( amrtf2 )
     6575        ENDIF ! nmrtbl > 0
     6576
     6577        IF ( rad_angular_discretization )  THEN
     6578#if defined( __parallel )
     6579!--        finalize MPI_RMA communication established to get global index of the surface from grid indices
     6580!--        flush all MPI window pending requests
     6581           CALL MPI_Win_flush_all(win_gridsurf, ierr)
     6582           IF ( ierr /= 0 ) THEN
     6583               WRITE(9,*) 'Error MPI_Win_flush_all1:', ierr, win_gridsurf
     6584               FLUSH(9)
     6585           ENDIF
     6586!--        unlock MPI window
     6587           CALL MPI_Win_unlock_all(win_gridsurf, ierr)
     6588           IF ( ierr /= 0 ) THEN
     6589               WRITE(9,*) 'Error MPI_Win_unlock_all1:', ierr, win_gridsurf
     6590               FLUSH(9)
     6591           ENDIF
     6592!--        free MPI window
     6593           CALL MPI_Win_free(win_gridsurf, ierr)
     6594           IF ( ierr /= 0 ) THEN
     6595               WRITE(9,*) 'Error MPI_Win_free1:', ierr, win_gridsurf
     6596               FLUSH(9)
     6597           ENDIF
     6598#else
     6599           DEALLOCATE ( gridsurf )
     6600#endif
     6601        ENDIF
     6602
     6603        CALL radiation_write_debug_log( 'End of calculation SVF' )
     6604        WRITE(msg, *) 'Raytracing skipped for maximum distance of ', &
     6605           max_raytracing_dist, ' m on ', ray_skip_maxdist, ' pairs.'
     6606        CALL radiation_write_debug_log( msg )
     6607        WRITE(msg, *) 'Raytracing skipped for minimum potential value of ', &
     6608           min_irrf_value , ' on ', ray_skip_minval, ' pairs.'
     6609        CALL radiation_write_debug_log( msg )
    59776610
    59786611        CALL location_message( '    waiting for completion of SVF and CSF calculation in all processes', .TRUE. )
     
    59836616!--         finalize mpi_rma communication and deallocate temporary arrays
    59846617#if defined( __parallel )
    5985             IF ( rma_lad_raytrace )  THEN
     6618            IF ( raytrace_mpi_rma )  THEN
    59866619                CALL MPI_Win_flush_all(win_lad, ierr)
     6620                IF ( ierr /= 0 ) THEN
     6621                    WRITE(9,*) 'Error MPI_Win_flush_all2:', ierr, win_lad
     6622                    FLUSH(9)
     6623                ENDIF
    59876624!--             unlock MPI window
    59886625                CALL MPI_Win_unlock_all(win_lad, ierr)
     6626                IF ( ierr /= 0 ) THEN
     6627                    WRITE(9,*) 'Error MPI_Win_unlock_all2:', ierr, win_lad
     6628                    FLUSH(9)
     6629                ENDIF
    59896630!--             free MPI window
    59906631                CALL MPI_Win_free(win_lad, ierr)
    5991                
     6632                IF ( ierr /= 0 ) THEN
     6633                    WRITE(9,*) 'Error MPI_Win_free2:', ierr, win_lad
     6634                    FLUSH(9)
     6635                ENDIF
    59926636!--             deallocate temporary arrays storing values for csf calculation during raytracing
    59936637                DEALLOCATE( lad_s_ray )
    5994 !--             sub_lad is the pointer to lad_s_rma in case of rma_lad_raytrace
     6638!--             sub_lad is the pointer to lad_s_rma in case of raytrace_mpi_rma
    59956639!--             and must not be deallocated here
    59966640            ELSE
     
    60096653        CALL location_message( '    calculation of the complete SVF array', .TRUE. )
    60106654
    6011 !         CALL radiation_write_debug_log( 'Start SVF sort' )
    6012 !--     sort svf ( a version of quicksort )
    6013         CALL quicksort_svf(asvf,1,nsvfl)
    6014 
    6015         !< load svf from the structure array to plain arrays
    6016 !         CALL radiation_write_debug_log( 'Load svf from the structure array to plain arrays' )
    6017         ALLOCATE( svf(ndsvf,nsvfl) )
    6018         ALLOCATE( svfsurf(idsvf,nsvfl) )
    6019         svfnorm_counts(:) = 0._wp
    6020         isurflt_prev = -1
    6021         ksvf = 1
    6022         svfsum = 0._wp
    6023         DO isvf = 1, nsvfl
    6024 !--         normalize svf per target face
    6025             IF ( asvf(ksvf)%isurflt /= isurflt_prev )  THEN
    6026                 IF ( isurflt_prev /= -1  .AND.  svfsum /= 0._wp )  THEN
    6027                     !< update histogram of logged svf normalization values
    6028                     i = searchsorted(svfnorm_report_thresh, svfsum / (1._wp-skyvf(isurflt_prev)))
    6029                     svfnorm_counts(i) = svfnorm_counts(i) + 1
    6030 
    6031                     svf(1, isvf_surflt:isvf-1) = svf(1, isvf_surflt:isvf-1) / svfsum * (1._wp-skyvf(isurflt_prev))
    6032                 ENDIF
    6033                 isurflt_prev = asvf(ksvf)%isurflt
    6034                 isvf_surflt = isvf
    6035                 svfsum = asvf(ksvf)%rsvf !?? / asvf(ksvf)%rtransp
    6036             ELSE
    6037                 svfsum = svfsum + asvf(ksvf)%rsvf !?? / asvf(ksvf)%rtransp
    6038             ENDIF
    6039 
    6040             svf(:, isvf) = (/ asvf(ksvf)%rsvf, asvf(ksvf)%rtransp /)
    6041             svfsurf(:, isvf) = (/ asvf(ksvf)%isurflt, asvf(ksvf)%isurfs /)
    6042 
    6043 !--         next element
    6044             ksvf = ksvf + 1
    6045         ENDDO
    6046 
    6047         IF ( isurflt_prev /= -1  .AND.  svfsum /= 0._wp )  THEN
    6048             i = searchsorted(svfnorm_report_thresh, svfsum / (1._wp-skyvf(isurflt_prev)))
    6049             svfnorm_counts(i) = svfnorm_counts(i) + 1
    6050 
    6051             svf(1, isvf_surflt:nsvfl) = svf(1, isvf_surflt:nsvfl) / svfsum * (1._wp-skyvf(isurflt_prev))
    6052         ENDIF
    6053         !TODO we should be able to deallocate skyvf, from now on we only need skyvft
     6655        IF ( rad_angular_discretization )  THEN
     6656           CALL radiation_write_debug_log( 'Load svf from the structure array to plain arrays' )
     6657           ALLOCATE( svf(ndsvf,nsvfl) )
     6658           ALLOCATE( svfsurf(idsvf,nsvfl) )
     6659
     6660           DO isvf = 1, nsvfl
     6661               svf(:, isvf) = (/ asvf(isvf)%rsvf, asvf(isvf)%rtransp /)
     6662               svfsurf(:, isvf) = (/ asvf(isvf)%isurflt, asvf(isvf)%isurfs /)
     6663           ENDDO
     6664        ELSE
     6665           CALL radiation_write_debug_log( 'Start SVF sort' )
     6666!--        sort svf ( a version of quicksort )
     6667           CALL quicksort_svf(asvf,1,nsvfl)
     6668
     6669           !< load svf from the structure array to plain arrays
     6670           CALL radiation_write_debug_log( 'Load svf from the structure array to plain arrays' )
     6671           ALLOCATE( svf(ndsvf,nsvfl) )
     6672           ALLOCATE( svfsurf(idsvf,nsvfl) )
     6673           svfnorm_counts(:) = 0._wp
     6674           isurflt_prev = -1
     6675           ksvf = 1
     6676           svfsum = 0._wp
     6677           DO isvf = 1, nsvfl
     6678!--            normalize svf per target face
     6679               IF ( asvf(ksvf)%isurflt /= isurflt_prev )  THEN
     6680                   IF ( isurflt_prev /= -1  .AND.  svfsum /= 0._wp )  THEN
     6681                       !< update histogram of logged svf normalization values
     6682                       i = searchsorted(svfnorm_report_thresh, svfsum / (1._wp-skyvf(isurflt_prev)))
     6683                       svfnorm_counts(i) = svfnorm_counts(i) + 1
     6684
     6685                       svf(1, isvf_surflt:isvf-1) = svf(1, isvf_surflt:isvf-1) / svfsum * (1._wp-skyvf(isurflt_prev))
     6686                   ENDIF
     6687                   isurflt_prev = asvf(ksvf)%isurflt
     6688                   isvf_surflt = isvf
     6689                   svfsum = asvf(ksvf)%rsvf !?? / asvf(ksvf)%rtransp
     6690               ELSE
     6691                   svfsum = svfsum + asvf(ksvf)%rsvf !?? / asvf(ksvf)%rtransp
     6692               ENDIF
     6693
     6694               svf(:, isvf) = (/ asvf(ksvf)%rsvf, asvf(ksvf)%rtransp /)
     6695               svfsurf(:, isvf) = (/ asvf(ksvf)%isurflt, asvf(ksvf)%isurfs /)
     6696
     6697!--            next element
     6698               ksvf = ksvf + 1
     6699           ENDDO
     6700
     6701           IF ( isurflt_prev /= -1  .AND.  svfsum /= 0._wp )  THEN
     6702               i = searchsorted(svfnorm_report_thresh, svfsum / (1._wp-skyvf(isurflt_prev)))
     6703               svfnorm_counts(i) = svfnorm_counts(i) + 1
     6704
     6705               svf(1, isvf_surflt:nsvfl) = svf(1, isvf_surflt:nsvfl) / svfsum * (1._wp-skyvf(isurflt_prev))
     6706           ENDIF
     6707           WRITE(9, *) 'SVF normalization histogram:', svfnorm_counts,  &
     6708                       'on thresholds:', svfnorm_report_thresh(1:svfnorm_report_num), '(val < thresh <= val)'
     6709           !TODO we should be able to deallocate skyvf, from now on we only need skyvft
     6710        ENDIF ! rad_angular_discretization
    60546711
    60556712!--     deallocate temporary asvf array
     
    60676724
    60686725            CALL location_message( '    calculation of the complete CSF array', .TRUE. )
    6069 !             CALL radiation_write_debug_log( 'Calculation of the complete CSF array' )
     6726            CALL radiation_write_debug_log( 'Calculation of the complete CSF array' )
    60706727!--         sort and merge csf for the last time, keeping the array size to minimum
    60716728            CALL merge_and_grow_csf(-1)
     
    60736730!--         aggregate csb among processors
    60746731!--         allocate necessary arrays
    6075             ALLOCATE( csflt(ndcsf,max(ncsfl,ndcsf)) )
    6076             ALLOCATE( kcsflt(kdcsf,max(ncsfl,kdcsf)) )
     6732            udim = max(ncsfl,1)
     6733            ALLOCATE( csflt_l(ndcsf*udim) )
     6734            csflt(1:ndcsf,1:udim) => csflt_l(1:ndcsf*udim)
     6735            ALLOCATE( kcsflt_l(kdcsf*udim) )
     6736            kcsflt(1:kdcsf,1:udim) => kcsflt_l(1:kdcsf*udim)
    60776737            ALLOCATE( icsflt(0:numprocs-1) )
    60786738            ALLOCATE( dcsflt(0:numprocs-1) )
     
    61086768!--             fill out real values of rsvf, rtransp
    61096769                csflt(1,kcsf) = acsf(kcsf)%rsvf
    6110                 csflt(2,kcsf) = acsf(kcsf)%rtransp
    61116770!--             fill out integer values of itz,ity,itx,isurfs
    61126771                kcsflt(1,kcsf) = acsf(kcsf)%itz
     
    61386797!--         scatter and gather the number of elements to and from all processor
    61396798!--         and calculate displacements
    6140 !             CALL radiation_write_debug_log( 'Scatter and gather the number of elements to and from all processor' )
     6799            CALL radiation_write_debug_log( 'Scatter and gather the number of elements to and from all processor' )
    61416800            CALL MPI_AlltoAll(icsflt,1,MPI_INTEGER,ipcsflt,1,MPI_INTEGER,comm2d, ierr)
    6142            
     6801            IF ( ierr /= 0 ) THEN
     6802                WRITE(9,*) 'Error MPI_AlltoAll1:', ierr, SIZE(icsflt), SIZE(ipcsflt)
     6803                FLUSH(9)
     6804            ENDIF
     6805
    61436806            npcsfl = SUM(ipcsflt)
    61446807            d = 0
     
    61496812
    61506813!--         exchange csf fields between processors
    6151 !             CALL radiation_write_debug_log( 'Exchange csf fields between processors' )
    6152             ALLOCATE( pcsflt(ndcsf,max(npcsfl,ndcsf)) )
    6153             ALLOCATE( kpcsflt(kdcsf,max(npcsfl,kdcsf)) )
    6154             CALL MPI_AlltoAllv(csflt, ndcsf*icsflt, ndcsf*dcsflt, MPI_REAL, &
    6155                 pcsflt, ndcsf*ipcsflt, ndcsf*dpcsflt, MPI_REAL, comm2d, ierr)
    6156             CALL MPI_AlltoAllv(kcsflt, kdcsf*icsflt, kdcsf*dcsflt, MPI_INTEGER, &
    6157                 kpcsflt, kdcsf*ipcsflt, kdcsf*dpcsflt, MPI_INTEGER, comm2d, ierr)
     6814            CALL radiation_write_debug_log( 'Exchange csf fields between processors' )
     6815            udim = max(npcsfl,1)
     6816            ALLOCATE( pcsflt_l(ndcsf*udim) )
     6817            pcsflt(1:ndcsf,1:udim) => pcsflt_l(1:ndcsf*udim)
     6818            ALLOCATE( kpcsflt_l(kdcsf*udim) )
     6819            kpcsflt(1:kdcsf,1:udim) => kpcsflt_l(1:kdcsf*udim)
     6820            CALL MPI_AlltoAllv(csflt_l, ndcsf*icsflt, ndcsf*dcsflt, MPI_REAL, &
     6821                pcsflt_l, ndcsf*ipcsflt, ndcsf*dpcsflt, MPI_REAL, comm2d, ierr)
     6822            IF ( ierr /= 0 ) THEN
     6823                WRITE(9,*) 'Error MPI_AlltoAllv1:', ierr, SIZE(ipcsflt), ndcsf*icsflt, &
     6824                            ndcsf*dcsflt, SIZE(pcsflt_l),ndcsf*ipcsflt, ndcsf*dpcsflt
     6825                FLUSH(9)
     6826            ENDIF
     6827
     6828            CALL MPI_AlltoAllv(kcsflt_l, kdcsf*icsflt, kdcsf*dcsflt, MPI_INTEGER, &
     6829                kpcsflt_l, kdcsf*ipcsflt, kdcsf*dpcsflt, MPI_INTEGER, comm2d, ierr)
     6830            IF ( ierr /= 0 ) THEN
     6831                WRITE(9,*) 'Error MPI_AlltoAllv2:', ierr, SIZE(kcsflt_l),kdcsf*icsflt, &
     6832                           kdcsf*dcsflt, SIZE(kpcsflt_l), kdcsf*ipcsflt, kdcsf*dpcsflt
     6833                FLUSH(9)
     6834            ENDIF
    61586835           
    61596836#else
     
    61666843
    61676844!--         deallocate temporary arrays
    6168             DEALLOCATE( csflt )
    6169             DEALLOCATE( kcsflt )
     6845            DEALLOCATE( csflt_l )
     6846            DEALLOCATE( kcsflt_l )
    61706847            DEALLOCATE( icsflt )
    61716848            DEALLOCATE( dcsflt )
     
    61746851
    61756852!--         sort csf ( a version of quicksort )
    6176 !             CALL radiation_write_debug_log( 'Sort csf' )
     6853            CALL radiation_write_debug_log( 'Sort csf' )
    61776854            CALL quicksort_csf2(kpcsflt, pcsflt, 1, npcsfl)
    61786855
    61796856!--         aggregate canopy sink factor records with identical box & source
    61806857!--         againg across all values from all processors
    6181 !             CALL radiation_write_debug_log( 'Aggregate canopy sink factor records with identical box' )
     6858            CALL radiation_write_debug_log( 'Aggregate canopy sink factor records with identical box' )
    61826859
    61836860            IF ( npcsfl > 0 )  THEN
     
    61906867                         kpcsflt(1,icsf) == kpcsflt(1,icsf+1)  .AND.  &
    61916868                         kpcsflt(4,icsf) == kpcsflt(4,icsf+1) )  THEN
    6192 !--                     We could simply take either first or second rtransp, both are valid. As a very simple heuristic about which ray
    6193 !--                     probably passes nearer the center of the target box, we choose DIF from the entry with greater CSF, since that
    6194 !--                     might mean that the traced beam passes longer through the canopy box.
    6195                         IF ( pcsflt(1,kcsf) < pcsflt(1,icsf+1) )  THEN
    6196                             pcsflt(2,kcsf) = pcsflt(2,icsf+1)
    6197                         ENDIF
     6869
    61986870                        pcsflt(1,kcsf) = pcsflt(1,kcsf) + pcsflt(1,icsf+1)
    61996871
     
    62246896           
    62256897!--         deallocation of temporary arrays
    6226             DEALLOCATE( pcsflt )
    6227             DEALLOCATE( kpcsflt )
    6228 !             CALL radiation_write_debug_log( 'End of aggregate csf' )
     6898            IF ( npcbl > 0 )  DEALLOCATE( gridpcbl )
     6899            DEALLOCATE( pcsflt_l )
     6900            DEALLOCATE( kpcsflt_l )
     6901            CALL radiation_write_debug_log( 'End of aggregate csf' )
    62296902           
    62306903        ENDIF
     
    62336906        CALL MPI_BARRIER( comm2d, ierr )
    62346907#endif
    6235 !         CALL radiation_write_debug_log( 'End of radiation_calc_svf (after mpi_barrier)' )
     6908        CALL radiation_write_debug_log( 'End of radiation_calc_svf (after mpi_barrier)' )
    62366909
    62376910        RETURN
     
    62616934!>    doesn't need to be the same in all three dimensions).
    62626935!------------------------------------------------------------------------------!
    6263     SUBROUTINE raytrace(src, targ, isrc, rirrf, atarg, create_csf, visible, transparency, win_lad)
     6936    SUBROUTINE raytrace(src, targ, isrc, rirrf, atarg, create_csf, visible, transparency)
    62646937        IMPLICIT NONE
    62656938
     
    62716944        LOGICAL, INTENT(out)                   :: visible
    62726945        REAL(wp), INTENT(out)                  :: transparency !< along whole path
    6273         INTEGER(iwp), INTENT(in)               :: win_lad
    62746946        INTEGER(iwp)                           :: i, k, d
    62756947        INTEGER(iwp)                           :: seldim       !< dimension to be incremented
     
    62946966        INTEGER(iwp)                           :: ig           !< 1D index of gridbox in global 2D array
    62956967        REAL(wp)                               :: lad_s_target !< recieved lad_s of particular grid box
    6296         REAL(wp), PARAMETER                    :: grow_factor = 1.5_wp !< factor of expansion of grow arrays
    62976968
    62986969!
     
    63857056        IF ( plant_canopy )  THEN
    63867057#if defined( __parallel )
    6387             IF ( rma_lad_raytrace )  THEN
     7058            IF ( raytrace_mpi_rma )  THEN
    63887059!--             send requests for lad_s to appropriate processor
    6389                 CALL cpu_log( log_point_s(77), 'rad_init_rma', 'start' )
     7060                CALL cpu_log( log_point_s(77), 'rad_rma_lad', 'start' )
    63907061                DO i = 1, ncsb
    63917062                    CALL MPI_Get(lad_s_ray(i), 1, MPI_REAL, lad_ip(i), lad_disp(i), &
    63927063                                 1, MPI_REAL, win_lad, ierr)
    63937064                    IF ( ierr /= 0 )  THEN
    6394                         WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Get'
    6395                         CALL message( 'raytrace', 'PA0519', 1, 2, 0, 6, 0 )
     7065                        WRITE(9,*) 'Error MPI_Get1:', ierr, lad_s_ray(i), &
     7066                                   lad_ip(i), lad_disp(i), win_lad
     7067                        FLUSH(9)
    63967068                    ENDIF
    63977069                ENDDO
     
    64007072                CALL MPI_Win_flush_local_all(win_lad, ierr)
    64017073                IF ( ierr /= 0 )  THEN
    6402                     WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Win_flush_local_all'
    6403                     CALL message( 'raytrace', 'PA0519', 1, 2, 0, 6, 0 )
     7074                    WRITE(9,*) 'Error MPI_Win_flush_local_all1:', ierr, win_lad
     7075                    FLUSH(9)
    64047076                ENDIF
    6405                 CALL cpu_log( log_point_s(77), 'rad_init_rma', 'stop' )
     7077                CALL cpu_log( log_point_s(77), 'rad_rma_lad', 'stop' )
    64067078               
    64077079            ENDIF
     
    64117083            DO i = 1, ncsb
    64127084#if defined( __parallel )
    6413                 IF ( rma_lad_raytrace )  THEN
     7085                IF ( raytrace_mpi_rma )  THEN
    64147086                    lad_s_target = lad_s_ray(i)
    64157087                ELSE
     
    64297101                    acsf(ncsfl)%itz = boxes(1,i)
    64307102                    acsf(ncsfl)%isurfs = isrc
    6431                     acsf(ncsfl)%rsvf = REAL(cursink*rirrf*atarg, wp) !-- we postpone multiplication by transparency
    6432                     acsf(ncsfl)%rtransp = REAL(transparency, wp)
     7103                    acsf(ncsfl)%rsvf = cursink*transparency*rirrf*atarg
    64337104                ENDIF  !< create_csf
    64347105
     
    64527123!> vertical_delta / horizontal_distance
    64537124!------------------------------------------------------------------------------!
    6454    SUBROUTINE raytrace_2d(origin, yxdir, zdirs, iorig, aorig, vffrac, &
    6455                               create_csf, skip_1st_pcb, win_lad, horizon, &
    6456                               transparency)
     7125   SUBROUTINE raytrace_2d(origin, yxdir, nrays, zdirs, iorig, aorig, vffrac, &
     7126                              calc_svf, create_csf, skip_1st_pcb,            &
     7127                              lowest_free_ray, transparency, itarget)
    64577128      IMPLICIT NONE
    64587129
    64597130      REAL(wp), DIMENSION(3), INTENT(IN)     ::  origin        !< z,y,x coordinates of ray origin
    64607131      REAL(wp), DIMENSION(2), INTENT(IN)     ::  yxdir         !< y,x *unit* vector of ray direction (in grid units)
    6461       REAL(wp), DIMENSION(:), INTENT(IN)     ::  zdirs         !< list of z directions to raytrace (z/hdist, in grid)
     7132      INTEGER(iwp)                           ::  nrays         !< number of rays (z directions) to raytrace
     7133      REAL(wp), DIMENSION(nrays), INTENT(IN) ::  zdirs         !< list of z directions to raytrace (z/hdist, grid, zenith->nadir)
    64627134      INTEGER(iwp), INTENT(in)               ::  iorig         !< index of origin face for csf
    64637135      REAL(wp), INTENT(in)                   ::  aorig         !< origin face area for csf
    6464       REAL(wp), DIMENSION(LBOUND(zdirs, 1):UBOUND(zdirs, 1)), INTENT(in) ::  vffrac !<
    6465                                                                !< view factor fractions of each ray for csf
    6466       LOGICAL, INTENT(in)                    ::  create_csf    !< whether to generate new CSFs during raytracing
     7136      REAL(wp), DIMENSION(nrays), INTENT(in) ::  vffrac        !< view factor fractions of each ray for csf
     7137      LOGICAL, INTENT(in)                    ::  calc_svf      !< whether to calculate SFV (identify obstacle surfaces)
     7138      LOGICAL, INTENT(in)                    ::  create_csf    !< whether to create canopy sink factors
    64677139      LOGICAL, INTENT(in)                    ::  skip_1st_pcb  !< whether to skip first plant canopy box during raytracing
    6468       INTEGER(iwp), INTENT(in)               ::  win_lad       !< leaf area density MPI window
    6469       REAL(wp), INTENT(OUT)                  ::  horizon       !< highest horizon found after raytracing (z/hdist)
    6470       REAL(wp), DIMENSION(LBOUND(zdirs, 1):UBOUND(zdirs, 1)), INTENT(OUT) ::  transparency !<
    6471                                                                 !< transparencies of zdirs paths
    6472       !--INTEGER(iwp), DIMENSION(3, LBOUND(zdirs, 1):UBOUND(zdirs, 1)), INTENT(OUT) :: itarget !<
    6473                                                                 !< (z,y,x) coordinates of target faces for zdirs
     7140      INTEGER(iwp), INTENT(out)              ::  lowest_free_ray !< index into zdirs
     7141      REAL(wp), DIMENSION(nrays), INTENT(OUT) ::  transparency !< transparencies of zdirs paths
     7142      INTEGER(iwp), DIMENSION(nrays), INTENT(OUT) ::  itarget  !< global indices of target faces for zdirs
     7143
     7144      INTEGER(iwp), DIMENSION(nrays)         ::  target_procs
     7145      REAL(wp)                               ::  horizon       !< highest horizon found after raytracing (z/hdist)
    64747146      INTEGER(iwp)                           ::  i, k, l, d
    64757147      INTEGER(iwp)                           ::  seldim       !< dimension to be incremented
     
    65067178      INTEGER(iwp)                           ::  iz
    65077179      INTEGER(iwp)                           ::  zsgn
    6508       REAL(wp), PARAMETER                    ::  grow_factor = 1.5_wp !< factor of expansion of grow arrays
     7180      INTEGER(iwp)                           ::  lowest_lad   !< lowest column cell for which we need LAD
     7181      INTEGER(iwp)                           ::  lastdir      !< wall direction before hitting this column
     7182      INTEGER(iwp), DIMENSION(2)             ::  lastcolumn
    65097183
    65107184#if defined( __parallel )
     
    65157189      transparency(:) = 1._wp !-- Pre-set the all rays to transparent before reducing
    65167190      horizon = -HUGE(1._wp)
    6517 
    6518       !--Determine distance to boundary (in 2D xy)
     7191      lowest_free_ray = nrays
     7192      IF ( rad_angular_discretization  .AND.  calc_svf )  THEN
     7193         ALLOCATE(target_surfl(nrays))
     7194         target_surfl(:) = -1
     7195         lastdir = -999
     7196         lastcolumn(:) = -999
     7197      ENDIF
     7198
     7199!--   Determine distance to boundary (in 2D xy)
    65197200      IF ( yxdir(1) > 0._wp )  THEN
    65207201         bdydim = ny + .5_wp !< north global boundary
     
    65487229!--   Since all face coordinates have values *.5 and we'd like to use
    65497230!--   integers, all these have .5 added
    6550       DO d = 1, 2
     7231      DO  d = 1, 2
    65517232          IF ( yxdir(d) == 0._wp )  THEN
    65527233              dimnext(d) = HUGE(1_iwp)
     
    65697250         seldim = minloc(dimnextdist, 1)
    65707251         nextdist = dimnextdist(seldim)
    6571          IF ( nextdist > distance ) nextdist = distance
     7252         IF ( nextdist > distance )  nextdist = distance
    65727253
    65737254         IF ( nextdist > lastdist )  THEN
    65747255            ntrack = ntrack + 1
    65757256            crmid = (lastdist + nextdist) * .5_wp
    6576             column = INT(yxorigin(:) + yxdir(:) * crmid, iwp)  !NINT(yxorigin(:) + yxdir(:) * crmid, iwp)
     7257            column = NINT(yxorigin(:) + yxdir(:) * crmid, iwp)
    65777258
    65787259!--         calculate index of the grid with global indices (column(1),column(2))
     
    65867267               horz_entry = -HUGE(1._wp)
    65877268            ELSE
    6588                horz_entry = (nzterr(ig) - origin(1)) / lastdist
     7269               horz_entry = (REAL(nzterr(ig), wp) + .5_wp - origin(1)) / lastdist
    65897270            ENDIF
    6590             horz_exit = (nzterr(ig) - origin(1)) / nextdist
     7271            horz_exit = (REAL(nzterr(ig), wp) + .5_wp - origin(1)) / nextdist
     7272
     7273            IF ( rad_angular_discretization  .AND.  calc_svf )  THEN
     7274!
     7275!--            Identify vertical obstacles hit by rays in current column
     7276               DO WHILE ( lowest_free_ray > 0 )
     7277                  IF ( zdirs(lowest_free_ray) > horz_entry )  EXIT
     7278!
     7279!--               This may only happen after 1st column, so lastdir and lastcolumn are valid
     7280                  CALL request_itarget(lastdir,                                         &
     7281                        CEILING(-0.5_wp + origin(1) + zdirs(lowest_free_ray)*lastdist), &
     7282                        lastcolumn(1), lastcolumn(2),                                   &
     7283                        target_surfl(lowest_free_ray), target_procs(lowest_free_ray))
     7284                  lowest_free_ray = lowest_free_ray - 1
     7285               ENDDO
     7286!
     7287!--            Identify horizontal obstacles hit by rays in current column
     7288               DO WHILE ( lowest_free_ray > 0 )
     7289                  IF ( zdirs(lowest_free_ray) > horz_exit )  EXIT
     7290                  CALL request_itarget(iup_u, nzterr(ig)+1, column(1), column(2), &
     7291                                       target_surfl(lowest_free_ray),           &
     7292                                       target_procs(lowest_free_ray))
     7293                  lowest_free_ray = lowest_free_ray - 1
     7294               ENDDO
     7295            ENDIF
     7296
    65917297            horizon = MAX(horizon, horz_entry, horz_exit)
    65927298
     
    65987304
    65997305         IF ( nextdist >= distance )  EXIT
     7306
     7307         IF ( rad_angular_discretization  .AND.  calc_svf )  THEN
     7308!
     7309!--         Save wall direction of coming building column (= this air column)
     7310            IF ( seldim == 1 )  THEN
     7311               IF ( dimdelta(seldim) == 1 )  THEN
     7312                  lastdir = isouth_u
     7313               ELSE
     7314                  lastdir = inorth_u
     7315               ENDIF
     7316            ELSE
     7317               IF ( dimdelta(seldim) == 1 )  THEN
     7318                  lastdir = iwest_u
     7319               ELSE
     7320                  lastdir = ieast_u
     7321               ENDIF
     7322            ENDIF
     7323            lastcolumn = column
     7324         ENDIF
    66007325         lastdist = nextdist
    66017326         dimnext(seldim) = dimnext(seldim) + dimdelta(seldim)
     
    66047329
    66057330      IF ( plant_canopy )  THEN
    6606          !--Request LAD WHERE applicable
    6607          !--
     7331!--      Request LAD WHERE applicable
     7332!--     
    66087333#if defined( __parallel )
    6609          IF ( rma_lad_raytrace )  THEN
     7334         IF ( raytrace_mpi_rma )  THEN
    66107335!--         send requests for lad_s to appropriate processor
    66117336            !CALL cpu_log( log_point_s(77), 'usm_init_rma', 'start' )
    6612             DO i = 1, ntrack
     7337            DO  i = 1, ntrack
    66137338               px = rt2_track(2,i)/nnx
    66147339               py = rt2_track(1,i)/nny
    66157340               ip = px*pdims(2)+py
    66167341               ig = ip*nnx*nny + (rt2_track(2,i)-px*nnx)*nny + rt2_track(1,i)-py*nny
    6617                IF ( plantt(ig) <= nzterr(ig) )  CYCLE
    6618                wdisp = (rt2_track(2,i)-px*nnx)*(nny*nzp) + (rt2_track(1,i)-py*nny)*nzp + nzterr(ig)+1-nzub
    6619                wcount = plantt(ig)-nzterr(ig)
     7342
     7343               IF ( rad_angular_discretization  .AND.  calc_svf )  THEN
     7344!
     7345!--               For fixed view resolution, we need plant canopy even for rays
     7346!--               to opposing surfaces
     7347                  lowest_lad = nzterr(ig) + 1
     7348               ELSE
     7349!
     7350!--               We only need LAD for rays directed above horizon (to sky)
     7351                  lowest_lad = CEILING( -0.5_wp + origin(1) +            &
     7352                                    MIN( horizon * rt2_track_dist(i-1),  & ! entry
     7353                                         horizon * rt2_track_dist(i)   ) ) ! exit
     7354               ENDIF
     7355!
     7356!--            Skip asking for LAD where all plant canopy is under requested level
     7357               IF ( plantt(ig) < lowest_lad )  CYCLE
     7358
     7359               wdisp = (rt2_track(2,i)-px*nnx)*(nny*nzp) + (rt2_track(1,i)-py*nny)*nzp + lowest_lad-nzub
     7360               wcount = plantt(ig)-lowest_lad+1
    66207361               ! TODO send request ASAP - even during raytracing
    6621                CALL MPI_Get(rt2_track_lad(nzterr(ig)+1:plantt(ig), i), wcount, MPI_REAL, ip,    &
     7362               CALL MPI_Get(rt2_track_lad(lowest_lad:plantt(ig), i), wcount, MPI_REAL, ip,    &
    66227363                            wdisp, wcount, MPI_REAL, win_lad, ierr)
    66237364               IF ( ierr /= 0 )  THEN
    6624                   WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Get'
    6625                   CALL message( 'raytrace_2d', 'PA0526', 1, 2, 0, 6, 0 )
     7365                  WRITE(9,*) 'Error MPI_Get2:', ierr, rt2_track_lad(lowest_lad:plantt(ig), i), &
     7366                             wcount, ip, wdisp, win_lad
     7367                  FLUSH(9)
    66267368               ENDIF
    66277369            ENDDO
     
    66317373            CALL MPI_Win_flush_local_all(win_lad, ierr)
    66327374            IF ( ierr /= 0 )  THEN
    6633                WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Win_flush_local_all'
    6634                CALL message( 'raytrace', 'PA0527', 1, 2, 0, 6, 0 )
     7375               WRITE(9,*) 'Error MPI_Win_flush_local_all2:', ierr, win_lad
     7376               FLUSH(9)
    66357377            ENDIF
    66367378            !CALL cpu_log( log_point_s(77), 'usm_init_rma', 'stop' )
    6637          ELSE ! rma_lad_raytrace
    6638             DO i = 1, ntrack
     7379
     7380         ELSE ! raytrace_mpi_rma = .F.
     7381            DO  i = 1, ntrack
    66397382               px = rt2_track(2,i)/nnx
    66407383               py = rt2_track(1,i)/nny
     
    66457388         ENDIF
    66467389#else
    6647          DO i = 1, ntrack
     7390         DO  i = 1, ntrack
    66487391            rt2_track_lad(nzub:plantt_max, i) = sub_lad(rt2_track(1,i), rt2_track(2,i), nzub:plantt_max)
    66497392         ENDDO
    66507393#endif
    6651 
    6652          !--Skip the PCB around origin if requested
    6653          !--
    6654          IF ( skip_1st_pcb )  THEN
     7394      ENDIF ! plant_canopy
     7395
     7396      IF ( rad_angular_discretization  .AND.  calc_svf )  THEN
     7397#if defined( __parallel )
     7398!--      wait for all gridsurf requests to complete
     7399         CALL MPI_Win_flush_local_all(win_gridsurf, ierr)
     7400         IF ( ierr /= 0 )  THEN
     7401            WRITE(9,*) 'Error MPI_Win_flush_local_all3:', ierr, win_gridsurf
     7402            FLUSH(9)
     7403         ENDIF
     7404#endif
     7405!
     7406!--      recalculate local surf indices into global ones
     7407         DO i = 1, nrays
     7408            IF ( target_surfl(i) == -1 )  THEN
     7409               itarget(i) = -1
     7410            ELSE
     7411               itarget(i) = target_surfl(i) + surfstart(target_procs(i))
     7412            ENDIF
     7413         ENDDO
     7414         
     7415         DEALLOCATE( target_surfl )
     7416         
     7417      ELSE
     7418         itarget(:) = -1
     7419      ENDIF ! rad_angular_discretization
     7420
     7421      IF ( plant_canopy )  THEN
     7422!--      Skip the PCB around origin if requested (for MRT, the PCB might not be there)
     7423!--     
     7424         IF ( skip_1st_pcb  .AND.  NINT(origin(1)) <= plantt_max )  THEN
    66557425            rt2_track_lad(NINT(origin(1), iwp), 1) = 0._wp
    66567426         ENDIF
    66577427
    6658          !--Assert that we have space allocated for CSFs
    6659          !--
    6660          maxboxes = (ntrack + MAX(origin(1) - nzub, nzut - origin(1))) * SIZE(zdirs, 1)
     7428!--      Assert that we have space allocated for CSFs
     7429!--     
     7430         maxboxes = (ntrack + MAX(CEILING(origin(1)-.5_wp) - nzub,          &
     7431                                  nzpt - CEILING(origin(1)-.5_wp))) * nrays
    66617432         IF ( ncsfl + maxboxes > ncsfla )  THEN
    66627433!--         use this code for growing by fixed exponential increments (equivalent to case where ncsfl always increases by 1)
     
    66687439         ENDIF
    66697440
    6670          !--Calculate transparencies and store new CSFs
    6671          !--
     7441!--      Calculate transparencies and store new CSFs
     7442!--     
    66727443         zbottom = REAL(nzub, wp) - .5_wp
    66737444         ztop = REAL(plantt_max, wp) + .5_wp
    66747445
    6675          !--Reverse direction of radiation (face->sky), only when create_csf
    6676          !--
    6677          IF ( create_csf )  THEN
    6678             DO i = 1, ntrack ! for each column
     7446!--      Reverse direction of radiation (face->sky), only when calc_svf
     7447!--     
     7448         IF ( calc_svf )  THEN
     7449            DO  i = 1, ntrack ! for each column
    66797450               dxxyy = ((dy*yxdir(1))**2 + (dx*yxdir(2))**2) * (rt2_track_dist(i)-rt2_track_dist(i-1))**2
    66807451               px = rt2_track(2,i)/nnx
     
    66827453               ip = px*pdims(2)+py
    66837454
    6684                DO k = LBOUND(zdirs, 1), UBOUND(zdirs, 1) ! for each ray
    6685                   IF ( zdirs(k) <= horizon )  THEN
    6686                      CYCLE
     7455               DO  k = 1, nrays ! for each ray
     7456!
     7457!--               NOTE 6778:
     7458!--               With traditional svf discretization, CSFs under the horizon
     7459!--               (i.e. for surface to surface radiation)  are created in
     7460!--               raytrace(). With rad_angular_discretization, we must create
     7461!--               CSFs under horizon only for one direction, otherwise we would
     7462!--               have duplicate amount of energy. Although we could choose
     7463!--               either of the two directions (they differ only by
     7464!--               discretization error with no bias), we choose the the backward
     7465!--               direction, because it tends to cumulate high canopy sink
     7466!--               factors closer to raytrace origin, i.e. it should potentially
     7467!--               cause less moiree.
     7468                  IF ( .NOT. rad_angular_discretization )  THEN
     7469                     IF ( zdirs(k) <= horizon )  CYCLE
    66877470                  ENDIF
    66887471
    6689                   zorig = REAL(origin(1), wp) + zdirs(k) * rt2_track_dist(i-1)
    6690                   IF ( zorig <= zbottom .OR. zorig >= ztop )  CYCLE
     7472                  zorig = origin(1) + zdirs(k) * rt2_track_dist(i-1)
     7473                  IF ( zorig <= zbottom  .OR. zorig >= ztop )  CYCLE
    66917474
    66927475                  zsgn = INT(SIGN(1._wp, zdirs(k)), iwp)
     
    66957478                     nz = 2
    66967479                     rt2_dist(nz) = SQRT(dxxyy)
    6697                      iz = NINT(zorig, iwp)
     7480                     iz = CEILING(-.5_wp + zorig, iwp)
    66987481                  ELSE
    6699                      zexit = MIN(MAX(REAL(origin(1), wp) + zdirs(k) * rt2_track_dist(i), zbottom), ztop)
     7482                     zexit = MIN(MAX(origin(1) + zdirs(k) * rt2_track_dist(i), zbottom), ztop)
    67007483
    67017484                     zb0 = FLOOR(  zorig * zsgn - .5_wp) + 1  ! because it must be greater than orig
     
    67087491                  ENDIF
    67097492
    6710                   DO l = 2, nz
     7493                  DO  l = 2, nz
    67117494                     IF ( rt2_track_lad(iz, i) > 0._wp )  THEN
    67127495                        curtrans = exp(-ext_coef * rt2_track_lad(iz, i) * (rt2_dist(l)-rt2_dist(l-1)))
    6713                         ncsfl = ncsfl + 1
    6714                         acsf(ncsfl)%ip = ip
    6715                         acsf(ncsfl)%itx = rt2_track(2,i)
    6716                         acsf(ncsfl)%ity = rt2_track(1,i)
    6717                         acsf(ncsfl)%itz = iz
    6718                         acsf(ncsfl)%isurfs = iorig
    6719                         acsf(ncsfl)%rsvf = REAL((1._wp - curtrans)*aorig*vffrac(k), wp) ! we postpone multiplication by transparency
    6720                         acsf(ncsfl)%rtransp = REAL(transparency(k), wp)
     7496
     7497                        IF ( create_csf )  THEN
     7498                           ncsfl = ncsfl + 1
     7499                           acsf(ncsfl)%ip = ip
     7500                           acsf(ncsfl)%itx = rt2_track(2,i)
     7501                           acsf(ncsfl)%ity = rt2_track(1,i)
     7502                           acsf(ncsfl)%itz = iz
     7503                           acsf(ncsfl)%isurfs = iorig
     7504                           acsf(ncsfl)%rsvf = (1._wp - curtrans)*transparency(k)*aorig*vffrac(k)
     7505                        ENDIF
    67217506
    67227507                        transparency(k) = transparency(k) * curtrans
     
    67247509                     iz = iz + zsgn
    67257510                  ENDDO ! l = 1, nz - 1
    6726                ENDDO ! k = LBOUND(zdirs, 1), UBOUND(zdirs, 1)
     7511               ENDDO ! k = 1, nrays
    67277512            ENDDO ! i = 1, ntrack
    67287513
    6729             transparency(:) = 1._wp !-- Reset all rays to transparent
     7514            transparency(1:lowest_free_ray) = 1._wp !-- Reset rays above horizon to transparent (see NOTE 6778)
    67307515         ENDIF
    67317516
    6732          !-- Forward direction of radiation (sky->face), always
    6733          !--
    6734          DO i = ntrack, 1, -1 ! for each column backwards
     7517!--      Forward direction of radiation (sky->face), always
     7518!--     
     7519         DO  i = ntrack, 1, -1 ! for each column backwards
    67357520            dxxyy = ((dy*yxdir(1))**2 + (dx*yxdir(2))**2) * (rt2_track_dist(i)-rt2_track_dist(i-1))**2
    67367521            px = rt2_track(2,i)/nnx
     
    67387523            ip = px*pdims(2)+py
    67397524
    6740             DO k = LBOUND(zdirs, 1), UBOUND(zdirs, 1) ! for each ray
    6741                IF ( zdirs(k) <= horizon )  THEN
    6742                   transparency(k) = 0._wp
    6743                   CYCLE
    6744                ENDIF
    6745 
    6746                zexit = REAL(origin(1), wp) + zdirs(k) * rt2_track_dist(i-1)
    6747                IF ( zexit <= zbottom .OR. zexit >= ztop )  CYCLE
     7525            DO  k = 1, nrays ! for each ray
     7526!
     7527!--            See NOTE 6778 above
     7528               IF ( zdirs(k) <= horizon )  CYCLE
     7529
     7530               zexit = origin(1) + zdirs(k) * rt2_track_dist(i-1)
     7531               IF ( zexit <= zbottom  .OR.  zexit >= ztop )  CYCLE
    67487532
    67497533               zsgn = -INT(SIGN(1._wp, zdirs(k)), iwp)
     
    67547538                  iz = NINT(zexit, iwp)
    67557539               ELSE
    6756                   zorig = MIN(MAX(REAL(origin(1), wp) + zdirs(k) * rt2_track_dist(i), zbottom), ztop)
     7540                  zorig = MIN(MAX(origin(1) + zdirs(k) * rt2_track_dist(i), zbottom), ztop)
    67577541
    67587542                  zb0 = FLOOR(  zorig * zsgn - .5_wp) + 1  ! because it must be greater than orig
     
    67657549               ENDIF
    67667550
    6767                DO l = 2, nz
     7551               DO  l = 2, nz
    67687552                  IF ( rt2_track_lad(iz, i) > 0._wp )  THEN
    67697553                     curtrans = exp(-ext_coef * rt2_track_lad(iz, i) * (rt2_dist(l)-rt2_dist(l-1)))
     
    67757559                        acsf(ncsfl)%ity = rt2_track(1,i)
    67767560                        acsf(ncsfl)%itz = iz
    6777                         acsf(ncsfl)%isurfs = -1 ! a special ID indicating sky
    6778                         acsf(ncsfl)%rsvf = REAL((1._wp - curtrans)*aorig*vffrac(k), wp) ! we postpone multiplication by transparency
    6779                         acsf(ncsfl)%rtransp = REAL(transparency(k), wp)
    6780                      ENDIF  !< create_csf
     7561                        acsf(ncsfl)%isurfs = itarget(k) !if above horizon, then itarget(k)==-1, which
     7562                                                        !is also a special ID indicating sky
     7563                        acsf(ncsfl)%rsvf = (1._wp - curtrans)*transparency(k)*aorig*vffrac(k)
     7564                     ENDIF  ! create_csf
    67817565
    67827566                     transparency(k) = transparency(k) * curtrans
     
    67847568                  iz = iz + zsgn
    67857569               ENDDO ! l = 1, nz - 1
    6786             ENDDO ! k = LBOUND(zdirs, 1), UBOUND(zdirs, 1)
     7570            ENDDO ! k = 1, nrays
    67877571         ENDDO ! i = 1, ntrack
    6788 
    6789       ELSE ! not plant_canopy
    6790          DO k = UBOUND(zdirs, 1), LBOUND(zdirs, 1), -1 ! TODO make more generic
    6791             IF ( zdirs(k) > horizon )  EXIT
    6792             transparency(k) = 0._wp
     7572      ENDIF ! plant_canopy
     7573
     7574      IF ( .NOT. (rad_angular_discretization  .AND.  calc_svf) )  THEN
     7575!
     7576!--      Just update lowest_free_ray according to horizon
     7577         DO WHILE ( lowest_free_ray > 0 )
     7578            IF ( zdirs(lowest_free_ray) > horizon )  EXIT
     7579            lowest_free_ray = lowest_free_ray - 1
    67937580         ENDDO
    67947581      ENDIF
     7582
     7583   CONTAINS
     7584      SUBROUTINE request_itarget(d, z, y, x, isurfl, iproc)
     7585         INTEGER(iwp), INTENT(in)            ::  d, z, y, x
     7586         INTEGER(iwp), TARGET, INTENT(out)   ::  isurfl
     7587         INTEGER(iwp), INTENT(out)           ::  iproc
     7588         INTEGER(kind=MPI_ADDRESS_KIND)      ::  target_displ  !< index of the grid in the local gridsurf array
     7589         INTEGER(iwp)                        ::  px, py        !< number of processors in x and y direction
     7590                                                               !< before the processor in the question
     7591
     7592!--      calculate target processor and index in the remote local target gridsurf array
     7593         px = x/nnx
     7594         py = y/nny
     7595         iproc = px*pdims(2)+py
     7596         target_displ = ((x-px*nnx)*nny + y-py*nny)*nzu*nsurf_type_u + (z-nzub)*nsurf_type_u + d
     7597
     7598#if defined( __parallel )
     7599!--      send MPI_Get request to obtain index target_surfl(i)
     7600         CALL MPI_Get(isurfl, 1, MPI_INTEGER, iproc, target_displ, &
     7601                        1, MPI_INTEGER, win_gridsurf, ierr)
     7602         IF ( ierr /= 0 )  THEN
     7603            WRITE(9,*) 'Error MPI_Get3:', ierr, isurfl, iproc, target_displ, win_gridsurf
     7604            FLUSH(9)
     7605         ENDIF
     7606#else
     7607!--      set index target_surfl(i)
     7608         isurfl = gridsurf(d,z,y,x)
     7609#endif
     7610      END SUBROUTINE request_itarget
    67957611
    67967612   END SUBROUTINE raytrace_2d
     
    68497665      ALLOCATE ( dsitrans(nsurfl, ndsidir) )
    68507666      ALLOCATE ( dsitransc(npcbl, ndsidir) )
     7667      IF ( nmrtbl > 0 )  ALLOCATE ( mrtdsit(nmrtbl, ndsidir) )
    68517668
    68527669      WRITE ( message_string, * ) 'Precalculated', ndsidir, ' solar positions', &
     
    69057722
    69067723!-- first check: are the two surfaces directed in the same direction
    6907         IF ( (d==iup_u  .OR.  d==iup_l  .OR.  d==iup_a )                             &
     7724        IF ( (d==iup_u  .OR.  d==iup_l )                             &
    69087725             .AND. (d2==iup_u  .OR. d2==iup_l) ) RETURN
    6909         IF ( (d==isouth_u  .OR.  d==isouth_l  .OR.  d==isouth_a ) &
     7726        IF ( (d==isouth_u  .OR.  d==isouth_l ) &
    69107727             .AND.  (d2==isouth_u  .OR.  d2==isouth_l) ) RETURN
    6911         IF ( (d==inorth_u  .OR.  d==inorth_l  .OR.  d==inorth_a ) &
     7728        IF ( (d==inorth_u  .OR.  d==inorth_l ) &
    69127729             .AND.  (d2==inorth_u  .OR.  d2==inorth_l) ) RETURN
    6913         IF ( (d==iwest_u  .OR.  d==iwest_l  .OR.  d==iwest_a )     &
     7730        IF ( (d==iwest_u  .OR.  d==iwest_l )     &
    69147731             .AND.  (d2==iwest_u  .OR.  d2==iwest_l ) ) RETURN
    6915         IF ( (d==ieast_u  .OR.  d==ieast_l  .OR.  d==ieast_a )     &
     7732        IF ( (d==ieast_u  .OR.  d==ieast_l )     &
    69167733             .AND.  (d2==ieast_u  .OR.  d2==ieast_l ) ) RETURN
    69177734
    69187735!-- second check: are surfaces facing away from each other
    69197736        SELECT CASE (d)
    6920             CASE (iup_u, iup_l, iup_a)              !< upward facing surfaces
     7737            CASE (iup_u, iup_l)                     !< upward facing surfaces
    69217738                IF ( z2 < z ) RETURN
    6922             CASE (idown_a)                          !< downward facing surfaces
    6923                 IF ( z2 > z ) RETURN
    6924             CASE (isouth_u, isouth_l, isouth_a)     !< southward facing surfaces
     7739            CASE (isouth_u, isouth_l)               !< southward facing surfaces
    69257740                IF ( y2 > y ) RETURN
    6926             CASE (inorth_u, inorth_l, inorth_a)     !< northward facing surfaces
     7741            CASE (inorth_u, inorth_l)               !< northward facing surfaces
    69277742                IF ( y2 < y ) RETURN
    6928             CASE (iwest_u, iwest_l, iwest_a)        !< westward facing surfaces
     7743            CASE (iwest_u, iwest_l)                 !< westward facing surfaces
    69297744                IF ( x2 > x ) RETURN
    6930             CASE (ieast_u, ieast_l, ieast_a)        !< eastward facing surfaces
     7745            CASE (ieast_u, ieast_l)                 !< eastward facing surfaces
    69317746                IF ( x2 < x ) RETURN
    69327747        END SELECT
     
    71437958!>    array for csf
    71447959!------------------------------------------------------------------------------!
     7960!-- quicksort.f -*-f90-*-
     7961!-- Author: t-nissie, adaptation J.Resler
     7962!-- License: GPLv3
     7963!-- Gist: https://gist.github.com/t-nissie/479f0f16966925fa29ea
     7964    RECURSIVE SUBROUTINE quicksort_itarget(itarget, vffrac, ztransp, first, last)
     7965        IMPLICIT NONE
     7966        INTEGER(iwp), DIMENSION(:), INTENT(INOUT)   :: itarget
     7967        REAL(wp), DIMENSION(:), INTENT(INOUT)       :: vffrac, ztransp
     7968        INTEGER(iwp), INTENT(IN)                    :: first, last
     7969        INTEGER(iwp)                                :: x, t
     7970        INTEGER(iwp)                                :: i, j
     7971        REAL(wp)                                    :: tr
     7972
     7973        IF ( first>=last ) RETURN
     7974        x = itarget((first+last)/2)
     7975        i = first
     7976        j = last
     7977        DO
     7978            DO WHILE ( itarget(i) < x )
     7979               i=i+1
     7980            ENDDO
     7981            DO WHILE ( x < itarget(j) )
     7982                j=j-1
     7983            ENDDO
     7984            IF ( i >= j ) EXIT
     7985            t = itarget(i);  itarget(i) = itarget(j);  itarget(j) = t
     7986            tr = vffrac(i);  vffrac(i) = vffrac(j);  vffrac(j) = tr
     7987            tr = ztransp(i);  ztransp(i) = ztransp(j);  ztransp(j) = tr
     7988            i=i+1
     7989            j=j-1
     7990        ENDDO
     7991        IF ( first < i-1 ) CALL quicksort_itarget(itarget, vffrac, ztransp, first, i-1)
     7992        IF ( j+1 < last )  CALL quicksort_itarget(itarget, vffrac, ztransp, j+1, last)
     7993    END SUBROUTINE quicksort_itarget
     7994
    71457995    PURE FUNCTION svf_lt(svf1,svf2) result (res)
    71467996      TYPE (t_svf), INTENT(in) :: svf1,svf2
     
    71538003      ENDIF
    71548004    END FUNCTION svf_lt
    7155    
    7156  
     8005
     8006
    71578007!-- quicksort.f -*-f90-*-
    71588008!-- Author: t-nissie, adaptation J.Resler
     
    71868036    END SUBROUTINE quicksort_svf
    71878037
    7188    
    71898038    PURE FUNCTION csf_lt(csf1,csf2) result (res)
    71908039      TYPE (t_csf), INTENT(in) :: csf1,csf2
     
    72418090        INTEGER(iwp)                            :: iread, iwrite
    72428091        TYPE(t_csf), DIMENSION(:), POINTER      :: acsfnew
    7243 !        CHARACTER(100)                          :: msg
     8092        CHARACTER(100)                          :: msg
    72448093
    72458094        IF ( newsize == -1 )  THEN
     
    72708119                         .AND.  acsfnew(iwrite)%itz == acsf(iread)%itz &
    72718120                         .AND.  acsfnew(iwrite)%isurfs == acsf(iread)%isurfs )  THEN
    7272 !--                 We could simply take either first or second rtransp, both are valid. As a very simple heuristic about which ray
    7273 !--                 probably passes nearer the center of the target box, we choose DIF from the entry with greater CSF, since that
    7274 !--                 might mean that the traced beam passes longer through the canopy box.
    7275                     IF ( acsfnew(iwrite)%rsvf < acsf(iread)%rsvf )  THEN
    7276                         acsfnew(iwrite)%rtransp = acsf(iread)%rtransp
    7277                     ENDIF
     8121
    72788122                    acsfnew(iwrite)%rsvf = acsfnew(iwrite)%rsvf + acsf(iread)%rsvf
    72798123!--                 advance reading index, keep writing index
     
    73108154        ncsfla = newsize
    73118155
    7312 !         WRITE(msg,'(A,2I12)') 'Grow acsf2:',ncsfl,ncsfla
    7313 !         CALL radiation_write_debug_log( msg )
     8156        WRITE(msg,'(A,2I12)') 'Grow acsf2:',ncsfl,ncsfla
     8157        CALL radiation_write_debug_log( msg )
    73148158
    73158159    END SUBROUTINE merge_and_grow_csf
     
    74438287               SELECT CASE ( d )
    74448288
    7445                CASE (iup_u,iup_l,iup_a)                    !- gridbox up_facing face
     8289               CASE (iup_u,iup_l)                          !- gridbox up_facing face
    74468290                  sw_gridbox(1) = surfinsw(l)
    74478291                  lw_gridbox(1) = surfinlw(l)
    74488292                  swd_gridbox(1) = surfinswdif(l)
    74498293
    7450                CASE (idown_a)                         !- gridbox down_facing face
    7451                   sw_gridbox(2) = surfinsw(l)
    7452                   lw_gridbox(2) = surfinlw(l)
    7453                   swd_gridbox(2) = surfinswdif(l)
    7454 
    7455                CASE (inorth_u,inorth_l,inorth_a)  !- gridbox north_facing face
     8294               CASE (inorth_u,inorth_l)                    !- gridbox north_facing face
    74568295                  sw_gridbox(3) = surfinsw(l)
    74578296                  lw_gridbox(3) = surfinlw(l)
    74588297                  swd_gridbox(3) = surfinswdif(l)
    74598298
    7460                CASE (isouth_u,isouth_l,isouth_a)  !- gridbox south_facing face
     8299               CASE (isouth_u,isouth_l)                    !- gridbox south_facing face
    74618300                  sw_gridbox(4) = surfinsw(l)
    74628301                  lw_gridbox(4) = surfinlw(l)
    74638302                  swd_gridbox(4) = surfinswdif(l)
    74648303
    7465                CASE (ieast_u,ieast_l,ieast_a)      !- gridbox east_facing face
     8304               CASE (ieast_u,ieast_l)                      !- gridbox east_facing face
    74668305                  sw_gridbox(5) = surfinsw(l)
    74678306                  lw_gridbox(5) = surfinlw(l)
    74688307                  swd_gridbox(5) = surfinswdif(l)
    74698308
    7470                CASE (iwest_u,iwest_l,iwest_a)      !- gridbox west_facing face
     8309               CASE (iwest_u,iwest_l)                      !- gridbox west_facing face
    74718310                  sw_gridbox(6) = surfinsw(l)
    74728311                  lw_gridbox(6) = surfinlw(l)
     
    75208359    INTEGER(iwp) ::  j !<
    75218360    INTEGER(iwp) ::  k !<
    7522     INTEGER(iwp) ::  m !< index of current surface element
     8361    INTEGER(iwp) ::  l, m !< index of current surface element
    75238362
    75248363    IF ( mode == 'allocate' )  THEN
     
    76048443                rad_sw_hr_av = 0.0_wp
    76058444
     8445             CASE ( 'rad_mrt_sw' )
     8446                IF ( .NOT. ALLOCATED( mrtinsw_av ) )  THEN
     8447                   ALLOCATE( mrtinsw_av(nmrtbl) )
     8448                ENDIF
     8449                mrtinsw_av = 0.0_wp
     8450
     8451             CASE ( 'rad_mrt_lw' )
     8452                IF ( .NOT. ALLOCATED( mrtinlw_av ) )  THEN
     8453                   ALLOCATE( mrtinlw_av(nmrtbl) )
     8454                ENDIF
     8455                mrtinlw_av = 0.0_wp
     8456
     8457             CASE ( 'rad_mrt' )
     8458                IF ( .NOT. ALLOCATED( mrt_av ) )  THEN
     8459                   ALLOCATE( mrt_av(nmrtbl) )
     8460                ENDIF
     8461                mrt_av = 0.0_wp
     8462
    76068463          CASE DEFAULT
    76078464             CONTINUE
     
    78198676             ENDIF
    78208677
     8678          CASE ( 'rad_mrt_sw' )
     8679             IF ( ALLOCATED( mrtinsw_av ) )  THEN
     8680                mrtinsw_av(:) = mrtinsw_av(:) + mrtinsw(:)
     8681             ENDIF
     8682
     8683          CASE ( 'rad_mrt_lw' )
     8684             IF ( ALLOCATED( mrtinlw_av ) )  THEN
     8685                mrtinlw_av(:) = mrtinlw_av(:) + mrtinlw(:)
     8686             ENDIF
     8687
     8688          CASE ( 'rad_mrt' )
     8689             IF ( ALLOCATED( mrt_av ) )  THEN
     8690                mrt_av(:) = mrt_av(:) + mrt(:)
     8691             ENDIF
     8692
    78218693          CASE DEFAULT
    78228694             CONTINUE
     
    78288700       SELECT CASE ( TRIM( variable ) )
    78298701
    7830          CASE ( 'rad_net*' )
     8702          CASE ( 'rad_net*' )
    78318703             IF ( ALLOCATED( rad_net_av ) ) THEN
    78328704                DO  i = nxlg, nxrg
     
    79748846             ENDIF
    79758847
     8848          CASE ( 'rad_mrt_sw' )
     8849             IF ( ALLOCATED( mrtinsw_av ) )  THEN
     8850                mrtinsw_av(:) = mrtinsw_av(:)  / REAL( average_count_3d, KIND=wp )
     8851             ENDIF
     8852
     8853          CASE ( 'rad_mrt_lw' )
     8854             IF ( ALLOCATED( mrtinlw_av ) )  THEN
     8855                mrtinlw_av(:) = mrtinlw_av(:) / REAL( average_count_3d, KIND=wp )
     8856             ENDIF
     8857
     8858          CASE ( 'rad_mrt' )
     8859             IF ( ALLOCATED( mrt_av ) )  THEN
     8860                mrt_av(:) = mrt_av(:) / REAL( average_count_3d, KIND=wp )
     8861             ENDIF
     8862
    79768863       END SELECT
    79778864
     
    80098896              'rad_sw_hr_xy', 'rad_lw_cs_hr_xz', 'rad_lw_hr_xz',               &
    80108897              'rad_sw_cs_hr_xz', 'rad_sw_hr_xz', 'rad_lw_cs_hr_yz',            &
    8011               'rad_lw_hr_yz', 'rad_sw_cs_hr_yz', 'rad_sw_hr_yz' )
     8898              'rad_lw_hr_yz', 'rad_sw_cs_hr_yz', 'rad_sw_hr_yz',               &
     8899              'rad_mrt', 'rad_mrt_sw', 'rad_mrt_lw' )
    80128900          grid_x = 'x'
    80138901          grid_y = 'y'
     
    80378925! Description:
    80388926! ------------
    8039 !> Subroutine defining 3D output variables
     8927!> Subroutine defining 2D output variables
    80408928!------------------------------------------------------------------------------!
    80418929 SUBROUTINE radiation_data_output_2d( av, variable, found, grid, mode,         &
     
    84569344    CHARACTER (LEN=*) ::  variable !<
    84579345
    8458     INTEGER(iwp) ::  av    !<
    8459     INTEGER(iwp) ::  i     !<
    8460     INTEGER(iwp) ::  j     !<
    8461     INTEGER(iwp) ::  k     !<
    8462     INTEGER(iwp) ::  nzb_do   !<
    8463     INTEGER(iwp) ::  nzt_do   !<
    8464 
    8465     LOGICAL      ::  found !<
    8466 
    8467     REAL(wp) ::  fill_value = -999.0_wp    !< value for the _FillValue attribute
    8468 
    8469     REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
    8470 
     9346    INTEGER(iwp) ::  av          !<
     9347    INTEGER(iwp) ::  i, j, k, l  !<
     9348    INTEGER(iwp) ::  nzb_do      !<
     9349    INTEGER(iwp) ::  nzt_do      !<
     9350
     9351    LOGICAL      ::  found       !<
     9352
     9353    REAL(wp)     ::  fill_value = -999.0_wp    !< value for the _FillValue attribute
     9354
     9355    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !<
    84719356
    84729357    found = .TRUE.
     
    86579542               ENDDO
    86589543            ENDDO
     9544         ENDIF
     9545
     9546      CASE ( 'rad_mrt_sw' )
     9547         local_pf = REAL( fill_value, KIND = wp )
     9548         IF ( av == 0 )  THEN
     9549            DO  l = 1, nmrtbl
     9550               i = mrtbl(ix,l)
     9551               j = mrtbl(iy,l)
     9552               k = mrtbl(iz,l)
     9553               local_pf(i,j,k) = mrtinsw(l)
     9554            ENDDO
     9555         ELSE
     9556            IF ( ALLOCATED( mrtinsw_av ) ) THEN
     9557               DO  l = 1, nmrtbl
     9558                  i = mrtbl(ix,l)
     9559                  j = mrtbl(iy,l)
     9560                  k = mrtbl(iz,l)
     9561                  local_pf(i,j,k) = mrtinsw_av(l)
     9562               ENDDO
     9563            ENDIF
     9564         ENDIF
     9565
     9566      CASE ( 'rad_mrt_lw' )
     9567         local_pf = REAL( fill_value, KIND = wp )
     9568         IF ( av == 0 )  THEN
     9569            DO  l = 1, nmrtbl
     9570               i = mrtbl(ix,l)
     9571               j = mrtbl(iy,l)
     9572               k = mrtbl(iz,l)
     9573               local_pf(i,j,k) = mrtinlw(l)
     9574            ENDDO
     9575         ELSE
     9576            IF ( ALLOCATED( mrtinlw_av ) ) THEN
     9577               DO  l = 1, nmrtbl
     9578                  i = mrtbl(ix,l)
     9579                  j = mrtbl(iy,l)
     9580                  k = mrtbl(iz,l)
     9581                  local_pf(i,j,k) = mrtinlw_av(l)
     9582               ENDDO
     9583            ENDIF
     9584         ENDIF
     9585
     9586      CASE ( 'rad_mrt' )
     9587         local_pf = REAL( fill_value, KIND = wp )
     9588         IF ( av == 0 )  THEN
     9589            DO  l = 1, nmrtbl
     9590               i = mrtbl(ix,l)
     9591               j = mrtbl(iy,l)
     9592               k = mrtbl(iz,l)
     9593               local_pf(i,j,k) = mrt(l)
     9594            ENDDO
     9595         ELSE
     9596            IF ( ALLOCATED( mrt_av ) ) THEN
     9597               DO  l = 1, nmrtbl
     9598                  i = mrtbl(ix,l)
     9599                  j = mrtbl(iy,l)
     9600                  k = mrtbl(iz,l)
     9601                  local_pf(i,j,k) = mrt_av(l)
     9602               ENDDO
     9603            ENDIF
    86599604         ENDIF
    86609605
     
    934210287!> Subroutine writes debug information
    934310288!------------------------------------------------------------------------------!
    9344 ! SUBROUTINE radiation_write_debug_log ( message )
     10289 SUBROUTINE radiation_write_debug_log ( message )
    934510290    !> it writes debug log with time stamp
    9346 !    CHARACTER(*)  :: message
    9347 !    CHARACTER(15) :: dtc
    9348 !    CHARACTER(8)  :: date
    9349 !    CHARACTER(10) :: time
    9350 !    CHARACTER(5)  :: zone
    9351 !    CALL date_and_time(date, time, zone)
    9352 !    dtc = date(7:8)//','//time(1:2)//':'//time(3:4)//':'//time(5:10)
    9353 !    WRITE(9,'(2A)') dtc, TRIM(message)
    9354 !    FLUSH(9)
    9355 ! END SUBROUTINE radiation_write_debug_log
     10291    CHARACTER(*)  :: message
     10292    CHARACTER(15) :: dtc
     10293    CHARACTER(8)  :: date
     10294    CHARACTER(10) :: time
     10295    CHARACTER(5)  :: zone
     10296    CALL date_and_time(date, time, zone)
     10297    dtc = date(7:8)//','//time(1:2)//':'//time(3:4)//':'//time(5:10)
     10298    WRITE(9,'(2A)') dtc, TRIM(message)
     10299    FLUSH(9)
     10300 END SUBROUTINE radiation_write_debug_log
    935610301
    935710302 END MODULE radiation_model_mod
  • palm/trunk/SOURCE/sum_up_3d_data.f90

    r3294 r3337  
    2525! -----------------
    2626! $Id$
     27! (from branch resler)
     28! Add biometeorology,
     29! fix chemistry output call,
     30! move usm calls
     31!
     32! 3294 2018-10-01 02:37:10Z raasch
    2733! changes concerning modularization of ocean option
    2834!
     
    260266    USE radiation_model_mod,                                                   &
    261267        ONLY:  radiation, radiation_3d_data_averaging
     268
     269    USE biometeorology_mod,                                                    &
     270        ONLY:  biometeorology_3d_data_averaging
    262271
    263272    USE surface_mod,                                                           &
     
    305314
    306315       DO  ii = 1, doav_n
    307 !
    308 !--       Temporary solution to account for data output within the new urban
    309 !--       surface model (urban_surface_mod.f90), see also SELECT CASE ( trimvar )
     316
    310317          trimvar = TRIM( doav(ii) )
    311           IF ( urban_surface  .AND.  trimvar(1:4) == 'usm_' )  THEN
    312              trimvar = 'usm_output'
    313           ENDIF
    314        
     318
    315319          SELECT CASE ( trimvar )
    316320
     
    495499                z0q_av = 0.0_wp
    496500
    497              CASE ( 'usm_output' )
    498 !             
    499 !--             Block of urban surface model outputs
    500                 CALL usm_average_3d_data( 'allocate', doav(ii) )
    501              
    502501
    503502             CASE DEFAULT
     
    505504!
    506505!--             Allocating and initializing data arrays for other modules
     506
     507                IF ( air_chemistry  .AND. &
     508                     (trimvar(1:3) == 'kc_' .OR. trimvar(1:3) == 'em_') )  THEN
     509                   CALL chem_3d_data_averaging( 'allocate', doav(ii) )
     510                ENDIF
     511
    507512                IF ( bulk_cloud_model )  THEN
    508513                   CALL bcm_3d_data_averaging( 'allocate', doav(ii) )
    509514                ENDIF
    510515
    511                 IF ( air_chemistry  .AND.  trimvar(1:3) == 'kc_')  THEN
    512                    CALL chem_3d_data_averaging( 'allocate', doav(ii) )
     516                IF ( radiation  .AND.  trimvar(1:4) == 'bio_' )  THEN
     517                   CALL biometeorology_3d_data_averaging( 'allocate', doav(ii) )
    513518                ENDIF
    514519
     
    529534                ENDIF
    530535
     536                CALL tcm_3d_data_averaging( 'allocate', doav(ii) )
     537
     538                IF ( urban_surface  .AND.  trimvar(1:4) == 'usm_' )  THEN
     539                   CALL usm_average_3d_data( 'allocate', doav(ii) )
     540                ENDIF
     541
    531542                IF ( uv_exposure  .AND.  trimvar(1:5) == 'uvem_')  THEN
    532543                   CALL uvem_3d_data_averaging( 'allocate', doav(ii) )
    533544                ENDIF
    534545
    535                 CALL tcm_3d_data_averaging( 'allocate', doav(ii) )
    536 
    537546!
    538547!--             User-defined quantities
     
    548557!-- Loop of all variables to be averaged.
    549558    DO  ii = 1, doav_n
    550 !
    551 !--       Temporary solution to account for data output within the new urban
    552 !--       surface model (urban_surface_mod.f90), see also SELECT CASE ( trimvar )
    553           trimvar = TRIM( doav(ii) )
    554           IF ( urban_surface  .AND.  trimvar(1:4) == 'usm_' )  THEN
    555              trimvar = 'usm_output'
    556           ENDIF
     559
     560       trimvar = TRIM( doav(ii) )
    557561!
    558562!--    Store the array chosen on the temporary array.
     
    11461150             ENDIF
    11471151
    1148           CASE ( 'usm_output' )
    1149 !             
    1150 !--          Block of urban surface model outputs.
     1152          CASE DEFAULT
     1153!
     1154!--          Summing up data from other modules
     1155             IF ( bulk_cloud_model )  THEN
     1156                CALL bcm_3d_data_averaging( 'sum', doav(ii) )
     1157             ENDIF
     1158
     1159             IF ( air_chemistry  .AND. &
     1160                  (trimvar(1:3) == 'kc_' .OR. trimvar(1:3) == 'em_') )  THEN
     1161                CALL chem_3d_data_averaging( 'sum',doav(ii) )
     1162             ENDIF
     1163
     1164             IF ( radiation   .AND.  trimvar(1:4) == 'bio_')  THEN
     1165                CALL biometeorology_3d_data_averaging( 'sum', doav(ii) )
     1166             ENDIF
     1167
     1168             IF ( gust_module_enabled )  THEN
     1169                CALL gust_3d_data_averaging( 'sum', doav(ii) )
     1170             ENDIF
     1171
     1172             IF ( land_surface )  THEN
     1173                CALL lsm_3d_data_averaging( 'sum', doav(ii) )
     1174             ENDIF
     1175
     1176             IF ( ocean_mode )  THEN
     1177                CALL ocean_3d_data_averaging( 'sum', doav(ii) )
     1178             ENDIF
     1179
     1180             IF ( radiation )  THEN
     1181                CALL radiation_3d_data_averaging( 'sum', doav(ii) )
     1182             ENDIF
     1183
     1184             CALL tcm_3d_data_averaging( 'sum', doav(ii) )
     1185
    11511186!--          In case of urban surface variables it should be always checked
    11521187!--          if respective arrays are allocated, at least in case of a restart
    11531188!--          run, as averaged usm arrays are not read from file at the moment.
    1154              CALL usm_average_3d_data( 'allocate', doav(ii) )
    1155              CALL usm_average_3d_data( 'sum', doav(ii) )
    1156 
    1157           CASE DEFAULT
    1158 !
    1159 !--          Summing up data from other modules
    1160              IF ( bulk_cloud_model )  THEN
    1161                 CALL bcm_3d_data_averaging( 'sum', doav(ii) )
    1162              ENDIF
    1163 
    1164              IF ( air_chemistry  .AND.  trimvar(1:3) == 'kc_')  THEN
    1165                 CALL chem_3d_data_averaging( 'sum',doav(ii) )
    1166              ENDIF
    1167 
    1168              IF ( gust_module_enabled )  THEN
    1169                 CALL gust_3d_data_averaging( 'sum', doav(ii) )
    1170              ENDIF
    1171 
    1172              IF ( land_surface )  THEN
    1173                 CALL lsm_3d_data_averaging( 'sum', doav(ii) )
    1174              ENDIF
    1175 
    1176              IF ( ocean_mode )  THEN
    1177                 CALL ocean_3d_data_averaging( 'sum', doav(ii) )
    1178              ENDIF
    1179 
    1180              IF ( radiation )  THEN
    1181                 CALL radiation_3d_data_averaging( 'sum', doav(ii) )
    1182              ENDIF
    1183 
    1184              CALL tcm_3d_data_averaging( 'sum', doav(ii) )
     1189             IF ( urban_surface  .AND.  trimvar(1:4) == 'usm_' )  THEN
     1190                CALL usm_average_3d_data( 'allocate', doav(ii) )
     1191                CALL usm_average_3d_data( 'sum', doav(ii) )
     1192             ENDIF
    11851193
    11861194             IF ( uv_exposure )  THEN
  • palm/trunk/SOURCE/time_integration_spinup.f90

    • Property svn:mergeinfo set to (toggle deleted branches)
      /palm/branches/chemistry/SOURCE/time_integration_spinup.f902047-3190,​3218-3297
      /palm/branches/palm4u/SOURCE/time_integration_spinup.f902540-2692
      /palm/branches/rans/SOURCE/time_integration_spinup.f902078-3128
      /palm/branches/resler/SOURCE/time_integration_spinup.f902023-3336
      /palm/branches/forwind/SOURCE/time_integration_spinup.f901564-1913
      /palm/branches/fricke/SOURCE/time_integration_spinup.f90942-977
      /palm/branches/hoffmann/SOURCE/time_integration_spinup.f90989-1052
      /palm/branches/letzel/masked_output/SOURCE/time_integration_spinup.f90296-409
      /palm/branches/suehring/time_integration_spinup.f90423-666
    r3274 r3337  
    2525! -----------------
    2626! $Id$
     27! (from branch resler)
     28! Add pt1 initialization
     29!
     30! 3274 2018-09-24 15:42:55Z knoop
    2731! Modularization of all bulk cloud physics code components
    2832!
     
    295299                k   = surf_usm_h%k(m)
    296300                pt(k,j,i) = pt_spinup
     301                !!!!!!!!!!!!!!!!HACK!!!!!!!!!!!!!
     302                surf_usm_h%pt1 = pt_spinup
     303                !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    297304             ENDDO
    298305
     
    303310                   k   = surf_usm_v(l)%k(m)
    304311                   pt(k,j,i) = pt_spinup
     312                   !!!!!!!!!!!!!!!!HACK!!!!!!!!!!!!!
     313                   surf_usm_v(l)%pt1 = pt_spinup
     314                   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    305315                ENDDO
    306316             ENDDO
  • palm/trunk/SOURCE/urban_surface_mod.f90

    r3274 r3337  
    2828! -----------------
    2929! $Id$
     30! Add output variables usm_rad_pc_inlw, usm_rad_pc_insw*
     31!
     32! 3274 2018-09-24 15:42:55Z knoop
    3033! Modularization of all bulk cloud physics code components
    3134!
     
    349352               spinup_pt_mean, spinup_time, time_do3d, dt_do3d,                &
    350353               average_count_3d, varnamelength, urban_surface,                 &
    351                plant_canopy
     354               plant_canopy, dz
    352355
    353356    USE cpulog,                                                                &
     
    381384               surfinl, surfinlwdif, rad_sw_in_dir, rad_sw_in_diff,            &
    382385               rad_lw_in_diff, surfouts, surfoutl, surfoutsl, surfoutll, surf, &
    383                surfl, nsurfl, nsurfs, surfstart, pcbinsw, pcbinlw,             &
    384                iup_u, inorth_u, isouth_u, ieast_u, iwest_u, iup_l,            &
     386               surfl, nsurfl, pcbinsw, pcbinlw, pcbinswdir,                    &
     387               pcbinswdif, iup_u, inorth_u, isouth_u, ieast_u, iwest_u, iup_l, &
    385388               inorth_l, isouth_l, ieast_l, iwest_l, id,                       &
    386                iz, iy, ix, idir, jdir, kdir,  nsurf_type, nsurf, idsvf, ndsvf, &
    387                iup_a, idown_a, inorth_a, isouth_a, ieast_a, iwest_a,           &
     389               iz, iy, ix,  nsurf, idsvf, ndsvf,                               &
    388390               idcsf, ndcsf, kdcsf, pct,                                       &
    389                startland, endland, startwall, endwall, skyvf, skyvft
     391               startland, endland, startwall, endwall, skyvf, skyvft, nzub,    &
     392               nzut, nzpt, npcbl, pcbl
    390393
    391394    USE statistics,                                                            &
     
    414417    INTEGER(iwp) ::  pedestrian_category = 2         !< default category for wall surface in pedestrian zone
    415418    INTEGER(iwp) ::  roof_category = 2               !< default category for root surface
     419    REAL(wp)     ::  roughness_concrete = 0.001_wp   !< roughness length of average concrete surface
    416420!
    417421!-- Indices of input attributes for (above) ground floor level
     
    605609    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfins_av       !< average of array of residua of sw radiation absorbed in surface after last reflection
    606610    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinl_av       !< average of array of residua of lw radiation absorbed in surface after last reflection
     611    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinlw_av       !< Average of pcbinlw
     612    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinsw_av       !< Average of pcbinsw
     613    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdir_av    !< Average of pcbinswdir
     614    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdif_av    !< Average of pcbinswdif
     615    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswref_av    !< Average of pcbinswref
    607616   
    608617
     
    12601269!--     find the real name of the variable
    12611270        ids = -1
     1271        l = -1
    12621272        var = TRIM(variable)
    12631273        DO i = 0, nd-1
    12641274            k = len(TRIM(var))
    12651275            j = len(TRIM(dirname(i)))
    1266             IF ( var(k-j+1:k) == dirname(i) )  THEN
     1276            IF ( TRIM(var(k-j+1:k)) == TRIM(dirname(i)) )  THEN
    12671277                ids = i
    12681278                idsint = dirint(ids)
     
    12711281            ENDIF
    12721282        ENDDO
     1283        l = idsint - 2  ! horisontal direction index - terible hack !
     1284        IF ( l < 0 .OR. l > 3 ) THEN
     1285           l = -1
     1286        END IF
    12731287        IF ( ids == -1 )  THEN
    12741288            var = TRIM(variable)
     
    13111325                CASE ( 'usm_rad_net' )
    13121326!--                 array of complete radiation balance
    1313                     IF ( .NOT.  ALLOCATED(surf_usm_h%rad_net_av) )  THEN
     1327                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%rad_net_av) )  THEN
    13141328                        ALLOCATE( surf_usm_h%rad_net_av(1:surf_usm_h%ns) )
    13151329                        surf_usm_h%rad_net_av = 0.0_wp
     1330                    ELSE
     1331                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%rad_net_av) )  THEN
     1332                            ALLOCATE( surf_usm_v(l)%rad_net_av(1:surf_usm_v(l)%ns) )
     1333                            surf_usm_v(l)%rad_net_av = 0.0_wp
     1334                        ENDIF
    13161335                    ENDIF
    1317                     DO  l = 0, 3
    1318                        IF ( .NOT.  ALLOCATED(surf_usm_v(l)%rad_net_av) )  THEN
    1319                            ALLOCATE( surf_usm_v(l)%rad_net_av(1:surf_usm_v(l)%ns) )
    1320                            surf_usm_v(l)%rad_net_av = 0.0_wp
    1321                        ENDIF
    1322                     ENDDO
    13231336                   
    13241337                CASE ( 'usm_rad_insw' )
     
    13981411                    ENDIF
    13991412                                   
     1413                CASE ( 'usm_rad_pc_inlw' )
     1414!--                 array of of lw radiation absorbed in plant canopy
     1415                    IF ( .NOT.  ALLOCATED(pcbinlw_av) )  THEN
     1416                        ALLOCATE( pcbinlw_av(1:npcbl) )
     1417                        pcbinlw_av = 0.0_wp
     1418                    ENDIF
     1419                                   
     1420                CASE ( 'usm_rad_pc_insw' )
     1421!--                 array of of sw radiation absorbed in plant canopy
     1422                    IF ( .NOT.  ALLOCATED(pcbinsw_av) )  THEN
     1423                        ALLOCATE( pcbinsw_av(1:npcbl) )
     1424                        pcbinsw_av = 0.0_wp
     1425                    ENDIF
     1426                                   
     1427                CASE ( 'usm_rad_pc_inswdir' )
     1428!--                 array of of direct sw radiation absorbed in plant canopy
     1429                    IF ( .NOT.  ALLOCATED(pcbinswdir_av) )  THEN
     1430                        ALLOCATE( pcbinswdir_av(1:npcbl) )
     1431                        pcbinswdir_av = 0.0_wp
     1432                    ENDIF
     1433                                   
     1434                CASE ( 'usm_rad_pc_inswdif' )
     1435!--                 array of of diffuse sw radiation absorbed in plant canopy
     1436                    IF ( .NOT.  ALLOCATED(pcbinswdif_av) )  THEN
     1437                        ALLOCATE( pcbinswdif_av(1:npcbl) )
     1438                        pcbinswdif_av = 0.0_wp
     1439                    ENDIF
     1440                                   
     1441                CASE ( 'usm_rad_pc_inswref' )
     1442!--                 array of of reflected sw radiation absorbed in plant canopy
     1443                    IF ( .NOT.  ALLOCATED(pcbinswref_av) )  THEN
     1444                        ALLOCATE( pcbinswref_av(1:npcbl) )
     1445                        pcbinswref_av = 0.0_wp
     1446                    ENDIF
     1447                                   
    14001448                CASE ( 'usm_rad_hf' )
    14011449!--                 array of heat flux from radiation for surfaces after i-th reflection
    1402                     IF ( .NOT.  ALLOCATED(surf_usm_h%surfhf_av) )  THEN
     1450                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%surfhf_av) )  THEN
    14031451                        ALLOCATE( surf_usm_h%surfhf_av(1:surf_usm_h%ns) )
    14041452                        surf_usm_h%surfhf_av = 0.0_wp
    1405                     ENDIF
    1406                     DO  l = 0, 3
     1453                    ELSE
    14071454                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%surfhf_av) )  THEN
    14081455                           ALLOCATE( surf_usm_v(l)%surfhf_av(1:surf_usm_v(l)%ns) )
    14091456                           surf_usm_v(l)%surfhf_av = 0.0_wp
    14101457                       ENDIF
    1411                     ENDDO
     1458                    ENDIF
    14121459
    14131460                CASE ( 'usm_wshf' )
    14141461!--                 array of sensible heat flux from surfaces
    14151462!--                 land surfaces
    1416                     IF ( .NOT.  ALLOCATED(surf_usm_h%wshf_eb_av) )  THEN
     1463                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%wshf_eb_av) )  THEN
    14171464                        ALLOCATE( surf_usm_h%wshf_eb_av(1:surf_usm_h%ns) )
    14181465                        surf_usm_h%wshf_eb_av = 0.0_wp
    1419                     ENDIF
    1420                     DO  l = 0, 3
     1466                    ELSE
    14211467                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wshf_eb_av) )  THEN
    14221468                           ALLOCATE( surf_usm_v(l)%wshf_eb_av(1:surf_usm_v(l)%ns) )
    14231469                           surf_usm_v(l)%wshf_eb_av = 0.0_wp
    14241470                       ENDIF
    1425                     ENDDO
     1471                    ENDIF
    14261472!
    14271473!--             Please note, the following output quantities belongs to the
     
    14311477                CASE ( 'usm_wghf' )
    14321478!--                 array of heat flux from ground (wall, roof, land)
    1433                     IF ( .NOT.  ALLOCATED(surf_usm_h%wghf_eb_av) )  THEN
     1479                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%wghf_eb_av) )  THEN
    14341480                        ALLOCATE( surf_usm_h%wghf_eb_av(1:surf_usm_h%ns) )
    14351481                        surf_usm_h%wghf_eb_av = 0.0_wp
    1436                     ENDIF
    1437                     DO  l = 0, 3
     1482                    ELSE
    14381483                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wghf_eb_av) )  THEN
    14391484                           ALLOCATE( surf_usm_v(l)%wghf_eb_av(1:surf_usm_v(l)%ns) )
    14401485                           surf_usm_v(l)%wghf_eb_av = 0.0_wp
    14411486                       ENDIF
    1442                     ENDDO
     1487                    ENDIF
    14431488
    14441489                CASE ( 'usm_wghf_window' )
    14451490!--                 array of heat flux from window ground (wall, roof, land)
    1446                     IF ( .NOT.  ALLOCATED(surf_usm_h%wghf_eb_window_av) )  THEN
     1491                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%wghf_eb_window_av) )  THEN
    14471492                        ALLOCATE( surf_usm_h%wghf_eb_window_av(1:surf_usm_h%ns) )
    14481493                        surf_usm_h%wghf_eb_window_av = 0.0_wp
    1449                     ENDIF
    1450                     DO  l = 0, 3
     1494                    ELSE
    14511495                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wghf_eb_window_av) )  THEN
    14521496                           ALLOCATE( surf_usm_v(l)%wghf_eb_window_av(1:surf_usm_v(l)%ns) )
    14531497                           surf_usm_v(l)%wghf_eb_window_av = 0.0_wp
    14541498                       ENDIF
    1455                     ENDDO
     1499                    ENDIF
    14561500
    14571501                CASE ( 'usm_wghf_green' )
    14581502!--                 array of heat flux from green ground (wall, roof, land)
    1459                     IF ( .NOT.  ALLOCATED(surf_usm_h%wghf_eb_green_av) )  THEN
     1503                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%wghf_eb_green_av) )  THEN
    14601504                        ALLOCATE( surf_usm_h%wghf_eb_green_av(1:surf_usm_h%ns) )
    14611505                        surf_usm_h%wghf_eb_green_av = 0.0_wp
    1462                     ENDIF
    1463                     DO  l = 0, 3
     1506                    ELSE
    14641507                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%wghf_eb_green_av) )  THEN
    14651508                           ALLOCATE( surf_usm_v(l)%wghf_eb_green_av(1:surf_usm_v(l)%ns) )
    14661509                           surf_usm_v(l)%wghf_eb_green_av = 0.0_wp
    14671510                       ENDIF
    1468                     ENDDO
     1511                    ENDIF
    14691512
    14701513                CASE ( 'usm_iwghf' )
    14711514!--                 array of heat flux from indoor ground (wall, roof, land)
    1472                     IF ( .NOT.  ALLOCATED(surf_usm_h%iwghf_eb_av) )  THEN
     1515                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%iwghf_eb_av) )  THEN
    14731516                        ALLOCATE( surf_usm_h%iwghf_eb_av(1:surf_usm_h%ns) )
    14741517                        surf_usm_h%iwghf_eb_av = 0.0_wp
    1475                     ENDIF
    1476                     DO  l = 0, 3
     1518                    ELSE
    14771519                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%iwghf_eb_av) )  THEN
    14781520                           ALLOCATE( surf_usm_v(l)%iwghf_eb_av(1:surf_usm_v(l)%ns) )
    14791521                           surf_usm_v(l)%iwghf_eb_av = 0.0_wp
    14801522                       ENDIF
    1481                     ENDDO
     1523                    ENDIF
    14821524
    14831525                CASE ( 'usm_iwghf_window' )
    14841526!--                 array of heat flux from indoor window ground (wall, roof, land)
    1485                     IF ( .NOT.  ALLOCATED(surf_usm_h%iwghf_eb_window_av) )  THEN
     1527                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%iwghf_eb_window_av) )  THEN
    14861528                        ALLOCATE( surf_usm_h%iwghf_eb_window_av(1:surf_usm_h%ns) )
    14871529                        surf_usm_h%iwghf_eb_window_av = 0.0_wp
    1488                     ENDIF
    1489                     DO  l = 0, 3
     1530                    ELSE
    14901531                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%iwghf_eb_window_av) )  THEN
    14911532                           ALLOCATE( surf_usm_v(l)%iwghf_eb_window_av(1:surf_usm_v(l)%ns) )
    14921533                           surf_usm_v(l)%iwghf_eb_window_av = 0.0_wp
    14931534                       ENDIF
    1494                     ENDDO
     1535                    ENDIF
    14951536                   
    14961537                CASE ( 'usm_t_surf' )
    14971538!--                 surface temperature for surfaces
    1498                     IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_av) )  THEN
     1539                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%t_surf_av) )  THEN
    14991540                        ALLOCATE( surf_usm_h%t_surf_av(1:surf_usm_h%ns) )
    15001541                        surf_usm_h%t_surf_av = 0.0_wp
    1501                     ENDIF
    1502                     DO  l = 0, 3
     1542                    ELSE
    15031543                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_av) )  THEN
    15041544                           ALLOCATE( surf_usm_v(l)%t_surf_av(1:surf_usm_v(l)%ns) )
    15051545                           surf_usm_v(l)%t_surf_av = 0.0_wp
    15061546                       ENDIF
    1507                     ENDDO
     1547                    ENDIF
    15081548
    15091549                CASE ( 'usm_t_surf_window' )
    15101550!--                 surface temperature for window surfaces
    1511                     IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_window_av) )  THEN
     1551                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%t_surf_window_av) )  THEN
    15121552                        ALLOCATE( surf_usm_h%t_surf_window_av(1:surf_usm_h%ns) )
    15131553                        surf_usm_h%t_surf_window_av = 0.0_wp
    1514                     ENDIF
    1515                     DO  l = 0, 3
     1554                    ELSE
    15161555                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_window_av) )  THEN
    15171556                           ALLOCATE( surf_usm_v(l)%t_surf_window_av(1:surf_usm_v(l)%ns) )
    15181557                           surf_usm_v(l)%t_surf_window_av = 0.0_wp
    15191558                       ENDIF
    1520                     ENDDO
     1559                    ENDIF
    15211560                   
    15221561                CASE ( 'usm_t_surf_green' )
    15231562!--                 surface temperature for green surfaces
    1524                     IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_green_av) )  THEN
     1563                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%t_surf_green_av) )  THEN
    15251564                        ALLOCATE( surf_usm_h%t_surf_green_av(1:surf_usm_h%ns) )
    15261565                        surf_usm_h%t_surf_green_av = 0.0_wp
    1527                     ENDIF
    1528                     DO  l = 0, 3
     1566                    ELSE
    15291567                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_green_av) )  THEN
    15301568                           ALLOCATE( surf_usm_v(l)%t_surf_green_av(1:surf_usm_v(l)%ns) )
    15311569                           surf_usm_v(l)%t_surf_green_av = 0.0_wp
    15321570                       ENDIF
    1533                     ENDDO
     1571                    ENDIF
    15341572               
    15351573                CASE ( 'usm_t_surf_10cm' )
    15361574!--                 near surface temperature for whole surfaces
    1537                     IF ( .NOT.  ALLOCATED(surf_usm_h%t_surf_10cm_av) )  THEN
     1575                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%t_surf_10cm_av) )  THEN
    15381576                        ALLOCATE( surf_usm_h%t_surf_10cm_av(1:surf_usm_h%ns) )
    15391577                        surf_usm_h%t_surf_10cm_av = 0.0_wp
    1540                     ENDIF
    1541                     DO  l = 0, 3
     1578                    ELSE
    15421579                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_surf_10cm_av) )  THEN
    15431580                           ALLOCATE( surf_usm_v(l)%t_surf_10cm_av(1:surf_usm_v(l)%ns) )
    15441581                           surf_usm_v(l)%t_surf_10cm_av = 0.0_wp
    15451582                       ENDIF
    1546                     ENDDO
     1583                    ENDIF
    15471584
    15481585                CASE ( 'usm_t_wall' )
    15491586!--                 wall temperature for iwl layer of walls and land
    1550                     IF ( .NOT.  ALLOCATED(surf_usm_h%t_wall_av) )  THEN
     1587                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%t_wall_av) )  THEN
    15511588                        ALLOCATE( surf_usm_h%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
    15521589                        surf_usm_h%t_wall_av = 0.0_wp
    1553                     ENDIF
    1554                     DO  l = 0, 3
     1590                    ELSE
    15551591                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_wall_av) )  THEN
    15561592                           ALLOCATE( surf_usm_v(l)%t_wall_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
    15571593                           surf_usm_v(l)%t_wall_av = 0.0_wp
    15581594                       ENDIF
    1559                     ENDDO
     1595                    ENDIF
    15601596
    15611597                CASE ( 'usm_t_window' )
    15621598!--                 window temperature for iwl layer of walls and land
    1563                     IF ( .NOT.  ALLOCATED(surf_usm_h%t_window_av) )  THEN
     1599                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%t_window_av) )  THEN
    15641600                        ALLOCATE( surf_usm_h%t_window_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
    15651601                        surf_usm_h%t_window_av = 0.0_wp
    1566                     ENDIF
    1567                     DO  l = 0, 3
     1602                    ELSE
    15681603                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_window_av) )  THEN
    15691604                           ALLOCATE( surf_usm_v(l)%t_window_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
    15701605                           surf_usm_v(l)%t_window_av = 0.0_wp
    15711606                       ENDIF
    1572                     ENDDO
     1607                    ENDIF
    15731608
    15741609                CASE ( 'usm_t_green' )
    15751610!--                 green temperature for iwl layer of walls and land
    1576                     IF ( .NOT.  ALLOCATED(surf_usm_h%t_green_av) )  THEN
     1611                    IF ( l == -1 .AND. .NOT.  ALLOCATED(surf_usm_h%t_green_av) )  THEN
    15771612                        ALLOCATE( surf_usm_h%t_green_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) )
    15781613                        surf_usm_h%t_green_av = 0.0_wp
    1579                     ENDIF
    1580                     DO  l = 0, 3
     1614                    ELSE
    15811615                       IF ( .NOT.  ALLOCATED(surf_usm_v(l)%t_green_av) )  THEN
    15821616                           ALLOCATE( surf_usm_v(l)%t_green_av(nzb_wall:nzt_wall,1:surf_usm_v(l)%ns) )
    15831617                           surf_usm_v(l)%t_green_av = 0.0_wp
    15841618                       ENDIF
    1585                     ENDDO
     1619                    ENDIF
    15861620
    15871621               CASE DEFAULT
     
    15961630                CASE ( 'usm_rad_net' )
    15971631!--                 array of complete radiation balance
    1598                     DO  m = 1, surf_usm_h%ns
    1599                        surf_usm_h%rad_net_av(m) =                              &
    1600                                           surf_usm_h%rad_net_av(m) +           &
    1601                                           surf_usm_h%rad_net_l(m)
    1602                     ENDDO
    1603                     DO  l = 0, 3
     1632                    IF ( l == -1 ) THEN
     1633                       DO  m = 1, surf_usm_h%ns
     1634                          surf_usm_h%rad_net_av(m) =                              &
     1635                                             surf_usm_h%rad_net_av(m) +           &
     1636                                             surf_usm_h%rad_net_l(m)
     1637                       ENDDO
     1638                    ELSE
    16041639                       DO  m = 1, surf_usm_v(l)%ns
    16051640                          surf_usm_v(l)%rad_net_av(m) =                        &
     
    16071642                                          surf_usm_v(l)%rad_net_l(m)
    16081643                       ENDDO
    1609                     ENDDO
    1610                    
     1644                    ENDIF
     1645
    16111646                CASE ( 'usm_rad_insw' )
    16121647!--                 array of sw radiation falling to surface after i-th reflection
     
    17001735                    ENDDO
    17011736                   
     1737                CASE ( 'usm_rad_pc_inlw' )
     1738                    pcbinlw_av(:) = pcbinlw_av(:) + pcbinlw(:)
     1739                   
     1740                CASE ( 'usm_rad_pc_insw' )
     1741                    pcbinsw_av(:) = pcbinsw_av(:) + pcbinsw(:)
     1742                   
     1743                CASE ( 'usm_rad_pc_inswdir' )
     1744                    pcbinswdir_av(:) = pcbinswdir_av(:) + pcbinswdir(:)
     1745                   
     1746                CASE ( 'usm_rad_pc_inswdif' )
     1747                    pcbinswdif_av(:) = pcbinswdif_av(:) + pcbinswdif(:)
     1748                   
     1749                CASE ( 'usm_rad_pc_inswref' )
     1750                    pcbinswref_av(:) = pcbinswref_av(:) + pcbinsw(:)     &
     1751                                                        - pcbinswdir(:)  &
     1752                                                        - pcbinswdif(:)
     1753                   
    17021754                CASE ( 'usm_rad_hf' )
    17031755!--                 array of heat flux from radiation for surfaces after i-th reflection
    1704                     DO  m = 1, surf_usm_h%ns
    1705                        surf_usm_h%surfhf_av(m) =                               &
    1706                                           surf_usm_h%surfhf_av(m) +            &
    1707                                           surf_usm_h%surfhf(m)
    1708                     ENDDO
    1709                     DO  l = 0, 3
     1756                    IF ( l == -1 ) THEN
     1757                       DO  m = 1, surf_usm_h%ns
     1758                          surf_usm_h%surfhf_av(m) =                               &
     1759                                             surf_usm_h%surfhf_av(m) +            &
     1760                                             surf_usm_h%surfhf(m)
     1761                       ENDDO
     1762                    ELSE
    17101763                       DO  m = 1, surf_usm_v(l)%ns
    17111764                          surf_usm_v(l)%surfhf_av(m) =                         &
     
    17131766                                          surf_usm_v(l)%surfhf(m)
    17141767                       ENDDO
    1715                     ENDDO
     1768                    ENDIF
    17161769                   
    17171770                CASE ( 'usm_wshf' )
    17181771!--                 array of sensible heat flux from surfaces (land, roof, wall)
    1719                     DO  m = 1, surf_usm_h%ns
    1720                        surf_usm_h%wshf_eb_av(m) =                              &
    1721                                           surf_usm_h%wshf_eb_av(m) +           &
    1722                                           surf_usm_h%wshf_eb(m)
    1723                     ENDDO
    1724                     DO  l = 0, 3
     1772                    IF ( l == -1 ) THEN
     1773                       DO  m = 1, surf_usm_h%ns
     1774                          surf_usm_h%wshf_eb_av(m) =                              &
     1775                                             surf_usm_h%wshf_eb_av(m) +           &
     1776                                             surf_usm_h%wshf_eb(m)
     1777                       ENDDO
     1778                    ELSE
    17251779                       DO  m = 1, surf_usm_v(l)%ns
    17261780                          surf_usm_v(l)%wshf_eb_av(m) =                        &
     
    17281782                                          surf_usm_v(l)%wshf_eb(m)
    17291783                       ENDDO
    1730                     ENDDO
     1784                    ENDIF
    17311785                   
    17321786                CASE ( 'usm_wghf' )
    17331787!--                 array of heat flux from ground (wall, roof, land)
    1734                     DO  m = 1, surf_usm_h%ns
    1735                        surf_usm_h%wghf_eb_av(m) =                              &
    1736                                           surf_usm_h%wghf_eb_av(m) +           &
    1737                                           surf_usm_h%wghf_eb(m)
    1738                     ENDDO
    1739                     DO  l = 0, 3
     1788                    IF ( l == -1 ) THEN
     1789                       DO  m = 1, surf_usm_h%ns
     1790                          surf_usm_h%wghf_eb_av(m) =                              &
     1791                                             surf_usm_h%wghf_eb_av(m) +           &
     1792                                             surf_usm_h%wghf_eb(m)
     1793                       ENDDO
     1794                    ELSE
    17401795                       DO  m = 1, surf_usm_v(l)%ns
    17411796                          surf_usm_v(l)%wghf_eb_av(m) =                        &
     
    17431798                                          surf_usm_v(l)%wghf_eb(m)
    17441799                       ENDDO
    1745                     ENDDO
     1800                    ENDIF
    17461801                   
    17471802                CASE ( 'usm_wghf_window' )
    17481803!--                 array of heat flux from window ground (wall, roof, land)
    1749                     DO  m = 1, surf_usm_h%ns
    1750                        surf_usm_h%wghf_eb_window_av(m) =                              &
    1751                                           surf_usm_h%wghf_eb_window_av(m) +           &
    1752                                           surf_usm_h%wghf_eb_window(m)
    1753                     ENDDO
    1754                     DO  l = 0, 3
     1804                    IF ( l == -1 ) THEN
     1805                       DO  m = 1, surf_usm_h%ns
     1806                          surf_usm_h%wghf_eb_window_av(m) =                              &
     1807                                             surf_usm_h%wghf_eb_window_av(m) +           &
     1808                                             surf_usm_h%wghf_eb_window(m)
     1809                       ENDDO
     1810                    ELSE
    17551811                       DO  m = 1, surf_usm_v(l)%ns
    17561812                          surf_usm_v(l)%wghf_eb_window_av(m) =                        &
     
    17581814                                          surf_usm_v(l)%wghf_eb_window(m)
    17591815                       ENDDO
    1760                     ENDDO
     1816                    ENDIF
    17611817
    17621818                CASE ( 'usm_wghf_green' )
    17631819!--                 array of heat flux from green ground (wall, roof, land)
    1764                     DO  m = 1, surf_usm_h%ns
    1765                        surf_usm_h%wghf_eb_green_av(m) =                              &
    1766                                           surf_usm_h%wghf_eb_green_av(m) +           &
    1767                                           surf_usm_h%wghf_eb_green(m)
    1768                     ENDDO
    1769                     DO  l = 0, 3
     1820                    IF ( l == -1 ) THEN
     1821                       DO  m = 1, surf_usm_h%ns
     1822                          surf_usm_h%wghf_eb_green_av(m) =                              &
     1823                                             surf_usm_h%wghf_eb_green_av(m) +           &
     1824                                             surf_usm_h%wghf_eb_green(m)
     1825                       ENDDO
     1826                    ELSE
    17701827                       DO  m = 1, surf_usm_v(l)%ns
    17711828                          surf_usm_v(l)%wghf_eb_green_av(m) =                        &
     
    17731830                                          surf_usm_v(l)%wghf_eb_green(m)
    17741831                       ENDDO
    1775                     ENDDO
     1832                    ENDIF
    17761833                   
    17771834                CASE ( 'usm_iwghf' )
    17781835!--                 array of heat flux from indoor ground (wall, roof, land)
    1779                     DO  m = 1, surf_usm_h%ns
    1780                        surf_usm_h%iwghf_eb_av(m) =                              &
    1781                                           surf_usm_h%iwghf_eb_av(m) +           &
    1782                                           surf_usm_h%iwghf_eb(m)
    1783                     ENDDO
    1784                     DO  l = 0, 3
     1836                    IF ( l == -1 ) THEN
     1837                       DO  m = 1, surf_usm_h%ns
     1838                          surf_usm_h%iwghf_eb_av(m) =                              &
     1839                                             surf_usm_h%iwghf_eb_av(m) +           &
     1840                                             surf_usm_h%iwghf_eb(m)
     1841                       ENDDO
     1842                    ELSE
    17851843                       DO  m = 1, surf_usm_v(l)%ns
    17861844                          surf_usm_v(l)%iwghf_eb_av(m) =                        &
     
    17881846                                          surf_usm_v(l)%iwghf_eb(m)
    17891847                       ENDDO
    1790                     ENDDO
     1848                    ENDIF
    17911849                   
    17921850                CASE ( 'usm_iwghf_window' )
    17931851!--                 array of heat flux from indoor window ground (wall, roof, land)
    1794                     DO  m = 1, surf_usm_h%ns
    1795                        surf_usm_h%iwghf_eb_window_av(m) =                              &
    1796                                           surf_usm_h%iwghf_eb_window_av(m) +           &
    1797                                           surf_usm_h%iwghf_eb_window(m)
    1798                     ENDDO
    1799                     DO  l = 0, 3
     1852                    IF ( l == -1 ) THEN
     1853                       DO  m = 1, surf_usm_h%ns
     1854                          surf_usm_h%iwghf_eb_window_av(m) =                              &
     1855                                             surf_usm_h%iwghf_eb_window_av(m) +           &
     1856                                             surf_usm_h%iwghf_eb_window(m)
     1857                       ENDDO
     1858                    ELSE
    18001859                       DO  m = 1, surf_usm_v(l)%ns
    18011860                          surf_usm_v(l)%iwghf_eb_window_av(m) =                        &
     
    18031862                                          surf_usm_v(l)%iwghf_eb_window(m)
    18041863                       ENDDO
    1805                     ENDDO
     1864                    ENDIF
    18061865                   
    18071866                CASE ( 'usm_t_surf' )
    18081867!--                 surface temperature for surfaces
    1809                     DO  m = 1, surf_usm_h%ns
    1810                        surf_usm_h%t_surf_av(m) =                               &
    1811                                           surf_usm_h%t_surf_av(m) +            &
    1812                                           t_surf_h(m)
    1813                     ENDDO
    1814                     DO  l = 0, 3
     1868                    IF ( l == -1 ) THEN
     1869                       DO  m = 1, surf_usm_h%ns
     1870                          surf_usm_h%t_surf_av(m) =                               &
     1871                                             surf_usm_h%t_surf_av(m) +            &
     1872                                             t_surf_h(m)
     1873                       ENDDO
     1874                    ELSE
    18151875                       DO  m = 1, surf_usm_v(l)%ns
    18161876                          surf_usm_v(l)%t_surf_av(m) =                         &
     
    18181878                                          t_surf_v(l)%t(m)
    18191879                       ENDDO
    1820                     ENDDO
     1880                    ENDIF
    18211881                   
    18221882                CASE ( 'usm_t_surf_window' )
    18231883!--                 surface temperature for window surfaces
    1824                     DO  m = 1, surf_usm_h%ns
    1825                        surf_usm_h%t_surf_window_av(m) =                               &
    1826                                           surf_usm_h%t_surf_window_av(m) +            &
    1827                                           t_surf_window_h(m)
    1828                     ENDDO
    1829                     DO  l = 0, 3
     1884                    IF ( l == -1 ) THEN
     1885                       DO  m = 1, surf_usm_h%ns
     1886                          surf_usm_h%t_surf_window_av(m) =                               &
     1887                                             surf_usm_h%t_surf_window_av(m) +            &
     1888                                             t_surf_window_h(m)
     1889                       ENDDO
     1890                    ELSE
    18301891                       DO  m = 1, surf_usm_v(l)%ns
    18311892                          surf_usm_v(l)%t_surf_window_av(m) =                         &
     
    18331894                                          t_surf_window_v(l)%t(m)
    18341895                       ENDDO
    1835                     ENDDO
     1896                    ENDIF
    18361897                   
    18371898                CASE ( 'usm_t_surf_green' )
    18381899!--                 surface temperature for green surfaces
    1839                     DO  m = 1, surf_usm_h%ns
    1840                        surf_usm_h%t_surf_green_av(m) =                               &
    1841                                           surf_usm_h%t_surf_green_av(m) +            &
    1842                                           t_surf_green_h(m)
    1843                     ENDDO
    1844                     DO  l = 0, 3
     1900                    IF ( l == -1 ) THEN
     1901                       DO  m = 1, surf_usm_h%ns
     1902                          surf_usm_h%t_surf_green_av(m) =                               &
     1903                                             surf_usm_h%t_surf_green_av(m) +            &
     1904                                             t_surf_green_h(m)
     1905                       ENDDO
     1906                    ELSE
    18451907                       DO  m = 1, surf_usm_v(l)%ns
    18461908                          surf_usm_v(l)%t_surf_green_av(m) =                         &
     
    18481910                                          t_surf_green_v(l)%t(m)
    18491911                       ENDDO
    1850                     ENDDO
     1912                    ENDIF
    18511913               
    18521914                CASE ( 'usm_t_surf_10cm' )
    18531915!--                 near surface temperature for whole surfaces
    1854                     DO  m = 1, surf_usm_h%ns
    1855                        surf_usm_h%t_surf_10cm_av(m) =                               &
    1856                                           surf_usm_h%t_surf_10cm_av(m) +            &
    1857                                           t_surf_10cm_h(m)
    1858                     ENDDO
    1859                     DO  l = 0, 3
     1916                    IF ( l == -1 ) THEN
     1917                       DO  m = 1, surf_usm_h%ns
     1918                          surf_usm_h%t_surf_10cm_av(m) =                               &
     1919                                             surf_usm_h%t_surf_10cm_av(m) +            &
     1920                                             t_surf_10cm_h(m)
     1921                       ENDDO
     1922                    ELSE
    18601923                       DO  m = 1, surf_usm_v(l)%ns
    18611924                          surf_usm_v(l)%t_surf_10cm_av(m) =                         &
     
    18631926                                          t_surf_10cm_v(l)%t(m)
    18641927                       ENDDO
    1865                     ENDDO
     1928                    ENDIF
    18661929
    18671930                   
    18681931                CASE ( 'usm_t_wall' )
    18691932!--                 wall temperature for  iwl layer of walls and land
    1870                     DO  m = 1, surf_usm_h%ns
    1871                        surf_usm_h%t_wall_av(iwl,m) =                           &
    1872                                           surf_usm_h%t_wall_av(iwl,m) +        &
    1873                                           t_wall_h(iwl,m)
    1874                     ENDDO
    1875                     DO  l = 0, 3
     1933                    IF ( l == -1 ) THEN
     1934                       DO  m = 1, surf_usm_h%ns
     1935                          surf_usm_h%t_wall_av(iwl,m) =                           &
     1936                                             surf_usm_h%t_wall_av(iwl,m) +        &
     1937                                             t_wall_h(iwl,m)
     1938                       ENDDO
     1939                    ELSE
    18761940                       DO  m = 1, surf_usm_v(l)%ns
    18771941                          surf_usm_v(l)%t_wall_av(iwl,m) =                     &
     
    18791943                                          t_wall_v(l)%t(iwl,m)
    18801944                       ENDDO
    1881                     ENDDO
     1945                    ENDIF
    18821946                   
    18831947                CASE ( 'usm_t_window' )
    18841948!--                 window temperature for  iwl layer of walls and land
    1885   &