Changeset 2920


Ignore:
Timestamp:
Mar 22, 2018 11:22:01 AM (4 years ago)
Author:
kanani
Message:

Optimize SVF calculation, clean-up, bugfixes

Location:
palm/trunk/SOURCE
Files:
5 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/init_3d_model.f90

    r2906 r2920  
    2525! -----------------
    2626! $Id$
     27! Add call for precalculating apparent solar positions (moh.hefny)
     28!
     29! 2906 2018-03-19 08:56:40Z Giersch
    2730! The variables read/write_svf_on_init have been removed. Instead ENVIRONMENT
    2831! variables read/write_svf have been introduced. Location_message has been
     
    494497               radiation_calc_svf, radiation_write_svf,                        &
    495498               radiation_interaction, radiation_interactions,                  &
    496                radiation_interaction_init, radiation_read_svf
     499               radiation_interaction_init, radiation_read_svf,                 &
     500               radiation_presimulate_solar_pos
    497501   
    498502    USE random_function_mod
     
    23492353!
    23502354!--    If required, initialize radiation interactions between surfaces
    2351 !--    via sky-view factors. This must be done befoe radiation is initialized.
     2355!--    via sky-view factors. This must be done before radiation is initialized.
    23522356       IF ( radiation_interactions )  CALL radiation_interaction_init
    23532357
     
    23572361       CALL radiation_init
    23582362       CALL location_message( 'finished', .TRUE. )
     2363
     2364!
     2365!--    Find all discretized apparent solar positions for radiation interaction.
     2366!--    This must be done after radiation_init.
     2367       IF ( radiation_interactions )  CALL radiation_presimulate_solar_pos
    23592368
    23602369!
  • palm/trunk/SOURCE/plant_canopy_model_mod.f90

    r2892 r2920  
    2525! -----------------
    2626! $Id$
     27! Move usm_lad_rma and prototype_lad to radiation_model (moh.hefny)
     28!
     29! 2892 2018-03-14 15:06:29Z suehring
    2730! Bugfix, read separate ASCII LAD files for parent and child model.
    2831!
     
    190193
    191194    LOGICAL ::  calc_beta_lad_profile = .FALSE. !< switch for calc. of lad from beta func.
    192     LOGICAL ::  usm_lad_rma = .TRUE.            !< use MPI RMA to access LAD for raytracing (instead of global array)
    193195
    194196    REAL(wp) ::  alpha_lad = 9999999.9_wp         !< coefficient for lad calculation
     
    205207    REAL(wp) ::  lsc = 0.0_wp                     !< leaf surface concentration
    206208    REAL(wp) ::  lsec = 0.0_wp                    !< leaf scalar exchange coeff.
    207     REAL(wp) ::  prototype_lad                    !< prototype leaf area density for computing effective optical depth
    208209
    209210    REAL(wp) ::  lad_vertical_gradient(10) = 0.0_wp              !< lad gradient
     
    230231!-- Public variables and constants
    231232    PUBLIC pc_heating_rate, canopy_mode, cthf, dt_plant_canopy, lad, lad_s,   &
    232            pch_index, prototype_lad, usm_lad_rma
     233           pch_index
    233234           
    234235
  • palm/trunk/SOURCE/radiation_model_mod.f90

    r2906 r2920  
    1515! PALM. If not, see <http://www.gnu.org/licenses/>.
    1616!
     17! Copyright 2015-2018 Czech Technical University in Prague
     18! Copyright 2015-2018 Institute of Computer Science of the
     19!                     Czech Academy of Sciences, Prague
    1720! Copyright 1997-2018 Leibniz Universitaet Hannover
    1821!------------------------------------------------------------------------------!
     
    2528! -----------------
    2629! $Id$
     30! - Bugfix: Initialize pcbl array (=-1)
     31! moh.hefny:
     32! - Use precalculated apparent solar positions for direct irradiance
     33! - Cumulative commit for radiation changes - merged RTM version 2.0:
     34! - New version of radiation interaction
     35! - Added new 2D raytracing process using whole vertical column at once (e.g.
     36!   memory efficiency & much less RMA)
     37! - Removed virtual surfaces
     38! - Separate processing of direct and diffuse solar radiation, new discreti
     39!   zation by azimuth and elevation angles
     40! - Diffuse radiation processed cumulatively using sky view factor
     41! - Enabled limiting of number of view factors between real surfaces, thus
     42!   greatly enhancing scalability
     43! - Minor bugfixes and enhancements
     44! - Fixing bugs from moving radiation interaction from urban_surface_mod
     45!
     46!
     47! 2906 2018-03-19 08:56:40Z Giersch
    2748! NAMELIST paramter read/write_svf_on_init have been removed, functions
    2849! check_open and close_file are used now for opening/closing files related to
     
    97118!
    98119! 2544 2017-10-13 18:09:32Z maronga
    99 ! Moved date and time quantitis to separate module date_and_time_mod 
     120! Moved date and time quantitis to separate module date_and_time_mod
    100121!
    101122! 2512 2017-10-04 08:26:59Z raasch
     
    262283               latitude, longitude, large_scale_forcing, lsf_surf,             &
    263284               message_string, microphysics_morrison, plant_canopy, pt_surface,&
    264                rho_surface, surface_pressure, time_since_reference_point
     285               rho_surface, surface_pressure, time_since_reference_point,      &
     286               urban_surface, land_surface, end_time, spinup_time, dt_spinup
    265287
    266288    USE cpulog,                                                                &
     
    272294    USE date_and_time_mod,                                                     &
    273295        ONLY:  calc_date_and_time, d_hours_day, d_seconds_hour, day_of_year,   &
    274                time_utc
     296               d_seconds_year, day_of_year_init, time_utc_init, time_utc
    275297
    276298    USE indices,                                                               &
     
    294316
    295317    USE plant_canopy_model_mod,                                                &
    296         ONLY:  pc_heating_rate, lad_s, usm_lad_rma
     318        ONLY:  pc_heating_rate, lad_s
    297319
    298320    USE pegrid
     
    380402                sun_direction = .FALSE.,              & !< flag parameter indicating whether solar direction shall be calculated
    381403                average_radiation = .FALSE.,          & !< flag to set the calculation of radiation averaging for the domain
    382                 atm_surfaces = .FALSE.,               & !< flag parameter indicating wheather surfaces of atmospheric cells will be considered in calculating SVF
    383                 radiation_interactions = .TRUE.,      & !< flag to control if radiation interactions via sky-view factors shall be considered
     404                radiation_interactions = .FALSE.,     & !< flag to control if radiation interactions via sky-view factors shall be considered
    384405                surf_reflections = .TRUE.               !< flag to switch the calculation of radiation interaction between surfaces.
    385406                                                        !< When it switched off, only the effect of buildings and trees shadow will
     
    581602    INTEGER(iwp), PARAMETER                        ::  ix = 4                             !< position of i-index in surfl and surf
    582603
    583     INTEGER(iwp), PARAMETER                        ::  nsurf_type = 21                    !< number of surf types incl. phys.(land+urban) & (atm.,sky,boundary) surfaces - 1
    584 
    585     INTEGER(iwp), PARAMETER                        ::  iup_u    = 0                       !< 0 - index of urban ubward surface (ground or roof)
     604    INTEGER(iwp), PARAMETER                        ::  nsurf_type = 16                    !< number of surf types incl. phys.(land+urban) & (atm.,sky,boundary) surfaces - 1
     605
     606    INTEGER(iwp), PARAMETER                        ::  iup_u    = 0                       !< 0 - index of urban upward surface (ground or roof)
    586607    INTEGER(iwp), PARAMETER                        ::  idown_u  = 1                       !< 1 - index of urban downward surface (overhanging)
    587608    INTEGER(iwp), PARAMETER                        ::  inorth_u = 2                       !< 2 - index of urban northward facing wall
     
    590611    INTEGER(iwp), PARAMETER                        ::  iwest_u  = 5                       !< 5 - index of urban westward facing wall
    591612
    592     INTEGER(iwp), PARAMETER                        ::  iup_l    = 6                       !< 6 - index of land ubward surface (ground or roof)
     613    INTEGER(iwp), PARAMETER                        ::  iup_l    = 6                       !< 6 - index of land upward surface (ground or roof)
    593614    INTEGER(iwp), PARAMETER                        ::  inorth_l = 7                       !< 7 - index of land northward facing wall
    594615    INTEGER(iwp), PARAMETER                        ::  isouth_l = 8                       !< 8 - index of land southward facing wall
     
    603624    INTEGER(iwp), PARAMETER                        ::  iwest_a  = 16                      !< 16- index of atm. cell westward facing virtual surface
    604625
    605     INTEGER(iwp), PARAMETER                        ::  isky     = 17                      !< 17 - index of top border of the urban surface layer ("urban sky")
    606     INTEGER(iwp), PARAMETER                        ::  inorth_b = 18                      !< 18 - index of free north border of the domain (south facing)
    607     INTEGER(iwp), PARAMETER                        ::  isouth_b = 19                      !< 19 - index of north south border of the domain (north facing)
    608     INTEGER(iwp), PARAMETER                        ::  ieast_b  = 20                      !< 20 - index of east border of the domain (west facing)
    609     INTEGER(iwp), PARAMETER                        ::  iwest_b  = 21                      !< 21 - index of wast border of the domain (east facing)
    610 
    611     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, 0, 0,0,-1,1/)   !< surface normal direction x indices
    612     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, 0,-1,1, 0,0/)   !< surface normal direction y indices
    613     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,-1, 0,0, 0,0/)   !< surface normal direction z indices
     626    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
     627    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
     628    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
    614629                                                                                          !< parameter but set in the code
    615630
    616631
    617632!-- indices and sizes of urban and land surface models
    618     INTEGER(iwp)                                   ::  nskys            !< number of sky surfaces in local processor
    619     INTEGER(iwp)                                   ::  startland        !< start index of block of land and roof surfaces!-- block variables needed for calculation of the plant canopy model inside the urban surface model
    620     INTEGER(iwp)                                   ::  endland          !< end index of block of land and roof surfaces    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  pct              !< top layer of the plant canopy
    621     INTEGER(iwp)                                   ::  nlands           !< number of land and roof surfaces in local processor    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  pch              !< heights of the plant canopy
    622     INTEGER(iwp)                                   ::  startwall        !< start index of block of wall surfaces    INTEGER(iwp)                                   ::  npcbl            !< number of the plant canopy gridboxes in local processor
    623     INTEGER(iwp)                                   ::  endwall          !< end index of block of wall surfaces    INTEGER(wp), DIMENSION(:,:), ALLOCATABLE       ::  pcbl             !< k,j,i coordinates of l-th local plant canopy box pcbl[:,l] = [k, j,
    624     INTEGER(iwp)                                   ::  nwalls           !< number of wall surfaces in local processor    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinsw          !< array of absorbed sw radiation for local plant canopy box
    625     INTEGER(iwp)                                   ::  nborder          !< number of border surfaces in local processor    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinlw          !< array of absorbed lw radiation for local plant canopy box
    626 
     633    INTEGER(iwp)                                   ::  startland        !< start index of block of land and roof surfaces
     634    INTEGER(iwp)                                   ::  endland          !< end index of block of land and roof surfaces
     635    INTEGER(iwp)                                   ::  nlands           !< number of land and roof surfaces in local processor
     636    INTEGER(iwp)                                   ::  startwall        !< start index of block of wall surfaces
     637    INTEGER(iwp)                                   ::  endwall          !< end index of block of wall surfaces
     638    INTEGER(iwp)                                   ::  nwalls           !< number of wall surfaces in local processor
    627639
    628640!-- indices and sizes of urban and land surface models
     
    631643    INTEGER(iwp)                                   ::  nsurfl           !< number of all surfaces in local processor
    632644    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  nsurfs           !< array of number of all surfaces in individual processors
    633     INTEGER(iwp)                                   ::  startsky         !< start index of block of sky
    634     INTEGER(iwp)                                   ::  endsky           !< end index of block of sky
    635     INTEGER(iwp)                                   ::  startenergy      !< start index of block of real surfaces (land, walls and roofs)
    636     INTEGER(iwp)                                   ::  endenergy        !< end index of block of real surfaces (land, walls and roofs)
    637     INTEGER(iwp)                                   ::  nenergy          !< number of real surfaces in local processor
    638645    INTEGER(iwp)                                   ::  nsurf            !< global number of surfaces in index array of surfaces (nsurf = proc nsurfs)
    639     INTEGER(iwp)                                   ::  startborder      !< start index of block of border
    640     INTEGER(iwp)                                   ::  endborder        !< end index of block of border
    641646    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  surfstart        !< starts of blocks of surfaces for individual processors in array surf
    642647                                                                        !< respective block for particular processor is surfstart[iproc]+1 : surfstart[iproc+1]
     
    648653    INTEGER(wp), DIMENSION(:,:), ALLOCATABLE       ::  pcbl             !< k,j,i coordinates of l-th local plant canopy box pcbl[:,l] = [k, j, i]
    649654    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinsw          !< array of absorbed sw radiation for local plant canopy box
     655    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdir       !< array of absorbed direct sw radiation for local plant canopy box
     656    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinswdif       !< array of absorbed diffusion sw radiation for local plant canopy box
    650657    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  pcbinlw          !< array of absorbed lw radiation for local plant canopy box
    651658
     
    653660    LOGICAL                                        ::  split_diffusion_radiation = .TRUE. !< split direct and diffusion dw radiation
    654661                                                                                          !< (.F. in case the radiation model already does it)   
    655     LOGICAL                                        ::  energy_balance_surf_h = .TRUE.     !< flag parameter indicating wheather the energy balance is calculated for horizontal surfaces
    656     LOGICAL                                        ::  energy_balance_surf_v = .TRUE.     !< flag parameter indicating wheather the energy balance is calculated for vertical surfaces
     662    LOGICAL                                        ::  rma_lad_raytrace = .FALSE.         !< use MPI RMA to access LAD for raytracing (instead of global array)
    657663    LOGICAL                                        ::  mrt_factors = .FALSE.              !< whether to generate MRT factor files during init
    658664    INTEGER(iwp)                                   ::  nrefsteps = 0                      !< number of reflection steps to perform
     
    660666    INTEGER(iwp), PARAMETER                        ::  svf_code_len = 15                  !< length of code for verification of the end of svf file
    661667    CHARACTER(svf_code_len), PARAMETER             ::  svf_code = '*** end svf ***'       !< code for verification of the end of svf file
    662     INTEGER(iwp), PARAMETER                        ::  usm_version_len = 10               !< length of identification string of usm version
    663     CHARACTER(usm_version_len), PARAMETER          ::  usm_version = 'USM v. 1.0'         !< identification of version of binary svf and restart files
     668    INTEGER(iwp), PARAMETER                        ::  rad_version_len = 10               !< length of identification string of rad version
     669    CHARACTER(rad_version_len), PARAMETER          ::  rad_version = 'RAD v. 1.0'         !< identification of version of binary svf and restart files
     670    INTEGER(iwp)                                   ::  raytrace_discrete_elevs = 40       !< number of discretization steps for elevation (nadir to zenith)
     671    INTEGER(iwp)                                   ::  raytrace_discrete_azims = 80       !< number of discretization steps for azimuth (out of 360 degrees)
     672    REAL(wp)                                       ::  max_raytracing_dist = -999.0_wp    !< maximum distance for raytracing (in metres)
     673    REAL(wp)                                       ::  min_irrf_value = 1e-6_wp           !< minimum potential irradiance factor value for raytracing
     674    REAL(wp), DIMENSION(1:30)                      ::  svfnorm_report_thresh = 1e21_wp    !< thresholds of SVF normalization values to report
     675    INTEGER(iwp)                                   ::  svfnorm_report_num                 !< number of SVF normalization thresholds to report
    664676
    665677!-- radiation related arrays to be used in radiation_interaction routine
     
    697709    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfins          !< array of sw radiation falling to local surface after i-th reflection
    698710    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinl          !< array of lw radiation for local surface after i-th reflection
    699    
    700                                                                         !< Inward radiation is also valid for virtual surfaces (radiation leaving domain)
     711
     712    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  skyvf            !< array of sky view factor for each local surface
     713    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  skyvft           !< array of sky view factor including transparency for each local surface
     714    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  dsitrans         !< dsidir[isvfl,i] = path transmittance of i-th
     715                                                                        !< direction of direct solar irradiance per target surface
     716    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  dsitransc        !< dtto per plant canopy box
     717    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  dsidir           !< dsidir[:,i] = unit vector of i-th
     718                                                                        !< direction of direct solar irradiance
     719    INTEGER(iwp)                                   ::  ndsidir          !< number of apparent solar directions used
     720    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  dsidir_rev       !< dsidir_rev[ielev,iazim] = i for dsidir or -1 if not present
     721
    701722    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinsw         !< array of sw radiation falling to local surface including radiation from reflections
    702723    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfinlw         !< array of lw radiation falling to local surface including radiation from reflections
     
    713734    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfoutlw        !< array of total lw radiation outgoing from nonvirtual surfaces surfaces after all reflection
    714735    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  surfhf           !< array of total radiation flux incoming to minus outgoing from local surface
    715     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  rad_net_l        !< local copy of rad_net (net radiation at surface)
    716736
    717737!-- block variables needed for calculation of the plant canopy model inside the urban surface model
    718738    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  csfsurf          !< csfsurf[:,icsf] = index of target surface and csf grid index for csf[icsf]
    719739    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  csf              !< array of plant canopy sink fators + direct irradiation factors (transparency)
    720     REAL(wp), DIMENSION(:,:,:), POINTER            ::  usm_lad          !< subset of lad_s within urban surface, transformed to plain Z coordinate
    721     REAL(wp), DIMENSION(:), POINTER                ::  usm_lad_g        !< usm_lad globalized (used to avoid MPI RMA calls in raytracing)
     740    REAL(wp), DIMENSION(:,:,:), POINTER            ::  sub_lad          !< subset of lad_s within urban surface, transformed to plain Z coordinate
     741    REAL(wp), DIMENSION(:), POINTER                ::  sub_lad_g        !< sub_lad globalized (used to avoid MPI RMA calls in raytracing)
     742    REAL(wp)                                       ::  prototype_lad    !< prototype leaf area density for computing effective optical depth
    722743    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  nzterr, plantt   !< temporary global arrays for raytracing
     744    INTEGER(iwp)                                   ::  plantt_max
    723745
    724746!-- arrays and variables for calculation of svf and csf
     
    736758    INTEGER(iwp)                                   ::  ncsfl            !< no. of csf in local processor
    737759                                                                        !< needed only during calc_svf but must be here because it is
    738                                                                         !< shared between subroutines usm_calc_svf and usm_raytrace
     760                                                                        !< shared between subroutines calc_svf and raytrace
    739761    INTEGER(iwp), DIMENSION(:,:,:), ALLOCATABLE    ::  gridpcbl         !< index of local pcb[k,j,i]
    740762
     
    748770                  DIMENSION(:), ALLOCATABLE        ::  lad_disp         !< array of displaycements of lad in local array of proc lad_ip
    749771#endif
    750     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  lad_s_ray        !< array of received lad_s for appropriate gridboxes crossed by ray
     772    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  lad_s_ray        !< array of received lad_s for appropriate gridboxes crossed by ray
     773    INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE      ::  rt2_track
     774    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  rt2_track_lad
     775    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  rt2_track_dist
     776    REAL(wp), DIMENSION(:), ALLOCATABLE            ::  rt2_dist
     777
    751778
    752779
     
    839866       MODULE PROCEDURE radiation_interaction_init
    840867    END INTERFACE radiation_interaction_init
     868 
     869    INTERFACE radiation_presimulate_solar_pos
     870       MODULE PROCEDURE radiation_presimulate_solar_pos
     871    END INTERFACE radiation_presimulate_solar_pos
    841872
    842873    INTERFACE radiation_radflux_gridbox
     
    872903           radiation_radflux_gridbox, radiation_calc_svf, radiation_write_svf, &
    873904           radiation_interaction, radiation_interaction_init,                  &
    874            radiation_read_svf
     905           radiation_read_svf, radiation_presimulate_solar_pos
    875906           
    876907
     
    888919           zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon,       &
    889920           split_diffusion_radiation,                                          &
    890            energy_balance_surf_h, energy_balance_surf_v,                       &
    891921           nrefsteps, mrt_factors, dist_max_svf, nsvfl, svf,                   &
    892922           svfsurf, surfinsw, surfinlw, surfins, surfinl, surfinswdir,         &
    893923           surfinswdif, surfoutsw, surfoutlw, surfinlwdif, rad_sw_in_dir,      &
    894924           rad_sw_in_diff, rad_lw_in_diff, surfouts, surfoutl, surfoutsl,      &
    895            surfoutll, idir, jdir, kdir, id, iz, iy, ix, isky, nenergy, nsurfs, &
    896            surfstart, surf, surfl, nsurfl, pcbinsw, pcbinlw, pcbl, npcbl,      &
    897            startenergy, endenergy, iup_u, inorth_u, isouth_u, ieast_u, iwest_u,&
    898            iup_l, inorth_l, isouth_l, ieast_l, iwest_l, startsky, endsky,      &
    899            startborder, endborder, nsurf_type, nzub, nzut, inorth_b,idown_a,   &
    900            isouth_b, ieast_b, iwest_b, nzu, pch, nsurf, iup_a, inorth_a,       &
    901            isouth_a, ieast_a, iwest_a, idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct, &
    902            radiation_interactions, startwall, startland, endland, endwall
    903 
    904 
     925           surfoutll, idir, jdir, kdir, id, iz, iy, ix, nsurfs, surfstart,     &
     926           surf, surfl, nsurfl, pcbinswdir, pcbinswdif, pcbinsw, pcbinlw,      &
     927           pcbl, npcbl, iup_u, inorth_u, isouth_u, ieast_u, iwest_u,           &
     928           iup_l, inorth_l, isouth_l, ieast_l, iwest_l,                        &
     929           nsurf_type, nzub, nzut, nzu, pch, nsurf,                            &
     930           iup_a, idown_a, inorth_a, isouth_a, ieast_a, iwest_a,               &
     931           idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct,                             &
     932           radiation_interactions, startwall, startland, endland, endwall,     &
     933           skyvf, skyvft
    905934
    906935#if defined ( __rrtmg )
     
    12461275!
    12471276!--    Radiation interactions
    1248        IF ( urban_surface .AND.  .NOT. radiation_interactions )  THEN
    1249           message_string = 'radiation_interactions = .T. is required '//       &
    1250                            'when using the urban surface model'
     1277       IF ( nrefsteps < 1  .AND.  radiation_interactions )  THEN
     1278          message_string = 'nrefsteps must be > 0 when using LSM/USM to' //    &
     1279                           'account for surface outgoing SW flux.'       //    &
     1280                           'You may set surf_reflections = .FALSE. to '  //    &
     1281                           'diable surface reflections instead.'
    12511282          CALL message( 'check_parameters', 'PA0999', 1, 2, 0, 6, 0 )
    12521283       ENDIF
     1284
     1285!
     1286!--    Incialize svf normalization reporting histogram
     1287       svfnorm_report_num = 1
     1288       DO WHILE ( svfnorm_report_thresh(svfnorm_report_num) < 1e20_wp          &
     1289                   .AND. svfnorm_report_num <= 30 )
     1290          svfnorm_report_num = svfnorm_report_num + 1
     1291       ENDDO
     1292       svfnorm_report_num = svfnorm_report_num - 1
     1293
    12531294
    12541295 
     
    14981539       IF ( radiation_scheme == 'clear-sky'  .OR.                              &
    14991540            radiation_scheme == 'constant')  THEN
     1541
     1542
     1543!
     1544!--       Allocate arrays for incoming/outgoing short/longwave radiation
     1545          IF ( .NOT. ALLOCATED ( rad_sw_in ) )  THEN
     1546             ALLOCATE ( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
     1547          ENDIF
     1548          IF ( .NOT. ALLOCATED ( rad_sw_out ) )  THEN
     1549             ALLOCATE ( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
     1550          ENDIF
     1551
     1552          IF ( .NOT. ALLOCATED ( rad_lw_in ) )  THEN
     1553             ALLOCATE ( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
     1554          ENDIF
     1555          IF ( .NOT. ALLOCATED ( rad_lw_out ) )  THEN
     1556             ALLOCATE ( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
     1557          ENDIF
     1558
    15001559!
    15011560!--       Allocate average arrays for incoming/outgoing short/longwave radiation
     
    15151574!
    15161575!--       Allocate arrays for broadband albedo, and level 1 initialization
    1517 !--       via namelist paramter.
    1518           IF ( .NOT. ALLOCATED(surf_def_h(0)%albedo) )                         &
     1576!--       via namelist paramter, unless already allocated.
     1577          IF ( .NOT. ALLOCATED(surf_def_h(0)%albedo) )  THEN
    15191578             ALLOCATE( surf_def_h(0)%albedo(0:0,1:surf_def_h(0)%ns) )
    1520           IF ( .NOT. ALLOCATED(surf_lsm_h%albedo) )                            &
     1579             surf_def_h(0)%albedo = albedo
     1580          ENDIF
     1581          IF ( .NOT. ALLOCATED(surf_lsm_h%albedo) )  THEN
    15211582             ALLOCATE( surf_lsm_h%albedo(0:2,1:surf_lsm_h%ns)     )
    1522           IF ( .NOT. ALLOCATED(surf_usm_h%albedo) )                            &
     1583             surf_lsm_h%albedo    = albedo
     1584          ENDIF
     1585          IF ( .NOT. ALLOCATED(surf_usm_h%albedo) )  THEN
    15231586             ALLOCATE( surf_usm_h%albedo(0:2,1:surf_usm_h%ns)     )
    1524 
    1525           surf_def_h(0)%albedo = albedo
    1526           surf_lsm_h%albedo    = albedo
    1527           surf_usm_h%albedo    = albedo
     1587             surf_usm_h%albedo    = albedo
     1588          ENDIF
     1589
    15281590          DO  l = 0, 3
    1529              IF ( .NOT. ALLOCATED( surf_def_v(l)%albedo ) )                    &
     1591             IF ( .NOT. ALLOCATED( surf_def_v(l)%albedo ) )  THEN
    15301592                ALLOCATE( surf_def_v(l)%albedo(0:0,1:surf_def_v(l)%ns) )
    1531              IF ( .NOT. ALLOCATED( surf_lsm_v(l)%albedo ) )                    &
     1593                surf_def_v(l)%albedo = albedo
     1594             ENDIF
     1595             IF ( .NOT. ALLOCATED( surf_lsm_v(l)%albedo ) )  THEN
    15321596                ALLOCATE( surf_lsm_v(l)%albedo(0:2,1:surf_lsm_v(l)%ns) )
    1533              IF ( .NOT. ALLOCATED( surf_usm_v(l)%albedo ) )                    &
     1597                surf_lsm_v(l)%albedo = albedo
     1598             ENDIF
     1599             IF ( .NOT. ALLOCATED( surf_usm_v(l)%albedo ) )  THEN
    15341600                ALLOCATE( surf_usm_v(l)%albedo(0:2,1:surf_usm_v(l)%ns) )
    1535 
    1536              surf_def_v(l)%albedo = albedo
    1537              surf_lsm_v(l)%albedo = albedo
    1538              surf_usm_v(l)%albedo = albedo
     1601                surf_usm_v(l)%albedo = albedo
     1602             ENDIF
    15391603          ENDDO
    15401604!
     
    25232587             ENDIF
    25242588
     2589!
     2590!--          Fill out values in radiation arrays
     2591             DO  m = 1, surf%ns
     2592                i = surf%i(m)
     2593                j = surf%j(m)
     2594                rad_sw_in(0,j,i) = surf%rad_sw_in(m)
     2595                rad_sw_out(0,j,i) = surf%rad_sw_out(m)
     2596                rad_lw_in(0,j,i) = surf%rad_lw_in(m)
     2597                rad_lw_out(0,j,i) = surf%rad_lw_out(m)
     2598             ENDDO
     2599
    25252600          END SUBROUTINE radiation_clearsky_surf
    25262601
     
    27112786
    27122787             ENDIF
     2788
     2789!
     2790!--          Fill out values in radiation arrays
     2791             DO  m = 1, surf%ns
     2792                i = surf%i(m)
     2793                j = surf%j(m)
     2794                rad_sw_in(0,j,i) = surf%rad_sw_in(m)
     2795                rad_sw_out(0,j,i) = surf%rad_sw_out(m)
     2796                rad_lw_in(0,j,i) = surf%rad_lw_in(m)
     2797                rad_lw_out(0,j,i) = surf%rad_lw_out(m)
     2798             ENDDO
    27132799
    27142800          END SUBROUTINE radiation_constant_surf
     
    28022888                                  sw_radiation, unscheduled_radiation_calls,   &
    28032889                                  split_diffusion_radiation,                   &
    2804                                   energy_balance_surf_h,                       &
    2805                                   energy_balance_surf_v,                       &
    2806                                   nrefsteps,                                   &
    2807                                   mrt_factors,                                 &
     2890                                  max_raytracing_dist, min_irrf_value,         &
     2891                                  nrefsteps, mrt_factors, rma_lad_raytrace,    &
    28082892                                  dist_max_svf,                                &
    28092893                                  average_radiation,                           &
    2810                                   radiation_interactions, atm_surfaces,        &
    2811                                   surf_reflections
     2894                                  surf_reflections, svfnorm_report_thresh
    28122895       
    28132896       line = ' '
     
    28292912!--    Set flag that indicates that the radiation model is switched on
    28302913       radiation = .TRUE.
     2914
     2915!--    Set radiation_interactions flag according to urban_ and land_surface flag
     2916       IF ( urban_surface  .OR.  land_surface ) radiation_interactions = .TRUE.
    28312917
    28322918 10    CONTINUE
     
    44724558 END SUBROUTINE radiation_tendency
    44734559
    4474 
    44754560!------------------------------------------------------------------------------!
    44764561! Description:
    44774562! ------------
    44784563!> This subroutine calculates interaction of the solar radiation
    4479 !> with urban and land surfaces and updates all surface heatfluxes, including
    4480 !> the vertual atmospheric cell faces. It calculates also the required parameters
    4481 !> for RRTMG lower BC.
    4482 !> 
     4564!> with urban and land surfaces and updates all surface heatfluxes.
     4565!> It calculates also the required parameters for RRTMG lower BC.
     4566!>
    44834567!> For more info. see Resler et al. 2017
    4484 !> 
     4568!>
     4569!> The new version 2.0 was radically rewriten, the discretization scheme
     4570!> has been changed. This new version significantly improves effectivity
     4571!> of the paralelization and the scalability of the model.
     4572!------------------------------------------------------------------------------!
     4573
     4574 SUBROUTINE radiation_interaction
     4575
     4576     IMPLICIT NONE
     4577
     4578     INTEGER(iwp)                      :: i, j, k, kk, is, js, d, ku, refstep, m, mm, l, ll
     4579     INTEGER(iwp)                      :: nzubl, nzutl, isurf, isurfsrc, isvf, icsf, ipcgb
     4580     INTEGER(iwp)                      :: isd                !< solar direction number
     4581     REAL(wp), DIMENSION(3,3)          :: mrot               !< grid rotation matrix (zyx)
     4582     REAL(wp), DIMENSION(3,0:nsurf_type):: vnorm             !< face direction normal vectors (zyx)
     4583     REAL(wp), DIMENSION(3)            :: sunorig            !< grid rotated solar direction unit vector (zyx)
     4584     REAL(wp), DIMENSION(3)            :: sunorig_grid       !< grid squashed solar direction unit vector (zyx)
     4585     REAL(wp), DIMENSION(0:nsurf_type) :: costheta           !< direct irradiance factor of solar angle
     4586     REAL(wp), DIMENSION(nzub:nzut)    :: pchf_prep          !< precalculated factor for canopy temp tendency
     4587     REAL(wp), PARAMETER               :: alpha = 0._wp      !< grid rotation (TODO: add to namelist or remove)
     4588     REAL(wp)                          :: pc_box_area, pc_abs_frac, pc_abs_eff
     4589     INTEGER(iwp)                      :: pc_box_dimshift    !< transform for best accuracy
     4590     INTEGER(iwp), DIMENSION(0:3)      :: reorder = (/ 1, 0, 3, 2 /)
     4591     REAL(wp), DIMENSION(0:nsurf_type) :: facearea
     4592     REAL(wp)                          :: pabsswl  = 0.0_wp  !< total absorbed SW radiation energy in local processor (W)
     4593     REAL(wp)                          :: pabssw   = 0.0_wp  !< total absorbed SW radiation energy in all processors (W)
     4594     REAL(wp)                          :: pabslwl  = 0.0_wp  !< total absorbed LW radiation energy in local processor (W)
     4595     REAL(wp)                          :: pabslw   = 0.0_wp  !< total absorbed LW radiation energy in all processors (W)
     4596     REAL(wp)                          :: pemitlwl = 0.0_wp  !< total emitted LW radiation energy in all processors (W)
     4597     REAL(wp)                          :: pemitlw  = 0.0_wp  !< total emitted LW radiation energy in all processors (W)
     4598     REAL(wp)                          :: pinswl   = 0.0_wp  !< total received SW radiation energy in local processor (W)
     4599     REAL(wp)                          :: pinsw    = 0.0_wp  !< total received SW radiation energy in all processor (W)
     4600     REAL(wp)                          :: pinlwl   = 0.0_wp  !< total received LW radiation energy in local processor (W)
     4601     REAL(wp)                          :: pinlw    = 0.0_wp  !< total received LW radiation energy in all processor (W)
     4602     REAL(wp)                          :: emiss_sum_surfl    !< sum of emissisivity of surfaces in local processor
     4603     REAL(wp)                          :: emiss_sum_surf     !< sum of emissisivity of surfaces in all processor
     4604     REAL(wp)                          :: area_surfl         !< total area of surfaces in local processor
     4605     REAL(wp)                          :: area_surf          !< total area of surfaces in all processor
     4606
     4607#if ! defined( __nopointer )
     4608     IF ( plant_canopy )  THEN
     4609         pchf_prep(:) = r_d * (hyp(nzub:nzut) / 100000.0_wp)**0.286_wp &
     4610                     / (cp * hyp(nzub:nzut) * dx*dy*dz) !< equals to 1 / (rho * c_p * Vbox * T)
     4611     ENDIF
     4612#endif
     4613     sun_direction = .TRUE.
     4614     CALL calc_zenith  !< required also for diffusion radiation
     4615
     4616!--     prepare rotated normal vectors and irradiance factor
     4617     vnorm(1,:) = kdir(:)
     4618     vnorm(2,:) = jdir(:)
     4619     vnorm(3,:) = idir(:)
     4620     mrot(1, :) = (/ 1._wp,  0._wp,      0._wp      /)
     4621     mrot(2, :) = (/ 0._wp,  COS(alpha), SIN(alpha) /)
     4622     mrot(3, :) = (/ 0._wp, -SIN(alpha), COS(alpha) /)
     4623     sunorig = (/ zenith(0), sun_dir_lat, sun_dir_lon /)
     4624     sunorig = MATMUL(mrot, sunorig)
     4625     DO d = 0, nsurf_type
     4626         costheta(d) = DOT_PRODUCT(sunorig, vnorm(:,d))
     4627     ENDDO
     4628
     4629     IF ( zenith(0) > 0 )  THEN
     4630!--         now we will "squash" the sunorig vector by grid box size in
     4631!--         each dimension, so that this new direction vector will allow us
     4632!--         to traverse the ray path within grid coordinates directly
     4633         sunorig_grid = (/ sunorig(1)/dz, sunorig(2)/dy, sunorig(3)/dx /)
     4634!--         sunorig_grid = sunorig_grid / norm2(sunorig_grid)
     4635         sunorig_grid = sunorig_grid / SQRT(SUM(sunorig_grid**2))
     4636
     4637         IF ( plant_canopy )  THEN
     4638!--            precompute effective box depth with prototype Leaf Area Density
     4639            pc_box_dimshift = MAXLOC(ABS(sunorig), 1) - 1
     4640            CALL box_absorb(CSHIFT((/dz,dy,dx/), pc_box_dimshift),      &
     4641                                60, prototype_lad,                          &
     4642                                CSHIFT(ABS(sunorig), pc_box_dimshift),      &
     4643                                pc_box_area, pc_abs_frac)
     4644            pc_box_area = pc_box_area * ABS(sunorig(pc_box_dimshift+1) / sunorig(1))
     4645            pc_abs_eff = LOG(1._wp - pc_abs_frac) / prototype_lad
     4646         ENDIF
     4647     ENDIF
     4648
     4649!--     split diffusion and direct part of the solar downward radiation
     4650!--     comming from radiation model and store it in 2D arrays
     4651!--     rad_sw_in_diff, rad_sw_in_dir and rad_lw_in_diff
     4652     IF ( split_diffusion_radiation )  THEN
     4653         CALL calc_diffusion_radiation
     4654     ELSE
     4655         rad_sw_in_diff = 0.0_wp
     4656         rad_sw_in_dir(:,:)  = rad_sw_in(0,:,:)
     4657         rad_lw_in_diff(:,:) = rad_lw_in(0,:,:)
     4658     ENDIF
     4659
     4660!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     4661!--     First pass: direct + diffuse irradiance + thermal
     4662!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     4663     surfinswdir   = 0._wp !nsurfl
     4664     surfins       = 0._wp !nsurfl
     4665     surfinl       = 0._wp !nsurfl
     4666     surfoutsl(:)  = 0.0_wp !start-end
     4667     surfoutll(:)  = 0.0_wp !start-end
     4668
     4669!--  Set up thermal radiation from surfaces
     4670!--  emiss_surf is defined only for surfaces for which energy balance is calculated
     4671!--  Workaround: reorder surface data type back on 1D array including all surfaces,
     4672!--  which implies to reorder horizontal and vertical surfaces
     4673!
     4674!--  Horizontal walls
     4675     mm = 1
     4676     DO  i = nxl, nxr
     4677        DO  j = nys, nyn
     4678!--           urban
     4679           DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
     4680              surfoutll(mm) = SUM ( surf_usm_h%frac(:,m) *                  &
     4681                                    surf_usm_h%emissivity(:,m) )            &
     4682                                  * sigma_sb                                &
     4683                                  * surf_usm_h%pt_surface(m)**4
     4684              albedo_surf(mm) = SUM ( surf_usm_h%frac(:,m) *                &
     4685                                      surf_usm_h%albedo(:,m) )
     4686              emiss_surf(mm)  = SUM ( surf_usm_h%frac(:,m) *                &
     4687                                      surf_usm_h%emissivity(:,m) )
     4688              mm = mm + 1
     4689           ENDDO
     4690!--           land
     4691           DO  m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
     4692              surfoutll(mm) = SUM ( surf_lsm_h%frac(:,m) *                  &
     4693                                    surf_lsm_h%emissivity(:,m) )            &
     4694                                  * sigma_sb                                &
     4695                                  * surf_lsm_h%pt_surface(m)**4
     4696              albedo_surf(mm) = SUM ( surf_lsm_h%frac(:,m) *                &
     4697                                      surf_lsm_h%albedo(:,m) )
     4698              emiss_surf(mm)  = SUM ( surf_lsm_h%frac(:,m) *                &
     4699                                      surf_lsm_h%emissivity(:,m) )
     4700              mm = mm + 1
     4701           ENDDO
     4702        ENDDO
     4703     ENDDO
     4704!
     4705!--     Vertical walls
     4706     DO  i = nxl, nxr
     4707        DO  j = nys, nyn
     4708           DO  ll = 0, 3
     4709              l = reorder(ll)
     4710!--              urban
     4711              DO  m = surf_usm_v(l)%start_index(j,i),                       &
     4712                      surf_usm_v(l)%end_index(j,i)
     4713                 surfoutll(mm) = SUM ( surf_usm_v(l)%frac(:,m) *            &
     4714                                       surf_usm_v(l)%emissivity(:,m) )      &
     4715                                  * sigma_sb                                &
     4716                                  * surf_usm_v(l)%pt_surface(m)**4
     4717                 albedo_surf(mm) = SUM ( surf_usm_v(l)%frac(:,m) *          &
     4718                                         surf_usm_v(l)%albedo(:,m) )
     4719                 emiss_surf(mm)  = SUM ( surf_usm_v(l)%frac(:,m) *          &
     4720                                         surf_usm_v(l)%emissivity(:,m) )
     4721                 mm = mm + 1
     4722              ENDDO
     4723!--              land
     4724              DO  m = surf_lsm_v(l)%start_index(j,i),                       &
     4725                      surf_lsm_v(l)%end_index(j,i)
     4726                 surfoutll(mm) = SUM ( surf_lsm_v(l)%frac(:,m) *            &
     4727                                       surf_lsm_v(l)%emissivity(:,m) )      &
     4728                                  * sigma_sb                                &
     4729                                  * surf_lsm_v(l)%pt_surface(m)**4
     4730                 albedo_surf(mm) = SUM ( surf_lsm_v(l)%frac(:,m) *          &
     4731                                         surf_lsm_v(l)%albedo(:,m) )
     4732                 emiss_surf(mm)  = SUM ( surf_lsm_v(l)%frac(:,m) *          &
     4733                                         surf_lsm_v(l)%emissivity(:,m) )
     4734                 mm = mm + 1
     4735              ENDDO
     4736           ENDDO
     4737        ENDDO
     4738     ENDDO
     4739
     4740#if defined( __parallel )
     4741!--     might be optimized and gather only values relevant for current processor
     4742     CALL MPI_AllGatherv(surfoutll, nsurfl, MPI_REAL, &
     4743                         surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr) !nsurf global
     4744#else
     4745     surfoutl(:) = surfoutll(:) !nsurf global
     4746#endif
     4747
     4748     DO isvf = 1, nsvfl
     4749         isurf = svfsurf(1, isvf)
     4750         k = surfl(iz, isurf)
     4751         j = surfl(iy, isurf)
     4752         i = surfl(ix, isurf)
     4753         isurfsrc = svfsurf(2, isvf)
     4754
     4755!--         for surface-to-surface factors we calculate thermal radiation in 1st pass
     4756         surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc)
     4757     ENDDO
     4758
     4759     !-- diffuse radiation using sky view factor, TODO: homogeneous rad_*w_in_diff because now it depends on no. of processors
     4760     surfinswdif(:) = rad_sw_in_diff(nyn,nxl) * skyvft(:)
     4761     surfinlwdif(:) = rad_lw_in_diff(nyn,nxl) * skyvf(:)
     4762
     4763     !-- direct radiation
     4764     IF ( zenith(0) > 0 )  THEN
     4765        !--Identify solar direction vector (discretized number) 1)
     4766        !--
     4767        j = FLOOR(ACOS(zenith(0)) / pi * raytrace_discrete_elevs)
     4768        i = MODULO(NINT(ATAN2(sun_dir_lon(0), sun_dir_lat(0))               &
     4769                        / (2._wp*pi) * raytrace_discrete_azims-.5_wp, iwp), &
     4770                   raytrace_discrete_azims)
     4771        isd = dsidir_rev(j, i)
     4772        DO isurf = 1, nsurfl
     4773           surfinswdir(isurf) = rad_sw_in_dir(nyn,nxl) * costheta(surfl(id, isurf)) * dsitrans(isurf, isd) / zenith(0)
     4774        ENDDO
     4775     ENDIF
     4776
     4777     IF ( plant_canopy )  THEN
     4778
     4779         pcbinswdir(:) = 0._wp
     4780         pcbinswdif(:) = 0._wp
     4781         pcbinlw(:) = 0._wp  !< will stay always 0 since we don't absorb lw anymore
     4782!
     4783!--         pcsf first pass
     4784         DO icsf = 1, ncsfl
     4785             ipcgb = csfsurf(1, icsf)
     4786             i = pcbl(ix,ipcgb)
     4787             j = pcbl(iy,ipcgb)
     4788             k = pcbl(iz,ipcgb)
     4789             isurfsrc = csfsurf(2, icsf)
     4790
     4791             IF ( isurfsrc == -1 )  THEN
     4792!--                 Diffuse rad from sky.
     4793                 pcbinswdif(ipcgb) = csf(1,icsf) * csf(2,icsf) * rad_sw_in_diff(j,i)
     4794
     4795                 !--Direct rad
     4796                 IF ( zenith(0) > 0 )  THEN
     4797                    !--Estimate directed box absorption
     4798                    pc_abs_frac = 1._wp - exp(pc_abs_eff * lad_s(k,j,i))
     4799
     4800                    !--isd has already been established, see 1)
     4801                    pcbinswdir(ipcgb) = rad_sw_in_dir(j, i) * pc_box_area &
     4802                                        * pc_abs_frac * dsitransc(ipcgb, isd)
     4803                 ENDIF
     4804
     4805                 EXIT ! only isurfsrc=-1 is processed here
     4806             ENDIF
     4807         ENDDO
     4808
     4809         pcbinsw(:) = pcbinswdir(:) + pcbinswdif(:)
     4810     ENDIF
     4811     surfins = surfinswdir + surfinswdif
     4812     surfinl = surfinl + surfinlwdif
     4813     surfinsw = surfins
     4814     surfinlw = surfinl
     4815     surfoutsw = 0.0_wp
     4816     surfoutlw = surfoutll
     4817!        surfhf = surfinsw + surfinlw - surfoutsw - surfoutlw
     4818
     4819!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     4820!--     Next passes - reflections
     4821!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     4822     DO refstep = 1, nrefsteps
     4823
     4824         surfoutsl = albedo_surf * surfins
     4825!--         for non-transparent surfaces, longwave albedo is 1 - emissivity
     4826         surfoutll = (1._wp - emiss_surf) * surfinl
     4827
     4828#if defined( __parallel )
     4829         CALL MPI_AllGatherv(surfoutsl, nsurfl, MPI_REAL, &
     4830             surfouts, nsurfs, surfstart, MPI_REAL, comm2d, ierr)
     4831         CALL MPI_AllGatherv(surfoutll, nsurfl, MPI_REAL, &
     4832             surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr)
     4833#else
     4834         surfouts = surfoutsl
     4835         surfoutl = surfoutll
     4836#endif
     4837
     4838!--         reset for next pass input
     4839         surfins = 0._wp
     4840         surfinl = 0._wp
     4841
     4842!--         reflected radiation
     4843         DO isvf = 1, nsvfl
     4844             isurf = svfsurf(1, isvf)
     4845             isurfsrc = svfsurf(2, isvf)
     4846             surfins(isurf) = surfins(isurf) + svf(1,isvf) * svf(2,isvf) * surfouts(isurfsrc)
     4847             surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc)
     4848         ENDDO
     4849
     4850!--         radiation absorbed by plant canopy
     4851         DO icsf = 1, ncsfl
     4852             ipcgb = csfsurf(1, icsf)
     4853             isurfsrc = csfsurf(2, icsf)
     4854             IF ( isurfsrc == -1 )  CYCLE ! sky->face only in 1st pass, not here
     4855
     4856             pcbinsw(ipcgb) = pcbinsw(ipcgb) + csf(1,icsf) * csf(2,icsf) * surfouts(isurfsrc)
     4857         ENDDO
     4858
     4859         surfinsw = surfinsw  + surfins
     4860         surfinlw = surfinlw  + surfinl
     4861         surfoutsw = surfoutsw + surfoutsl
     4862         surfoutlw = surfoutlw + surfoutll
     4863!            surfhf = surfinsw + surfinlw - surfoutsw - surfoutlw
     4864
     4865     ENDDO
     4866
     4867!--  push heat flux absorbed by plant canopy to respective 3D arrays
     4868     IF ( plant_canopy )  THEN
     4869         pc_heating_rate(:,:,:) = 0._wp
     4870         DO ipcgb = 1, npcbl
     4871                 
     4872             j = pcbl(iy, ipcgb)
     4873             i = pcbl(ix, ipcgb)
     4874             k = pcbl(iz, ipcgb)
     4875!
     4876!--             Following expression equals former kk = k - nzb_s_inner(j,i)
     4877             kk = k - get_topography_top_index_ji( j, i, 's' )  !- lad arrays are defined flat
     4878             pc_heating_rate(kk, j, i) = (pcbinsw(ipcgb) + pcbinlw(ipcgb)) &
     4879                 * pchf_prep(k) * pt(k, j, i) !-- = dT/dt
     4880         ENDDO
     4881     ENDIF
     4882!
     4883!--     Transfer radiation arrays required for energy balance to the respective data types
     4884     DO  i = 1, nsurfl
     4885        m  = surfl(5,i)
     4886!
     4887!--     (1) Urban surfaces
     4888!--     upward-facing
     4889        IF ( surfl(1,i) == iup_u )  THEN
     4890           surf_usm_h%rad_sw_in(m)  = surfinsw(i)
     4891           surf_usm_h%rad_sw_out(m) = surfoutsw(i)
     4892           surf_usm_h%rad_lw_in(m)  = surfinlw(i)
     4893           surf_usm_h%rad_lw_out(m) = surfoutlw(i)
     4894           surf_usm_h%rad_net(m)    = surfinsw(i) - surfoutsw(i) +          &
     4895                                      surfinlw(i) - surfoutlw(i)
     4896!
     4897!--     northward-facding
     4898        ELSEIF ( surfl(1,i) == inorth_u )  THEN
     4899           surf_usm_v(0)%rad_sw_in(m)  = surfinsw(i)
     4900           surf_usm_v(0)%rad_sw_out(m) = surfoutsw(i)
     4901           surf_usm_v(0)%rad_lw_in(m)  = surfinlw(i)
     4902           surf_usm_v(0)%rad_lw_out(m) = surfoutlw(i)
     4903           surf_usm_v(0)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
     4904                                         surfinlw(i) - surfoutlw(i)
     4905!
     4906!--     southward-facding
     4907        ELSEIF ( surfl(1,i) == isouth_u )  THEN
     4908           surf_usm_v(1)%rad_sw_in(m)  = surfinsw(i)
     4909           surf_usm_v(1)%rad_sw_out(m) = surfoutsw(i)
     4910           surf_usm_v(1)%rad_lw_in(m)  = surfinlw(i)
     4911           surf_usm_v(1)%rad_lw_out(m) = surfoutlw(i)
     4912           surf_usm_v(1)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
     4913                                         surfinlw(i) - surfoutlw(i)
     4914!
     4915!--     eastward-facing
     4916        ELSEIF ( surfl(1,i) == ieast_u )  THEN
     4917           surf_usm_v(2)%rad_sw_in(m)  = surfinsw(i)
     4918           surf_usm_v(2)%rad_sw_out(m) = surfoutsw(i)
     4919           surf_usm_v(2)%rad_lw_in(m)  = surfinlw(i)
     4920           surf_usm_v(2)%rad_lw_out(m) = surfoutlw(i)
     4921           surf_usm_v(2)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
     4922                                         surfinlw(i) - surfoutlw(i)
     4923!
     4924!--     westward-facding
     4925        ELSEIF ( surfl(1,i) == iwest_u )  THEN
     4926           surf_usm_v(3)%rad_sw_in(m)  = surfinsw(i)
     4927           surf_usm_v(3)%rad_sw_out(m) = surfoutsw(i)
     4928           surf_usm_v(3)%rad_lw_in(m)  = surfinlw(i)
     4929           surf_usm_v(3)%rad_lw_out(m) = surfoutlw(i)
     4930           surf_usm_v(3)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
     4931                                         surfinlw(i) - surfoutlw(i)
     4932!
     4933!--     (2) land surfaces
     4934!--     upward-facing
     4935        ELSEIF ( surfl(1,i) == iup_l )  THEN
     4936           surf_lsm_h%rad_sw_in(m)  = surfinsw(i)
     4937           surf_lsm_h%rad_sw_out(m) = surfoutsw(i)
     4938           surf_lsm_h%rad_lw_in(m)  = surfinlw(i)
     4939           surf_lsm_h%rad_lw_out(m) = surfoutlw(i)
     4940           surf_lsm_h%rad_net(m)    = surfinsw(i) - surfoutsw(i) +          &
     4941                                      surfinlw(i) - surfoutlw(i)
     4942!
     4943!--     northward-facding
     4944        ELSEIF ( surfl(1,i) == inorth_l )  THEN
     4945           surf_lsm_v(0)%rad_sw_in(m)  = surfinsw(i)
     4946           surf_lsm_v(0)%rad_sw_out(m) = surfoutsw(i)
     4947           surf_lsm_v(0)%rad_lw_in(m)  = surfinlw(i)
     4948           surf_lsm_v(0)%rad_lw_out(m) = surfoutlw(i)
     4949           surf_lsm_v(0)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
     4950                                         surfinlw(i) - surfoutlw(i)
     4951!
     4952!--     southward-facding
     4953        ELSEIF ( surfl(1,i) == isouth_l )  THEN
     4954           surf_lsm_v(1)%rad_sw_in(m)  = surfinsw(i)
     4955           surf_lsm_v(1)%rad_sw_out(m) = surfoutsw(i)
     4956           surf_lsm_v(1)%rad_lw_in(m)  = surfinlw(i)
     4957           surf_lsm_v(1)%rad_lw_out(m) = surfoutlw(i)
     4958           surf_lsm_v(1)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
     4959                                         surfinlw(i) - surfoutlw(i)
     4960!
     4961!--     eastward-facing
     4962        ELSEIF ( surfl(1,i) == ieast_l )  THEN
     4963           surf_lsm_v(2)%rad_sw_in(m)  = surfinsw(i)
     4964           surf_lsm_v(2)%rad_sw_out(m) = surfoutsw(i)
     4965           surf_lsm_v(2)%rad_lw_in(m)  = surfinlw(i)
     4966           surf_lsm_v(2)%rad_lw_out(m) = surfoutlw(i)
     4967           surf_lsm_v(2)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
     4968                                         surfinlw(i) - surfoutlw(i)
     4969!
     4970!--     westward-facing
     4971        ELSEIF ( surfl(1,i) == iwest_l )  THEN
     4972           surf_lsm_v(3)%rad_sw_in(m)  = surfinsw(i)
     4973           surf_lsm_v(3)%rad_sw_out(m) = surfoutsw(i)
     4974           surf_lsm_v(3)%rad_lw_in(m)  = surfinlw(i)
     4975           surf_lsm_v(3)%rad_lw_out(m) = surfoutlw(i)
     4976           surf_lsm_v(3)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
     4977                                         surfinlw(i) - surfoutlw(i)
     4978        ENDIF
     4979
     4980     ENDDO
     4981
     4982     DO  m = 1, surf_usm_h%ns
     4983        surf_usm_h%surfhf(m) = surf_usm_h%rad_sw_in(m)  +                   &
     4984                               surf_usm_h%rad_lw_in(m)  -                   &
     4985                               surf_usm_h%rad_sw_out(m) -                   &
     4986                               surf_usm_h%rad_lw_out(m)
     4987     ENDDO
     4988     DO  m = 1, surf_lsm_h%ns
     4989        surf_lsm_h%surfhf(m) = surf_lsm_h%rad_sw_in(m)  +                   &
     4990                               surf_lsm_h%rad_lw_in(m)  -                   &
     4991                               surf_lsm_h%rad_sw_out(m) -                   &
     4992                               surf_lsm_h%rad_lw_out(m)
     4993     ENDDO
     4994
     4995     DO  l = 0, 3
     4996!--     urban
     4997        DO  m = 1, surf_usm_v(l)%ns
     4998           surf_usm_v(l)%surfhf(m) = surf_usm_v(l)%rad_sw_in(m)  +          &
     4999                                     surf_usm_v(l)%rad_lw_in(m)  -          &
     5000                                     surf_usm_v(l)%rad_sw_out(m) -          &
     5001                                     surf_usm_v(l)%rad_lw_out(m)
     5002        ENDDO
     5003!--     land
     5004        DO  m = 1, surf_lsm_v(l)%ns
     5005           surf_lsm_v(l)%surfhf(m) = surf_lsm_v(l)%rad_sw_in(m)  +          &
     5006                                     surf_lsm_v(l)%rad_lw_in(m)  -          &
     5007                                     surf_lsm_v(l)%rad_sw_out(m) -          &
     5008                                     surf_lsm_v(l)%rad_lw_out(m)
     5009
     5010        ENDDO
     5011     ENDDO
     5012!
     5013!--  Calculate the average temperature, albedo, and emissivity for urban/land
     5014!--  domain when using average_radiation in the respective radiation model
     5015
     5016     IF ( average_radiation )  THEN
     5017
     5018!--     precalculate face areas for all face directions using normal vector
     5019        DO d = 0, nsurf_type
     5020            facearea(d) = 1._wp
     5021            IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx
     5022            IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy
     5023            IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz
     5024        ENDDO
     5025!
     5026!--     absorbed/received SW & LW and emitted LW energy of all physical
     5027!--     surfaces (land and urban) in local processor
     5028        pinswl = 0._wp
     5029        pinlwl = 0._wp
     5030        pabsswl = 0._wp
     5031        pabslwl = 0._wp
     5032        pemitlwl = 0._wp
     5033        emiss_sum_surfl = 0._wp
     5034        area_surfl = 0._wp
     5035        DO  i = 1, nsurfl
     5036           d = surfl(id, i)
     5037!--        received SW & LW
     5038           pinswl = pinswl + surfinsw(i) * facearea(d)
     5039           pinlwl = pinlwl + surfinlw(i) * facearea(d)
     5040!--        absorbed SW & LW
     5041           pabsswl = pabsswl + (1._wp - albedo_surf(i)) *                   &
     5042                                                  surfinsw(i) * facearea(d)
     5043           pabslwl = pabslwl + emiss_surf(i) * surfinlw(i) * facearea(d)
     5044!--        emitted LW
     5045           pemitlwl = pemitlwl + surfoutlw(i) * facearea(d)
     5046!--        emissivity and area sum
     5047           emiss_sum_surfl = emiss_sum_surfl + emiss_surf(i) * facearea(d)
     5048           area_surfl = area_surfl + facearea(d)
     5049        END DO
     5050!
     5051!--     add the absorbed SW energy by plant canopy
     5052        IF ( plant_canopy )  THEN
     5053           pabsswl = pabsswl + SUM(pcbinsw)
     5054           pabslwl = pabslwl + SUM(pcbinlw)
     5055        ENDIF
     5056!
     5057!--     gather all rad flux energy in all processors
     5058#if defined( __parallel )
     5059        CALL MPI_ALLREDUCE( pinswl, pinsw, 1, MPI_REAL, MPI_SUM, comm2d, ierr)
     5060        CALL MPI_ALLREDUCE( pinlwl, pinlw, 1, MPI_REAL, MPI_SUM, comm2d, ierr)
     5061        CALL MPI_ALLREDUCE( pabsswl, pabssw, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
     5062        CALL MPI_ALLREDUCE( pabslwl, pabslw, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
     5063        CALL MPI_ALLREDUCE( pemitlwl, pemitlw, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
     5064        CALL MPI_ALLREDUCE( emiss_sum_surfl, emiss_sum_surf, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
     5065        CALL MPI_ALLREDUCE( area_surfl, area_surf, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
     5066#else
     5067        pinsw = pinswl
     5068        pinlw = pinlwl
     5069        pabssw = pabsswl
     5070        pabslwl = pabslw
     5071        pemitlwl = pemitlw
     5072        emiss_sum_surf = emiss_sum_surfl
     5073        area_surf = area_surfl
     5074#endif
     5075
     5076!--     (1) albedo
     5077        IF ( pinsw /= 0.0_wp )  albedo_urb = 1._wp - pabssw / pinsw
     5078
     5079!--     (2) average emmsivity
     5080        IF ( area_surf /= 0.0_wp ) emissivity_urb = emiss_sum_surf / area_surf
     5081
     5082!--     (3) temperature
     5083        t_rad_urb = ((pemitlw - pabslw + emissivity_urb*pinlw)/(emissivity_urb*sigma_sb*area_surf))**0.25_wp
     5084
     5085     ENDIF
     5086
     5087    CONTAINS
     5088
     5089!------------------------------------------------------------------------------!
     5090!> Calculates radiation absorbed by box with given size and LAD.
     5091!>
     5092!> Simulates resol**2 rays (by equally spacing a bounding horizontal square
     5093!> conatining all possible rays that would cross the box) and calculates
     5094!> average transparency per ray. Returns fraction of absorbed radiation flux
     5095!> and area for which this fraction is effective.
     5096!------------------------------------------------------------------------------!
     5097    PURE SUBROUTINE box_absorb(boxsize, resol, dens, uvec, area, absorb)
     5098       IMPLICIT NONE
     5099
     5100       REAL(wp), DIMENSION(3), INTENT(in) :: &
     5101            boxsize, &      !< z, y, x size of box in m
     5102            uvec            !< z, y, x unit vector of incoming flux
     5103       INTEGER(iwp), INTENT(in) :: &
     5104            resol           !< No. of rays in x and y dimensions
     5105       REAL(wp), INTENT(in) :: &
     5106            dens            !< box density (e.g. Leaf Area Density)
     5107       REAL(wp), INTENT(out) :: &
     5108            area, &         !< horizontal area for flux absorbtion
     5109            absorb          !< fraction of absorbed flux
     5110       REAL(wp) :: &
     5111            xshift, yshift, &
     5112            xmin, xmax, ymin, ymax, &
     5113            xorig, yorig, &
     5114            dx1, dy1, dz1, dx2, dy2, dz2, &
     5115            crdist, &
     5116            transp
     5117       INTEGER(iwp) :: &
     5118            i, j
     5119
     5120       xshift = uvec(3) / uvec(1) * boxsize(1)
     5121       xmin = min(0._wp, -xshift)
     5122       xmax = boxsize(3) + max(0._wp, -xshift)
     5123       yshift = uvec(2) / uvec(1) * boxsize(1)
     5124       ymin = min(0._wp, -yshift)
     5125       ymax = boxsize(2) + max(0._wp, -yshift)
     5126
     5127       transp = 0._wp
     5128       DO i = 1, resol
     5129          xorig = xmin + (xmax-xmin) * (i-.5_wp) / resol
     5130          DO j = 1, resol
     5131             yorig = ymin + (ymax-ymin) * (j-.5_wp) / resol
     5132
     5133             dz1 = 0._wp
     5134             dz2 = boxsize(1)/uvec(1)
     5135
     5136             IF ( uvec(2) > 0._wp )  THEN
     5137                dy1 = -yorig             / uvec(2) !< crossing with y=0
     5138                dy2 = (boxsize(2)-yorig) / uvec(2) !< crossing with y=boxsize(2)
     5139             ELSE !uvec(2)==0
     5140                dy1 = -huge(1._wp)
     5141                dy2 = huge(1._wp)
     5142             ENDIF
     5143
     5144             IF ( uvec(3) > 0._wp )  THEN
     5145                dx1 = -xorig             / uvec(3) !< crossing with x=0
     5146                dx2 = (boxsize(3)-xorig) / uvec(3) !< crossing with x=boxsize(3)
     5147             ELSE !uvec(3)==0
     5148                dx1 = -huge(1._wp)
     5149                dx2 = huge(1._wp)
     5150             ENDIF
     5151
     5152             crdist = max(0._wp, (min(dz2, dy2, dx2) - max(dz1, dy1, dx1)))
     5153             transp = transp + exp(-ext_coef * dens * crdist)
     5154          ENDDO
     5155       ENDDO
     5156       transp = transp / resol**2
     5157       area = (boxsize(3)+xshift)*(boxsize(2)+yshift)
     5158       absorb = 1._wp - transp
     5159
     5160    END SUBROUTINE box_absorb
     5161
     5162!------------------------------------------------------------------------------!
     5163! Description:
     5164! ------------
     5165!> This subroutine splits direct and diffusion dw radiation
     5166!> It sould not be called in case the radiation model already does it
     5167!> It follows <CITATION>
     5168!------------------------------------------------------------------------------!
     5169    SUBROUTINE calc_diffusion_radiation
     5170   
     5171        REAL(wp), PARAMETER                          :: lowest_solarUp = 0.1_wp  !< limit the sun elevation to protect stability of the calculation
     5172        INTEGER(iwp)                                 :: i, j
     5173        REAL(wp)                                     ::  year_angle              !< angle
     5174        REAL(wp)                                     ::  etr                     !< extraterestrial radiation
     5175        REAL(wp)                                     ::  corrected_solarUp       !< corrected solar up radiation
     5176        REAL(wp)                                     ::  horizontalETR           !< horizontal extraterestrial radiation
     5177        REAL(wp)                                     ::  clearnessIndex          !< clearness index
     5178        REAL(wp)                                     ::  diff_frac               !< diffusion fraction of the radiation
     5179
     5180       
     5181!--     Calculate current day and time based on the initial values and simulation time
     5182        year_angle = ( (day_of_year_init * 86400) + time_utc_init              &
     5183                        + time_since_reference_point )  * d_seconds_year       &
     5184                        * 2.0_wp * pi
     5185       
     5186        etr = solar_constant * (1.00011_wp +                                   &
     5187                          0.034221_wp * cos(year_angle) +                      &
     5188                          0.001280_wp * sin(year_angle) +                      &
     5189                          0.000719_wp * cos(2.0_wp * year_angle) +             &
     5190                          0.000077_wp * sin(2.0_wp * year_angle))
     5191       
     5192!--   
     5193!--     Under a very low angle, we keep extraterestrial radiation at
     5194!--     the last small value, therefore the clearness index will be pushed
     5195!--     towards 0 while keeping full continuity.
     5196!--   
     5197        IF ( zenith(0) <= lowest_solarUp )  THEN
     5198            corrected_solarUp = lowest_solarUp
     5199        ELSE
     5200            corrected_solarUp = zenith(0)
     5201        ENDIF
     5202       
     5203        horizontalETR = etr * corrected_solarUp
     5204       
     5205        DO i = nxl, nxr
     5206            DO j = nys, nyn
     5207                clearnessIndex = rad_sw_in(0,j,i) / horizontalETR
     5208                diff_frac = 1.0_wp / (1.0_wp + exp(-5.0033_wp + 8.6025_wp * clearnessIndex))
     5209                rad_sw_in_diff(j,i) = rad_sw_in(0,j,i) * diff_frac
     5210                rad_sw_in_dir(j,i)  = rad_sw_in(0,j,i) * (1.0_wp - diff_frac)
     5211                rad_lw_in_diff(j,i) = rad_lw_in(0,j,i)
     5212            ENDDO
     5213        ENDDO
     5214       
     5215    END SUBROUTINE calc_diffusion_radiation
     5216
     5217
     5218 END SUBROUTINE radiation_interaction
     5219   
     5220!------------------------------------------------------------------------------!
     5221! Description:
     5222! ------------
     5223!> This subroutine initializes structures needed for radiative transfer
     5224!> model. This model calculates transformation processes of the
     5225!> radiation inside urban and land canopy layer. The module includes also
     5226!> the interaction of the radiation with the resolved plant canopy.
     5227!>
     5228!> For more info. see Resler et al. 2017
     5229!>
     5230!> The new version 2.0 was radically rewriten, the discretization scheme
     5231!> has been changed. This new version significantly improves effectivity
     5232!> of the paralelization and the scalability of the model.
     5233!>
    44855234!------------------------------------------------------------------------------!
    44865235    SUBROUTINE radiation_interaction_init
    4487    
     5236
    44885237       USE netcdf_data_input_mod,                                              &
    44895238           ONLY:  leaf_area_density_f
    44905239
    4491        USE plant_canopy_model_mod,                                             &     
    4492            ONLY:  pch_index, pc_heating_rate, lad_s, prototype_lad, usm_lad_rma       
    4493        
     5240       USE plant_canopy_model_mod,                                             &
     5241           ONLY:  pch_index, pc_heating_rate, lad_s
     5242
    44945243       IMPLICIT NONE
    44955244
     
    44995248       INTEGER(iwp) :: nzubl, nzutl, isurf, ipcgb
    45005249       INTEGER(iwp) :: procid
    4501 
    4502        INTEGER(iwp), DIMENSION(1:4,inorth_b:iwest_b)  ::  ijdb                               !< start and end of the local domain border coordinates (set in code)
    4503        LOGICAL, DIMENSION(inorth_b:iwest_b)           ::  isborder                           !< is PE on the border of the domain in four corresponding directions
    4504 
    4505 !
    4506 !--    Find nzub, nzut, nzu via wall_flag_0 array (nzb_s_inner will be
     5250       REAL(wp)     :: mrl
     5251
     5252
     5253       !INTEGER(iwp), DIMENSION(1:4,inorth_b:iwest_b)  ::  ijdb                               !< start and end of the local domain border coordinates (set in code)
     5254       !LOGICAL, DIMENSION(inorth_b:iwest_b)           ::  isborder                           !< is PE on the border of the domain in four corresponding directions
     5255
     5256!
     5257!--    Find nzub, nzut, nzu via wall_flag_0 array (nzb_s_inner will be
    45075258!--    removed later). The following contruct finds the lowest / largest index
    4508 !--    for any upward-facing wall (see bit 12). 
     5259!--    for any upward-facing wall (see bit 12).
    45095260       nzubl = MINVAL( get_topography_top_index( 's' ) )
    45105261       nzutl = MAXVAL( get_topography_top_index( 's' ) )
     
    45385289               ENDDO
    45395290           ENDDO
    4540            
     5291
    45415292           nzutl = MAX( nzutl, MAXVAL( pct ) )
    45425293!--        code of plant canopy model uses parameter pch_index
     
    45445295!--        (pch_index, lad_s and other arrays in PCM are defined flat)
    45455296           pch_index = MERGE( leaf_area_density_f%nz - 1, MAXVAL( pch ),       &
    4546                               leaf_area_density_f%from_file ) 
     5297                              leaf_area_density_f%from_file )
    45475298
    45485299           prototype_lad = MAXVAL( lad_s ) * .9_wp  !< better be *1.0 if lad is either 0 or maxval(lad) everywhere
     
    45525303           !CALL message('usm_init_urban_surface', 'PA0520', 0, 0, -1, 6, 0)
    45535304       ENDIF
    4554        
     5305
    45555306       nzutl = MIN( nzutl + nzut_free, nzt )
    45565307
     
    45665317       nzu = nzut - nzub + 1
    45675318!
     5319!--    check max_raytracing_dist relative to urban surface layer height
     5320       mrl = 2.0_wp * nzu * dz
     5321       IF ( max_raytracing_dist <= mrl ) THEN
     5322          IF ( max_raytracing_dist /= -999.0_wp ) THEN
     5323!--          max_raytracing_dist too low
     5324             WRITE(message_string, '(a,f6.1)') 'Max_raytracing_dist too low, ' &
     5325                   // 'override to value ', mrl
     5326             CALL message('radiation_interaction_init', 'PA0521', 0, 0, -1, 6, 0)
     5327          ENDIF
     5328          max_raytracing_dist = mrl
     5329       ENDIF
     5330!
    45685331!--    allocate urban surfaces grid
    45695332!--    calc number of surfaces in local proc
     
    45795342
    45805343!
    4581 !--    Number of vertical surfaces in both USM and LSM. Note that all vertical surface elements are 
     5344!--    Number of vertical surfaces in both USM and LSM. Note that all vertical surface elements are
    45825345!--    already counted in surface_mod.
    45835346       startwall = nsurfl+1
     
    45875350       endwall = nsurfl
    45885351       nwalls  = endwall - startwall + 1
    4589        
    4590 !--    range of energy balance surfaces  ! will be treated separately by surf_usm_h and surf_usm_v
    4591 !--    Do we really need usm_energy_balance_land??!!
    4592 !--    !!! Attention: if usm_energy_balance_land = false then only vertical surfaces will be considered here
    4593        nenergy = 0
    4594        IF ( energy_balance_surf_h )  THEN
    4595            startenergy = startland
    4596            nenergy = nenergy + nlands
    4597        ELSE
    4598            startenergy = startwall
    4599        ENDIF
    4600        IF ( energy_balance_surf_v )  THEN
    4601            endenergy = endwall
    4602            nenergy = nenergy + nwalls
    4603        ELSE
    4604            endenergy = endland
    4605        ENDIF
    4606 
    4607 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    4608 !--    block of virtual surfaces
    4609 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    4610 !--    calculate sky surfaces  ! not used so far!
    4611        startsky = nsurfl+1
    4612        nsurfl = nsurfl+(nxr-nxl+1)*(nyn-nys+1)
    4613        endsky = nsurfl
    4614        nskys = endsky-startsky+1
    4615        
    4616 !--    border flags
    4617 #if defined( __parallel )
    4618        isborder = (/ north_border_pe, south_border_pe, right_border_pe, left_border_pe /)
    4619 #else
    4620        isborder = (/.TRUE.,.TRUE.,.TRUE.,.TRUE./)
    4621 #endif
    4622 !--    fill array of the limits of the local domain borders
    4623        ijdb = RESHAPE( (/ nxl,nxr,nyn,nyn,nxl,nxr,nys,nys,nxr,nxr,nys,nyn,nxl,nxl,nys,nyn /), (/4, 4/) )
    4624 !--    calulation of the free borders of the domain
    4625        startborder = nsurfl + 1
    4626        DO  ids = inorth_b,iwest_b
    4627           IF ( isborder(ids) )  THEN
    4628 !--          free border of the domain in direction ids
    4629              DO  i = ijdb(1,ids), ijdb(2,ids)
    4630                 DO  j = ijdb(3,ids), ijdb(4,ids)
    4631 
    4632                    k_topo  = get_topography_top_index_ji( j, i, 's' )
    4633                    k_topo2 = get_topography_top_index_ji( j-jdir(ids), i-idir(ids), 's' )
    4634 
    4635 
    4636                    k = nzut - MAX( k_topo, k_topo2 )
    4637                    nsurfl = nsurfl + k
    4638                 ENDDO
    4639              ENDDO
    4640           ENDIF
    4641        ENDDO
    4642        endborder = nsurfl
    4643        nborder = endborder - startborder + 1
    4644 
    4645 !--    calulation of the atmospheric virtual surfaces
    4646 !--    each atmospheric cell has 6 faces
    4647        IF ( atm_surfaces ) THEN
    4648           DO i = nxl, nxr
    4649              DO j = nys, nyn
    4650 !--              Find topography top index
    4651                  k_topo = get_topography_top_index_ji( j, i, 's' )
    4652                  k = nzut - k_topo
    4653                  nsurfl = nsurfl + 6 * k
    4654              ENDDO
    4655           ENDDO
    4656 !--       exclude the local physical surfaces
    4657           nsurfl = nsurfl - nlands - nwalls
    4658 !--       exclude the local virtual surfaces
    4659           nsurfl = nsurfl - nskys - nborder
    4660        ENDIF
    46615352
    46625353!--    fill gridpcbl and pcbl
     
    46645355           ALLOCATE( pcbl(iz:ix, 1:npcbl) )
    46655356           ALLOCATE( gridpcbl(nzub:nzut,nys:nyn,nxl:nxr) )
     5357           pcbl = -1
    46665358           gridpcbl(:,:,:) = 0
    46675359           ipcgb = 0
     
    46795371               ENDDO
    46805372           ENDDO
    4681 
    46825373           ALLOCATE( pcbinsw( 1:npcbl ) )
     5374           ALLOCATE( pcbinswdir( 1:npcbl ) )
     5375           ALLOCATE( pcbinswdif( 1:npcbl ) )
    46835376           ALLOCATE( pcbinlw( 1:npcbl ) )
    46845377       ENDIF
    46855378
    4686 !--    fill surfl
    4687        ALLOCATE(surfl(5,nsurfl))  ! is it mecessary to allocate it with (5,nsurfl)?       
     5379!--    fill surfl (the ordering of local surfaces given by the following
     5380!--    cycles must not be altered, certain file input routines may depend
     5381!--    on it)
     5382       ALLOCATE(surfl(5,nsurfl))  ! is it mecessary to allocate it with (5,nsurfl)?
    46885383       isurf = 0
    4689        
     5384
    46905385!--    add horizontal surface elements (land and urban surfaces)
    46915386!--    TODO: add urban overhanging surfaces (idown_u)
     
    47055400                 surfl(:,isurf) = (/iup_l,k,j,i,m/)
    47065401              ENDDO
    4707              
     5402
    47085403           ENDDO
    47095404       ENDDO
    47105405
    47115406!--    add vertical surface elements (land and urban surfaces)
    4712 !--    TODO: remove the hard coding of l = 0 to l = idirection       
     5407!--    TODO: remove the hard coding of l = 0 to l = idirection
    47135408       DO i = nxl, nxr
    47145409           DO j = nys, nyn
     
    47715466       ENDDO
    47725467
    4773 !--    add sky
    4774        DO i = nxl, nxr
    4775            DO j = nys, nyn
    4776                isurf = isurf + 1
    4777                k = nzut
    4778                surfl(:,isurf) = (/isky,k,j,i,-1/)
    4779            ENDDO
    4780        ENDDO
    4781        
    4782 !--    calulation of the free borders of the domain
    4783        DO ids = inorth_b,iwest_b
    4784            IF ( isborder(ids) )  THEN
    4785 !--            free border of the domain in direction ids
    4786                DO i = ijdb(1,ids), ijdb(2,ids)
    4787                    DO j = ijdb(3,ids), ijdb(4,ids)
    4788                        k_topo  = get_topography_top_index_ji( j, i, 's' )
    4789                        k_topo2 = get_topography_top_index_ji( j-jdir(ids), i-idir(ids), 's' )
    4790 
    4791                        DO k = MAX(k_topo,k_topo2)+1, nzut
    4792                            isurf = isurf + 1
    4793                            surfl(:,isurf) = (/ids,k,j,i,-1/)
    4794                        ENDDO
    4795                    ENDDO
    4796                ENDDO
    4797            ENDIF
    4798        ENDDO
    4799 
    4800 !--    adding the atmospheric virtual surfaces
    4801        IF ( atm_surfaces ) THEN
    4802 !-- TODO: use flags to identfy atmospheric cells and its coresponding surfaces           
    4803 !--    add horizontal surface
    4804           DO i = nxl, nxr
    4805              DO j = nys, nyn
    4806                 k_topo = get_topography_top_index_ji( j, i, 's' )
    4807 
    4808 !--             add upward surface
    4809                 DO k = (k_topo+1), nzut-1
    4810                    isurf = isurf + 1
    4811                    surfl(:,isurf) = (/iup_a,k+1,j,i,-1/)
    4812                 ENDDO
    4813 
    4814 !--             add downward surface
    4815                 DO k = (k_topo+1), nzut-1
    4816                    isurf = isurf + 1
    4817                    surfl(:,isurf) = (/idown_a,k,j,i,-1/)
    4818                 ENDDO
    4819              ENDDO
    4820           ENDDO
    4821 
    4822 !--       add vertical surfaces
    4823           DO i = nxl, nxr
    4824              DO j = nys, nyn
    4825                 k_topo = get_topography_top_index_ji( j, i, 's' )
    4826 !--             north
    4827                 IF ( j /= ny ) THEN
    4828                    ids = inorth_a
    4829                    jr = min(max(j-jdir(ids),0),ny)
    4830                    ir = min(max(i-idir(ids),0),nx)
    4831                    k_topo2 = get_topography_top_index_ji( jr, ir, 's' )
    4832                    DO k = MAX(k_topo,k_topo2)+1, nzut
    4833                       isurf = isurf + 1
    4834                       surfl(:,isurf) = (/inorth_a,k,j,i,-1/)
    4835                    ENDDO
    4836                 END IF
    4837 !--             south
    4838                 IF ( j /= 0 ) THEN
    4839                    ids = isouth_a
    4840                    jr = min(max(j-jdir(ids),0),ny)
    4841                    ir = min(max(i-idir(ids),0),nx)
    4842                    k_topo2 = get_topography_top_index_ji( jr, ir, 's' )
    4843 
    4844                    DO k = MAX(k_topo,k_topo2)+1, nzut
    4845                       isurf = isurf + 1
    4846                       surfl(:,isurf) = (/isouth_a,k,j,i,-1/)
    4847                    ENDDO
    4848                 END IF
    4849 !--             east
    4850                 IF ( i /= nx ) THEN
    4851                    ids = ieast_a
    4852                    jr = min(max(j-jdir(ids),0),ny)
    4853                    ir = min(max(i-idir(ids),0),nx)
    4854                    k_topo2 = get_topography_top_index_ji( jr, ir, 's' )
    4855 
    4856                    DO k = MAX(k_topo,k_topo2)+1, nzut
    4857                       isurf = isurf + 1
    4858                       surfl(:,isurf) = (/ieast_a,k,j,i,-1/)
    4859                    ENDDO
    4860                 END IF
    4861 !--             west
    4862                 IF ( i /= 0 ) THEN
    4863                    ids = iwest_a
    4864                    jr = min(max(j-jdir(ids),0),ny)
    4865                    ir = min(max(i-idir(ids),0),nx)
    4866                    k_topo2 = get_topography_top_index_ji( jr, ir, 's' )
    4867 
    4868                    DO k = MAX(k_topo,k_topo2)+1, nzut
    4869                       isurf = isurf + 1
    4870                       surfl(:,isurf) = (/iwest_a,k,j,i,-1/)
    4871                    ENDDO
    4872                 END IF
    4873              ENDDO
    4874           ENDDO
    4875 
    4876        ENDIF
    4877 
    4878 !
    4879 !--     broadband albedo of the land, roof and wall surface
    4880 !--     for domain border and sky set artifically to 1.0
    4881 !--     what allows us to calculate heat flux leaving over
    4882 !--     side and top borders of the domain
    4883         ALLOCATE ( albedo_surf(nsurfl) )
    4884         albedo_surf = 1.0_wp
    4885 !
    4886 !--     Also allocate further array for emissivity with identical order of
    4887 !--     surface elements as radiation arrays.
    4888 !--     MS: Why startenergy:endenergy and albedo surf from 1:nsurfl ? 
    4889         ALLOCATE ( emiss_surf(startenergy:endenergy)  )
     5468!
     5469!--    broadband albedo of the land, roof and wall surface
     5470!--    for domain border and sky set artifically to 1.0
     5471!--    what allows us to calculate heat flux leaving over
     5472!--    side and top borders of the domain
     5473       ALLOCATE ( albedo_surf(nsurfl) )
     5474       albedo_surf = 1.0_wp
     5475!
     5476!--    Also allocate further array for emissivity with identical order of
     5477!--    surface elements as radiation arrays.
     5478       ALLOCATE ( emiss_surf(nsurfl)  )
    48905479
    48915480
     
    48935482!--    global array surf of indices of surfaces and displacement index array surfstart
    48945483       ALLOCATE(nsurfs(0:numprocs-1))
    4895        
     5484
    48965485#if defined( __parallel )
    48975486       CALL MPI_Allgather(nsurfl,1,MPI_INTEGER,nsurfs,1,MPI_INTEGER,comm2d,ierr)
     
    49085497       nsurf = k
    49095498       ALLOCATE(surf(5,nsurf))
    4910        
     5499
    49115500#if defined( __parallel )
    4912        CALL MPI_AllGatherv(surfl, nsurfl*5, MPI_INTEGER, surf, nsurfs*5, surfstart*5, MPI_INTEGER, comm2d, ierr)
     5501       CALL MPI_AllGatherv(surfl, nsurfl*5, MPI_INTEGER, surf, nsurfs*5, &
     5502           surfstart(0:numprocs-1)*5, MPI_INTEGER, comm2d, ierr)
    49135503#else
    49145504       surf = surfl
     
    49205510!--    rad_sw_in, rad_lw_in are computed in radiation model,
    49215511!--    splitting of direct and diffusion part is done
    4922 !--    in usm_calc_diffusion_radiation for now
     5512!--    in calc_diffusion_radiation for now
    49235513
    49245514       ALLOCATE( rad_sw_in_dir(nysg:nyng,nxlg:nxrg) )
     
    49275517       rad_sw_in_dir  = 0.0_wp
    49285518       rad_sw_in_diff = 0.0_wp
    4929        rad_lw_in_diff = 0.0_wp 
    4930        
     5519       rad_lw_in_diff = 0.0_wp
     5520
    49315521!--    allocate radiation arrays
    49325522       ALLOCATE( surfins(nsurfl) )
     
    49375527       ALLOCATE( surfinswdif(nsurfl) )
    49385528       ALLOCATE( surfinlwdif(nsurfl) )
    4939        ALLOCATE( surfoutsl(startenergy:endenergy) )
    4940        ALLOCATE( surfoutll(startenergy:endenergy) )
    4941        ALLOCATE( surfoutsw(startenergy:endenergy) )
    4942        ALLOCATE( surfoutlw(startenergy:endenergy) )
    4943        ALLOCATE( surfouts(nsurf) ) !TODO: global surfaces without virtual
    4944        ALLOCATE( surfoutl(nsurf) ) !TODO: global surfaces without virtual
    4945 
    4946 !
    4947 !--    @Mohamed
     5529       ALLOCATE( surfoutsl(nsurfl) )
     5530       ALLOCATE( surfoutll(nsurfl) )
     5531       ALLOCATE( surfoutsw(nsurfl) )
     5532       ALLOCATE( surfoutlw(nsurfl) )
     5533       ALLOCATE( surfouts(nsurf) )
     5534       ALLOCATE( surfoutl(nsurf) )
     5535       ALLOCATE( skyvf(nsurfl) )
     5536       ALLOCATE( skyvft(nsurfl) )
     5537
     5538!
    49485539!--    In case of average_radiation, aggregated surface albedo and emissivity,
    4949 !--    also set initial value of t_rad_urb.
    4950 !--    For the moment set an arbitrary initial value.
     5540!--    also set initial value for t_rad_urb.
     5541!--    For now set an arbitrary initial value.
    49515542       IF ( average_radiation )  THEN
    49525543          albedo_urb = 0.5_wp
    49535544          emissivity_urb = 0.5_wp
    4954           t_rad_urb = pt_surface   
    4955        ENDIF 
     5545          t_rad_urb = pt_surface
     5546       ENDIF
    49565547
    49575548    END SUBROUTINE radiation_interaction_init
     5549
    49585550!------------------------------------------------------------------------------!
    49595551! Description:
    49605552! ------------
    4961 !> This subroutine calculates interaction of the solar radiation
    4962 !> with urban and land surfaces and updates all surface heatfluxes, including
    4963 !> the vertual atmospheric cell faces. It calculates also the required parameters
    4964 !> for RRTMG lower BC.
    4965 !> 
    4966 !> For more info. see Resler et al. 2017
    4967 !> 
     5553!> Calculates shape view factors (SVF), plant sink canopy factors (PCSF),
     5554!> sky-view factors, discretized path for direct solar radiation, MRT factors
     5555!> and other preprocessed data needed for radiation_interaction.
    49685556!------------------------------------------------------------------------------!
    4969     SUBROUTINE radiation_interaction
    4970    
    4971      
    4972       USE control_parameters
    4973 
    4974       USE plant_canopy_model_mod,                                                &
    4975            ONLY: prototype_lad
     5557    SUBROUTINE radiation_calc_svf
    49765558   
    49775559        IMPLICIT NONE
    4978        
    4979         INTEGER(iwp)               :: i, j, k, kk, is, js, d, ku, refstep, m, mm, l, ll
    4980         INTEGER(iwp)               :: ii, jj !< running indices
    4981         INTEGER(iwp)               :: nzubl, nzutl, isurf, isurfsrc, isurf1, isvf, icsf, ipcgb
    4982         INTEGER(iwp), DIMENSION(4) :: bdycross
    4983         REAL(wp), DIMENSION(3,3)   :: mrot            !< grid rotation matrix (xyz)
    4984         REAL(wp), DIMENSION(3,0:nsurf_type) :: vnorm  !< face direction normal vectors (xyz)
    4985         REAL(wp), DIMENSION(3)     :: sunorig         !< grid rotated solar direction unit vector (xyz)
    4986         REAL(wp), DIMENSION(3)     :: sunorig_grid    !< grid squashed solar direction unit vector (zyx)
    4987         REAL(wp), DIMENSION(0:nsurf_type)  :: costheta        !< direct irradiance factor of solar angle
    4988         REAL(wp), DIMENSION(nzub:nzut) :: pchf_prep   !< precalculated factor for canopy temp tendency
    4989         REAL(wp), PARAMETER        :: alpha = 0._wp   !< grid rotation (TODO: add to namelist or remove)
    4990         REAL(wp)                   :: rx, ry, rz
    4991         REAL(wp)                   :: pc_box_area, pc_abs_frac, pc_abs_eff
    4992         INTEGER(iwp)               :: pc_box_dimshift !< transform for best accuracy
    4993         INTEGER(iwp), DIMENSION(0:3) :: reorder = (/ 1, 0, 3, 2 /)
    4994         REAL(wp),     DIMENSION(0:nsurf_type)       :: facearea
    4995         REAL(wp)                   :: pabsswl  = 0.0_wp  !< total absorbed SW radiation energy in local processor (W)
    4996         REAL(wp)                   :: pabssw   = 0.0_wp  !< total absorbed SW radiation energy in all processors (W)
    4997         REAL(wp)                   :: pabslwl  = 0.0_wp  !< total absorbed LW radiation energy in local processor (W)
    4998         REAL(wp)                   :: pabslw   = 0.0_wp  !< total absorbed LW radiation energy in all processors (W)
    4999         REAL(wp)                   :: pemitlwl = 0.0_wp  !< total emitted LW radiation energy in all processors (W)
    5000         REAL(wp)                   :: pemitlw  = 0.0_wp  !< total emitted LW radiation energy in all processors (W)
    5001         REAL(wp)                   :: pinswl   = 0.0_wp  !< total received SW radiation energy in local processor (W)
    5002         REAL(wp)                   :: pinsw    = 0.0_wp  !< total received SW radiation energy in all processor (W)
    5003         REAL(wp)                   :: pinlwl   = 0.0_wp  !< total received LW radiation energy in local processor (W)
    5004         REAL(wp)                   :: pinlw    = 0.0_wp  !< total received LW radiation energy in all processor (W)
    5005         REAL(wp)                   :: emiss_sum_surfl    !< sum of emissisivity of surfaces in local processor
    5006         REAL(wp)                   :: emiss_sum_surf     !< sum of emissisivity of surfaces in all processor
    5007         REAL(wp)                   :: area_surfl         !< total area of surfaces in local processor
    5008         REAL(wp)                   :: area_surf          !< total area of surfaces in all processor
    50095560       
    5010         IF ( plant_canopy )  THEN
    5011             pchf_prep(:) = r_d * (hyp(nzub:nzut) / 100000.0_wp)**0.286_wp &
    5012                         / (cp * hyp(nzub:nzut) * dx*dy*dz) !< equals to 1 / (rho * c_p * Vbox * T)
    5013         ENDIF
    5014 
    5015         sun_direction = .TRUE.
    5016         CALL calc_zenith  !< required also for diffusion radiation
    5017 
    5018 !--     prepare rotated normal vectors and irradiance factor
    5019         vnorm(1,:) = idir(:)
    5020         vnorm(2,:) = jdir(:)
    5021         vnorm(3,:) = kdir(:)
    5022         mrot(1, :) = (/ cos(alpha), -sin(alpha), 0._wp /)
    5023         mrot(2, :) = (/ sin(alpha),  cos(alpha), 0._wp /)
    5024         mrot(3, :) = (/ 0._wp,       0._wp,      1._wp /)
    5025         sunorig = (/ sun_dir_lon, sun_dir_lat, zenith(0) /)
    5026         sunorig = matmul(mrot, sunorig)
    5027         DO d = 0, nsurf_type
    5028             costheta(d) = dot_product(sunorig, vnorm(:,d))
    5029         ENDDO
    5030        
    5031         IF ( zenith(0) > 0 )  THEN
    5032 !--         now we will "squash" the sunorig vector by grid box size in
    5033 !--         each dimension, so that this new direction vector will allow us
    5034 !--         to traverse the ray path within grid coordinates directly
    5035             sunorig_grid = (/ sunorig(3)/dz, sunorig(2)/dy, sunorig(1)/dx /)
    5036 !--         sunorig_grid = sunorig_grid / norm2(sunorig_grid)
    5037             sunorig_grid = sunorig_grid / SQRT(SUM(sunorig_grid**2))
    5038 
    5039             IF ( plant_canopy )  THEN
    5040 !--            precompute effective box depth with prototype Leaf Area Density
    5041                pc_box_dimshift = maxloc(sunorig, 1) - 1
    5042                CALL box_absorb(cshift((/dx,dy,dz/), pc_box_dimshift),          &
    5043                                    60, prototype_lad,                          &
    5044                                    cshift(sunorig, pc_box_dimshift),           &
    5045                                    pc_box_area, pc_abs_frac)
    5046                pc_box_area = pc_box_area * sunorig(pc_box_dimshift+1) / sunorig(3)
    5047                pc_abs_eff = log(1._wp - pc_abs_frac) / prototype_lad
    5048             ENDIF
    5049         ENDIF
    5050        
    5051 !--     split diffusion and direct part of the solar downward radiation
    5052 !--     comming from radiation model and store it in 2D arrays
    5053 !--     rad_sw_in_diff, rad_sw_in_dir and rad_lw_in_diff
    5054         IF ( split_diffusion_radiation )  THEN
    5055             CALL calc_diffusion_radiation
    5056         ELSE
    5057            DO  i = nxl, nxr
    5058               DO  j = nys, nyn
    5059                  DO  m = surf_def_h(0)%start_index(j,i),                       &
    5060                          surf_def_h(0)%end_index(j,i)
    5061                     rad_sw_in_diff(j,i) = 0.0_wp
    5062                     rad_sw_in_dir(j,i)  = surf_def_h(0)%rad_sw_in(m)
    5063                     rad_lw_in_diff(j,i) = surf_def_h(0)%rad_lw_in(m)
    5064                  ENDDO
    5065                  DO  m = surf_lsm_h%start_index(j,i),                          &
    5066                          surf_lsm_h%end_index(j,i)
    5067                     rad_sw_in_diff(j,i) = 0.0_wp
    5068                     rad_sw_in_dir(j,i)  = surf_lsm_h%rad_sw_in(m)
    5069                     rad_lw_in_diff(j,i) = surf_lsm_h%rad_lw_in(m)
    5070                  ENDDO
    5071                  DO  m = surf_usm_h%start_index(j,i),                          &
    5072                          surf_usm_h%end_index(j,i)
    5073                     rad_sw_in_diff(j,i) = 0.0_wp
    5074                     rad_sw_in_dir(j,i)  = surf_usm_h%rad_sw_in(m)
    5075                     rad_lw_in_diff(j,i) = surf_usm_h%rad_lw_in(m)
    5076                  ENDDO
    5077               ENDDO
    5078            ENDDO
    5079         ENDIF
    5080 
    5081 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    5082 !--     First pass: direct + diffuse irradiance
    5083 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    5084         surfinswdir   = 0._wp !nsurfl
    5085         surfinswdif   = 0._wp !nsurfl
    5086         surfinlwdif   = 0._wp !nsurfl
    5087         surfins       = 0._wp !nsurfl
    5088         surfinl       = 0._wp !nsurfl
    5089         surfoutsl(:)  = 0.0_wp !start-end
    5090         surfoutll(:)  = 0.0_wp !start-end
    5091        
    5092 !--     Set up thermal radiation from surfaces
    5093 !--     emiss_surf is defined only for surfaces for which energy balance is calculated
    5094 !--     Workaround: reorder surface data type back on 1D array including all surfaces,
    5095 !--     which implies to reorder horizontal and vertical surfaces
    5096 !
    5097 !--     Horizontal walls
    5098         mm = 1
    5099         DO  i = nxl, nxr
    5100            DO  j = nys, nyn
    5101 !--           urban
    5102               DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
    5103                  surfoutll(mm) = SUM ( surf_usm_h%frac(:,m) *                  &
    5104                                        surf_usm_h%emissivity(:,m) )            &
    5105                                      * sigma_sb                                &
    5106                                      * surf_usm_h%pt_surface(m)**4
    5107                  albedo_surf(mm) = SUM ( surf_usm_h%frac(:,m) *                &
    5108                                          surf_usm_h%albedo(:,m) )       
    5109                  emiss_surf(mm)  = SUM ( surf_usm_h%frac(:,m) *                &
    5110                                          surf_usm_h%emissivity(:,m) ) 
    5111                  mm = mm + 1
    5112               ENDDO
    5113 !--           land
    5114               DO  m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
    5115                  surfoutll(mm) = SUM ( surf_lsm_h%frac(:,m) *                  &
    5116                                        surf_lsm_h%emissivity(:,m) )            &
    5117                                      * sigma_sb                                &
    5118                                      * surf_lsm_h%pt_surface(m)**4
    5119                  albedo_surf(mm) = SUM ( surf_lsm_h%frac(:,m) *                &
    5120                                          surf_lsm_h%albedo(:,m) )       
    5121                  emiss_surf(mm)  = SUM ( surf_lsm_h%frac(:,m) *                &
    5122                                          surf_lsm_h%emissivity(:,m) )   
    5123                  mm = mm + 1
    5124               ENDDO
    5125            ENDDO
    5126         ENDDO
    5127 !
    5128 !--     Vertical walls
    5129         DO  i = nxl, nxr
    5130            DO  j = nys, nyn
    5131               DO  ll = 0, 3
    5132                  l = reorder(ll)
    5133 !--              urban
    5134                  DO  m = surf_usm_v(l)%start_index(j,i),                       &
    5135                          surf_usm_v(l)%end_index(j,i)
    5136                     surfoutll(mm) = SUM ( surf_usm_v(l)%frac(:,m) *            &
    5137                                           surf_usm_v(l)%emissivity(:,m) )      &
    5138                                      * sigma_sb                                &
    5139                                      * surf_usm_v(l)%pt_surface(m)**4
    5140                     albedo_surf(mm) = SUM ( surf_usm_v(l)%frac(:,m) *          &
    5141                                             surf_usm_v(l)%albedo(:,m) )   
    5142                     emiss_surf(mm)  = SUM ( surf_usm_v(l)%frac(:,m) *          &
    5143                                             surf_usm_v(l)%emissivity(:,m) ) 
    5144                     mm = mm + 1
    5145                  ENDDO
    5146 !--              land
    5147                  DO  m = surf_lsm_v(l)%start_index(j,i),                       &
    5148                          surf_lsm_v(l)%end_index(j,i)
    5149                     surfoutll(mm) = SUM ( surf_lsm_v(l)%frac(:,m) *            &
    5150                                           surf_lsm_v(l)%emissivity(:,m) )      &
    5151                                      * sigma_sb                                &
    5152                                      * surf_lsm_v(l)%pt_surface(m)**4
    5153                     albedo_surf(mm) = SUM ( surf_lsm_v(l)%frac(:,m) *          &
    5154                                             surf_lsm_v(l)%albedo(:,m) )   
    5155                     emiss_surf(mm)  = SUM ( surf_lsm_v(l)%frac(:,m) *          &
    5156                                             surf_lsm_v(l)%emissivity(:,m) )
    5157                     mm = mm + 1
    5158                  ENDDO
    5159               ENDDO
    5160            ENDDO
    5161         ENDDO
    5162 
     5561        INTEGER(iwp)                                  :: i, j, k, l, d, ip, jp
     5562        INTEGER(iwp)                                  :: isvf, ksvf, icsf, kcsf, npcsfl, isvf_surflt, imrtt, imrtf, ipcgb
     5563        INTEGER(iwp)                                  :: sd, td, ioln, iproc
     5564        INTEGER(iwp)                                  :: iaz, izn      !< azimuth, zenith counters
     5565        INTEGER(iwp)                                  :: naz, nzn      !< azimuth, zenith num of steps
     5566        REAL(wp)                                      :: az0, zn0      !< starting azimuth/zenith
     5567        REAL(wp)                                      :: azs, zns      !< azimuth/zenith cycle step
     5568        REAL(wp)                                      :: az1, az2      !< relative azimuth of section borders
     5569        REAL(wp)                                      :: azmid         !< ray (center) azimuth
     5570        REAL(wp)                                      :: horizon       !< computed horizon height (tangent of elevation)
     5571        REAL(wp)                                      :: azen          !< zenith angle
     5572        REAL(wp), DIMENSION(:), ALLOCATABLE           :: zdirs         !< directions in z (tangent of elevation)
     5573        REAL(wp), DIMENSION(:), ALLOCATABLE           :: zbdry         !< zenith angle boundaries
     5574        REAL(wp), DIMENSION(:), ALLOCATABLE           :: vffrac        !< view factor fractions for individual rays
     5575        REAL(wp), DIMENSION(:), ALLOCATABLE           :: ztransp       !< array of transparency in z steps
     5576        REAL(wp),     DIMENSION(0:nsurf_type)         :: facearea
     5577        INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE     :: nzterrl
     5578        REAL(wp),     DIMENSION(:,:), ALLOCATABLE     :: csflt, pcsflt
     5579        INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE     :: kcsflt,kpcsflt
     5580        INTEGER(iwp), DIMENSION(:), ALLOCATABLE       :: icsflt,dcsflt,ipcsflt,dpcsflt
     5581        REAL(wp), DIMENSION(3)                        :: uv
     5582        LOGICAL                                       :: visible
     5583        REAL(wp), DIMENSION(3)                        :: sa, ta          !< real coordinates z,y,x of source and target
     5584        REAL(wp)                                      :: transparency, rirrf, sqdist, svfsum
     5585        INTEGER(iwp)                                  :: isurflt, isurfs, isurflt_prev
     5586        INTEGER(iwp)                                  :: itx, ity, itz
     5587        INTEGER(idp)                                  :: ray_skip_maxdist, ray_skip_minval !< skipped raytracing counts
     5588        INTEGER(iwp)                                  :: max_track_len !< maximum 2d track length
     5589        CHARACTER(len=7)                              :: pid_char = ''
     5590        INTEGER(iwp)                                  :: win_lad, minfo
     5591        REAL(wp), DIMENSION(:,:,:), POINTER           :: lad_s_rma       !< fortran pointer, but lower bounds are 1
     5592        TYPE(c_ptr)                                   :: lad_s_rma_p     !< allocated c pointer
    51635593#if defined( __parallel )
    5164 !--     might be optimized and gather only values relevant for current processor
    5165        
    5166         CALL MPI_AllGatherv(surfoutll, nenergy, MPI_REAL, &
    5167                             surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr) !nsurf global
    5168 #else
    5169         surfoutl(:) = surfoutll(:) !nsurf global
     5594        INTEGER(kind=MPI_ADDRESS_KIND)                :: size_lad_rma
    51705595#endif
    5171        
    5172         isurf1 = -1   !< previous processed surface
    5173         DO isvf = 1, nsvfl
    5174             isurf = svfsurf(1, isvf)
    5175             k = surfl(iz, isurf)
    5176             j = surfl(iy, isurf)
    5177             i = surfl(ix, isurf)
    5178             isurfsrc = svfsurf(2, isvf)
    5179             IF ( zenith(0) > 0  .AND.  isurf /= isurf1 )  THEN
    5180 !--             locate the virtual surface where the direct solar ray crosses domain boundary
    5181 !--             (once per target surface)
    5182                 d = surfl(id, isurf)
    5183                 rz = REAL(k, wp) - 0.5_wp * kdir(d)
    5184                 ry = REAL(j, wp) - 0.5_wp * jdir(d)
    5185                 rx = REAL(i, wp) - 0.5_wp * idir(d)
    5186                
    5187                 CALL find_boundary_face( (/ rz, ry, rx /), sunorig_grid, bdycross)
    5188                
    5189                 isurf1 = isurf
    5190             ENDIF
    5191 
    5192             IF ( surf(id, isurfsrc) >= isky )  THEN
    5193 !--             diffuse rad from boundary surfaces. Since it is a simply
    5194 !--             calculated value, it is not assigned to surfref(s/l),
    5195 !--             instead it is used directly here
    5196 !--             we consider the radiation from the radiation model falling on surface
    5197 !--             as the radiation falling on the top of urban layer into the place of the source surface
    5198 !--             we consider it as a very reasonable simplification which allow as avoid
    5199 !--             necessity of other global range arrays and some all to all mpi communication
    5200                 surfinswdif(isurf) = surfinswdif(isurf) + rad_sw_in_diff(j,i) * svf(1,isvf) * svf(2,isvf)
    5201                                                                 !< canopy shading is applied only to shortwave
    5202                 surfinlwdif(isurf) = surfinlwdif(isurf) + rad_lw_in_diff(j,i) * svf(1,isvf)
    5203             ELSE
    5204 !--             for surface-to-surface factors we calculate thermal radiation in 1st pass
    5205                 surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc)
    5206             ENDIF
    5207 
    5208             IF ( zenith(0) > 0  .AND.  all( surf(1:4,isurfsrc) == bdycross ) )  THEN
    5209 !--             found svf between model boundary and the face => face isn't shaded
    5210                 surfinswdir(isurf) = rad_sw_in_dir(j,i) &
    5211                     * costheta(surfl(id, isurf)) * svf(2,isvf) / zenith(0)
    5212 
    5213             ENDIF
    5214         ENDDO
    5215 
    5216         IF ( plant_canopy )  THEN
    5217        
    5218             pcbinsw(:) = 0._wp
    5219             pcbinlw(:) = 0._wp  !< will stay always 0 since we don't absorb lw anymore
    5220             !
    5221 !--         pcsf first pass
    5222             isurf1 = -1  !< previous processed pcgb
    5223             DO icsf = 1, ncsfl
    5224                 ipcgb = csfsurf(1, icsf)
    5225                 i = pcbl(ix,ipcgb)
    5226                 j = pcbl(iy,ipcgb)
    5227                 k = pcbl(iz,ipcgb)
    5228                 isurfsrc = csfsurf(2, icsf)
    5229 
    5230                 IF ( zenith(0) > 0  .AND.  ipcgb /= isurf1 )  THEN
    5231 !--                 locate the virtual surface where the direct solar ray crosses domain boundary
    5232 !--                 (once per target PC gridbox)
    5233                     rz = REAL(k, wp)
    5234                     ry = REAL(j, wp)
    5235                     rx = REAL(i, wp)
    5236                     CALL find_boundary_face( (/ rz, ry, rx /), &
    5237                         sunorig_grid, bdycross)
    5238 
    5239                     isurf1 = ipcgb
    5240                 ENDIF
    5241 
    5242                 IF ( surf(id, isurfsrc) >= isky )  THEN
    5243 !--                 Diffuse rad from boundary surfaces. See comments for svf above.
    5244                     pcbinsw(ipcgb) = pcbinsw(ipcgb) + csf(1,icsf) * csf(2,icsf) * rad_sw_in_diff(j,i)
    5245 !--                 canopy shading is applied only to shortwave, therefore no absorbtion for lw
    5246 !--                 pcbinlw(ipcgb) = pcbinlw(ipcgb) + svf(1,isvf) * rad_lw_in_diff(j,i)
    5247                 !ELSE
    5248 !--                 Thermal radiation in 1st pass
    5249 !--                 pcbinlw(ipcgb) = pcbinlw(ipcgb) + svf(1,isvf) * surfoutl(isurfsrc)
    5250                 ENDIF
    5251 
    5252                 IF ( zenith(0) > 0  .AND.  ALL( surf(1:4,isurfsrc) == bdycross ) )  THEN
    5253 !--                 found svf between model boundary and the pcgb => pcgb isn't shaded
    5254                     pc_abs_frac = 1._wp - EXP(pc_abs_eff * lad_s(k,j,i))
    5255                     pcbinsw(ipcgb) = pcbinsw(ipcgb) &
    5256                         + rad_sw_in_dir(j, i) * pc_box_area * csf(2,icsf) * pc_abs_frac
    5257                 ENDIF
    5258             ENDDO
    5259         ENDIF
    5260 
    5261         surfins(startenergy:endenergy) = surfinswdir(startenergy:endenergy) + surfinswdif(startenergy:endenergy)
    5262         surfinl(startenergy:endenergy) = surfinl(startenergy:endenergy) + surfinlwdif(startenergy:endenergy)
    5263         surfinsw(:) = surfins(:)
    5264         surfinlw(:) = surfinl(:)
    5265         surfoutsw(:) = 0.0_wp
    5266         surfoutlw(:) = surfoutll(:)
    5267 !         surfhf(startenergy:endenergy) = surfinsw(startenergy:endenergy) + surfinlw(startenergy:endenergy) &
    5268 !                                       - surfoutsw(startenergy:endenergy) - surfoutlw(startenergy:endenergy)
    5269        
    5270 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    5271 !--     Next passes - reflections
    5272 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    5273         DO refstep = 1, nrefsteps
    5274        
    5275             surfoutsl(startenergy:endenergy) = albedo_surf(startenergy:endenergy) * surfins(startenergy:endenergy)
    5276 !--         for non-transparent surfaces, longwave albedo is 1 - emissivity
    5277             surfoutll(startenergy:endenergy) = (1._wp - emiss_surf(startenergy:endenergy)) * surfinl(startenergy:endenergy)
    5278 
    5279 #if defined( __parallel )
    5280             CALL MPI_AllGatherv(surfoutsl, nsurfl, MPI_REAL, &
    5281                 surfouts, nsurfs, surfstart, MPI_REAL, comm2d, ierr)
    5282             CALL MPI_AllGatherv(surfoutll, nsurfl, MPI_REAL, &
    5283                 surfoutl, nsurfs, surfstart, MPI_REAL, comm2d, ierr)
    5284 #else
    5285             surfouts(:) = surfoutsl(:)
    5286             surfoutl(:) = surfoutll(:)
    5287 #endif
    5288 
    5289 !--         reset for next pass input
    5290             surfins(:) = 0._wp
    5291             surfinl(:) = 0._wp
    5292            
    5293 !--         reflected radiation
    5294             DO isvf = 1, nsvfl
    5295                 isurf = svfsurf(1, isvf)
    5296                 isurfsrc = svfsurf(2, isvf)
    5297 
    5298 !--             TODO: to remove if, use start+end for isvf
    5299                 IF ( surf(id, isurfsrc) < isky )  THEN
    5300                     surfins(isurf) = surfins(isurf) + svf(1,isvf) * svf(2,isvf) * surfouts(isurfsrc)
    5301                     surfinl(isurf) = surfinl(isurf) + svf(1,isvf) * surfoutl(isurfsrc)
    5302                 ENDIF
    5303             ENDDO
    5304 
    5305 !--         radiation absorbed by plant canopy
    5306             DO icsf = 1, ncsfl
    5307                 ipcgb = csfsurf(1, icsf)
    5308                 isurfsrc = csfsurf(2, icsf)
    5309 
    5310                 IF ( surf(id, isurfsrc) < isky )  THEN
    5311                     pcbinsw(ipcgb) = pcbinsw(ipcgb) + csf(1,icsf) * csf(2,icsf) * surfouts(isurfsrc)
    5312 !--                 pcbinlw(ipcgb) = pcbinlw(ipcgb) + csf(1,icsf) * surfoutl(isurfsrc)
    5313                 ENDIF
    5314             ENDDO
    5315            
    5316             surfinsw(:) = surfinsw(:)  + surfins(:)
    5317             surfinlw(:) = surfinlw(:)  + surfinl(:)
    5318             surfoutsw(startenergy:endenergy) = surfoutsw(startenergy:endenergy) + surfoutsl(startenergy:endenergy)
    5319             surfoutlw(startenergy:endenergy) = surfoutlw(startenergy:endenergy) + surfoutll(startenergy:endenergy)
    5320 !             surfhf(startenergy:endenergy) = surfinsw(startenergy:endenergy) + surfinlw(startenergy:endenergy) &
    5321 !                                           - surfoutsw(startenergy:endenergy) - surfoutlw(startenergy:endenergy)
    5322        
    5323         ENDDO
    5324 
    5325 !--     push heat flux absorbed by plant canopy to respective 3D arrays
    5326         IF ( plant_canopy )  THEN
    5327             pc_heating_rate(:,:,:) = 0._wp
    5328             DO ipcgb = 1, npcbl
    5329                 j = pcbl(iy, ipcgb)
    5330                 i = pcbl(ix, ipcgb)
    5331                 k = pcbl(iz, ipcgb)
    5332 !
    5333 !--             Following expression equals former kk = k - nzb_s_inner(j,i)
    5334                 kk = k - get_topography_top_index_ji( j, i, 's' )  !- lad arrays are defined flat
    5335                 pc_heating_rate(kk, j, i) = (pcbinsw(ipcgb) + pcbinlw(ipcgb)) &
    5336                     * pchf_prep(k) * pt(k, j, i) !-- = dT/dt
    5337             ENDDO
    5338         ENDIF
    5339 !
    5340 !--     Transfer radiation arrays required for energy balance to the respective data types
    5341         DO  i = startenergy, endenergy
    5342            m  = surfl(5,i)         
    5343 !
    5344 !--        (1) Urban surfaces
    5345 !--        upward-facing
    5346            IF ( surfl(1,i) == iup_u )  THEN
    5347               surf_usm_h%rad_sw_in(m)  = surfinsw(i)
    5348               surf_usm_h%rad_sw_out(m) = surfoutsw(i)
    5349               surf_usm_h%rad_lw_in(m)  = surfinlw(i)
    5350               surf_usm_h%rad_lw_out(m) = surfoutlw(i)
    5351               surf_usm_h%rad_net(m)    = surfinsw(i) - surfoutsw(i) +          &
    5352                                          surfinlw(i) - surfoutlw(i)
    5353 !
    5354 !--        northward-facding
    5355            ELSEIF ( surfl(1,i) == inorth_u )  THEN
    5356               surf_usm_v(0)%rad_sw_in(m)  = surfinsw(i)
    5357               surf_usm_v(0)%rad_sw_out(m) = surfoutsw(i)
    5358               surf_usm_v(0)%rad_lw_in(m)  = surfinlw(i)
    5359               surf_usm_v(0)%rad_lw_out(m) = surfoutlw(i)
    5360               surf_usm_v(0)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
    5361                                             surfinlw(i) - surfoutlw(i)
    5362 !
    5363 !--        southward-facding
    5364            ELSEIF ( surfl(1,i) == isouth_u )  THEN
    5365               surf_usm_v(1)%rad_sw_in(m)  = surfinsw(i)
    5366               surf_usm_v(1)%rad_sw_out(m) = surfoutsw(i)
    5367               surf_usm_v(1)%rad_lw_in(m)  = surfinlw(i)
    5368               surf_usm_v(1)%rad_lw_out(m) = surfoutlw(i)
    5369               surf_usm_v(1)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
    5370                                             surfinlw(i) - surfoutlw(i)
    5371 !
    5372 !--        eastward-facing
    5373            ELSEIF ( surfl(1,i) == ieast_u )  THEN
    5374               surf_usm_v(2)%rad_sw_in(m)  = surfinsw(i)
    5375               surf_usm_v(2)%rad_sw_out(m) = surfoutsw(i)
    5376               surf_usm_v(2)%rad_lw_in(m)  = surfinlw(i)
    5377               surf_usm_v(2)%rad_lw_out(m) = surfoutlw(i)
    5378               surf_usm_v(2)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
    5379                                             surfinlw(i) - surfoutlw(i)
    5380 !
    5381 !--        westward-facding
    5382            ELSEIF ( surfl(1,i) == iwest_u )  THEN
    5383               surf_usm_v(3)%rad_sw_in(m)  = surfinsw(i)
    5384               surf_usm_v(3)%rad_sw_out(m) = surfoutsw(i)
    5385               surf_usm_v(3)%rad_lw_in(m)  = surfinlw(i)
    5386               surf_usm_v(3)%rad_lw_out(m) = surfoutlw(i)
    5387               surf_usm_v(3)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
    5388                                             surfinlw(i) - surfoutlw(i)
    5389 !
    5390 !--        (2) land surfaces
    5391 !--        upward-facing
    5392            ELSEIF ( surfl(1,i) == iup_l )  THEN
    5393               surf_lsm_h%rad_sw_in(m)  = surfinsw(i)
    5394               surf_lsm_h%rad_sw_out(m) = surfoutsw(i)
    5395               surf_lsm_h%rad_lw_in(m)  = surfinlw(i)
    5396               surf_lsm_h%rad_lw_out(m) = surfoutlw(i)
    5397               surf_lsm_h%rad_net(m)    = surfinsw(i) - surfoutsw(i) +          &
    5398                                          surfinlw(i) - surfoutlw(i)
    5399 !
    5400 !--        northward-facding
    5401            ELSEIF ( surfl(1,i) == inorth_l )  THEN
    5402               surf_lsm_v(0)%rad_sw_in(m)  = surfinsw(i)
    5403               surf_lsm_v(0)%rad_sw_out(m) = surfoutsw(i)
    5404               surf_lsm_v(0)%rad_lw_in(m)  = surfinlw(i)
    5405               surf_lsm_v(0)%rad_lw_out(m) = surfoutlw(i)
    5406               surf_lsm_v(0)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
    5407                                             surfinlw(i) - surfoutlw(i)
    5408 !
    5409 !--        southward-facding
    5410            ELSEIF ( surfl(1,i) == isouth_l )  THEN
    5411               surf_lsm_v(1)%rad_sw_in(m)  = surfinsw(i)
    5412               surf_lsm_v(1)%rad_sw_out(m) = surfoutsw(i)
    5413               surf_lsm_v(1)%rad_lw_in(m)  = surfinlw(i)
    5414               surf_lsm_v(1)%rad_lw_out(m) = surfoutlw(i)
    5415               surf_lsm_v(1)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
    5416                                             surfinlw(i) - surfoutlw(i)
    5417 !
    5418 !--        eastward-facing
    5419            ELSEIF ( surfl(1,i) == ieast_l )  THEN
    5420               surf_lsm_v(2)%rad_sw_in(m)  = surfinsw(i)
    5421               surf_lsm_v(2)%rad_sw_out(m) = surfoutsw(i)
    5422               surf_lsm_v(2)%rad_lw_in(m)  = surfinlw(i)
    5423               surf_lsm_v(2)%rad_lw_out(m) = surfoutlw(i)
    5424               surf_lsm_v(2)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
    5425                                             surfinlw(i) - surfoutlw(i)
    5426 !
    5427 !--        westward-facing
    5428            ELSEIF ( surfl(1,i) == iwest_l )  THEN
    5429               surf_lsm_v(3)%rad_sw_in(m)  = surfinsw(i)
    5430               surf_lsm_v(3)%rad_sw_out(m) = surfoutsw(i)
    5431               surf_lsm_v(3)%rad_lw_in(m)  = surfinlw(i)
    5432               surf_lsm_v(3)%rad_lw_out(m) = surfoutlw(i)
    5433               surf_lsm_v(3)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
    5434                                             surfinlw(i) - surfoutlw(i)
    5435            ENDIF
    5436 
    5437         ENDDO
    5438 
    5439         DO  m = 1, surf_usm_h%ns
    5440            surf_usm_h%surfhf(m) = surf_usm_h%rad_sw_in(m)  +                   &
    5441                                   surf_usm_h%rad_lw_in(m)  -                   &
    5442                                   surf_usm_h%rad_sw_out(m) -                   &
    5443                                   surf_usm_h%rad_lw_out(m)
    5444         ENDDO
    5445         DO  m = 1, surf_lsm_h%ns
    5446            surf_lsm_h%surfhf(m) = surf_lsm_h%rad_sw_in(m)  +                   &
    5447                                   surf_lsm_h%rad_lw_in(m)  -                   &
    5448                                   surf_lsm_h%rad_sw_out(m) -                   &
    5449                                   surf_lsm_h%rad_lw_out(m)
    5450         ENDDO
    5451 
    5452         DO  l = 0, 3
    5453 !--        urban
    5454            DO  m = 1, surf_usm_v(l)%ns
    5455               surf_usm_v(l)%surfhf(m) = surf_usm_v(l)%rad_sw_in(m)  +          &
    5456                                         surf_usm_v(l)%rad_lw_in(m)  -          &
    5457                                         surf_usm_v(l)%rad_sw_out(m) -          &
    5458                                         surf_usm_v(l)%rad_lw_out(m)
    5459            ENDDO
    5460 !--        land
    5461            DO  m = 1, surf_lsm_v(l)%ns
    5462               surf_lsm_v(l)%surfhf(m) = surf_lsm_v(l)%rad_sw_in(m)  +          &
    5463                                         surf_lsm_v(l)%rad_lw_in(m)  -          &
    5464                                         surf_lsm_v(l)%rad_sw_out(m) -          &
    5465                                         surf_lsm_v(l)%rad_lw_out(m)
    5466 
    5467            ENDDO
    5468         ENDDO
    5469 !
    5470 !--     Calculate the average temperature, albedo, and emissivity for urban/land domain
    5471 !--     in case of using average_radiation in the respective radiation model
    5472         IF ( average_radiation )  THEN
    5473 
    5474 !--
    5475 !--        precalculate face areas for different face directions using normal vector
    5476 !--        TODO: make facearea a globale variable because it is used in more than one subroutine
    5477            DO d = 0, nsurf_type
    5478                facearea(d) = 1._wp
    5479                IF ( idir(d) == 0 ) facearea(d) = facearea(d) * dx
    5480                IF ( jdir(d) == 0 ) facearea(d) = facearea(d) * dy
    5481                IF ( kdir(d) == 0 ) facearea(d) = facearea(d) * dz
    5482            ENDDO
    5483 !
    5484 !--        total absorbed SW & LW and emitted LW energy by all physical surfaces (land and urban) in local processor
    5485            pabsswl = 0._wp
    5486            pabslwl = 0._wp
    5487            pemitlwl = 0._wp
    5488            emiss_sum_surfl = 0._wp
    5489            area_surfl = 0._wp
    5490            DO  i = startenergy, endenergy
    5491               d = surfl(id, i)
    5492               pabsswl = pabsswl + (1._wp - albedo_surf(i)) * surfinsw(i) * facearea(d)
    5493               pabslwl = pabslwl + emiss_surf(i) * surfinlw(i) * facearea(d)
    5494               pemitlwl = pemitlwl + surfoutlw(i) * facearea(d)
    5495               emiss_sum_surfl = emiss_sum_surfl + emiss_surf(i) * facearea(d)
    5496               area_surfl = area_surfl + facearea(d)
    5497            END DO
    5498 !
    5499 !--        add the absorbed SW energy by plant canopy
    5500            IF ( plant_canopy )  THEN
    5501               pabsswl = pabsswl + SUM(pcbinsw)
    5502               pabslwl = pabslwl + SUM(pcbinlw)
    5503            ENDIF
    5504 !
    5505 !--        gather all absorbed SW energy in all processors
    5506 #if defined( __parallel )
    5507            CALL MPI_ALLREDUCE( pabsswl, pabssw, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
    5508            CALL MPI_ALLREDUCE( pabslwl, pabslw, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
    5509            CALL MPI_ALLREDUCE( pemitlwl, pemitlw, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
    5510            CALL MPI_ALLREDUCE( emiss_sum_surfl, emiss_sum_surf, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
    5511            CALL MPI_ALLREDUCE( area_surfl, area_surf, 1, MPI_REAL, MPI_SUM, comm2d, ierr )
    5512 #else
    5513            pabssw = pabsswl
    5514            pabslwl = pabslw
    5515            pemitlwl = pemitlw
    5516            emiss_sum_surf = emiss_sum_surfl
    5517            area_surf = area_surfl
    5518 #endif
    5519 !
    5520 !--        total received SW energy in local processor !!!!!! cos??!!!!
    5521            pinswl = 0._wp
    5522            pinlwl = 0._wp
    5523 !-- sky
    5524            DO  i = startsky, endsky
    5525               d = surfl(id, i)
    5526               ii = surfl(ix, i)
    5527               jj = surfl(iy, i)
    5528               pinswl = pinswl + (rad_sw_in_dir(jj,ii) + rad_sw_in_diff(jj,ii)) * facearea(d)
    5529               pinlwl = pinlwl + rad_lw_in_diff(jj,ii) * facearea(d)
    5530            ENDDO
    5531 !-- boundary
    5532            DO  i = startborder, endborder
    5533               d = surfl(id, i)
    5534               ii = surfl(ix, i)
    5535               jj = surfl(iy, i)
    5536               pinswl = pinswl + (rad_sw_in_dir(jj,ii) + rad_sw_in_diff(jj,ii)) * facearea(d)
    5537               pinlwl = pinlwl + rad_lw_in_diff(jj,ii) * facearea(d)
    5538            ENDDO
    5539 !--        gather all received SW energy in all processors
    5540 #if defined( __parallel )
    5541            CALL MPI_ALLREDUCE( pinswl, pinsw, 1, MPI_REAL, MPI_SUM, comm2d, ierr)
    5542            CALL MPI_ALLREDUCE( pinlwl, pinlw, 1, MPI_REAL, MPI_SUM, comm2d, ierr)
    5543 #else
    5544            pinsw = pinswl
    5545            pinlw = pinlwl
    5546 #endif
    5547 !--        (1) albedo
    5548            IF ( pinsw /= 0.0_wp )  albedo_urb = 1._wp - pabssw / pinsw
    5549        
    5550 !--        (2) average emmsivity
    5551            emissivity_urb = emiss_sum_surf / area_surf
    5552 
    5553 !--        (3) temerature
    5554            t_rad_urb = ((pemitlw - pabslw + emissivity_urb*pinlw)/(emissivity_urb*sigma_sb*area_surf))**0.25_wp
    5555 
    5556         ENDIF
    5557        
    5558 !--     return surface radiation to horizontal surfaces
    5559 !--     to rad_sw_in, rad_lw_in and rad_net for outputs
    5560         !!!!!!!!!!
    5561 !--     we need the original radiation on urban top layer
    5562 !--     for calculation of MRT so we can't do adjustment here for now
    5563         !!!!!!!!!!
    5564         !!!DO isurf = 1, nsurfl
    5565         !!!    i = surfl(ix,isurf)
    5566         !!!    j = surfl(iy,isurf)
    5567         !!!    k = surfl(iz,isurf)
    5568         !!!    d = surfl(id,isurf)
    5569         !!!    IF ( d==iroof )  THEN
    5570         !!!        rad_sw_in(:,j,i) = surfinsw(isurf)
    5571         !!!        rad_lw_in(:,j,i) = surfinlw(isurf)
    5572         !!!        rad_net(j,i) = rad_sw_in(k,j,i) - rad_sw_out(k,j,i) + rad_lw_in(k,j,i) - rad_lw_out(k,j,i)
    5573         !!!    ENDIF
    5574         !!!ENDDO
    5575 
    5576       CONTAINS
    5577 
    5578 !------------------------------------------------------------------------------!
    5579 ! Description:
    5580 ! ------------
    5581 !> This subroutine splits direct and diffusion dw radiation
    5582 !> It sould not be called in case the radiation model already does it
    5583 !> It follows <CITATION>
    5584 !------------------------------------------------------------------------------!
    5585         SUBROUTINE calc_diffusion_radiation
    5586 
    5587           USE date_and_time_mod,                                               &
    5588               ONLY:  day_of_year_init, time_utc_init
    5589          
    5590           REAL(wp), PARAMETER                          ::  sol_const = 1367.0_wp   !< solar conbstant
    5591           REAL(wp), PARAMETER                          ::  lowest_solarUp = 0.1_wp  !< limit the sun elevation to protect stability of the calculation
    5592           INTEGER(iwp)                                 ::  i, j
    5593           REAL(wp), PARAMETER                          ::  year_seconds = 86400._wp * 365._wp
    5594           REAL(wp)                                     ::  year_angle              !< angle
    5595           REAL(wp)                                     ::  etr                     !< extraterestrial radiation
    5596           REAL(wp)                                     ::  corrected_solarUp       !< corrected solar up radiation
    5597           REAL(wp)                                     ::  horizontalETR           !< horizontal extraterestrial radiation
    5598           REAL(wp)                                     ::  clearnessIndex          !< clearness index
    5599           REAL(wp)                                     ::  diff_frac               !< diffusion fraction of the radiation
    5600          
    5601        
    5602 !--     Calculate current day and time based on the initial values and simulation time
    5603           year_angle = ((day_of_year_init*86400)                               &
    5604                                   +  time_utc_init+time_since_reference_point) &
    5605                                   /  year_seconds * 2.0_wp * pi
    5606          
    5607           etr = sol_const * (1.00011_wp +                               &
    5608                0.034221_wp * cos(year_angle) +                          &
    5609                0.001280_wp * sin(year_angle) +                          &
    5610                0.000719_wp * cos(2.0_wp * year_angle) +                 &
    5611                0.000077_wp * sin(2.0_wp * year_angle))
    5612          
    5613 !--   
    5614 !--     Under a very low angle, we keep extraterestrial radiation at
    5615 !--     the last small value, therefore the clearness index will be pushed
    5616 !--     towards 0 while keeping full continuity.
    5617 !--   
    5618           IF ( zenith(0) <= lowest_solarUp )  THEN
    5619              corrected_solarUp = lowest_solarUp
    5620           ELSE
    5621              corrected_solarUp = zenith(0)
    5622           ENDIF
    5623          
    5624           horizontalETR = etr * corrected_solarUp
    5625          
    5626           DO i = nxl, nxr
    5627              DO j = nys, nyn
    5628 
    5629                 DO  m = surf_def_h(0)%start_index(j,i),                        &
    5630                         surf_def_h(0)%end_index(j,i)
    5631                    clearnessIndex = surf_def_h(0)%rad_sw_in(m) / horizontalETR
    5632                    diff_frac      = 1.0_wp /                                   &
    5633                         (1.0_wp + exp(-5.0033_wp + 8.6025_wp * clearnessIndex))
    5634                    rad_sw_in_diff(j,i) = surf_def_h(0)%rad_sw_in(m) * diff_frac
    5635                    rad_sw_in_dir(j,i)  = surf_def_h(0)%rad_sw_in(m) *          &
    5636                                             (1.0_wp - diff_frac)
    5637                    rad_lw_in_diff(j,i) = surf_def_h(0)%rad_lw_in(m)
    5638                 ENDDO
    5639                 DO  m = surf_lsm_h%start_index(j,i),                           &
    5640                         surf_lsm_h%end_index(j,i)
    5641                    clearnessIndex = surf_lsm_h%rad_sw_in(m) / horizontalETR
    5642                    diff_frac      = 1.0_wp /                                   &
    5643                         (1.0_wp + exp(-5.0033_wp + 8.6025_wp * clearnessIndex))
    5644                    rad_sw_in_diff(j,i) = surf_lsm_h%rad_sw_in(m) * diff_frac
    5645                    rad_sw_in_dir(j,i)  = surf_lsm_h%rad_sw_in(m) *             &
    5646                                             (1.0_wp - diff_frac)
    5647                    rad_lw_in_diff(j,i) = surf_lsm_h%rad_lw_in(m)
    5648                 ENDDO
    5649                 DO  m = surf_usm_h%start_index(j,i),                           &
    5650                         surf_usm_h%end_index(j,i)
    5651                    clearnessIndex = surf_usm_h%rad_sw_in(m) / horizontalETR
    5652                    diff_frac      = 1.0_wp /                                   &
    5653                         (1.0_wp + exp(-5.0033_wp + 8.6025_wp * clearnessIndex))
    5654                    rad_sw_in_diff(j,i) = surf_usm_h%rad_sw_in(m) * diff_frac
    5655                    rad_sw_in_dir(j,i)  = surf_usm_h%rad_sw_in(m) *             &
    5656                                             (1.0_wp - diff_frac)
    5657                    rad_lw_in_diff(j,i) = surf_usm_h%rad_lw_in(m)
    5658                 ENDDO
    5659              ENDDO
    5660           ENDDO
    5661          
    5662         END SUBROUTINE calc_diffusion_radiation
    5663 
    5664 !------------------------------------------------------------------------------!
    5665 !> Finds first model boundary crossed by a ray
    5666 !------------------------------------------------------------------------------!
    5667         PURE SUBROUTINE find_boundary_face(origin, uvect, bdycross)
    5668          
    5669           IMPLICIT NONE
    5670          
    5671           INTEGER(iwp) ::  d       !<
    5672           INTEGER(iwp) ::  seldim  !< found fist crossing index
    5673          
    5674           INTEGER(iwp), DIMENSION(3)              ::  bdyd      !< boundary direction       
    5675           INTEGER(iwp), DIMENSION(4), INTENT(out) ::  bdycross  !< found boundary crossing (d, z, y, x)
    5676          
    5677           REAL(wp)                                ::  bdydim  !<
    5678           REAL(wp)                                ::  dist    !<
    5679          
    5680           REAL(wp), DIMENSION(3)             ::  crossdist  !< crossing distance
    5681           REAL(wp), DIMENSION(3), INTENT(in) ::  origin     !< ray origin
    5682           REAL(wp), DIMENSION(3), INTENT(in) ::  uvect      !< ray unit vector
    5683          
    5684          
    5685           bdydim       = nzut + .5_wp  !< top boundary
    5686           bdyd(1)      = isky
    5687           crossdist(1) = ( bdydim - origin(1) ) / uvect(1)  !< subroutine called only when uvect(1)>0
    5688          
    5689           IF ( uvect(2) == 0._wp )  THEN
    5690              crossdist(2) = huge(1._wp)
    5691           ELSE
    5692              IF ( uvect(2) >= 0._wp )  THEN
    5693                 bdydim  = ny + .5_wp  !< north global boundary
    5694                 bdyd(2) = inorth_b
    5695              ELSE
    5696                 bdydim  = -.5_wp  !< south global boundary
    5697                 bdyd(2) = isouth_b
    5698              ENDIF
    5699              crossdist(2) = ( bdydim - origin(2) ) / uvect(2)
    5700           ENDIF
    5701          
    5702           IF ( uvect(3) == 0._wp )  THEN
    5703              crossdist(3) = huge(1._wp)
    5704           ELSE
    5705              IF ( uvect(3) >= 0._wp )  THEN
    5706                 bdydim  = nx + .5_wp  !< east global boundary
    5707                 bdyd(3) = ieast_b
    5708              ELSE
    5709                 bdydim  = -.5_wp  !< west global boundary
    5710                 bdyd(3) = iwest_b
    5711              ENDIF
    5712              crossdist(3) = ( bdydim - origin(3) ) / uvect(3)
    5713           ENDIF
    5714          
    5715           seldim = minloc(crossdist, 1)
    5716           dist   = crossdist(seldim)
    5717           d      = bdyd(seldim)
    5718          
    5719           bdycross(1)   = d
    5720           bdycross(2:4) = NINT( origin(:) + uvect(:) * dist &
    5721                + .5_wp * (/ kdir(d), jdir(d), idir(d) /) )
    5722          
    5723         END SUBROUTINE find_boundary_face
    5724 !------------------------------------------------------------------------------!
    5725 !> Calculates radiation absorbed by box with given size and LAD.
    5726 !>
    5727 !> Simulates resol**2 rays (by equally spacing a bounding horizontal square
    5728 !> conatining all possible rays that would cross the box) and calculates
    5729 !> average transparency per ray. Returns fraction of absorbed radiation flux
    5730 !> and area for which this fraction is effective.
    5731 !------------------------------------------------------------------------------!
    5732         PURE SUBROUTINE box_absorb(boxsize, resol, dens, uvec, area, absorb)
    5733           IMPLICIT NONE
    5734          
    5735           REAL(wp), DIMENSION(3), INTENT(in) :: &
    5736                boxsize, &      !< z, y, x size of box in m
    5737                uvec            !< z, y, x unit vector of incoming flux
    5738           INTEGER(iwp), INTENT(in) :: &
    5739                resol           !< No. of rays in x and y dimensions
    5740           REAL(wp), INTENT(in) :: &
    5741                dens            !< box density (e.g. Leaf Area Density)
    5742           REAL(wp), INTENT(out) :: &
    5743                area, &         !< horizontal area for flux absorbtion
    5744                absorb          !< fraction of absorbed flux
    5745           REAL(wp) :: &
    5746                xshift, yshift, &
    5747                xmin, xmax, ymin, ymax, &
    5748                xorig, yorig, &
    5749                dx1, dy1, dz1, dx2, dy2, dz2, &
    5750                crdist, &
    5751                transp
    5752           INTEGER(iwp) :: &
    5753                i, j
    5754          
    5755           xshift = uvec(3) / uvec(1) * boxsize(1)
    5756           xmin = min(0._wp, -xshift)
    5757           xmax = boxsize(3) + max(0._wp, -xshift)
    5758           yshift = uvec(2) / uvec(1) * boxsize(1)
    5759           ymin = min(0._wp, -yshift)
    5760           ymax = boxsize(2) + max(0._wp, -yshift)
    5761          
    5762           transp = 0._wp
    5763           DO i = 1, resol
    5764              xorig = xmin + (xmax-xmin) * (i-.5_wp) / resol
    5765              DO j = 1, resol
    5766                 yorig = ymin + (ymax-ymin) * (j-.5_wp) / resol
    5767                
    5768                 dz1 = 0._wp
    5769                 dz2 = boxsize(1)/uvec(1)
    5770                
    5771                 IF ( uvec(2) > 0._wp )  THEN
    5772                    dy1 = -yorig             / uvec(2) !< crossing with y=0
    5773                    dy2 = (boxsize(2)-yorig) / uvec(2) !< crossing with y=boxsize(2)
    5774                 ELSE IF ( uvec(2) < 0._wp )  THEN
    5775                    dy1 = (boxsize(2)-yorig) / uvec(2) !< crossing with y=boxsize(2)
    5776                    dy2 = -yorig             / uvec(2) !< crossing with y=0
    5777                 ELSE !uvec(2)==0
    5778                    dy1 = -huge(1._wp)
    5779                    dy2 = huge(1._wp)
    5780                 ENDIF
    5781                
    5782                 IF ( uvec(3) > 0._wp )  THEN
    5783                    dx1 = -xorig             / uvec(3) !< crossing with x=0
    5784                    dx2 = (boxsize(3)-xorig) / uvec(3) !< crossing with x=boxsize(3)
    5785                 ELSE IF ( uvec(3) < 0._wp )  THEN
    5786                    dx1 = (boxsize(3)-xorig) / uvec(3) !< crossing with x=boxsize(3)
    5787                    dx2 = -xorig             / uvec(3) !< crossing with x=0
    5788                 ELSE !uvec(1)==0
    5789                    dx1 = -huge(1._wp)
    5790                    dx2 = huge(1._wp)
    5791                 ENDIF
    5792                
    5793                 crdist = max(0._wp, (min(dz2, dy2, dx2) - max(dz1, dy1, dx1)))
    5794                 transp = transp + exp(-ext_coef * dens * crdist)
    5795              ENDDO
    5796           ENDDO
    5797           transp = transp / resol**2
    5798           area = (boxsize(3)+xshift)*(boxsize(2)+yshift)
    5799           absorb = 1._wp - transp
    5800          
    5801         END SUBROUTINE box_absorb
    5802 
    5803        
    5804     END SUBROUTINE radiation_interaction
    5805 
    5806 
    5807 !------------------------------------------------------------------------------!
    5808 ! Description:
    5809 ! ------------
    5810 !> Calculates shape view factors SVF and plant sink canopy factors PSCF
    5811 !> !!!!!DESCRIPTION!!!!!!!!!!
    5812 !------------------------------------------------------------------------------!
    5813     SUBROUTINE radiation_calc_svf
    5814 
    5815         IMPLICIT NONE
    5816        
    5817         INTEGER(iwp)                                :: i, j, k, l, d, ip, jp
    5818         INTEGER(iwp)                                :: isvf, ksvf, icsf, kcsf, npcsfl, isvf_surflt, imrtt, imrtf
    5819         INTEGER(iwp)                                :: sd, td, ioln, iproc
    5820         REAL(wp),     DIMENSION(0:nsurf_type)       :: facearea
    5821         INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE   :: nzterrl, planthl
    5822         REAL(wp),     DIMENSION(:,:), ALLOCATABLE   :: csflt, pcsflt
    5823         INTEGER(iwp), DIMENSION(:,:), ALLOCATABLE   :: kcsflt,kpcsflt
    5824         INTEGER(iwp), DIMENSION(:), ALLOCATABLE     :: icsflt,dcsflt,ipcsflt,dpcsflt
    5825         REAL(wp), DIMENSION(3)                      :: uv
    5826         LOGICAL                                     :: visible
    5827         REAL(wp), DIMENSION(3)                      :: sa, ta          !< real coordinates z,y,x of source and target
    5828         REAL(wp)                                    :: transparency, rirrf, sqdist, svfsum
    5829         INTEGER(iwp)                                :: isurflt, isurfs, isurflt_prev
    5830         INTEGER(iwp)                                :: itx, ity, itz
    5831         CHARACTER(len=7)                            :: pid_char = ''
    5832         INTEGER(iwp)                                :: win_lad, minfo
    5833         REAL(wp), DIMENSION(:,:,:), POINTER         :: lad_s_rma       !< fortran pointer, but lower bounds are 1
    5834         TYPE(c_ptr)                                 :: lad_s_rma_p     !< allocated c pointer
    5835 #if defined( __parallel )
    5836         INTEGER(kind=MPI_ADDRESS_KIND)              :: size_lad_rma
    5837 #endif
    5838         REAL(wp), DIMENSION(0:nsurf_type)           :: svf_threshold   !< threshold to ignore very small svf between far surfaces
    5839        
    58405596!   
     5597        INTEGER(iwp), DIMENSION(0:svfnorm_report_num) :: svfnorm_counts
     5598        CHARACTER(200)                                :: msg
     5599
    58415600!--     calculation of the SVF
    58425601        CALL location_message( '    calculation of SVF and CSF', .TRUE. )
    5843         CALL cpu_log( log_point_s(79), 'radiation_calc_svf', 'start' )
    5844 !
     5602!         CALL radiation_write_debug_log('Start calculation of SVF and CSF')
     5603
    58455604!--     precalculate face areas for different face directions using normal vector
    58465605        DO d = 0, nsurf_type
     
    58515610        ENDDO
    58525611
    5853 !--     calculate the svf threshold
    5854         svf_threshold = 0._wp
    5855         IF ( dist_max_svf > 0._wp ) THEN
    5856             DO d = 0, nsurf_type
    5857                sqdist = dist_max_svf * dist_max_svf
    5858                svf_threshold(d) = 1._wp / (pi * sqdist) * facearea(d)
    5859             ENDDO
    5860          ENDIF
    5861          
    58625612!--     initialize variables and temporary arrays for calculation of svf and csf
    58635613        nsvfl  = 0
     
    58735623            acsf => acsf1
    58745624        ENDIF
     5625        ray_skip_maxdist = 0
     5626        ray_skip_minval = 0
    58755627       
    58765628!--     initialize temporary terrain and plant canopy height arrays (global 2D array!)
     
    58885640            ALLOCATE( plantt(0:(nx+1)*(ny+1)-1) )
    58895641            maxboxesg = nx + ny + nzu + 1
     5642            max_track_len = nx + ny + 1
    58905643!--         temporary arrays storing values for csf calculation during raytracing
    58915644            ALLOCATE( boxes(3, maxboxesg) )
     
    58935646
    58945647#if defined( __parallel )
    5895             ALLOCATE( planthl(nys:nyn,nxl:nxr) )
    5896             planthl = pch(nys:nyn,nxl:nxr)
    5897        
    5898             CALL MPI_AllGather( planthl, nnx*nny, MPI_INTEGER, &
     5648            CALL MPI_AllGather( pct, nnx*nny, MPI_INTEGER, &
    58995649                                plantt, nnx*nny, MPI_INTEGER, comm2d, ierr )
    5900             DEALLOCATE( planthl )
    59015650           
    59025651!--         temporary arrays storing values for csf calculation during raytracing
     
    59045653            ALLOCATE( lad_disp(maxboxesg) )
    59055654
    5906             IF ( usm_lad_rma )  THEN
     5655            IF ( rma_lad_raytrace )  THEN
    59075656                ALLOCATE( lad_s_ray(maxboxesg) )
    59085657               
     
    59225671                                        lad_s_rma_p, win_lad, ierr)
    59235672                CALL c_f_pointer(lad_s_rma_p, lad_s_rma, (/ nzu, nny, nnx /))
    5924                 usm_lad(nzub:, nys:, nxl:) => lad_s_rma(:,:,:)
     5673                sub_lad(nzub:, nys:, nxl:) => lad_s_rma(:,:,:)
    59255674            ELSE
    5926                 ALLOCATE(usm_lad(nzub:nzut, nys:nyn, nxl:nxr))
     5675                ALLOCATE(sub_lad(nzub:nzut, nys:nyn, nxl:nxr))
    59275676            ENDIF
    59285677#else
    59295678            plantt = RESHAPE( pct(nys:nyn,nxl:nxr), (/(nx+1)*(ny+1)/) )
    5930             ALLOCATE(usm_lad(nzub:nzut, nys:nyn, nxl:nxr))
     5679            ALLOCATE(sub_lad(nzub:nzut, nys:nyn, nxl:nxr))
    59315680#endif
    5932             usm_lad(:,:,:) = 0._wp
     5681            plantt_max = MAXVAL(plantt)
     5682            ALLOCATE( rt2_track(2, max_track_len), rt2_track_lad(nzub:plantt_max, max_track_len), &
     5683                      rt2_track_dist(0:max_track_len), rt2_dist(plantt_max-nzub+2) )
     5684
     5685            sub_lad(:,:,:) = 0._wp
    59335686            DO i = nxl, nxr
    59345687                DO j = nys, nyn
    59355688                    k = get_topography_top_index_ji( j, i, 's' )
    59365689
    5937                     usm_lad(k:nzut, j, i) = lad_s(0:nzut-k, j, i)
     5690                    sub_lad(k:nzut, j, i) = lad_s(0:nzut-k, j, i)
    59385691                ENDDO
    59395692            ENDDO
    59405693
    59415694#if defined( __parallel )
    5942             IF ( usm_lad_rma )  THEN
     5695            IF ( rma_lad_raytrace )  THEN
    59435696                CALL MPI_Info_free(minfo, ierr)
    59445697                CALL MPI_Win_lock_all(0, win_lad, ierr)
    59455698            ELSE
    5946                 ALLOCATE( usm_lad_g(0:(nx+1)*(ny+1)*nzu-1) )
    5947                 CALL MPI_AllGather( usm_lad, nnx*nny*nzu, MPI_REAL, &
    5948                                     usm_lad_g, nnx*nny*nzu, MPI_REAL, comm2d, ierr )
     5699                ALLOCATE( sub_lad_g(0:(nx+1)*(ny+1)*nzu-1) )
     5700                CALL MPI_AllGather( sub_lad, nnx*nny*nzu, MPI_REAL, &
     5701                                    sub_lad_g, nnx*nny*nzu, MPI_REAL, comm2d, ierr )
    59495702            ENDIF
    59505703#endif
     
    60115764            CLOSE(154)
    60125765        ENDIF  !< mrt_factors
    6013 
     5766       
     5767        !--Directions opposite to face normals are not even calculated,
     5768        !--they must be preset to 0
     5769        !--
     5770        dsitrans(:,:) = 0._wp
    60145771       
    60155772        DO isurflt = 1, nsurfl
    60165773!--         determine face centers
    60175774            td = surfl(id, isurflt)
    6018             IF ( td >= isky  .AND.  .NOT.  plant_canopy ) CYCLE
    60195775            ta = (/ REAL(surfl(iz, isurflt), wp) - 0.5_wp * kdir(td),  &
    60205776                      REAL(surfl(iy, isurflt), wp) - 0.5_wp * jdir(td),  &
    60215777                      REAL(surfl(ix, isurflt), wp) - 0.5_wp * idir(td)  /)
     5778
     5779            !--Calculate sky view factor and raytrace DSI paths
     5780            skyvf(isurflt) = 0._wp
     5781            skyvft(isurflt) = 0._wp
     5782
     5783            !--Select a proper half-sphere for 2D raytracing
     5784            SELECT CASE ( td )
     5785               CASE ( iup_u, iup_l )
     5786                  az0 = 0._wp
     5787                  naz = raytrace_discrete_azims
     5788                  azs = 2._wp * pi / REAL(naz, wp)
     5789                  zn0 = 0._wp
     5790                  nzn = raytrace_discrete_elevs / 2
     5791                  zns = pi / 2._wp / REAL(nzn, wp)
     5792               CASE ( isouth_u, isouth_l )
     5793                  az0 = pi / 2._wp
     5794                  naz = raytrace_discrete_azims / 2
     5795                  azs = pi / REAL(naz, wp)
     5796                  zn0 = 0._wp
     5797                  nzn = raytrace_discrete_elevs
     5798                  zns = pi / REAL(nzn, wp)
     5799               CASE ( inorth_u, inorth_l )
     5800                  az0 = - pi / 2._wp
     5801                  naz = raytrace_discrete_azims / 2
     5802                  azs = pi / REAL(naz, wp)
     5803                  zn0 = 0._wp
     5804                  nzn = raytrace_discrete_elevs
     5805                  zns = pi / REAL(nzn, wp)
     5806               CASE ( iwest_u, iwest_l )
     5807                  az0 = pi
     5808                  naz = raytrace_discrete_azims / 2
     5809                  azs = pi / REAL(naz, wp)
     5810                  zn0 = 0._wp
     5811                  nzn = raytrace_discrete_elevs
     5812                  zns = pi / REAL(nzn, wp)
     5813               CASE ( ieast_u, ieast_l )
     5814                  az0 = 0._wp
     5815                  naz = raytrace_discrete_azims / 2
     5816                  azs = pi / REAL(naz, wp)
     5817                  zn0 = 0._wp
     5818                  nzn = raytrace_discrete_elevs
     5819                  zns = pi / REAL(nzn, wp)
     5820               CASE DEFAULT
     5821                  WRITE(message_string, *) 'ERROR: the surface type ',td , ' is not supported for calculating SVF'
     5822                  CALL message( 'radiation_calc_svf', 'PA0XXX', 1, 2, 0, 6, 0 )
     5823            END SELECT
     5824
     5825            ALLOCATE ( zdirs(1:nzn), zbdry(0:nzn), vffrac(1:nzn), ztransp(1:nzn) )
     5826            zdirs(:) = (/( TAN(pi/2 - (zn0+(REAL(izn,wp)-.5_wp)*zns)), izn=1, nzn )/)
     5827            zbdry(:) = (/( zn0+REAL(izn,wp)*zns, izn=0, nzn )/)
     5828            IF ( td == iup_u  .OR.  td == iup_l )  THEN
     5829               !-- For horizontal target, vf fractions are constant per azimuth
     5830               vffrac(:) = (COS(2 * zbdry(0:nzn-1)) - COS(2 * zbdry(1:nzn))) / 2._wp / REAL(naz, wp)
     5831               !--sum of vffrac for all iaz equals 1, verified
     5832            ENDIF
     5833
     5834            !--Calculate sky-view factor and direct solar visibility using 2D raytracing
     5835            DO iaz = 1, naz
     5836               azmid = az0 + (REAL(iaz, wp) - .5_wp) * azs
     5837               IF ( td /= iup_u  .AND.  td /= iup_l )  THEN
     5838                  az2 = REAL(iaz, wp) * azs - pi/2._wp
     5839                  az1 = az2 - azs
     5840                  !TODO precalculate after 1st line
     5841                  vffrac(:) = (SIN(az2) - SIN(az1))                           &
     5842                              * (zbdry(1:nzn) - zbdry(0:nzn-1)                &
     5843                                 + SIN(zbdry(0:nzn-1))*COS(zbdry(0:nzn-1))    &
     5844                                 - SIN(zbdry(1:nzn))*COS(zbdry(1:nzn)))       &
     5845                              / (2._wp * pi)
     5846                  !--sum of vffrac for all iaz equals 1, verified
     5847               ENDIF
     5848               CALL raytrace_2d(ta, (/ COS(azmid), SIN(azmid) /), zdirs,      &
     5849                                    surfstart(myid) + isurflt, facearea(td),  &
     5850                                    vffrac, .TRUE., .FALSE., win_lad, horizon,&
     5851                                    ztransp) !FIXME unit vect in grid units + zdirs
     5852
     5853               azen = pi/2 - ATAN(horizon)
     5854               IF ( td == iup_u  .OR.  td == iup_l )  THEN
     5855                  azen = MIN(azen, pi/2) !only above horizontal direction
     5856                  skyvf(isurflt) = skyvf(isurflt) + (1._wp - COS(2*azen)) /   &
     5857                     (2._wp * raytrace_discrete_azims)
     5858               ELSE
     5859                  skyvf(isurflt) = skyvf(isurflt) + (SIN(az2) - SIN(az1)) *   &
     5860                              (azen - SIN(azen)*COS(azen)) / (2._wp*pi)
     5861               ENDIF
     5862               skyvft(isurflt) = skyvft(isurflt) + SUM(ztransp(:) * vffrac(:))
     5863 
     5864               !--Save direct solar transparency
     5865               j = MODULO(NINT(azmid/                                          &
     5866                               (2._wp*pi)*raytrace_discrete_azims-.5_wp, iwp), &
     5867                          raytrace_discrete_azims)
     5868
     5869               DO k = 1, raytrace_discrete_elevs/2
     5870                  i = dsidir_rev(k-1, j)
     5871                  IF ( i /= -1 )  dsitrans(isurflt, i) = ztransp(k)
     5872               ENDDO
     5873            ENDDO
     5874
     5875            DEALLOCATE ( zdirs, zbdry, vffrac, ztransp )
     5876
    60225877            DO isurfs = 1, nsurf
    6023 !--             cycle for atmospheric surfaces since they are not source surfaces
    6024                 sd = surf(id, isurfs)
    6025                 IF ( sd > iwest_l  .AND.  sd < isky ) CYCLE
    6026 !--             if reflections between target surfaces (urban and land) are neglected (surf_reflection set to
    6027 !--             FALSE) cycle. This will reduce the number of SVFs and keep SVFs between only ertual surfaces to
    6028 !--             physical surfaces
    6029                 IF ( .NOT.  surf_reflections  .AND. sd < isky ) CYCLE
    6030 !--             cycle if the target and the source surfaces are not facing each other
    60315878                IF ( .NOT.  surface_facing(surfl(ix, isurflt), surfl(iy, isurflt), &
    60325879                    surfl(iz, isurflt), surfl(id, isurflt), &
     
    60365883                ENDIF
    60375884                 
     5885                sd = surf(id, isurfs)
    60385886                sa = (/ REAL(surf(iz, isurfs), wp) - 0.5_wp * kdir(sd),  &
    60395887                        REAL(surf(iy, isurfs), wp) - 0.5_wp * jdir(sd),  &
     
    60445892                sqdist = SUM(uv(:)**2)
    60455893                uv = uv / SQRT(sqdist)
     5894
     5895!--             reject raytracing above max distance
     5896                IF ( SQRT(sqdist) > max_raytracing_dist ) THEN
     5897                    ray_skip_maxdist = ray_skip_maxdist + 1
     5898                    CYCLE
     5899                ENDIF
    60465900               
    60475901!--             irradiance factor (our unshaded shape view factor) = view factor per differential target area * source area
     
    60515905                    * facearea(sd)
    60525906
    6053 !--             skip svf less than svf_threshold
    6054                 IF ( rirrf < svf_threshold(sd) .AND.  sd < isky ) CYCLE
     5907!--             reject raytracing for potentially too small view factor values
     5908                IF ( rirrf < min_irrf_value ) THEN
     5909                    ray_skip_minval = ray_skip_minval + 1
     5910                    CYCLE
     5911                ENDIF
    60555912
    60565913!--             raytrace + process plant canopy sinks within
    60575914                CALL raytrace(sa, ta, isurfs, rirrf, facearea(td), .TRUE., &
    60585915                        visible, transparency, win_lad)
    6059                
     5916
    60605917                IF ( .NOT.  visible ) CYCLE
    6061                 IF ( td >= isky ) CYCLE !< we calculated these only for raytracing
    6062                                         !< to find plant canopy sinks, we don't need svf for them
     5918                ! rsvf = rirrf * transparency
    60635919
    60645920!--             write to the svf array
     
    60805936                        DEALLOCATE( asvf1 )
    60815937                    ENDIF
     5938
     5939!                     WRITE(msg,'(A,3I12)') 'Grow asvf:',nsvfl,nsvfla,k
     5940!                     CALL radiation_write_debug_log( msg )
     5941                   
    60825942                    nsvfla = k
    60835943                ENDIF
     
    60905950        ENDDO
    60915951
     5952        !--Raytrace to canopy boxes to fill dsitransc TODO optimize
     5953        !--
     5954        dsitransc(:,:) = -999._wp !FIXME
     5955        az0 = 0._wp
     5956        naz = raytrace_discrete_azims
     5957        azs = 2._wp * pi / REAL(naz, wp)
     5958        zn0 = 0._wp
     5959        nzn = raytrace_discrete_elevs / 2
     5960        zns = pi / 2._wp / REAL(nzn, wp)
     5961        ALLOCATE ( zdirs(1:nzn), vffrac(1:nzn), ztransp(1:nzn) )
     5962        zdirs(:) = (/( TAN(pi/2 - (zn0+(REAL(izn,wp)-.5_wp)*zns)), izn=1, nzn )/)
     5963        vffrac(:) = 0._wp
     5964
     5965        DO ipcgb = 1, npcbl
     5966           ta = (/ REAL(pcbl(iz, ipcgb), wp),  &
     5967                   REAL(pcbl(iy, ipcgb), wp),  &
     5968                   REAL(pcbl(ix, ipcgb), wp) /)
     5969           !--Calculate sky-view factor and direct solar visibility using 2D raytracing
     5970           DO iaz = 1, naz
     5971              azmid = az0 + (REAL(iaz, wp) - .5_wp) * azs
     5972              CALL raytrace_2d(ta, (/ COS(azmid), SIN(azmid) /), zdirs,     &
     5973                                   -999, -999._wp, vffrac, .FALSE., .TRUE., &
     5974                                   win_lad, horizon, ztransp) !FIXME unit vect in grid units + zdirs
     5975
     5976              !--Save direct solar transparency
     5977              j = MODULO(NINT(azmid/                                         &
     5978                             (2._wp*pi)*raytrace_discrete_azims-.5_wp, iwp), &
     5979                         raytrace_discrete_azims)
     5980              DO k = 1, raytrace_discrete_elevs/2
     5981                 i = dsidir_rev(k-1, j)
     5982                 IF ( i /= -1 )  dsitransc(ipcgb, i) = ztransp(k)
     5983              ENDDO
     5984           ENDDO
     5985        ENDDO
     5986        DEALLOCATE ( zdirs, vffrac, ztransp )
     5987
     5988!         CALL radiation_write_debug_log( 'End of calculation SVF' )
     5989!         WRITE(msg, *) 'Raytracing skipped for maximum distance of ', &
     5990!             max_raytracing_dist, ' m on ', ray_skip_maxdist, ' pairs.'
     5991!         CALL radiation_write_debug_log( msg )
     5992!         WRITE(msg, *) 'Raytracing skipped for minimum potential value of ', &
     5993!             min_irrf_value , ' on ', ray_skip_minval, ' pairs.'
     5994!         CALL radiation_write_debug_log( msg )
     5995
    60925996        CALL location_message( '    waiting for completion of SVF and CSF calculation in all processes', .TRUE. )
    60935997!--     deallocate temporary global arrays
     
    60976001!--         finalize mpi_rma communication and deallocate temporary arrays
    60986002#if defined( __parallel )
    6099             IF ( usm_lad_rma )  THEN
     6003            IF ( rma_lad_raytrace )  THEN
    61006004                CALL MPI_Win_flush_all(win_lad, ierr)
    61016005!--             unlock MPI window
     
    61066010!--             deallocate temporary arrays storing values for csf calculation during raytracing
    61076011                DEALLOCATE( lad_s_ray )
    6108 !--             usm_lad is the pointer to lad_s_rma in case of usm_lad_rma
     6012!--             sub_lad is the pointer to lad_s_rma in case of rma_lad_raytrace
    61096013!--             and must not be deallocated here
    61106014            ELSE
    6111                 DEALLOCATE(usm_lad)
    6112                 DEALLOCATE(usm_lad_g)
     6015                DEALLOCATE(sub_lad)
     6016                DEALLOCATE(sub_lad_g)
    61136017            ENDIF
    61146018#else
    6115             DEALLOCATE(usm_lad)
     6019            DEALLOCATE(sub_lad)
    61166020#endif
    61176021            DEALLOCATE( boxes )
    61186022            DEALLOCATE( crlens )
    61196023            DEALLOCATE( plantt )
     6024            DEALLOCATE( rt2_track, rt2_track_lad, rt2_track_dist, rt2_dist )
    61206025        ENDIF
    61216026
    61226027        CALL location_message( '    calculation of the complete SVF array', .TRUE. )
    61236028
     6029!         CALL radiation_write_debug_log( 'Start SVF sort' )
    61246030!--     sort svf ( a version of quicksort )
    61256031        CALL quicksort_svf(asvf,1,nsvfl)
    61266032
     6033        !< load svf from the structure array to plain arrays
     6034!         CALL radiation_write_debug_log( 'Load svf from the structure array to plain arrays' )
    61276035        ALLOCATE( svf(ndsvf,nsvfl) )
    61286036        ALLOCATE( svfsurf(idsvf,nsvfl) )
    6129 
    6130         !< load svf from the structure array to plain arrays
     6037        svfnorm_counts(:) = 0._wp
    61316038        isurflt_prev = -1
    61326039        ksvf = 1
     
    61366043            IF ( asvf(ksvf)%isurflt /= isurflt_prev )  THEN
    61376044                IF ( isurflt_prev /= -1  .AND.  svfsum /= 0._wp )  THEN
    6138 !--                 TODO detect and log when normalization differs too much from 1
    6139                     svf(1, isvf_surflt:isvf-1) = svf(1, isvf_surflt:isvf-1) / svfsum
     6045                    !< update histogram of logged svf normalization values
     6046                    i = searchsorted(svfnorm_report_thresh, svfsum / (1._wp-skyvf(isurflt_prev)))
     6047                    svfnorm_counts(i) = svfnorm_counts(i) + 1
     6048
     6049                    svf(1, isvf_surflt:isvf-1) = svf(1, isvf_surflt:isvf-1) / svfsum * (1._wp-skyvf(isurflt_prev))
    61406050                ENDIF
    61416051                isurflt_prev = asvf(ksvf)%isurflt
     
    61546064
    61556065        IF ( isurflt_prev /= -1  .AND.  svfsum /= 0._wp )  THEN
    6156 !--         TODO detect and log when normalization differs too much from 1
    6157             svf(1, isvf_surflt:nsvfl) = svf(1, isvf_surflt:nsvfl) / svfsum
     6066            i = searchsorted(svfnorm_report_thresh, svfsum / (1._wp-skyvf(isurflt_prev)))
     6067            svfnorm_counts(i) = svfnorm_counts(i) + 1
     6068
     6069            svf(1, isvf_surflt:nsvfl) = svf(1, isvf_surflt:nsvfl) / svfsum * (1._wp-skyvf(isurflt_prev))
    61586070        ENDIF
     6071        !TODO we should be able to deallocate skyvf, from now on we only need skyvft
    61596072
    61606073!--     deallocate temporary asvf array
     
    61726085
    61736086            CALL location_message( '    calculation of the complete CSF array', .TRUE. )
    6174 
     6087!             CALL radiation_write_debug_log( 'Calculation of the complete CSF array' )
    61756088!--         sort and merge csf for the last time, keeping the array size to minimum
    61766089            CALL merge_and_grow_csf(-1)
     
    62436156!--         scatter and gather the number of elements to and from all processor
    62446157!--         and calculate displacements
     6158!             CALL radiation_write_debug_log( 'Scatter and gather the number of elements to and from all processor' )
    62456159            CALL MPI_AlltoAll(icsflt,1,MPI_INTEGER,ipcsflt,1,MPI_INTEGER,comm2d, ierr)
    62466160           
     
    62516165                d = d + ipcsflt(i)
    62526166            ENDDO
    6253        
     6167
    62546168!--         exchange csf fields between processors
     6169!             CALL radiation_write_debug_log( 'Exchange csf fields between processors' )
    62556170            ALLOCATE( pcsflt(ndcsf,max(npcsfl,ndcsf)) )
    62566171            ALLOCATE( kpcsflt(kdcsf,max(npcsfl,kdcsf)) )
     
    62776192
    62786193!--         sort csf ( a version of quicksort )
     6194!             CALL radiation_write_debug_log( 'Sort csf' )
    62796195            CALL quicksort_csf2(kpcsflt, pcsflt, 1, npcsfl)
    62806196
    62816197!--         aggregate canopy sink factor records with identical box & source
    62826198!--         againg across all values from all processors
     6199!             CALL radiation_write_debug_log( 'Aggregate canopy sink factor records with identical box' )
     6200
    62836201            IF ( npcsfl > 0 )  THEN
    62846202                icsf = 1 !< reading index
     
    63266244            DEALLOCATE( pcsflt )
    63276245            DEALLOCATE( kpcsflt )
    6328             IF ( ALLOCATED( gridpcbl ) )  DEALLOCATE( gridpcbl )
     6246!             CALL radiation_write_debug_log( 'End of aggregate csf' )
    63296247           
    63306248        ENDIF
    6331        
     6249
     6250        CALL MPI_BARRIER( comm2d, ierr )
     6251!         CALL radiation_write_debug_log( 'End of radiation_calc_svf (after mpi_barrier)' )
     6252
    63326253        RETURN
    63336254       
     
    63366257            'plant canopy sink factors / direct irradiance factors.'
    63376258        CALL message( 'init_urban_surface', 'PA0502', 2, 2, 0, 6, 0 )
    6338 
    6339         CALL cpu_log( log_point_s(79), 'radiation_calc_svf', 'stop' )
    6340 
    6341 
     6259       
    63426260    END SUBROUTINE radiation_calc_svf
    63436261
    6344 
     6262   
    63456263!------------------------------------------------------------------------------!
    63466264! Description:
     
    63976315!--     Maximum number of gridboxes visited equals to maximum number of boundaries crossed in each dimension plus one. That's also
    63986316!--     the maximum number of plant canopy boxes written. We grow the acsf array accordingly using exponential factor.
    6399         maxboxes = SUM(ABS(NINT(targ) - NINT(src))) + 1
     6317        maxboxes = SUM(ABS(NINT(targ, iwp) - NINT(src, iwp))) + 1
    64006318        IF ( plant_canopy  .AND.  ncsfl + maxboxes > ncsfla )  THEN
    64016319!--         use this code for growing by fixed exponential increments (equivalent to case where ncsfl always increases by 1)
     
    64496367            IF ( crlen > .001_wp )  THEN
    64506368                crmid = (lastdist + nextdist) * .5_wp
    6451                 box = NINT(src(:) + uvect(:) * crmid)
     6369                box = NINT(src(:) + uvect(:) * crmid, iwp)
    64526370
    64536371!--             calculate index of the grid with global indices (box(2),box(3))
     
    64836401        IF ( plant_canopy )  THEN
    64846402#if defined( __parallel )
    6485             IF ( usm_lad_rma )  THEN
     6403            IF ( rma_lad_raytrace )  THEN
    64866404!--             send requests for lad_s to appropriate processor
    6487                 CALL cpu_log( log_point_s(77), 'usm_init_rma', 'start' )
     6405                CALL cpu_log( log_point_s(77), 'rad_init_rma', 'start' )
    64886406                DO i = 1, ncsb
    64896407                    CALL MPI_Get(lad_s_ray(i), 1, MPI_REAL, lad_ip(i), lad_disp(i), &
     
    65016419                    CALL message( 'raytrace', 'PA0519', 1, 2, 0, 6, 0 )
    65026420                ENDIF
    6503                 CALL cpu_log( log_point_s(77), 'usm_init_rma', 'stop' )
     6421                CALL cpu_log( log_point_s(77), 'rad_init_rma', 'stop' )
    65046422               
    65056423            ENDIF
     
    65096427            DO i = 1, ncsb
    65106428#if defined( __parallel )
    6511                 IF ( usm_lad_rma )  THEN
     6429                IF ( rma_lad_raytrace )  THEN
    65126430                    lad_s_target = lad_s_ray(i)
    65136431                ELSE
    6514                     lad_s_target = usm_lad_g(lad_ip(i)*nnx*nny*nzu + lad_disp(i))
     6432                    lad_s_target = sub_lad_g(lad_ip(i)*nnx*nny*nzu + lad_disp(i))
    65156433                ENDIF
    65166434#else
    6517                 lad_s_target = usm_lad(boxes(1,i),boxes(2,i),boxes(3,i))
     6435                lad_s_target = sub_lad(boxes(1,i),boxes(2,i),boxes(3,i))
    65186436#endif
    65196437                cursink = 1._wp - exp(-ext_coef * lad_s_target * crlens(i)*realdist)
     
    65396457
    65406458    END SUBROUTINE raytrace
    6541 
    6542 
     6459   
     6460 
    65436461!------------------------------------------------------------------------------!
    65446462! Description:
    65456463! ------------
    6546 !> Determines whether two faces are oriented towards each other. Since the
    6547 !> surfaces follow the gird box surfaces, it checks first whether the two surfaces
    6548 !> are directed in the same direction, then it checks if the two surfaces are     
     6464!> A new, more efficient version of ray tracing algorithm that processes a whole
     6465!> arc instead of a single ray.
     6466!>
     6467!> In all comments, horizon means tangent of horizon angle, i.e.
     6468!> vertical_delta / horizontal_distance
     6469!------------------------------------------------------------------------------!
     6470   SUBROUTINE raytrace_2d(origin, yxdir, zdirs, iorig, aorig, vffrac, &
     6471                              create_csf, skip_1st_pcb, win_lad, horizon, &
     6472                              transparency)
     6473      IMPLICIT NONE
     6474
     6475      REAL(wp), DIMENSION(3), INTENT(IN)     ::  origin        !< z,y,x coordinates of ray origin
     6476      REAL(wp), DIMENSION(2), INTENT(IN)     ::  yxdir         !< y,x *unit* vector of ray direction (in grid units)
     6477      REAL(wp), DIMENSION(:), INTENT(IN)     ::  zdirs         !< list of z directions to raytrace (z/hdist, in grid)
     6478      INTEGER(iwp), INTENT(in)               ::  iorig         !< index of origin face for csf
     6479      REAL(wp), INTENT(in)                   ::  aorig         !< origin face area for csf
     6480      REAL(wp), DIMENSION(LBOUND(zdirs, 1):UBOUND(zdirs, 1)), INTENT(in) ::  vffrac !<
     6481                                                               !< view factor fractions of each ray for csf
     6482      LOGICAL, INTENT(in)                    ::  create_csf    !< whether to generate new CSFs during raytracing
     6483      LOGICAL, INTENT(in)                    ::  skip_1st_pcb  !< whether to skip first plant canopy box during raytracing
     6484      INTEGER(iwp), INTENT(in)               ::  win_lad       !< leaf area density MPI window
     6485      REAL(wp), INTENT(OUT)                  ::  horizon       !< highest horizon found after raytracing (z/hdist)
     6486      REAL(wp), DIMENSION(LBOUND(zdirs, 1):UBOUND(zdirs, 1)), INTENT(OUT) ::  transparency !<
     6487                                                                !< transparencies of zdirs paths
     6488      !--INTEGER(iwp), DIMENSION(3, LBOUND(zdirs, 1):UBOUND(zdirs, 1)), INTENT(OUT) :: itarget !<
     6489                                                                !< (z,y,x) coordinates of target faces for zdirs
     6490      INTEGER(iwp)                           ::  i, k, l, d
     6491      INTEGER(iwp)                           ::  seldim       !< dimension to be incremented
     6492      REAL(wp), DIMENSION(2)                 ::  yxorigin     !< horizontal copy of origin (y,x)
     6493      REAL(wp)                               ::  distance     !< euclidean along path
     6494      REAL(wp)                               ::  lastdist     !< beginning of current crossing
     6495      REAL(wp)                               ::  nextdist     !< end of current crossing
     6496      REAL(wp)                               ::  crmid        !< midpoint of crossing
     6497      REAL(wp)                               ::  horz_entry   !< horizon at entry to column
     6498      REAL(wp)                               ::  horz_exit    !< horizon at exit from column
     6499      REAL(wp)                               ::  bdydim       !< boundary for current dimension
     6500      REAL(wp), DIMENSION(2)                 ::  crossdist    !< distances to boundary for dimensions
     6501      REAL(wp), DIMENSION(2)                 ::  dimnextdist  !< distance for each dimension increments
     6502      INTEGER(iwp), DIMENSION(2)             ::  column       !< grid column being crossed
     6503      INTEGER(iwp), DIMENSION(2)             ::  dimnext      !< next dimension increments along path
     6504      INTEGER(iwp), DIMENSION(2)             ::  dimdelta     !< dimension direction = +- 1
     6505      INTEGER(iwp)                           ::  px, py       !< number of processors in x and y dir before
     6506                                                              !< the processor in the question
     6507      INTEGER(iwp)                           ::  ip           !< number of processor where gridbox reside
     6508      INTEGER(iwp)                           ::  ig           !< 1D index of gridbox in global 2D array
     6509      INTEGER(MPI_ADDRESS_KIND)              ::  wdisp        !< RMA window displacement
     6510      INTEGER(iwp)                           ::  wcount       !< RMA window item count
     6511      INTEGER(iwp)                           ::  maxboxes     !< max no of CSF created
     6512      INTEGER(iwp)                           ::  nly          !< maximum  plant canopy height
     6513      INTEGER(iwp)                           ::  ntrack
     6514      REAL(wp)                               ::  zbottom, ztop !< urban surface boundary in real numbers
     6515      REAL(wp)                               ::  zorig        !< z coordinate of ray column entry
     6516      REAL(wp)                               ::  zexit        !< z coordinate of ray column exit
     6517      REAL(wp)                               ::  qdist        !< ratio of real distance to z coord difference
     6518      REAL(wp)                               ::  dxxyy        !< square of real horizontal distance
     6519      REAL(wp)                               ::  curtrans     !< transparency of current PC box crossing
     6520      INTEGER(iwp)                           ::  zb0
     6521      INTEGER(iwp)                           ::  zb1
     6522      INTEGER(iwp)                           ::  nz
     6523      INTEGER(iwp)                           ::  iz
     6524      INTEGER(iwp)                           ::  zsgn
     6525      REAL(wp), PARAMETER                    ::  grow_factor = 1.5_wp !< factor of expansion of grow arrays
     6526
     6527     
     6528      yxorigin(:) = origin(2:3)
     6529      transparency(:) = 1._wp !-- Pre-set the all rays to transparent before reducing
     6530      horizon = -HUGE(1._wp)
     6531
     6532      !--Determine distance to boundary (in 2D xy)
     6533      IF ( yxdir(1) > 0._wp )  THEN
     6534         bdydim = ny + .5_wp !< north global boundary
     6535         crossdist(1) = (bdydim - yxorigin(1)) / yxdir(1)
     6536      ELSEIF ( yxdir(1) == 0._wp )  THEN
     6537         crossdist(1) = HUGE(1._wp)
     6538      ELSE
     6539          bdydim = -.5_wp !< south global boundary
     6540          crossdist(1) = (bdydim - yxorigin(1)) / yxdir(1)
     6541      ENDIF
     6542
     6543      IF ( yxdir(2) >= 0._wp )  THEN
     6544          bdydim = nx + .5_wp !< east global boundary
     6545          crossdist(2) = (bdydim - yxorigin(2)) / yxdir(2)
     6546      ELSEIF ( yxdir(2) == 0._wp )  THEN
     6547         crossdist(2) = HUGE(1._wp)
     6548      ELSE
     6549          bdydim = -.5_wp !< west global boundary
     6550          crossdist(2) = (bdydim - yxorigin(2)) / yxdir(2)
     6551      ENDIF
     6552      distance = minval(crossdist, 1)
     6553
     6554      IF ( plant_canopy )  THEN
     6555         rt2_track_dist(0) = 0._wp
     6556         rt2_track_lad(:,:) = 0._wp
     6557         nly = plantt_max - nzub + 1
     6558      ENDIF
     6559
     6560      lastdist = 0._wp
     6561
     6562!--   Since all face coordinates have values *.5 and we'd like to use
     6563!--   integers, all these have .5 added
     6564      DO d = 1, 2
     6565          IF ( yxdir(d) == 0._wp )  THEN
     6566              dimnext(d) = HUGE(1_iwp)
     6567              dimdelta(d) = HUGE(1_iwp)
     6568              dimnextdist(d) = HUGE(1._wp)
     6569          ELSE IF ( yxdir(d) > 0._wp )  THEN
     6570              dimnext(d) = FLOOR(yxorigin(d) + .5_wp) + 1
     6571              dimdelta(d) = 1
     6572              dimnextdist(d) = (dimnext(d) - .5_wp - yxorigin(d)) / yxdir(d)
     6573          ELSE
     6574              dimnext(d) = CEILING(yxorigin(d) + .5_wp) - 1
     6575              dimdelta(d) = -1
     6576              dimnextdist(d) = (dimnext(d) - .5_wp - yxorigin(d)) / yxdir(d)
     6577          ENDIF
     6578      ENDDO
     6579
     6580      ntrack = 0
     6581      DO
     6582!--      along what dimension will the next wall crossing be?
     6583         seldim = minloc(dimnextdist, 1)
     6584         nextdist = dimnextdist(seldim)
     6585         IF ( nextdist > distance ) nextdist = distance
     6586
     6587         IF ( nextdist > lastdist )  THEN
     6588            ntrack = ntrack + 1
     6589            crmid = (lastdist + nextdist) * .5_wp
     6590            column = NINT(yxorigin(:) + yxdir(:) * crmid, iwp)
     6591
     6592!--         calculate index of the grid with global indices (column(1),column(2))
     6593!--         in the array nzterr and plantt and id of the coresponding processor
     6594            px = column(2)/nnx
     6595            py = column(1)/nny
     6596            ip = px*pdims(2)+py
     6597            ig = ip*nnx*nny + (column(2)-px*nnx)*nny + column(1)-py*nny
     6598
     6599            IF ( lastdist == 0._wp )  THEN
     6600               horz_entry = -HUGE(1._wp)
     6601            ELSE
     6602               horz_entry = (nzterr(ig) - origin(1)) / lastdist
     6603            ENDIF
     6604            horz_exit = (nzterr(ig) - origin(1)) / nextdist
     6605            horizon = MAX(horizon, horz_entry, horz_exit)
     6606
     6607            IF ( plant_canopy )  THEN
     6608               rt2_track(:, ntrack) = column(:)
     6609               rt2_track_dist(ntrack) = nextdist
     6610            ENDIF
     6611         ENDIF
     6612
     6613         IF ( nextdist >= distance )  EXIT
     6614         lastdist = nextdist
     6615         dimnext(seldim) = dimnext(seldim) + dimdelta(seldim)
     6616         dimnextdist(seldim) = (dimnext(seldim) - .5_wp - yxorigin(seldim)) / yxdir(seldim)
     6617      ENDDO
     6618
     6619      IF ( plant_canopy )  THEN
     6620         !--Request LAD WHERE applicable
     6621         !--
     6622#if defined( __parallel )
     6623         IF ( rma_lad_raytrace )  THEN
     6624!--         send requests for lad_s to appropriate processor
     6625            !CALL cpu_log( log_point_s(77), 'usm_init_rma', 'start' )
     6626            DO i = 1, ntrack
     6627               px = rt2_track(2,i)/nnx
     6628               py = rt2_track(1,i)/nny
     6629               ip = px*pdims(2)+py
     6630               ig = ip*nnx*nny + (rt2_track(2,i)-px*nnx)*nny + rt2_track(1,i)-py*nny
     6631               IF ( plantt(ig) <= nzterr(ig) )  CYCLE
     6632               wdisp = (rt2_track(2,i)-px*nnx)*(nny*nzu) + (rt2_track(1,i)-py*nny)*nzu + nzterr(ig)+1-nzub
     6633               wcount = plantt(ig)-nzterr(ig)
     6634               ! TODO send request ASAP - even during raytracing
     6635               CALL MPI_Get(rt2_track_lad(nzterr(ig)+1:plantt(ig), i), wcount, MPI_REAL, ip,    &
     6636                            wdisp, wcount, MPI_REAL, win_lad, ierr)
     6637               IF ( ierr /= 0 )  THEN
     6638                  WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Get'
     6639                  CALL message( 'raytrace_2d', 'PA0526', 1, 2, 0, 6, 0 )
     6640               ENDIF
     6641            ENDDO
     6642
     6643!--         wait for all pending local requests complete
     6644            ! TODO WAIT selectively for each column later when needed
     6645            CALL MPI_Win_flush_local_all(win_lad, ierr)
     6646            IF ( ierr /= 0 )  THEN
     6647               WRITE(message_string, *) 'MPI error ', ierr, ' at MPI_Win_flush_local_all'
     6648               CALL message( 'raytrace', 'PA0527', 1, 2, 0, 6, 0 )
     6649            ENDIF
     6650            !CALL cpu_log( log_point_s(77), 'usm_init_rma', 'stop' )
     6651         ELSE ! rma_lad_raytrace
     6652            DO i = 1, ntrack
     6653               px = rt2_track(2,i)/nnx
     6654               py = rt2_track(1,i)/nny
     6655               ip = px*pdims(2)+py
     6656               ig = ip*nnx*nny*nzu + (rt2_track(2,i)-px*nnx)*(nny*nzu) + (rt2_track(1,i)-py*nny)*nzu
     6657               rt2_track_lad(nzub:plantt_max, i) = sub_lad_g(ig:ig+nly-1)
     6658            ENDDO
     6659         ENDIF
     6660#else
     6661         DO i = 1, ntrack
     6662            rt2_track_lad(nzub:plantt_max, i) = sub_lad(rt2_track(1,i), rt2_track(2,i), nzub:plantt_max)
     6663         ENDDO
     6664#endif
     6665
     6666         !--Skip the PCB around origin if requested
     6667         !--
     6668         IF ( skip_1st_pcb )  THEN
     6669            rt2_track_lad(NINT(origin(1), iwp), 1) = 0._wp
     6670         ENDIF
     6671
     6672         !--Assert that we have space allocated for CSFs
     6673         !--
     6674         maxboxes = (ntrack + MAX(origin(1) - nzub, nzut - origin(1))) * SIZE(zdirs, 1)
     6675         IF ( ncsfl + maxboxes > ncsfla )  THEN
     6676!--         use this code for growing by fixed exponential increments (equivalent to case where ncsfl always increases by 1)
     6677!--         k = CEILING(grow_factor ** real(CEILING(log(real(ncsfl + maxboxes, kind=wp)) &
     6678!--                                                / log(grow_factor)), kind=wp))
     6679!--         or use this code to simply always keep some extra space after growing
     6680            k = CEILING(REAL(ncsfl + maxboxes, kind=wp) * grow_factor)
     6681            CALL merge_and_grow_csf(k)
     6682         ENDIF
     6683
     6684         !--Calculate transparencies and store new CSFs
     6685         !--
     6686         zbottom = REAL(nzub, wp) - .5_wp
     6687         ztop = REAL(plantt_max, wp) + .5_wp
     6688
     6689         !--Reverse direction of radiation (face->sky), only when create_csf
     6690         !--
     6691         IF ( create_csf )  THEN
     6692            DO i = 1, ntrack ! for each column
     6693               dxxyy = ((dy*yxdir(1))**2 + (dx*yxdir(2))**2) * (rt2_track_dist(i)-rt2_track_dist(i-1))**2
     6694               px = rt2_track(2,i)/nnx
     6695               py = rt2_track(1,i)/nny
     6696               ip = px*pdims(2)+py
     6697
     6698               DO k = LBOUND(zdirs, 1), UBOUND(zdirs, 1) ! for each ray
     6699                  IF ( zdirs(k) <= horizon )  THEN
     6700                     CYCLE
     6701                  ENDIF
     6702
     6703                  zorig = REAL(origin(1), wp) + zdirs(k) * rt2_track_dist(i-1)
     6704                  IF ( zorig <= zbottom .OR. zorig >= ztop )  CYCLE
     6705
     6706                  zsgn = INT(SIGN(1._wp, zdirs(k)), iwp)
     6707                  rt2_dist(1) = 0._wp
     6708                  IF ( zdirs(k) == 0._wp )  THEN ! ray is exactly horizontal
     6709                     nz = 2
     6710                     rt2_dist(nz) = SQRT(dxxyy)
     6711                     iz = NINT(zorig, iwp)
     6712                  ELSE
     6713                     zexit = MIN(MAX(REAL(origin(1), wp) + zdirs(k) * rt2_track_dist(i), zbottom), ztop)
     6714
     6715                     zb0 = FLOOR(  zorig * zsgn - .5_wp) + 1  ! because it must be greater than orig
     6716                     zb1 = CEILING(zexit * zsgn - .5_wp) - 1  ! because it must be smaller than exit
     6717                     nz = MAX(zb1 - zb0 + 3, 2)
     6718                     rt2_dist(nz) = SQRT(((zexit-zorig)*dz)**2 + dxxyy)
     6719                     qdist = rt2_dist(nz) / (zexit-zorig)
     6720                     rt2_dist(2:nz-1) = (/( ((REAL(l, wp) + .5_wp) * zsgn - zorig) * qdist , l = zb0, zb1 )/)
     6721                     iz = zb0 * zsgn
     6722                  ENDIF
     6723
     6724                  DO l = 2, nz
     6725                     IF ( rt2_track_lad(iz, i) > 0._wp )  THEN
     6726                        curtrans = exp(-ext_coef * rt2_track_lad(iz, i) * (rt2_dist(l)-rt2_dist(l-1)))
     6727
     6728                        ncsfl = ncsfl + 1
     6729                        acsf(ncsfl)%ip = ip
     6730                        acsf(ncsfl)%itx = rt2_track(2,i)
     6731                        acsf(ncsfl)%ity = rt2_track(1,i)
     6732                        acsf(ncsfl)%itz = iz
     6733                        acsf(ncsfl)%isurfs = iorig
     6734                        acsf(ncsfl)%rsvf = REAL((1._wp - curtrans)*aorig*vffrac(k), wp) ! we postpone multiplication by transparency
     6735                        acsf(ncsfl)%rtransp = REAL(transparency(k), wp)
     6736
     6737                        transparency(k) = transparency(k) * curtrans
     6738                     ENDIF
     6739                     iz = iz + zsgn
     6740                  ENDDO ! l = 1, nz - 1
     6741               ENDDO ! k = LBOUND(zdirs, 1), UBOUND(zdirs, 1)
     6742            ENDDO ! i = 1, ntrack
     6743
     6744            transparency(:) = 1._wp !-- Reset all rays to transparent
     6745         ENDIF
     6746
     6747         !-- Forward direction of radiation (sky->face), always
     6748         !--
     6749         DO i = ntrack, 1, -1 ! for each column backwards
     6750            dxxyy = ((dy*yxdir(1))**2 + (dx*yxdir(2))**2) * (rt2_track_dist(i)-rt2_track_dist(i-1))**2
     6751            px = rt2_track(2,i)/nnx
     6752            py = rt2_track(1,i)/nny
     6753            ip = px*pdims(2)+py
     6754
     6755            DO k = LBOUND(zdirs, 1), UBOUND(zdirs, 1) ! for each ray
     6756               IF ( zdirs(k) <= horizon )  THEN
     6757                  transparency(k) = 0._wp
     6758                  CYCLE
     6759               ENDIF
     6760
     6761               zexit = REAL(origin(1), wp) + zdirs(k) * rt2_track_dist(i-1)
     6762               IF ( zexit <= zbottom .OR. zexit >= ztop )  CYCLE
     6763
     6764               zsgn = -INT(SIGN(1._wp, zdirs(k)), iwp)
     6765               rt2_dist(1) = 0._wp
     6766               IF ( zdirs(k) == 0._wp )  THEN ! ray is exactly horizontal
     6767                  nz = 2
     6768                  rt2_dist(nz) = SQRT(dxxyy)
     6769                  iz = NINT(zexit, iwp)
     6770               ELSE
     6771                  zorig = MIN(MAX(REAL(origin(1), wp) + zdirs(k) * rt2_track_dist(i), zbottom), ztop)
     6772
     6773                  zb0 = FLOOR(  zorig * zsgn - .5_wp) + 1  ! because it must be greater than orig
     6774                  zb1 = CEILING(zexit * zsgn - .5_wp) - 1  ! because it must be smaller than exit
     6775                  nz = MAX(zb1 - zb0 + 3, 2)
     6776                  rt2_dist(nz) = SQRT(((zexit-zorig)*dz)**2 + dxxyy)
     6777                  qdist = rt2_dist(nz) / (zexit-zorig)
     6778                  rt2_dist(2:nz-1) = (/( ((REAL(l, wp) + .5_wp) * zsgn - zorig) * qdist , l = zb0, zb1 )/)
     6779                  iz = zb0 * zsgn
     6780               ENDIF
     6781
     6782               DO l = 2, nz
     6783                  IF ( rt2_track_lad(iz, i) > 0._wp )  THEN
     6784                     curtrans = exp(-ext_coef * rt2_track_lad(iz, i) * (rt2_dist(l)-rt2_dist(l-1)))
     6785
     6786                     IF ( create_csf )  THEN
     6787                        ncsfl = ncsfl + 1
     6788                        acsf(ncsfl)%ip = ip
     6789                        acsf(ncsfl)%itx = rt2_track(2,i)
     6790                        acsf(ncsfl)%ity = rt2_track(1,i)
     6791                        acsf(ncsfl)%itz = iz
     6792                        acsf(ncsfl)%isurfs = -1 ! a special ID indicating sky
     6793                        acsf(ncsfl)%rsvf = REAL((1._wp - curtrans)*aorig*vffrac(k), wp) ! we postpone multiplication by transparency
     6794                        acsf(ncsfl)%rtransp = REAL(transparency(k), wp)
     6795                     ENDIF  !< create_csf
     6796
     6797                     transparency(k) = transparency(k) * curtrans
     6798                  ENDIF
     6799                  iz = iz + zsgn
     6800               ENDDO ! l = 1, nz - 1
     6801            ENDDO ! k = LBOUND(zdirs, 1), UBOUND(zdirs, 1)
     6802         ENDDO ! i = 1, ntrack
     6803
     6804      ELSE ! not plant_canopy
     6805         DO k = UBOUND(zdirs, 1), LBOUND(zdirs, 1), -1 ! TODO make more generic
     6806            IF ( zdirs(k) > horizon )  EXIT
     6807            transparency(k) = 0._wp
     6808         ENDDO
     6809      ENDIF
     6810
     6811   END SUBROUTINE raytrace_2d
     6812 
     6813
     6814!------------------------------------------------------------------------------!
     6815!
     6816! Description:
     6817! ------------
     6818!> Calculates apparent solar positions for all timesteps and stores discretized
     6819!> positions.
     6820!------------------------------------------------------------------------------!
     6821   SUBROUTINE radiation_presimulate_solar_pos
     6822      IMPLICIT NONE
     6823
     6824      INTEGER(iwp)                              ::  it, i, j
     6825      REAL(wp)                                  ::  tsrp_prev
     6826      REAL(wp), DIMENSION(:,:), ALLOCATABLE     ::  dsidir_tmp       !< dsidir_tmp[:,i] = unit vector of i-th
     6827                                                                     !< appreant solar direction
     6828
     6829      ALLOCATE ( dsidir_rev(0:raytrace_discrete_elevs/2-1,                 &
     6830                            0:raytrace_discrete_azims-1) )
     6831      dsidir_rev(:,:) = -1
     6832      ALLOCATE ( dsidir_tmp(3,                                             &
     6833                     raytrace_discrete_elevs/2*raytrace_discrete_azims) )
     6834      ndsidir = 0
     6835
     6836!
     6837!--   We will artificialy update time_since_reference_point and return to
     6838!--   true value later
     6839      tsrp_prev = time_since_reference_point
     6840      sun_direction = .TRUE.
     6841
     6842!
     6843!--   Process spinup time if configured
     6844      IF ( spinup_time > 0._wp )  THEN
     6845         DO  it = 0, CEILING(spinup_time / dt_spinup)
     6846            time_since_reference_point = -spinup_time + REAL(it, wp) * dt_spinup
     6847            CALL simulate_pos
     6848         ENDDO
     6849      ENDIF
     6850!
     6851!--   Process simulation time
     6852      DO  it = 0, CEILING(end_time / dt_radiation)
     6853         time_since_reference_point = REAL(it, wp) * dt_radiation
     6854         CALL simulate_pos
     6855      ENDDO
     6856
     6857      time_since_reference_point = tsrp_prev
     6858
     6859!--   Allocate global vars which depend on ndsidir
     6860      ALLOCATE ( dsidir ( 3, ndsidir ) )
     6861      dsidir(:,:) = dsidir_tmp(:, 1:ndsidir)
     6862      DEALLOCATE ( dsidir_tmp )
     6863      ALLOCATE ( dsitrans(nsurfl, ndsidir) )
     6864      ALLOCATE ( dsitransc(npcbl, ndsidir) )
     6865
     6866      WRITE ( message_string, * ) 'Precalculated', ndsidir, ' solar positions', &
     6867                                  'from', it, ' timesteps.'
     6868      CALL message( 'radiation_presimulate_solar_pos', 'UI0013', 0, 0, 0, 6, 0 )
     6869
     6870      CONTAINS
     6871
     6872      !------------------------------------------------------------------------!
     6873      ! Description:
     6874      ! ------------
     6875      !> Simuates a single position
     6876      !------------------------------------------------------------------------!
     6877      SUBROUTINE simulate_pos
     6878         IMPLICIT NONE
     6879!
     6880!--      Update apparent solar position based on modified t_s_r_p
     6881         CALL calc_zenith
     6882         IF ( zenith(0) > 0 )  THEN
     6883!--         
     6884!--         Identify solar direction vector (discretized number) 1)
     6885            i = MODULO(NINT(ATAN2(sun_dir_lon(0), sun_dir_lat(0))               &
     6886                            / (2._wp*pi) * raytrace_discrete_azims-.5_wp, iwp), &
     6887                       raytrace_discrete_azims)
     6888            j = FLOOR(ACOS(zenith(0)) / pi * raytrace_discrete_elevs)
     6889            IF ( dsidir_rev(j, i) == -1 )  THEN
     6890               ndsidir = ndsidir + 1
     6891               dsidir_tmp(:, ndsidir) =                                              &
     6892                     (/ COS((REAL(j,wp)+.5_wp) * pi      / raytrace_discrete_elevs), &
     6893                        SIN((REAL(j,wp)+.5_wp) * pi      / raytrace_discrete_elevs)  &
     6894                      * COS((REAL(i,wp)+.5_wp) * 2_wp*pi / raytrace_discrete_azims), &
     6895                        SIN((REAL(j,wp)+.5_wp) * pi      / raytrace_discrete_elevs)  &
     6896                      * SIN((REAL(i,wp)+.5_wp) * 2_wp*pi / raytrace_discrete_azims) /)
     6897               dsidir_rev(j, i) = ndsidir
     6898            ENDIF
     6899         ENDIF
     6900      END SUBROUTINE simulate_pos
     6901
     6902   END SUBROUTINE radiation_presimulate_solar_pos
     6903
     6904
     6905
     6906!------------------------------------------------------------------------------!
     6907! Description:
     6908! ------------
     6909!> Determines whether two faces are oriented towards each other. Since the
     6910!> surfaces follow the gird box surfaces, it checks first whether the two surfaces
     6911!> are directed in the same direction, then it checks if the two surfaces are
    65496912!> located in confronted direction but facing away from each other, e.g. <--| |-->
    65506913!------------------------------------------------------------------------------!
     
    65586921        IF ( (d==iup_u  .OR.  d==iup_l  .OR.  d==iup_a )                             &
    65596922             .AND. (d2==iup_u  .OR. d2==iup_l) ) RETURN
    6560         IF ( (d==isky  .OR.  d==idown_a)  .AND.  d2==isky ) RETURN
    6561         IF ( (d==isouth_u  .OR.  d==isouth_l  .OR.  d==isouth_a  .OR.  d==inorth_b ) &
    6562              .AND.  (d2==isouth_u  .OR.  d2==isouth_l  .OR.  d2==inorth_b) ) RETURN
    6563         IF ( (d==inorth_u  .OR.  d==inorth_l  .OR.  d==inorth_a  .OR.  d==isouth_b ) &
    6564              .AND.  (d2==inorth_u  .OR.  d2==inorth_l  .OR.  d2==isouth_b) ) RETURN
    6565         IF ( (d==iwest_u  .OR.  d==iwest_l  .OR.  d==iwest_a  .OR.  d==ieast_b )     &
    6566              .AND.  (d2==iwest_u  .OR.  d2==iwest_l  .OR.  d2==ieast_b ) ) RETURN
    6567         IF ( (d==ieast_u  .OR.  d==ieast_l  .OR.  d==ieast_a  .OR.  d==iwest_b )     &
    6568              .AND.  (d2==ieast_u  .OR.  d2==ieast_l  .OR.  d2==iwest_b ) ) RETURN
     6923        IF ( (d==isouth_u  .OR.  d==isouth_l  .OR.  d==isouth_a ) &
     6924             .AND.  (d2==isouth_u  .OR.  d2==isouth_l) ) RETURN
     6925        IF ( (d==inorth_u  .OR.  d==inorth_l  .OR.  d==inorth_a ) &
     6926             .AND.  (d2==inorth_u  .OR.  d2==inorth_l) ) RETURN
     6927        IF ( (d==iwest_u  .OR.  d==iwest_l  .OR.  d==iwest_a )     &
     6928             .AND.  (d2==iwest_u  .OR.  d2==iwest_l ) ) RETURN
     6929        IF ( (d==ieast_u  .OR.  d==ieast_l  .OR.  d==ieast_a )     &
     6930             .AND.  (d2==ieast_u  .OR.  d2==ieast_l ) ) RETURN
    65696931
    65706932!-- second check: are surfaces facing away from each other
    65716933        SELECT CASE (d)
    6572             CASE (iup_u, iup_l, iup_a)                    !< upward facing surfaces
     6934            CASE (iup_u, iup_l, iup_a)              !< upward facing surfaces
    65736935                IF ( z2 < z ) RETURN
    6574             CASE (isky, idown_a)                          !< downward facing surfaces
     6936            CASE (idown_a)                          !< downward facing surfaces
    65756937                IF ( z2 > z ) RETURN
    6576             CASE (isouth_u, isouth_l, isouth_a, inorth_b) !< southward facing surfaces
     6938            CASE (isouth_u, isouth_l, isouth_a)    !< southward facing surfaces
    65776939                IF ( y2 > y ) RETURN
    6578             CASE (inorth_u, inorth_l, inorth_a, isouth_b) !< northward facing surfaces
     6940            CASE (inorth_u, inorth_l, inorth_a)    !< northward facing surfaces
    65796941                IF ( y2 < y ) RETURN
    6580             CASE (iwest_u, iwest_l, iwest_a, ieast_b)     !< westward facing surfaces
     6942            CASE (iwest_u, iwest_l, iwest_a)        !< westward facing surfaces
    65816943                IF ( x2 > x ) RETURN
    6582             CASE (ieast_u, ieast_l, ieast_a, iwest_b)     !< eastward facing surfaces
     6944            CASE (ieast_u, ieast_l, ieast_a)        !< eastward facing surfaces
    65836945                IF ( x2 < x ) RETURN
    65846946        END SELECT
    65856947
    65866948        SELECT CASE (d2)
    6587             CASE (iup_u)                        !< ground, roof
     6949            CASE (iup_u)                            !< ground, roof
    65886950                IF ( z < z2 ) RETURN
    6589             CASE (isky)                         !< sky
    6590                 IF ( z > z2 ) RETURN
    6591             CASE (isouth_u, isouth_l, inorth_b) !< south facing
     6951            CASE (isouth_u, isouth_l)               !< south facing
    65926952                IF ( y > y2 ) RETURN
    6593             CASE (inorth_u, inorth_l, isouth_b) !< north facing
     6953            CASE (inorth_u, inorth_l)              !< north facing
    65946954                IF ( y < y2 ) RETURN
    6595             CASE (iwest_u, iwest_l, ieast_b)    !< west facing
     6955            CASE (iwest_u, iwest_l)                 !< west facing
    65966956                IF ( x > x2 ) RETURN
    6597             CASE (ieast_u, ieast_l, iwest_b)    !< east facing
     6957            CASE (ieast_u, ieast_l)                 !< east facing
    65986958                IF ( x < x2 ) RETURN
    65996959            CASE (-1)
     
    66046964       
    66056965    END FUNCTION surface_facing
     6966
    66066967
    66076968!------------------------------------------------------------------------------!
     
    66136974    SUBROUTINE radiation_read_svf
    66146975
    6615        IMPLICIT NONE
    6616        INTEGER(iwp)                 :: fsvf = 88
    6617        INTEGER(iwp)                 :: i
    6618        CHARACTER(usm_version_len)   :: usm_version_field
    6619        CHARACTER(svf_code_len)      :: svf_code_field
    6620 
    6621        DO  i = 0, io_blocks-1
    6622           IF ( i == io_group )  THEN
     6976        IMPLICIT NONE
     6977        INTEGER(iwp)                 :: fsvf = 88
     6978        INTEGER(iwp)                 :: i
     6979        CHARACTER(rad_version_len)   :: rad_version_field
     6980        CHARACTER(svf_code_len)      :: svf_code_field
     6981
     6982        DO  i = 0, io_blocks-1
     6983            IF ( i == io_group )  THEN
    66236984
    66246985!
     
    66266987             CALL check_open( fsvf )
    66276988
    6628 
    6629 !--          read and check version
    6630              READ ( fsvf ) usm_version_field
    6631              IF ( TRIM(usm_version_field) /= TRIM(usm_version) )  THEN
    6632                  WRITE( message_string, * ) 'Version of binary SVF file "',           &
    6633                                          TRIM(usm_version_field), '" does not match ',            &
    6634                                          'the version of model "', TRIM(usm_version), '"'
    6635                  CALL message( 'radiation_read_svf', 'PA0482', 1, 2, 0, 6, 0 )
    6636              ENDIF
     6989!--             read and check version
     6990                READ ( fsvf ) rad_version_field
     6991                IF ( TRIM(rad_version_field) /= TRIM(rad_version) )  THEN
     6992                    WRITE( message_string, * ) 'Version of binary SVF file "',           &
     6993                                TRIM(rad_version_field), '" does not match ',            &
     6994                                'the version of model "', TRIM(rad_version), '"'
     6995                    CALL message( 'radiation_read_svf', 'PA0482', 1, 2, 0, 6, 0 )
     6996                ENDIF
    66376997               
    6638 !--          read nsvfl, ncsfl
    6639              READ ( fsvf ) nsvfl, ncsfl
    6640              IF ( nsvfl <= 0  .OR.  ncsfl < 0 )  THEN
    6641                 WRITE( message_string, * ) 'Wrong number of SVF or CSF'
    6642                 CALL message( 'radiation_read_svf', 'PA0483', 1, 2, 0, 6, 0 )
    6643              ELSE
    6644                 WRITE(message_string,*) '    Number of SVF and CSF to read', nsvfl, ncsfl
    6645                 CALL location_message( message_string, .TRUE. )
    6646              ENDIF
     6998!--             read nsvfl, ncsfl
     6999                READ ( fsvf ) nsvfl, ncsfl, nsurfl
     7000                IF ( nsvfl <= 0  .OR.  ncsfl < 0 )  THEN
     7001                    WRITE( message_string, * ) 'Wrong number of SVF or CSF'
     7002                    CALL message( 'radiation_read_svf', 'PA0483', 1, 2, 0, 6, 0 )
     7003                ELSE
     7004                    WRITE(message_string,*) '    Number of SVF, CSF, and nsurfl to read '&
     7005                         , nsvfl, ncsfl, nsurfl
     7006                    CALL location_message( message_string, .TRUE. )
     7007                ENDIF
    66477008               
    6648              ALLOCATE(svf(ndsvf,nsvfl))
    6649              ALLOCATE(svfsurf(idsvf,nsvfl))
    6650 
    6651              READ(fsvf) svf
    6652              READ(fsvf) svfsurf
    6653 
    6654              IF ( plant_canopy )  THEN
    6655                  ALLOCATE(csf(ndcsf,ncsfl))
    6656                  ALLOCATE(csfsurf(idcsf,ncsfl))
    6657                  READ(fsvf) csf
    6658                  READ(fsvf) csfsurf
    6659              ENDIF
    6660 
    6661              READ ( fsvf ) svf_code_field
    6662              IF ( TRIM(svf_code_field) /= TRIM(svf_code) )  THEN
    6663                 WRITE( message_string, * ) 'Wrong structure of binary svf file'
    6664                 CALL message( 'radiation_read_svf', 'PA0484', 1, 2, 0, 6, 0 )
    6665              ENDIF
     7009                ALLOCATE(skyvf(nsurfl))
     7010                ALLOCATE(skyvft(nsurfl))
     7011                ALLOCATE(svf(ndsvf,nsvfl))
     7012                ALLOCATE(svfsurf(idsvf,nsvfl))
     7013                READ(fsvf) skyvf
     7014                READ(fsvf) skyvft
     7015                READ(fsvf) svf
     7016                READ(fsvf) svfsurf
     7017                IF ( plant_canopy )  THEN
     7018                    ALLOCATE(csf(ndcsf,ncsfl))
     7019                    ALLOCATE(csfsurf(idcsf,ncsfl))
     7020                    READ(fsvf) csf
     7021                    READ(fsvf) csfsurf
     7022                ENDIF
     7023                READ ( fsvf ) svf_code_field
    66667024               
     7025                IF ( TRIM(svf_code_field) /= TRIM(svf_code) )  THEN
     7026                    WRITE( message_string, * ) 'Wrong structure of binary svf file'
     7027                    CALL message( 'radiation_read_svf', 'PA0484', 1, 2, 0, 6, 0 )
     7028                ENDIF       
    66677029!
    66687030!--          Close binary file                 
    66697031             CALL close_file( fsvf )
    66707032               
    6671           ENDIF
     7033            ENDIF
    66727034#if defined( __parallel )
    6673           CALL MPI_BARRIER( comm2d, ierr )
     7035            CALL MPI_BARRIER( comm2d, ierr )
    66747036#endif
    6675        ENDDO
     7037        ENDDO
    66767038
    66777039    END SUBROUTINE radiation_read_svf
     
    66867048    SUBROUTINE radiation_write_svf
    66877049
    6688 
    6689        IMPLICIT NONE
    6690 
    6691        INTEGER(iwp)        :: fsvf = 89
    6692        INTEGER(iwp)        :: i
    6693 
    6694 
    6695        DO  i = 0, io_blocks-1
    6696           IF ( i == io_group )  THEN
    6697            
     7050        IMPLICIT NONE
     7051        INTEGER(iwp)        :: fsvf = 89
     7052        INTEGER(iwp)        :: i
     7053
     7054        DO  i = 0, io_blocks-1
     7055            IF ( i == io_group )  THEN
    66987056!
    66997057!--          Open binary file
    67007058             CALL check_open( fsvf )
    67017059
    6702              WRITE ( fsvf )  usm_version
    6703              WRITE ( fsvf )  nsvfl, ncsfl
    6704              WRITE ( fsvf )  svf
    6705              WRITE ( fsvf )  svfsurf
    6706              IF ( plant_canopy )  THEN
    6707                 WRITE ( fsvf )  csf
    6708                 WRITE ( fsvf )  csfsurf
    6709              ENDIF
    6710              WRITE ( fsvf )  TRIM(svf_code)
    6711 
     7060                WRITE ( fsvf )  rad_version
     7061                WRITE ( fsvf )  nsvfl, ncsfl, nsurfl
     7062                WRITE ( fsvf )  skyvf
     7063                WRITE ( fsvf )  skyvft
     7064                WRITE ( fsvf )  svf
     7065                WRITE ( fsvf )  svfsurf
     7066                IF ( plant_canopy )  THEN
     7067                    WRITE ( fsvf )  csf
     7068                    WRITE ( fsvf )  csfsurf
     7069                ENDIF
     7070                WRITE ( fsvf )  TRIM(svf_code)
    67127071!
    67137072!--          Close binary file                 
     
    67167075          ENDIF
    67177076#if defined( __parallel )
    6718           CALL MPI_BARRIER( comm2d, ierr )
     7077                CALL MPI_BARRIER( comm2d, ierr )
    67197078#endif
    6720        ENDDO
    6721 
     7079        ENDDO
    67227080    END SUBROUTINE radiation_write_svf
    67237081
    6724 
     7082!------------------------------------------------------------------------------!
     7083!
     7084! Description:
     7085! ------------
     7086!> Block of auxiliary subroutines:
     7087!> 1. quicksort and corresponding comparison
     7088!> 2. merge_and_grow_csf for implementation of "dynamical growing"
     7089!>    array for csf
     7090!------------------------------------------------------------------------------!
     7091    PURE FUNCTION svf_lt(svf1,svf2) result (res)
     7092      TYPE (t_svf), INTENT(in) :: svf1,svf2
     7093      LOGICAL                  :: res
     7094      IF ( svf1%isurflt < svf2%isurflt  .OR.    &
     7095          (svf1%isurflt == svf2%isurflt  .AND.  svf1%isurfs < svf2%isurfs) )  THEN
     7096          res = .TRUE.
     7097      ELSE
     7098          res = .FALSE.
     7099      ENDIF
     7100    END FUNCTION svf_lt
     7101   
     7102 
     7103!-- quicksort.f -*-f90-*-
     7104!-- Author: t-nissie, adaptation J.Resler
     7105!-- License: GPLv3
     7106!-- Gist: https://gist.github.com/t-nissie/479f0f16966925fa29ea
     7107    RECURSIVE SUBROUTINE quicksort_svf(svfl, first, last)
     7108        IMPLICIT NONE
     7109        TYPE(t_svf), DIMENSION(:), INTENT(INOUT)  :: svfl
     7110        INTEGER(iwp), INTENT(IN)                  :: first, last
     7111        TYPE(t_svf)                               :: x, t
     7112        INTEGER(iwp)                              :: i, j
     7113
     7114        IF ( first>=last ) RETURN
     7115        x = svfl( (first+last) / 2 )
     7116        i = first
     7117        j = last
     7118        DO
     7119            DO while ( svf_lt(svfl(i),x) )
     7120               i=i+1
     7121            ENDDO
     7122            DO while ( svf_lt(x,svfl(j)) )
     7123                j=j-1
     7124            ENDDO
     7125            IF ( i >= j ) EXIT
     7126            t = svfl(i);  svfl(i) = svfl(j);  svfl(j) = t
     7127            i=i+1
     7128            j=j-1
     7129        ENDDO
     7130        IF ( first < i-1 ) CALL quicksort_svf(svfl, first, i-1)
     7131        IF ( j+1 < last )  CALL quicksort_svf(svfl, j+1, last)
     7132    END SUBROUTINE quicksort_svf
     7133
     7134   
     7135    PURE FUNCTION csf_lt(csf1,csf2) result (res)
     7136      TYPE (t_csf), INTENT(in) :: csf1,csf2
     7137      LOGICAL                  :: res
     7138      IF ( csf1%ip < csf2%ip  .OR.    &
     7139           (csf1%ip == csf2%ip  .AND.  csf1%itx < csf2%itx)  .OR.  &
     7140           (csf1%ip == csf2%ip  .AND.  csf1%itx == csf2%itx  .AND.  csf1%ity < csf2%ity)  .OR.  &
     7141           (csf1%ip == csf2%ip  .AND.  csf1%itx == csf2%itx  .AND.  csf1%ity == csf2%ity  .AND.   &
     7142            csf1%itz < csf2%itz)  .OR.  &
     7143           (csf1%ip == csf2%ip  .AND.  csf1%itx == csf2%itx  .AND.  csf1%ity == csf2%ity  .AND.   &
     7144            csf1%itz == csf2%itz  .AND.  csf1%isurfs < csf2%isurfs) )  THEN
     7145          res = .TRUE.
     7146      ELSE
     7147          res = .FALSE.
     7148      ENDIF
     7149    END FUNCTION csf_lt
     7150
     7151
     7152!-- quicksort.f -*-f90-*-
     7153!-- Author: t-nissie, adaptation J.Resler
     7154!-- License: GPLv3
     7155!-- Gist: https://gist.github.com/t-nissie/479f0f16966925fa29ea
     7156    RECURSIVE SUBROUTINE quicksort_csf(csfl, first, last)
     7157        IMPLICIT NONE
     7158        TYPE(t_csf), DIMENSION(:), INTENT(INOUT)  :: csfl
     7159        INTEGER(iwp), INTENT(IN)                  :: first, last
     7160        TYPE(t_csf)                               :: x, t
     7161        INTEGER(iwp)                              :: i, j
     7162
     7163        IF ( first>=last ) RETURN
     7164        x = csfl( (first+last)/2 )
     7165        i = first
     7166        j = last
     7167        DO
     7168            DO while ( csf_lt(csfl(i),x) )
     7169                i=i+1
     7170            ENDDO
     7171            DO while ( csf_lt(x,csfl(j)) )
     7172                j=j-1
     7173            ENDDO
     7174            IF ( i >= j ) EXIT
     7175            t = csfl(i);  csfl(i) = csfl(j);  csfl(j) = t
     7176            i=i+1
     7177            j=j-1
     7178        ENDDO
     7179        IF ( first < i-1 ) CALL quicksort_csf(csfl, first, i-1)
     7180        IF ( j+1 < last )  CALL quicksort_csf(csfl, j+1, last)
     7181    END SUBROUTINE quicksort_csf
     7182
     7183   
     7184    SUBROUTINE merge_and_grow_csf(newsize)
     7185        INTEGER(iwp), INTENT(in)                :: newsize  !< new array size after grow, must be >= ncsfl
     7186                                                            !< or -1 to shrink to minimum
     7187        INTEGER(iwp)                            :: iread, iwrite
     7188        TYPE(t_csf), DIMENSION(:), POINTER      :: acsfnew
     7189        CHARACTER(100)                          :: msg
     7190
     7191        IF ( newsize == -1 )  THEN
     7192!--         merge in-place
     7193            acsfnew => acsf
     7194        ELSE
     7195!--         allocate new array
     7196            IF ( mcsf == 0 )  THEN
     7197                ALLOCATE( acsf1(newsize) )
     7198                acsfnew => acsf1
     7199            ELSE
     7200                ALLOCATE( acsf2(newsize) )
     7201                acsfnew => acsf2
     7202            ENDIF
     7203        ENDIF
     7204
     7205        IF ( ncsfl >= 1 )  THEN
     7206!--         sort csf in place (quicksort)
     7207            CALL quicksort_csf(acsf,1,ncsfl)
     7208
     7209!--         while moving to a new array, aggregate canopy sink factor records with identical box & source
     7210            acsfnew(1) = acsf(1)
     7211            iwrite = 1
     7212            DO iread = 2, ncsfl
     7213!--             here acsf(kcsf) already has values from acsf(icsf)
     7214                IF ( acsfnew(iwrite)%itx == acsf(iread)%itx &
     7215                         .AND.  acsfnew(iwrite)%ity == acsf(iread)%ity &
     7216                         .AND.  acsfnew(iwrite)%itz == acsf(iread)%itz &
     7217                         .AND.  acsfnew(iwrite)%isurfs == acsf(iread)%isurfs )  THEN
     7218!--                 We could simply take either first or second rtransp, both are valid. As a very simple heuristic about which ray
     7219!--                 probably passes nearer the center of the target box, we choose DIF from the entry with greater CSF, since that
     7220!--                 might mean that the traced beam passes longer through the canopy box.
     7221                    IF ( acsfnew(iwrite)%rsvf < acsf(iread)%rsvf )  THEN
     7222                        acsfnew(iwrite)%rtransp = acsf(iread)%rtransp
     7223                    ENDIF
     7224                    acsfnew(iwrite)%rsvf = acsfnew(iwrite)%rsvf + acsf(iread)%rsvf
     7225!--                 advance reading index, keep writing index
     7226                ELSE
     7227!--                 not identical, just advance and copy
     7228                    iwrite = iwrite + 1
     7229                    acsfnew(iwrite) = acsf(iread)
     7230                ENDIF
     7231            ENDDO
     7232            ncsfl = iwrite
     7233        ENDIF
     7234
     7235        IF ( newsize == -1 )  THEN
     7236!--         allocate new array and copy shrinked data
     7237            IF ( mcsf == 0 )  THEN
     7238                ALLOCATE( acsf1(ncsfl) )
     7239                acsf1(1:ncsfl) = acsf2(1:ncsfl)
     7240            ELSE
     7241                ALLOCATE( acsf2(ncsfl) )
     7242                acsf2(1:ncsfl) = acsf1(1:ncsfl)
     7243            ENDIF
     7244        ENDIF
     7245
     7246!--     deallocate old array
     7247        IF ( mcsf == 0 )  THEN
     7248            mcsf = 1
     7249            acsf => acsf1
     7250            DEALLOCATE( acsf2 )
     7251        ELSE
     7252            mcsf = 0
     7253            acsf => acsf2
     7254            DEALLOCATE( acsf1 )
     7255        ENDIF
     7256        ncsfla = newsize
     7257
     7258!         WRITE(msg,'(A,2I12)') 'Grow acsf2:',ncsfl,ncsfla
     7259!         CALL radiation_write_debug_log( msg )
     7260
     7261    END SUBROUTINE merge_and_grow_csf
     7262
     7263   
     7264!-- quicksort.f -*-f90-*-
     7265!-- Author: t-nissie, adaptation J.Resler
     7266!-- License: GPLv3
     7267!-- Gist: https://gist.github.com/t-nissie/479f0f16966925fa29ea
     7268    RECURSIVE SUBROUTINE quicksort_csf2(kpcsflt, pcsflt, first, last)
     7269        IMPLICIT NONE
     7270        INTEGER(iwp), DIMENSION(:,:), INTENT(INOUT)  :: kpcsflt
     7271        REAL(wp), DIMENSION(:,:), INTENT(INOUT)      :: pcsflt
     7272        INTEGER(iwp), INTENT(IN)                     :: first, last
     7273        REAL(wp), DIMENSION(ndcsf)                   :: t2
     7274        INTEGER(iwp), DIMENSION(kdcsf)               :: x, t1
     7275        INTEGER(iwp)                                 :: i, j
     7276
     7277        IF ( first>=last ) RETURN
     7278        x = kpcsflt(:, (first+last)/2 )
     7279        i = first
     7280        j = last
     7281        DO
     7282            DO while ( csf_lt2(kpcsflt(:,i),x) )
     7283                i=i+1
     7284            ENDDO
     7285            DO while ( csf_lt2(x,kpcsflt(:,j)) )
     7286                j=j-1
     7287            ENDDO
     7288            IF ( i >= j ) EXIT
     7289            t1 = kpcsflt(:,i);  kpcsflt(:,i) = kpcsflt(:,j);  kpcsflt(:,j) = t1
     7290            t2 = pcsflt(:,i);  pcsflt(:,i) = pcsflt(:,j);  pcsflt(:,j) = t2
     7291            i=i+1
     7292            j=j-1
     7293        ENDDO
     7294        IF ( first < i-1 ) CALL quicksort_csf2(kpcsflt, pcsflt, first, i-1)
     7295        IF ( j+1 < last )  CALL quicksort_csf2(kpcsflt, pcsflt, j+1, last)
     7296    END SUBROUTINE quicksort_csf2
     7297   
     7298
     7299    PURE FUNCTION csf_lt2(item1, item2) result(res)
     7300        INTEGER(iwp), DIMENSION(kdcsf), INTENT(in)  :: item1, item2
     7301        LOGICAL                                     :: res
     7302        res = ( (item1(3) < item2(3))                                                        &
     7303             .OR.  (item1(3) == item2(3)  .AND.  item1(2) < item2(2))                            &
     7304             .OR.  (item1(3) == item2(3)  .AND.  item1(2) == item2(2)  .AND.  item1(1) < item2(1)) &
     7305             .OR.  (item1(3) == item2(3)  .AND.  item1(2) == item2(2)  .AND.  item1(1) == item2(1) &
     7306                 .AND.  item1(4) < item2(4)) )
     7307    END FUNCTION csf_lt2
     7308
     7309    PURE FUNCTION searchsorted(athresh, val) result(ind)
     7310        REAL(wp), DIMENSION(:), INTENT(IN)  :: athresh
     7311        REAL(wp), INTENT(IN)                :: val
     7312        INTEGER(iwp)                        :: ind
     7313        INTEGER(iwp)                        :: i
     7314
     7315        DO i = LBOUND(athresh, 1), UBOUND(athresh, 1)
     7316            IF ( val < athresh(i) ) THEN
     7317                ind = i - 1
     7318                RETURN
     7319            ENDIF
     7320        ENDDO
     7321        ind = UBOUND(athresh, 1)
     7322    END FUNCTION searchsorted
    67257323
    67267324!------------------------------------------------------------------------------!
     
    67967394                  swd_gridbox(1) = surfinswdif(l)
    67977395
    6798                CASE (isky,idown_a)                         !- gridbox down_facing face
     7396               CASE (idown_a)                         !- gridbox down_facing face
    67997397                  sw_gridbox(2) = surfinsw(l)
    68007398                  lw_gridbox(2) = surfinlw(l)
    68017399                  swd_gridbox(2) = surfinswdif(l)
    68027400
    6803                CASE (inorth_u,inorth_l,inorth_a,isouth_b)  !- gridbox north_facing face
     7401               CASE (inorth_u,inorth_l,inorth_a)  !- gridbox north_facing face
    68047402                  sw_gridbox(3) = surfinsw(l)
    68057403                  lw_gridbox(3) = surfinlw(l)
    68067404                  swd_gridbox(3) = surfinswdif(l)
    68077405
    6808                CASE (isouth_u,isouth_l,isouth_a,inorth_b)  !- gridbox south_facing face
     7406               CASE (isouth_u,isouth_l,isouth_a)  !- gridbox south_facing face
    68097407                  sw_gridbox(4) = surfinsw(l)
    68107408                  lw_gridbox(4) = surfinlw(l)
    68117409                  swd_gridbox(4) = surfinswdif(l)
    68127410
    6813                CASE (ieast_u,ieast_l,ieast_a,iwest_b)      !- gridbox east_facing face
     7411               CASE (ieast_u,ieast_l,ieast_a)      !- gridbox east_facing face
    68147412                  sw_gridbox(5) = surfinsw(l)
    68157413                  lw_gridbox(5) = surfinlw(l)
    68167414                  swd_gridbox(5) = surfinswdif(l)
    68177415
    6818                CASE (iwest_u,iwest_l,iwest_a,ieast_b)      !- gridbox west_facing face
     7416               CASE (iwest_u,iwest_l,iwest_a)      !- gridbox west_facing face
    68197417                  sw_gridbox(6) = surfinsw(l)
    68207418                  lw_gridbox(6) = surfinlw(l)
     
    68417439       
    68427440    END SUBROUTINE radiation_radflux_gridbox
    6843 
    6844 
    6845 !------------------------------------------------------------------------------!
    6846 !
    6847 ! Description:
    6848 ! ------------
    6849 !> Block of auxiliary subroutines:
    6850 !> 1. quicksort and corresponding comparison
    6851 !> 2. merge_and_grow_csf for implementation of "dynamical growing"
    6852 !>    array for csf
    6853 !------------------------------------------------------------------------------!   
    6854     PURE FUNCTION svf_lt(svf1,svf2) result (res)
    6855       TYPE (t_svf), INTENT(in) :: svf1,svf2
    6856       LOGICAL                  :: res
    6857       IF ( svf1%isurflt < svf2%isurflt  .OR.    &
    6858           (svf1%isurflt == svf2%isurflt  .AND.  svf1%isurfs < svf2%isurfs) )  THEN
    6859           res = .TRUE.
    6860       ELSE
    6861           res = .FALSE.
    6862       ENDIF
    6863     END FUNCTION svf_lt
    6864    
    6865  
    6866 !-- quicksort.f -*-f90-*-
    6867 !-- Author: t-nissie, adaptation J.Resler
    6868 !-- License: GPLv3
    6869 !-- Gist: https://gist.github.com/t-nissie/479f0f16966925fa29ea
    6870     RECURSIVE SUBROUTINE quicksort_svf(svfl, first, last)
    6871         IMPLICIT NONE
    6872         TYPE(t_svf), DIMENSION(:), INTENT(INOUT)  :: svfl
    6873         INTEGER(iwp), INTENT(IN)                  :: first, last
    6874         TYPE(t_svf)                               :: x, t
    6875         INTEGER(iwp)                              :: i, j
    6876 
    6877         IF ( first>=last ) RETURN
    6878         x = svfl( (first+last) / 2 )
    6879         i = first
    6880         j = last
    6881         DO
    6882             DO while ( svf_lt(svfl(i),x) )
    6883                 i=i+1
    6884             ENDDO
    6885             DO while ( svf_lt(x,svfl(j)) )
    6886                 j=j-1
    6887             ENDDO
    6888             IF ( i >= j ) EXIT
    6889             t = svfl(i);  svfl(i) = svfl(j);  svfl(j) = t
    6890             i=i+1
    6891             j=j-1
    6892         ENDDO
    6893         IF ( first < i-1 ) CALL quicksort_svf(svfl, first, i-1)
    6894         IF ( j+1 < last )  CALL quicksort_svf(svfl, j+1, last)
    6895     END SUBROUTINE quicksort_svf
    6896 
    6897    
    6898     PURE FUNCTION csf_lt(csf1,csf2) result (res)
    6899       TYPE (t_csf), INTENT(in) :: csf1,csf2
    6900       LOGICAL                  :: res
    6901       IF ( csf1%ip < csf2%ip  .OR.    &
    6902            (csf1%ip == csf2%ip  .AND.  csf1%itx < csf2%itx)  .OR.  &
    6903            (csf1%ip == csf2%ip  .AND.  csf1%itx == csf2%itx  .AND.  csf1%ity < csf2%ity)  .OR.  &
    6904            (csf1%ip == csf2%ip  .AND.  csf1%itx == csf2%itx  .AND.  csf1%ity == csf2%ity  .AND.   &
    6905             csf1%itz < csf2%itz)  .OR.  &
    6906            (csf1%ip == csf2%ip  .AND.  csf1%itx == csf2%itx  .AND.  csf1%ity == csf2%ity  .AND.   &
    6907             csf1%itz == csf2%itz  .AND.  csf1%isurfs < csf2%isurfs) )  THEN
    6908           res = .TRUE.
    6909       ELSE
    6910           res = .FALSE.
    6911       ENDIF
    6912     END FUNCTION csf_lt
    6913 
    6914 
    6915 !-- quicksort.f -*-f90-*-
    6916 !-- Author: t-nissie, adaptation J.Resler
    6917 !-- License: GPLv3
    6918 !-- Gist: https://gist.github.com/t-nissie/479f0f16966925fa29ea
    6919     RECURSIVE SUBROUTINE quicksort_csf(csfl, first, last)
    6920         IMPLICIT NONE
    6921         TYPE(t_csf), DIMENSION(:), INTENT(INOUT)  :: csfl
    6922         INTEGER(iwp), INTENT(IN)                  :: first, last
    6923         TYPE(t_csf)                               :: x, t
    6924         INTEGER(iwp)                              :: i, j
    6925 
    6926         IF ( first>=last ) RETURN
    6927         x = csfl( (first+last)/2 )
    6928         i = first
    6929         j = last
    6930         DO
    6931             DO while ( csf_lt(csfl(i),x) )
    6932                 i=i+1
    6933             ENDDO
    6934             DO while ( csf_lt(x,csfl(j)) )
    6935                 j=j-1
    6936             ENDDO
    6937             IF ( i >= j ) EXIT
    6938             t = csfl(i);  csfl(i) = csfl(j);  csfl(j) = t
    6939             i=i+1
    6940             j=j-1
    6941         ENDDO
    6942         IF ( first < i-1 ) CALL quicksort_csf(csfl, first, i-1)
    6943         IF ( j+1 < last )  CALL quicksort_csf(csfl, j+1, last)
    6944     END SUBROUTINE quicksort_csf
    6945 
    6946    
    6947     SUBROUTINE merge_and_grow_csf(newsize)
    6948         INTEGER(iwp), INTENT(in)                :: newsize  !< new array size after grow, must be >= ncsfl
    6949                                                             !< or -1 to shrink to minimum
    6950         INTEGER(iwp)                            :: iread, iwrite
    6951         TYPE(t_csf), DIMENSION(:), POINTER      :: acsfnew
    6952 
    6953         IF ( newsize == -1 )  THEN
    6954 !--         merge in-place
    6955             acsfnew => acsf
    6956         ELSE
    6957 !--         allocate new array
    6958             IF ( mcsf == 0 )  THEN
    6959                 ALLOCATE( acsf1(newsize) )
    6960                 acsfnew => acsf1
    6961             ELSE
    6962                 ALLOCATE( acsf2(newsize) )
    6963                 acsfnew => acsf2
    6964             ENDIF
    6965         ENDIF
    6966 
    6967         IF ( ncsfl >= 1 )  THEN
    6968 !--         sort csf in place (quicksort)
    6969             CALL quicksort_csf(acsf,1,ncsfl)
    6970 
    6971 !--         while moving to a new array, aggregate canopy sink factor records with identical box & source
    6972             acsfnew(1) = acsf(1)
    6973             iwrite = 1
    6974             DO iread = 2, ncsfl
    6975 !--             here acsf(kcsf) already has values from acsf(icsf)
    6976                 IF ( acsfnew(iwrite)%itx == acsf(iread)%itx &
    6977                          .AND.  acsfnew(iwrite)%ity == acsf(iread)%ity &
    6978                          .AND.  acsfnew(iwrite)%itz == acsf(iread)%itz &
    6979                          .AND.  acsfnew(iwrite)%isurfs == acsf(iread)%isurfs )  THEN
    6980 !--                 We could simply take either first or second rtransp, both are valid. As a very simple heuristic about which ray
    6981 !--                 probably passes nearer the center of the target box, we choose DIF from the entry with greater CSF, since that
    6982 !--                 might mean that the traced beam passes longer through the canopy box.
    6983                     IF ( acsfnew(iwrite)%rsvf < acsf(iread)%rsvf )  THEN
    6984                         acsfnew(iwrite)%rtransp = acsf(iread)%rtransp
    6985                     ENDIF
    6986                     acsfnew(iwrite)%rsvf = acsfnew(iwrite)%rsvf + acsf(iread)%rsvf
    6987 !--                 advance reading index, keep writing index
    6988                 ELSE
    6989 !--                 not identical, just advance and copy
    6990                     iwrite = iwrite + 1
    6991                     acsfnew(iwrite) = acsf(iread)
    6992                 ENDIF
    6993             ENDDO
    6994             ncsfl = iwrite
    6995         ENDIF
    6996 
    6997         IF ( newsize == -1 )  THEN
    6998 !--         allocate new array and copy shrinked data
    6999             IF ( mcsf == 0 )  THEN
    7000                 ALLOCATE( acsf1(ncsfl) )
    7001                 acsf1(1:ncsfl) = acsf2(1:ncsfl)
    7002             ELSE
    7003                 ALLOCATE( acsf2(ncsfl) )
    7004                 acsf2(1:ncsfl) = acsf1(1:ncsfl)
    7005             ENDIF
    7006         ENDIF
    7007 
    7008 !--     deallocate old array
    7009         IF ( mcsf == 0 )  THEN
    7010             mcsf = 1
    7011             acsf => acsf1
    7012             DEALLOCATE( acsf2 )
    7013         ELSE
    7014             mcsf = 0
    7015             acsf => acsf2
    7016             DEALLOCATE( acsf1 )
    7017         ENDIF
    7018         ncsfla = newsize
    7019     END SUBROUTINE merge_and_grow_csf
    7020 
    7021    
    7022 !-- quicksort.f -*-f90-*-
    7023 !-- Author: t-nissie, adaptation J.Resler
    7024 !-- License: GPLv3
    7025 !-- Gist: https://gist.github.com/t-nissie/479f0f16966925fa29ea
    7026     RECURSIVE SUBROUTINE quicksort_csf2(kpcsflt, pcsflt, first, last)
    7027         IMPLICIT NONE
    7028         INTEGER(iwp), DIMENSION(:,:), INTENT(INOUT)  :: kpcsflt
    7029         REAL(wp), DIMENSION(:,:), INTENT(INOUT)      :: pcsflt
    7030         INTEGER(iwp), INTENT(IN)                     :: first, last
    7031         REAL(wp), DIMENSION(ndcsf)                   :: t2
    7032         INTEGER(iwp), DIMENSION(kdcsf)               :: x, t1
    7033         INTEGER(iwp)                                 :: i, j
    7034 
    7035         IF ( first>=last ) RETURN
    7036         x = kpcsflt(:, (first+last)/2 )
    7037         i = first
    7038         j = last
    7039         DO
    7040             DO while ( csf_lt2(kpcsflt(:,i),x) )
    7041                 i=i+1
    7042             ENDDO
    7043             DO while ( csf_lt2(x,kpcsflt(:,j)) )
    7044                 j=j-1
    7045             ENDDO
    7046             IF ( i >= j ) EXIT
    7047             t1 = kpcsflt(:,i);  kpcsflt(:,i) = kpcsflt(:,j);  kpcsflt(:,j) = t1
    7048             t2 = pcsflt(:,i);  pcsflt(:,i) = pcsflt(:,j);  pcsflt(:,j) = t2
    7049             i=i+1
    7050             j=j-1
    7051         ENDDO
    7052         IF ( first < i-1 ) CALL quicksort_csf2(kpcsflt, pcsflt, first, i-1)
    7053         IF ( j+1 < last )  CALL quicksort_csf2(kpcsflt, pcsflt, j+1, last)
    7054     END SUBROUTINE quicksort_csf2
    7055    
    7056 
    7057     PURE FUNCTION csf_lt2(item1, item2) result(res)
    7058         INTEGER(iwp), DIMENSION(kdcsf), INTENT(in)  :: item1, item2
    7059         LOGICAL                                     :: res
    7060         res = ( (item1(3) < item2(3))                                                        &
    7061              .OR.  (item1(3) == item2(3)  .AND.  item1(2) < item2(2))                            &
    7062              .OR.  (item1(3) == item2(3)  .AND.  item1(2) == item2(2)  .AND.  item1(1) < item2(1)) &
    7063              .OR.  (item1(3) == item2(3)  .AND.  item1(2) == item2(2)  .AND.  item1(1) == item2(1) &
    7064                  .AND.  item1(4) < item2(4)) )
    7065     END FUNCTION csf_lt2
    70667441
    70677442!------------------------------------------------------------------------------!
     
    80558430
    80568431!------------------------------------------------------------------------------!
    8057 !
    80588432! Description:
    80598433! ------------
    8060 !> Subroutine writes the respective restart data
     8434!> Subroutine writes local (subdomain) restart data
    80618435!------------------------------------------------------------------------------!
    80628436 SUBROUTINE radiation_wrd_local
     
    81548528 END SUBROUTINE radiation_wrd_local
    81558529
    8156 
    8157 SUBROUTINE radiation_rrd_local( i, k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,     &
     8530!------------------------------------------------------------------------------!
     8531! Description:
     8532! ------------
     8533!> Subroutine reads local (subdomain) restart data
     8534!------------------------------------------------------------------------------!
     8535 SUBROUTINE radiation_rrd_local( i, k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,    &
    81588536                                nxr_on_file, nynf, nync, nyn_on_file, nysf,    &
    81598537                                nysc, nys_on_file, tmp_2d, tmp_3d, found )
     
    81988576
    81998577
    8200        SELECT CASE ( restart_string(1:length) )
    8201 
    8202            CASE ( 'rad_net_av' )
    8203               IF ( .NOT. ALLOCATED( rad_net_av ) )  THEN
    8204                  ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
    8205               ENDIF 
    8206               IF ( k == 1 )  READ ( 13 )  tmp_2d
    8207               rad_net_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =           &
    8208                             tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    8209            CASE ( 'rad_lw_in' )
    8210               IF ( .NOT. ALLOCATED( rad_lw_in ) )  THEN
    8211                  IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    8212                       radiation_scheme == 'constant')  THEN
    8213                     ALLOCATE( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
    8214                  ELSE
    8215                     ALLOCATE( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    8216                  ENDIF
    8217               ENDIF 
    8218               IF ( k == 1 )  THEN
    8219                  IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    8220                       radiation_scheme == 'constant')  THEN
    8221                     READ ( 13 )  tmp_3d2
    8222                     rad_lw_in(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =   &
    8223                        tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    8224                  ELSE
    8225                     READ ( 13 )  tmp_3d
    8226                     rad_lw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =     &
    8227                         tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    8228                  ENDIF
    8229               ENDIF
    8230 
    8231            CASE ( 'rad_lw_in_av' )
    8232               IF ( .NOT. ALLOCATED( rad_lw_in_av ) )  THEN
    8233                  IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    8234                       radiation_scheme == 'constant')  THEN
    8235                     ALLOCATE( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
    8236                  ELSE
    8237                     ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    8238                  ENDIF
    8239               ENDIF 
    8240               IF ( k == 1 )  THEN
    8241                  IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    8242                       radiation_scheme == 'constant')  THEN
    8243                     READ ( 13 )  tmp_3d2
    8244                     rad_lw_in_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
    8245                         tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    8246                  ELSE
    8247                     READ ( 13 )  tmp_3d
    8248                     rad_lw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =  &
    8249                         tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    8250                  ENDIF
    8251               ENDIF
    8252 
    8253            CASE ( 'rad_lw_out' )
    8254               IF ( .NOT. ALLOCATED( rad_lw_out ) )  THEN
    8255                  IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    8256                       radiation_scheme == 'constant')  THEN
    8257                     ALLOCATE( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
    8258                  ELSE
    8259                     ALLOCATE( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    8260                  ENDIF
    8261               ENDIF 
    8262               IF ( k == 1 )  THEN
    8263                  IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    8264                       radiation_scheme == 'constant')  THEN
    8265                     READ ( 13 )  tmp_3d2
    8266                     rad_lw_out(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =  &
    8267                         tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    8268                  ELSE
    8269                     READ ( 13 )  tmp_3d
    8270                     rad_lw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =    &
    8271                         tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    8272                  ENDIF
    8273               ENDIF
    8274 
    8275            CASE ( 'rad_lw_out_av' )
    8276               IF ( .NOT. ALLOCATED( rad_lw_out_av ) )  THEN
    8277                  IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    8278                       radiation_scheme == 'constant')  THEN
    8279                     ALLOCATE( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
    8280                  ELSE
    8281                     ALLOCATE( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    8282                  ENDIF
    8283               ENDIF 
    8284               IF ( k == 1 )  THEN
    8285                  IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    8286                       radiation_scheme == 'constant')  THEN
    8287                     READ ( 13 )  tmp_3d2
    8288                     rad_lw_out_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) &
    8289                        = tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    8290                  ELSE
    8291                     READ ( 13 )  tmp_3d
    8292                     rad_lw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
    8293                         tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    8294                  ENDIF
    8295               ENDIF
    8296 
    8297            CASE ( 'rad_lw_cs_hr' )
    8298               IF ( .NOT. ALLOCATED( rad_lw_cs_hr ) )  THEN
    8299                  ALLOCATE( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    8300               ENDIF
    8301               IF ( k == 1 )  READ ( 13 )  tmp_3d
    8302               rad_lw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =        &
    8303                       tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    8304 
    8305            CASE ( 'rad_lw_cs_hr_av' )
    8306               IF ( .NOT. ALLOCATED( rad_lw_cs_hr_av ) )  THEN
    8307                  ALLOCATE( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    8308               ENDIF
    8309               IF ( k == 1 )  READ ( 13 )  tmp_3d
    8310               rad_lw_cs_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =     &
    8311                       tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    8312 
    8313            CASE ( 'rad_lw_hr' )
    8314               IF ( .NOT. ALLOCATED( rad_lw_hr ) )  THEN
    8315                  ALLOCATE( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    8316               ENDIF
    8317               IF ( k == 1 )  READ ( 13 )  tmp_3d
    8318               rad_lw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =           &
    8319                       tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    8320 
    8321            CASE ( 'rad_lw_hr_av' )
    8322               IF ( .NOT. ALLOCATED( rad_lw_hr_av ) )  THEN
    8323                  ALLOCATE( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    8324               ENDIF
    8325               IF ( k == 1 )  READ ( 13 )  tmp_3d
    8326               rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =        &
    8327                       tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    8328 
    8329            CASE ( 'rad_sw_in' )
    8330               IF ( .NOT. ALLOCATED( rad_sw_in ) )  THEN
    8331                  IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    8332                       radiation_scheme == 'constant')  THEN
    8333                     ALLOCATE( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
    8334                  ELSE
    8335                     ALLOCATE( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    8336                  ENDIF
    8337               ENDIF 
    8338               IF ( k == 1 )  THEN
    8339                  IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    8340                       radiation_scheme == 'constant')  THEN
    8341                     READ ( 13 )  tmp_3d2
    8342                     rad_sw_in(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =   &
    8343                         tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    8344                  ELSE
    8345                     READ ( 13 )  tmp_3d
    8346                     rad_sw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =     &
    8347                         tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    8348                  ENDIF
    8349               ENDIF
    8350 
    8351            CASE ( 'rad_sw_in_av' )
    8352               IF ( .NOT. ALLOCATED( rad_sw_in_av ) )  THEN
    8353                  IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    8354                       radiation_scheme == 'constant')  THEN
    8355                     ALLOCATE( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
    8356                  ELSE
    8357                     ALLOCATE( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    8358                  ENDIF
    8359               ENDIF 
    8360               IF ( k == 1 )  THEN
    8361                  IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    8362                       radiation_scheme == 'constant')  THEN
    8363                     READ ( 13 )  tmp_3d2
    8364                     rad_sw_in_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
    8365                         tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    8366                  ELSE
    8367                     READ ( 13 )  tmp_3d
    8368                     rad_sw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =  &
    8369                         tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    8370                  ENDIF
    8371               ENDIF
    8372 
    8373            CASE ( 'rad_sw_out' )
    8374               IF ( .NOT. ALLOCATED( rad_sw_out ) )  THEN
    8375                  IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    8376                       radiation_scheme == 'constant')  THEN
    8377                     ALLOCATE( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
    8378                  ELSE
    8379                     ALLOCATE( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    8380                  ENDIF
    8381               ENDIF 
    8382               IF ( k == 1 )  THEN
    8383                  IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    8384                       radiation_scheme == 'constant')  THEN
    8385                     READ ( 13 )  tmp_3d2
    8386                     rad_sw_out(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =  &
    8387                         tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    8388                  ELSE
    8389                     READ ( 13 )  tmp_3d
    8390                     rad_sw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =    &
    8391                         tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    8392                  ENDIF
    8393               ENDIF
    8394 
    8395            CASE ( 'rad_sw_out_av' )
    8396               IF ( .NOT. ALLOCATED( rad_sw_out_av ) )  THEN
    8397                  IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    8398                       radiation_scheme == 'constant')  THEN
    8399                     ALLOCATE( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
    8400                  ELSE
    8401                     ALLOCATE( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    8402                  ENDIF
    8403               ENDIF 
    8404               IF ( k == 1 )  THEN
    8405                  IF ( radiation_scheme == 'clear-sky'  .OR.                    &
    8406                       radiation_scheme == 'constant')  THEN
    8407                     READ ( 13 )  tmp_3d2
    8408                     rad_sw_out_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) &
    8409                        = tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    8410                  ELSE
    8411                     READ ( 13 )  tmp_3d
    8412                     rad_sw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
    8413                         tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    8414                  ENDIF
    8415               ENDIF
    8416 
    8417            CASE ( 'rad_sw_cs_hr' )
    8418               IF ( .NOT. ALLOCATED( rad_sw_cs_hr ) )  THEN
    8419                  ALLOCATE( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    8420               ENDIF
    8421               IF ( k == 1 )  READ ( 13 )  tmp_3d
    8422               rad_sw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =        &
    8423                       tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    8424 
    8425            CASE ( 'rad_sw_cs_hr_av' )
    8426               IF ( .NOT. ALLOCATED( rad_sw_cs_hr_av ) )  THEN
    8427                  ALLOCATE( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    8428               ENDIF
    8429               IF ( k == 1 )  READ ( 13 )  tmp_3d
    8430               rad_sw_cs_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =     &
    8431                       tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    8432 
    8433            CASE ( 'rad_sw_hr' )
    8434               IF ( .NOT. ALLOCATED( rad_sw_hr ) )  THEN
    8435                  ALLOCATE( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    8436               ENDIF
    8437               IF ( k == 1 )  READ ( 13 )  tmp_3d
    8438               rad_sw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =           &
    8439                       tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    8440 
    8441            CASE ( 'rad_sw_hr_av' )
    8442               IF ( .NOT. ALLOCATED( rad_sw_hr_av ) )  THEN
    8443                  ALLOCATE( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    8444               ENDIF
    8445               IF ( k == 1 )  READ ( 13 )  tmp_3d
    8446               rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =        &
    8447                       tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
    8448 
    8449            CASE DEFAULT
    8450 
    8451               found = .FALSE.
    8452 
    8453        END SELECT
    8454 
     8578    SELECT CASE ( restart_string(1:length) )
     8579
     8580       CASE ( 'rad_net_av' )
     8581          IF ( .NOT. ALLOCATED( rad_net_av ) )  THEN
     8582             ALLOCATE( rad_net_av(nysg:nyng,nxlg:nxrg) )
     8583          ENDIF 
     8584          IF ( k == 1 )  READ ( 13 )  tmp_2d
     8585          rad_net_av(nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp)  =           &
     8586                        tmp_2d(nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     8587       CASE ( 'rad_lw_in' )
     8588          IF ( .NOT. ALLOCATED( rad_lw_in ) )  THEN
     8589             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
     8590                  radiation_scheme == 'constant')  THEN
     8591                ALLOCATE( rad_lw_in(0:0,nysg:nyng,nxlg:nxrg) )
     8592             ELSE
     8593                ALLOCATE( rad_lw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     8594             ENDIF
     8595          ENDIF 
     8596          IF ( k == 1 )  THEN
     8597             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
     8598                  radiation_scheme == 'constant')  THEN
     8599                READ ( 13 )  tmp_3d2
     8600                rad_lw_in(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =   &
     8601                   tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     8602             ELSE
     8603                READ ( 13 )  tmp_3d
     8604                rad_lw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =     &
     8605                    tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     8606             ENDIF
     8607          ENDIF
     8608
     8609       CASE ( 'rad_lw_in_av' )
     8610          IF ( .NOT. ALLOCATED( rad_lw_in_av ) )  THEN
     8611             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
     8612                  radiation_scheme == 'constant')  THEN
     8613                ALLOCATE( rad_lw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
     8614             ELSE
     8615                ALLOCATE( rad_lw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     8616             ENDIF
     8617          ENDIF 
     8618          IF ( k == 1 )  THEN
     8619             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
     8620                  radiation_scheme == 'constant')  THEN
     8621                READ ( 13 )  tmp_3d2
     8622                rad_lw_in_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
     8623                    tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     8624             ELSE
     8625                READ ( 13 )  tmp_3d
     8626                rad_lw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =  &
     8627                    tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     8628             ENDIF
     8629          ENDIF
     8630
     8631       CASE ( 'rad_lw_out' )
     8632          IF ( .NOT. ALLOCATED( rad_lw_out ) )  THEN
     8633             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
     8634                  radiation_scheme == 'constant')  THEN
     8635                ALLOCATE( rad_lw_out(0:0,nysg:nyng,nxlg:nxrg) )
     8636             ELSE
     8637                ALLOCATE( rad_lw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     8638             ENDIF
     8639          ENDIF 
     8640          IF ( k == 1 )  THEN
     8641             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
     8642                  radiation_scheme == 'constant')  THEN
     8643                READ ( 13 )  tmp_3d2
     8644                rad_lw_out(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =  &
     8645                    tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     8646             ELSE
     8647                READ ( 13 )  tmp_3d
     8648                rad_lw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =    &
     8649                    tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     8650             ENDIF
     8651          ENDIF
     8652
     8653       CASE ( 'rad_lw_out_av' )
     8654          IF ( .NOT. ALLOCATED( rad_lw_out_av ) )  THEN
     8655             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
     8656                  radiation_scheme == 'constant')  THEN
     8657                ALLOCATE( rad_lw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
     8658             ELSE
     8659                ALLOCATE( rad_lw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     8660             ENDIF
     8661          ENDIF 
     8662          IF ( k == 1 )  THEN
     8663             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
     8664                  radiation_scheme == 'constant')  THEN
     8665                READ ( 13 )  tmp_3d2
     8666                rad_lw_out_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) &
     8667                   = tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     8668             ELSE
     8669                READ ( 13 )  tmp_3d
     8670                rad_lw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     8671                    tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     8672             ENDIF
     8673          ENDIF
     8674
     8675       CASE ( 'rad_lw_cs_hr' )
     8676          IF ( .NOT. ALLOCATED( rad_lw_cs_hr ) )  THEN
     8677             ALLOCATE( rad_lw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     8678          ENDIF
     8679          IF ( k == 1 )  READ ( 13 )  tmp_3d
     8680          rad_lw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =        &
     8681                  tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     8682
     8683       CASE ( 'rad_lw_cs_hr_av' )
     8684          IF ( .NOT. ALLOCATED( rad_lw_cs_hr_av ) )  THEN
     8685             ALLOCATE( rad_lw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     8686          ENDIF
     8687          IF ( k == 1 )  READ ( 13 )  tmp_3d
     8688          rad_lw_cs_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =     &
     8689                  tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     8690
     8691       CASE ( 'rad_lw_hr' )
     8692          IF ( .NOT. ALLOCATED( rad_lw_hr ) )  THEN
     8693             ALLOCATE( rad_lw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     8694          ENDIF
     8695          IF ( k == 1 )  READ ( 13 )  tmp_3d
     8696          rad_lw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =           &
     8697                  tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     8698
     8699       CASE ( 'rad_lw_hr_av' )
     8700          IF ( .NOT. ALLOCATED( rad_lw_hr_av ) )  THEN
     8701             ALLOCATE( rad_lw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     8702          ENDIF
     8703          IF ( k == 1 )  READ ( 13 )  tmp_3d
     8704          rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =        &
     8705                  tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     8706
     8707       CASE ( 'rad_sw_in' )
     8708          IF ( .NOT. ALLOCATED( rad_sw_in ) )  THEN
     8709             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
     8710                  radiation_scheme == 'constant')  THEN
     8711                ALLOCATE( rad_sw_in(0:0,nysg:nyng,nxlg:nxrg) )
     8712             ELSE
     8713                ALLOCATE( rad_sw_in(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     8714             ENDIF
     8715          ENDIF 
     8716          IF ( k == 1 )  THEN
     8717             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
     8718                  radiation_scheme == 'constant')  THEN
     8719                READ ( 13 )  tmp_3d2
     8720                rad_sw_in(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =   &
     8721                    tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     8722             ELSE
     8723                READ ( 13 )  tmp_3d
     8724                rad_sw_in(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =     &
     8725                    tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     8726             ENDIF
     8727          ENDIF
     8728
     8729       CASE ( 'rad_sw_in_av' )
     8730          IF ( .NOT. ALLOCATED( rad_sw_in_av ) )  THEN
     8731             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
     8732                  radiation_scheme == 'constant')  THEN
     8733                ALLOCATE( rad_sw_in_av(0:0,nysg:nyng,nxlg:nxrg) )
     8734             ELSE
     8735                ALLOCATE( rad_sw_in_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     8736             ENDIF
     8737          ENDIF 
     8738          IF ( k == 1 )  THEN
     8739             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
     8740                  radiation_scheme == 'constant')  THEN
     8741                READ ( 13 )  tmp_3d2
     8742                rad_sw_in_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =&
     8743                    tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     8744             ELSE
     8745                READ ( 13 )  tmp_3d
     8746                rad_sw_in_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =  &
     8747                    tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     8748             ENDIF
     8749          ENDIF
     8750
     8751       CASE ( 'rad_sw_out' )
     8752          IF ( .NOT. ALLOCATED( rad_sw_out ) )  THEN
     8753             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
     8754                  radiation_scheme == 'constant')  THEN
     8755                ALLOCATE( rad_sw_out(0:0,nysg:nyng,nxlg:nxrg) )
     8756             ELSE
     8757                ALLOCATE( rad_sw_out(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     8758             ENDIF
     8759          ENDIF 
     8760          IF ( k == 1 )  THEN
     8761             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
     8762                  radiation_scheme == 'constant')  THEN
     8763                READ ( 13 )  tmp_3d2
     8764                rad_sw_out(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =  &
     8765                    tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     8766             ELSE
     8767                READ ( 13 )  tmp_3d
     8768                rad_sw_out(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =    &
     8769                    tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     8770             ENDIF
     8771          ENDIF
     8772
     8773       CASE ( 'rad_sw_out_av' )
     8774          IF ( .NOT. ALLOCATED( rad_sw_out_av ) )  THEN
     8775             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
     8776                  radiation_scheme == 'constant')  THEN
     8777                ALLOCATE( rad_sw_out_av(0:0,nysg:nyng,nxlg:nxrg) )
     8778             ELSE
     8779                ALLOCATE( rad_sw_out_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     8780             ENDIF
     8781          ENDIF 
     8782          IF ( k == 1 )  THEN
     8783             IF ( radiation_scheme == 'clear-sky'  .OR.                    &
     8784                  radiation_scheme == 'constant')  THEN
     8785                READ ( 13 )  tmp_3d2
     8786                rad_sw_out_av(0:0,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) &
     8787                   = tmp_3d2(0:0,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     8788             ELSE
     8789                READ ( 13 )  tmp_3d
     8790                rad_sw_out_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) = &
     8791                    tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     8792             ENDIF
     8793          ENDIF
     8794
     8795       CASE ( 'rad_sw_cs_hr' )
     8796          IF ( .NOT. ALLOCATED( rad_sw_cs_hr ) )  THEN
     8797             ALLOCATE( rad_sw_cs_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     8798          ENDIF
     8799          IF ( k == 1 )  READ ( 13 )  tmp_3d
     8800          rad_sw_cs_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =        &
     8801                  tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     8802
     8803       CASE ( 'rad_sw_cs_hr_av' )
     8804          IF ( .NOT. ALLOCATED( rad_sw_cs_hr_av ) )  THEN
     8805             ALLOCATE( rad_sw_cs_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     8806          ENDIF
     8807          IF ( k == 1 )  READ ( 13 )  tmp_3d
     8808          rad_sw_cs_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =     &
     8809                  tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     8810
     8811       CASE ( 'rad_sw_hr' )
     8812          IF ( .NOT. ALLOCATED( rad_sw_hr ) )  THEN
     8813             ALLOCATE( rad_sw_hr(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     8814          ENDIF
     8815          IF ( k == 1 )  READ ( 13 )  tmp_3d
     8816          rad_sw_hr(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =           &
     8817                  tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     8818
     8819       CASE ( 'rad_sw_hr_av' )
     8820          IF ( .NOT. ALLOCATED( rad_sw_hr_av ) )  THEN
     8821             ALLOCATE( rad_sw_hr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     8822          ENDIF
     8823          IF ( k == 1 )  READ ( 13 )  tmp_3d
     8824          rad_lw_hr_av(:,nysc-nbgp:nync+nbgp,nxlc-nbgp:nxrc+nbgp) =        &
     8825                  tmp_3d(:,nysf-nbgp:nynf+nbgp,nxlf-nbgp:nxrf+nbgp)
     8826
     8827       CASE DEFAULT
     8828
     8829          found = .FALSE.
     8830
     8831    END SELECT
    84558832
    84568833 END SUBROUTINE radiation_rrd_local
    84578834
     8835!------------------------------------------------------------------------------!
     8836! Description:
     8837! ------------
     8838!> Subroutine writes debug information
     8839!------------------------------------------------------------------------------!
     8840 SUBROUTINE radiation_write_debug_log ( message )
     8841    !> it writes debug log with time stamp
     8842    CHARACTER(*)  :: message
     8843    CHARACTER(15) :: dtc
     8844    CHARACTER(8)  :: date
     8845    CHARACTER(10) :: time
     8846    CHARACTER(5)  :: zone
     8847    CALL date_and_time(date, time, zone)
     8848    dtc = date(7:8)//','//time(1:2)//':'//time(3:4)//':'//time(5:10)
     8849    WRITE(9,'(2A)') dtc, TRIM(message)
     8850    FLUSH(9)
     8851 END SUBROUTINE radiation_write_debug_log
    84588852
    84598853 END MODULE radiation_model_mod
  • palm/trunk/SOURCE/surface_mod.f90

    r2894 r2920  
    2626! -----------------
    2727! $Id$
     28! Correct comment for surface directions
     29!
     30! 2894 2018-03-15 09:17:58Z Giersch
    2831! Calculations of the index range of the subdomain on file which overlaps with
    2932! the current subdomain are already done in read_restart_data_mod,
     
    420423
    421424    TYPE (surf_type), DIMENSION(0:2), TARGET ::  surf_def_h  !< horizontal default surfaces (Up, Down, and Top)
    422     TYPE (surf_type), DIMENSION(0:3), TARGET ::  surf_def_v  !< vertical default surfaces (North, South, West, East)
     425    TYPE (surf_type), DIMENSION(0:3), TARGET ::  surf_def_v  !< vertical default surfaces (North, South, East, West)
    423426    TYPE (surf_type)                , TARGET ::  surf_lsm_h  !< horizontal natural land surfaces, so far only upward-facing
    424     TYPE (surf_type), DIMENSION(0:3), TARGET ::  surf_lsm_v  !< vertical land surfaces (North, South, West, East)
     427    TYPE (surf_type), DIMENSION(0:3), TARGET ::  surf_lsm_v  !< vertical land surfaces (North, South, East, West)
    425428    TYPE (surf_type)                , TARGET ::  surf_usm_h  !< horizontal urban surfaces, so far only upward-facing
    426     TYPE (surf_type), DIMENSION(0:3), TARGET ::  surf_usm_v  !< vertical urban surfaces (North, South, West, East)
     429    TYPE (surf_type), DIMENSION(0:3), TARGET ::  surf_usm_v  !< vertical urban surfaces (North, South, East, West)
    427430
    428431    INTEGER(iwp) ::  ns_h_on_file(0:2)                       !< total number of horizontal surfaces with the same facing, required for writing restart data
  • palm/trunk/SOURCE/urban_surface_mod.f90

    r2906 r2920  
    1616!
    1717! Copyright 2015-2018 Czech Technical University in Prague
     18! Copyright 2015-2018 Institute of Computer Science of the
     19!                     Czech Academy of Sciences, Prague
    1820! Copyright 1997-2018 Leibniz Universitaet Hannover
    1921!------------------------------------------------------------------------------!
     
    2628! -----------------
    2729! $Id$
     30! Remove unused pcbl, npcbl from ONLY list
     31! moh.hefny:
     32! Fixed bugs introduced by new structures and by moving radiation interaction
     33! into radiation_model_mod.f90.
     34! Bugfix: usm data output 3D didn't respect directions
     35!
     36! 2906 2018-03-19 08:56:40Z Giersch
    2837! Local variable ids has to be initialized with a value of -1 in
    2938! usm_average_3d_data
     
    217226!> Further work:
    218227!> -------------
    219 !> 1. Reduce number of shape view factors by merging factors for distant surfaces
    220 !>    under shallow angles. Idea: Iteratively select the smallest shape view
    221 !>    factor by value (among all sources and targets) which has a similarly
    222 !>    oriented source neighbor (or near enough) SVF and merge them by adding
    223 !>    value of the smaller SVF to the larger one and deleting the smaller one.
    224 !>    This will allow for better scaling at higher resolutions.
    225 !>
    226 !> 2. Remove global arrays surfouts, surfoutl and only keep track of radiosity
     228!> 1. Remove global arrays surfouts, surfoutl and only keep track of radiosity
    227229!>    from surfaces that are visible from local surfaces (i.e. there is a SVF
    228230!>    where target is local). To do that, radiosity will be exchanged after each
    229231!>    reflection step using MPI_Alltoall instead of current MPI_Allgather.
    230232!>
    231 !> 3. Temporarily large values of surface heat flux can be observed, up to
     233!> 2. Temporarily large values of surface heat flux can be observed, up to
    232234!>    1.2 Km/s, which seem to be not realistic.
    233235!>
     
    248250#if ! defined( __nopointer )
    249251    USE arrays_3d,                                                             &
    250         ONLY:  zu, pt, pt_1, pt_2, p, u, v, w, hyp, tend
     252        ONLY:  hyp, zu, pt, pt_1, pt_2, p, u, v, w, hyp, tend
    251253#endif
    252254
     
    272274
    273275    USE date_and_time_mod,                                                     &
    274         ONLY:  d_seconds_year, day_of_year_init, time_utc_init
     276        ONLY:  time_utc_init
    275277
    276278    USE grid_variables,                                                        &
     
    288290   
    289291    USE plant_canopy_model_mod,                                                &
    290         ONLY:  pc_heating_rate, usm_lad_rma
     292        ONLY:  pc_heating_rate
    291293   
    292294    USE radiation_model_mod,                                                   &
    293         ONLY:  albedo_type, radiation, calc_zenith, zenith,                    &
     295        ONLY:  albedo_type, radiation_interaction, calc_zenith, zenith,                    &
    294296               rad_sw_in, rad_lw_in, rad_sw_out, rad_lw_out,                   &
    295297               sigma_sb, solar_constant, sun_direction, sun_dir_lat,           &
     
    297299               force_radiation_call, surfinsw, surfinlw, surfinswdir,          &
    298300               surfinswdif, surfoutsw, surfoutlw, surfins,nsvfl, svf, svfsurf, &
    299                surfinl, surfinlwdif, energy_balance_surf_h,                    &
    300                energy_balance_surf_v, rad_sw_in_dir, rad_sw_in_diff,           &
     301               surfinl, surfinlwdif, rad_sw_in_dir, rad_sw_in_diff,            &
    301302               rad_lw_in_diff, surfouts, surfoutl, surfoutsl, surfoutll, surf, &
    302                surfl, nsurfl, nsurfs, surfstart, pcbinsw, pcbinlw, pcbl, npcbl, startenergy,      &
    303                endenergy, nenergy, iup_u, inorth_u, isouth_u, ieast_u, iwest_u, iup_l,  &
    304                inorth_l, isouth_l, ieast_l, iwest_l, startsky, endsky,id,      &
    305                iz, iy, ix, idir, jdir, kdir, startborder, endborder, nsurf_type, nzub, nzut,   &
    306                isky, inorth_b,idown_a, isouth_b, ieast_b, iwest_b, nzu, pch, nsurf,  &
    307                iup_a, inorth_a, isouth_a, ieast_a, iwest_a, idsvf, ndsvf,      &
    308                idcsf, ndcsf, kdcsf, pct, startland, endland, startwall, endwall
     303               surfl, nsurfl, nsurfs, surfstart, pcbinsw, pcbinlw,             &
     304               iup_u, inorth_u, isouth_u, ieast_u, iwest_u, iup_l,             &
     305               inorth_l, isouth_l, ieast_l, iwest_l, id,                       &
     306               iz, iy, ix, idir, jdir, kdir,  nsurf_type, nzub, nzut,          &
     307               nzu, pch, nsurf, idsvf, ndsvf,                                  &
     308               iup_a, idown_a, inorth_a, isouth_a, ieast_a, iwest_a,           &
     309               idcsf, ndcsf, kdcsf, pct,                                       &
     310               startland, endland, startwall, endwall, skyvf, skyvft
    309311
    310312    USE statistics,                                                            &
     
    322324    LOGICAL                                        ::  force_radiation_call_l = .FALSE.   !< flag parameter for unscheduled radiation model calls
    323325    LOGICAL                                        ::  indoor_model = .FALSE.              !< whether to use the indoor model
    324 
    325    
     326    LOGICAL                                        ::  read_wall_temp_3d = .FALSE.
     327
     328
    326329    INTEGER(iwp)                                   ::  building_type = 1                  !< default building type (preleminary setting)
    327330    INTEGER(iwp)                                   ::  land_category = 2                  !< default category for land surface
    328331    INTEGER(iwp)                                   ::  wall_category = 2                  !< default category for wall surface over pedestrian zone
    329     INTEGER(iwp)                                   ::  pedestrant_category = 2            !< default category for wall surface in pedestrian zone
     332    INTEGER(iwp)                                   ::  pedestrian_category = 2            !< default category for wall surface in pedestrian zone
    330333    INTEGER(iwp)                                   ::  roof_category = 2                  !< default category for root surface
    331 
     334    REAL(wp)                                       ::  roughness_concrete = 0.001_wp      !< roughness length of average concrete surface
    332335!
    333336!-- Indices of input attributes for (above) ground floor level
     
    380383    REAL(wp)  ::  roof_height_limit = 4._wp          !< height for distinguish between land surfaces and roofs
    381384    REAL(wp)  ::  ground_floor_level = 4.0_wp        !< default ground floor level
    382     REAL(wp)  ::  ra_horiz_coef = 5.0_wp             !< mysterious coefficient for correction of overestimation
    383                                                                                           !< of r_a for horizontal surfaces -> TODO
     385
    384386
    385387    CHARACTER(37), DIMENSION(0:6), PARAMETER :: building_type_name = (/     &
     
    524526!-- anthropogenic heat sources
    525527!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    526     REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  aheat             !< daily average of anthropogenic heat (W/m2)
    527     REAL(wp), DIMENSION(:), ALLOCATABLE            ::  aheatprof         !< diurnal profile of anthropogenic heat
     528    REAL(wp), DIMENSION(:,:,:), ALLOCATABLE        ::  aheat             !< daily average of anthropogenic heat (W/m2)
     529    REAL(wp), DIMENSION(:,:), ALLOCATABLE          ::  aheatprof         !< diurnal profiles of anthropogenic heat for particular layers
     530    INTEGER(wp)                                    ::  naheatlayers = 1  !< number of layers of anthropogenic heat
    528531
    529532!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     
    547550!    REAL(wp), DIMENSION(nzb_wall:nzt_wall)         :: zwn_default_green = (/0.33_wp, 0.66_wp, 1.0_wp /)
    548551
    549                                                                        
    550     REAL(wp)                                       ::   wall_inner_temperature = 296.0_wp    !< temperature of the inner wall surface (~23 degrees C) (K)
    551     REAL(wp)                                       ::   roof_inner_temperature = 296.0_wp    !< temperature of the inner roof surface (~23 degrees C) (K)
    552     REAL(wp)                                       ::   soil_inner_temperature = 283.0_wp    !< temperature of the deep soil (~10 degrees C) (K)
    553     REAL(wp)                                       ::   window_inner_temperature = 296.0_wp  !< temperature of the inner window surface (~23 degrees C) (K)
     552
     553    REAL(wp)                                       :: wall_inner_temperature = 295.0_wp    !< temperature of the inner wall surface (~22 degrees C) (K)
     554    REAL(wp)                                       :: roof_inner_temperature = 295.0_wp    !< temperature of the inner roof surface (~22 degrees C) (K)
     555    REAL(wp)                                       :: soil_inner_temperature = 288.0_wp    !< temperature of the deep soil (~15 degrees C) (K)
     556    REAL(wp)                                       :: window_inner_temperature = 295.0_wp  !< temperature of the inner window surface (~22 degrees C) (K)
    554557
    555558!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
     
    674677!-- albedo, emissivity, lambda_surf, roughness, thickness, volumetric heat capacity, thermal conductivity
    675678    INTEGER(iwp)                                   :: n_surface_types      !< number of the wall type categories
    676     INTEGER(iwp), PARAMETER                        :: n_surface_params = 8 !< number of parameters for each type of the wall
     679    INTEGER(iwp), PARAMETER                        :: n_surface_params = 9 !< number of parameters for each type of the wall
    677680    INTEGER(iwp), PARAMETER                        :: ialbedo  = 1         !< albedo of the surface
    678681    INTEGER(iwp), PARAMETER                        :: iemiss   = 2         !< emissivity of the surface
    679     INTEGER(iwp), PARAMETER