Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

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

reintegrate branch resler to trunk

Location:
palm/trunk/SOURCE
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE

  • palm/trunk/SOURCE/radiation_model_mod.f90

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