Changeset 4653


Ignore:
Timestamp:
Aug 27, 2020 8:54:43 AM (4 years ago)
Author:
pavelkrc
Message:

Radiation Transfer Model version 4.0

Location:
palm/trunk/SOURCE
Files:
3 edited

Legend:

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

    r4648 r4653  
    24112411!-- Pre-calculate topography top indices (former get_topography_top_index
    24122412!-- function)
    2413     ALLOCATE( topo_top_ind(nysg:nyng,nxlg:nxrg,0:4) )
     2413    ALLOCATE( topo_top_ind(nysg:nyng,nxlg:nxrg,0:5) )
    24142414!
    24152415!-- Uppermost topography index on scalar grid
     
    24372437    topo_top_ind(:,:,4) = MAXLOC( MERGE( 1, 0, BTEST( wall_flags_total_0(:,:,:), ibit ) ), DIM=1 ) &
    24382438                          - 1
     2439!
     2440!-- Uppermost topography index including full-3D geometry
     2441    ibit = 12
     2442    DO  k = nzb, nzt+1
     2443       WHERE( BTEST( wall_flags_total_0(k,:,:), ibit ) )  topo_top_ind(:,:,5) = k
     2444    ENDDO
    24392445
    24402446 END SUBROUTINE set_topo_flags
  • palm/trunk/SOURCE/radiation_model_mod.f90

    r4644 r4653  
    2828! -----------------
    2929! $Id$
     30! RTM version 4.0:
     31! - Support full-3D geometry:
     32!   * Downward faces can only be default type (temporary PALM limitation)
     33!   * rad_angular_discretization = .FALSE. is not supported for full-3D cases
     34! - Order faces in surf(l) by column
     35! - Remove separate direction indices for urban and land faces
     36! - Bugfix: skyvft output
     37!
     38! 4644 2020-08-22 15:05:46Z pavelkrc
    3039! Bugfix: remove invalid reordering of vertical surfaces in RTM
    3140!
     
    426435
    427436    USE surface_mod,                                                           &
    428         ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win,                       &
     437        ONLY:  ind_pav_green, ind_veg_wall, ind_wat_win, surf_def_h,           &
    429438               surf_lsm_h, surf_lsm_v, surf_type, surf_usm_h, surf_usm_v,      &
    430439               vertical_surfaces_exist
     
    711720    INTEGER(iwp), PARAMETER                        ::  iy = 3                             !< position of j-index in surfl and surf
    712721    INTEGER(iwp), PARAMETER                        ::  ix = 4                             !< position of i-index in surfl and surf
    713     INTEGER(iwp), PARAMETER                        ::  im = 5                             !< position of surface m-index in surfl and surf
    714     INTEGER(iwp), PARAMETER                        ::  nidx_surf = 5                      !< number of indices in surfl and surf
    715 
    716     INTEGER(iwp), PARAMETER                        ::  nsurf_type = 10                    !< number of surf types incl. phys.(land+urban) & (atm.,sky,boundary) surfaces - 1
    717 
    718     INTEGER(iwp), PARAMETER                        ::  iup_u    = 0                       !< 0 - index of urban upward surface (ground or roof)
    719     INTEGER(iwp), PARAMETER                        ::  idown_u  = 1                       !< 1 - index of urban downward surface (overhanging)
    720     INTEGER(iwp), PARAMETER                        ::  inorth_u = 2                       !< 2 - index of urban northward facing wall
    721     INTEGER(iwp), PARAMETER                        ::  isouth_u = 3                       !< 3 - index of urban southward facing wall
    722     INTEGER(iwp), PARAMETER                        ::  ieast_u  = 4                       !< 4 - index of urban eastward facing wall
    723     INTEGER(iwp), PARAMETER                        ::  iwest_u  = 5                       !< 5 - index of urban westward facing wall
    724 
    725     INTEGER(iwp), PARAMETER                        ::  iup_l    = 6                       !< 6 - index of land upward surface (ground or roof)
    726     INTEGER(iwp), PARAMETER                        ::  inorth_l = 7                       !< 7 - index of land northward facing wall
    727     INTEGER(iwp), PARAMETER                        ::  isouth_l = 8                       !< 8 - index of land southward facing wall
    728     INTEGER(iwp), PARAMETER                        ::  ieast_l  = 9                       !< 9 - index of land eastward facing wall
    729     INTEGER(iwp), PARAMETER                        ::  iwest_l  = 10                      !< 10- index of land westward facing wall
    730 
    731     INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  idir = (/0, 0,0, 0,1,-1,0,0, 0,1,-1/)   !< surface normal direction x indices
    732     INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  jdir = (/0, 0,1,-1,0, 0,0,1,-1,0, 0/)   !< surface normal direction y indices
    733     INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  kdir = (/1,-1,0, 0,0, 0,1,0, 0,0, 0/)   !< surface normal direction z indices
     722    INTEGER(iwp), PARAMETER                        ::  nidx_surf = 4                      !< number of indices in surfl and surf
     723
     724    INTEGER(iwp), PARAMETER                        ::  nsurf_type = 5                     !< number of surf types = surface directions
     725
     726    INTEGER(iwp), PARAMETER                        ::  iup    = 0                         !< 0 - index of upward surface (ground or roof)
     727    INTEGER(iwp), PARAMETER                        ::  idown  = 1                         !< 1 - index of downward surface (overhanging)
     728    INTEGER(iwp), PARAMETER                        ::  inorth = 2                         !< 2 - index of northward facing wall
     729    INTEGER(iwp), PARAMETER                        ::  isouth = 3                         !< 3 - index of southward facing wall
     730    INTEGER(iwp), PARAMETER                        ::  ieast  = 4                         !< 4 - index of eastward facing wall
     731    INTEGER(iwp), PARAMETER                        ::  iwest  = 5                         !< 5 - index of westward facing wall
     732
     733    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  idir = (/0, 0, 0, 0, 1,-1/)      !< surface normal direction x indices
     734    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  jdir = (/0, 0, 1,-1, 0, 0/)      !< surface normal direction y indices
     735    INTEGER(iwp), DIMENSION(0:nsurf_type), PARAMETER ::  kdir = (/1,-1, 0, 0, 0, 0/)      !< surface normal direction z indices
    734736    REAL(wp),     DIMENSION(0:nsurf_type)          :: facearea                            !< area of single face in respective
    735737                                                                                          !< direction (will be calc'd)
    736738
    737 
    738 !-- indices and sizes of urban and land surface models
    739     INTEGER(iwp)                                   ::  startland        !< start index of block of land and roof surfaces
    740     INTEGER(iwp)                                   ::  endland          !< end index of block of land and roof surfaces
    741     INTEGER(iwp)                                   ::  nlands           !< number of land and roof surfaces in local processor
    742     INTEGER(iwp)                                   ::  startwall        !< start index of block of wall surfaces
    743     INTEGER(iwp)                                   ::  endwall          !< end index of block of wall surfaces
    744     INTEGER(iwp)                                   ::  nwalls           !< number of wall surfaces in local processor
    745739
    746740!-- indices needed for RTM netcdf output subroutines
    747741    INTEGER(iwp), PARAMETER                        :: nd = 5
    748742    CHARACTER(LEN=6), DIMENSION(0:nd-1), PARAMETER :: dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
    749     INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     :: dirint_u = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /)
    750     INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     :: dirint_l = (/ iup_l, isouth_l, inorth_l, iwest_l, ieast_l /)
    751     INTEGER(iwp), DIMENSION(0:nd-1)                :: dirstart
    752     INTEGER(iwp), DIMENSION(0:nd-1)                :: dirend
     743    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER     :: dirint = (/ iup, isouth, inorth, iwest, ieast /)
    753744
    754745!-- indices and sizes of urban and land surface models
     
    762753    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfstart        !< starts of blocks of surfaces for individual processors in array surf (indexed from 1)
    763754                                                                        !< respective block for particular processor is surfstart[iproc+1]+1 : surfstart[iproc+1]+nsurfs[iproc+1]
     755    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfl_col_start  !< start of surfaces in surfl for each x,y column (local surfaces)
     756    INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET ::  surfg_col_start  !< start of surfaces in surfl for each x,y column (all surfaces)
    764757
    765758!-- block variables needed for calculation of the plant canopy model inside the urban surface model
     
    786779    REAL(wp), PARAMETER                            ::  ext_coef = 0.6_wp                  !< extinction coefficient (a.k.a. alpha)
    787780    INTEGER(iwp), PARAMETER                        ::  rad_version_len = 10               !< length of identification string of rad version
    788     CHARACTER(rad_version_len), PARAMETER          ::  rad_version = 'RAD v. 3.0'         !< identification of version of binary svf and restart files
     781    CHARACTER(rad_version_len), PARAMETER          ::  rad_version = 'RAD v. 4.0'         !< identification of version of binary svf and restart files
    789782    INTEGER(iwp)                                   ::  raytrace_discrete_elevs = 40       !< number of discretization steps for elevation (nadir to zenith)
    790783    INTEGER(iwp)                                   ::  raytrace_discrete_azims = 80       !< number of discretization steps for azimuth (out of 360 degrees)
     
    879872#endif
    880873    REAL(wp)                                       ::  prototype_lad    !< prototype leaf area density for computing effective optical depth
    881     INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  nzterr, plantt   !< temporary global arrays for raytracing
     874    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  nzterrt          !< temporary global arrays for raytracing
     875    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  nzterrb          !< temporary global arrays for raytracing
     876    INTEGER(iwp), DIMENSION(:), ALLOCATABLE        ::  plantt           !< temporary global arrays for raytracing
    882877    INTEGER(iwp)                                   ::  plantt_max
    883878
     
    10931088           cos_zenith, calc_zenith, sun_direction, sun_dir_lat, sun_dir_lon,   &
    10941089           idir, jdir, kdir, id, iz, iy, ix,                                   &
    1095            iup_u, inorth_u, isouth_u, ieast_u, iwest_u,                        &
    1096            iup_l, inorth_l, isouth_l, ieast_l, iwest_l,                        &
    1097            nsurf_type, nz_urban_b, nz_urban_t, nz_urban, pch, nsurf,                 &
     1090           iup, inorth, isouth, ieast, iwest,                                  &
     1091           nsurf_type, nz_urban_b, nz_urban_t, nz_urban, pch, nsurf,           &
    10981092           idsvf, ndsvf, idcsf, ndcsf, kdcsf, pct,                             &
    1099            radiation_interactions, startwall, startland, endland, endwall,     &
    1100            skyvf, skyvft, radiation_interactions_on, average_radiation,        &
    1101            rad_sw_in_diff, rad_sw_in_dir
     1093           radiation_interactions, skyvf, skyvft, radiation_interactions_on,   &
     1094           average_radiation, rad_sw_in_diff, rad_sw_in_dir
    11021095
    11031096
     
    59665959!
    59675960!--  Set up thermal radiation from surfaces
    5968 !--  Horizontal walls
    59695961     mm = 1
    59705962     DO  i = nxl, nxr
    59715963        DO  j = nys, nyn
    59725964!
     5965!--        Horizontal walls
    59735966!--        urban
    59745967           DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
     
    59965989              mm = mm + 1
    59975990           ENDDO
    5998         ENDDO
    5999      ENDDO
    6000 !
    6001 !--  Vertical walls
    6002      DO  i = nxl, nxr
    6003         DO  j = nys, nyn
     5991!
     5992!--        Downward facing default surfaces
     5993!--        TODO: Downward-facing surfaces are only available as default surfaces which lack
     5994!--        temperature and heat balance, therefore they have to be simulated as neutral to
     5995!--        radiation, that means having sw/lw reflectivity=1, so they neither absorb nor emit
     5996!--        radiation and they have zero net heat flux.
     5997           DO  m = surf_def_h(1)%start_index(j,i), surf_def_h(1)%end_index(j,i)
     5998              surfoutll(mm) = 0.0_wp
     5999              albedo_surf(mm) = 1.0_wp
     6000              emiss_surf(mm)  = 0.0_wp
     6001              mm = mm + 1
     6002           ENDDO
     6003!
     6004!--        Vertical walls
    60046005           DO  l = 0, 3
    60056006!
     
    64006401!--  Transfer radiation arrays required for energy balance to the respective data types
    64016402!    and claculate relevant radiation model-RTM coupling terms
     6403     mm = 1
     6404     DO  i = nxl, nxr
     6405        DO  j = nys, nyn
     6406!
     6407!--        Horizontal walls
     6408!--        urban
     6409           DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
     6410              surf_usm_h%rad_sw_in(m)  = surfinsw(mm)
     6411              surf_usm_h%rad_sw_out(m) = surfoutsw(mm)
     6412              surf_usm_h%rad_sw_dir(m) = surfinswdir(mm)
     6413              surf_usm_h%rad_sw_dif(m) = surfinswdif(mm)
     6414              surf_usm_h%rad_sw_ref(m) = surfinsw(mm) - surfinswdir(mm) -        &
     6415                                         surfinswdif(mm)
     6416              surf_usm_h%rad_sw_res(m) = surfins(mm)
     6417              surf_usm_h%rad_lw_in(m)  = surfinlw(mm)
     6418              surf_usm_h%rad_lw_out(m) = surfoutlw(mm)
     6419              surf_usm_h%rad_net(m)    = surfinsw(mm) - surfoutsw(mm) +          &
     6420                                         surfinlw(mm) - surfoutlw(mm)
     6421              surf_usm_h%rad_net_l(m)  = surf_usm_h%rad_net(m)
     6422              surf_usm_h%rad_lw_dif(m) = surfinlwdif(mm)
     6423              surf_usm_h%rad_lw_ref(m) = surfinlw(mm) - surfinlwdif(mm)
     6424              surf_usm_h%rad_lw_res(m) = surfinl(mm)
     6425              mm = mm + 1
     6426           ENDDO
     6427!
     6428!--        land
     6429           DO  m = surf_lsm_h%start_index(j,i), surf_lsm_h%end_index(j,i)
     6430              surf_lsm_h%rad_sw_in(m)  = surfinsw(mm)
     6431              surf_lsm_h%rad_sw_out(m) = surfoutsw(mm)
     6432              surf_lsm_h%rad_sw_dir(m) = surfinswdir(mm)
     6433              surf_lsm_h%rad_sw_dif(m) = surfinswdif(mm)
     6434              surf_lsm_h%rad_sw_ref(m) = surfinsw(mm) - surfinswdir(mm) -        &
     6435                                            surfinswdif(mm)
     6436              surf_lsm_h%rad_sw_res(m) = surfins(mm)
     6437              surf_lsm_h%rad_lw_in(m)  = surfinlw(mm)
     6438              surf_lsm_h%rad_lw_out(m) = surfoutlw(mm)
     6439              surf_lsm_h%rad_net(m)    = surfinsw(mm) - surfoutsw(mm) +          &
     6440                                         surfinlw(mm) - surfoutlw(mm)
     6441              surf_lsm_h%rad_lw_dif(m) = surfinlwdif(mm)
     6442              surf_lsm_h%rad_lw_ref(m) = surfinlw(mm) - surfinlwdif(mm)
     6443              surf_lsm_h%rad_lw_res(m) = surfinl(mm)
     6444              mm = mm + 1
     6445           ENDDO
     6446!
     6447!--        Downward facing default surfaces
     6448!--        TODO: Downward-facing surfaces are only available as default surfaces which lack
     6449!--        temperature and heat balance, therefore they have to be simulated as neutral to
     6450!--        radiation, that means having sw/lw reflectivity=1, so they neither absorb nor emit
     6451!--        radiation and they have zero net heat flux.
     6452           DO  m = surf_def_h(1)%start_index(j,i), surf_def_h(1)%end_index(j,i)
     6453              mm = mm + 1
     6454           ENDDO
     6455!
     6456!--        Vertical walls
     6457           DO  l = 0, 3
     6458!
     6459!--           urban
     6460              DO  m = surf_usm_v(l)%start_index(j,i),                       &
     6461                      surf_usm_v(l)%end_index(j,i)
     6462                 surf_usm_v(l)%rad_sw_in(m)  = surfinsw(mm)
     6463                 surf_usm_v(l)%rad_sw_out(m) = surfoutsw(mm)
     6464                 surf_usm_v(l)%rad_sw_dir(m) = surfinswdir(mm)
     6465                 surf_usm_v(l)%rad_sw_dif(m) = surfinswdif(mm)
     6466                 surf_usm_v(l)%rad_sw_ref(m) = surfinsw(mm) - surfinswdir(mm) -     &
     6467                                               surfinswdif(mm)
     6468                 surf_usm_v(l)%rad_sw_res(m) = surfins(mm)
     6469                 surf_usm_v(l)%rad_lw_in(m)  = surfinlw(mm)
     6470                 surf_usm_v(l)%rad_lw_out(m) = surfoutlw(mm)
     6471                 surf_usm_v(l)%rad_net(m)    = surfinsw(mm) - surfoutsw(mm) +       &
     6472                                               surfinlw(mm) - surfoutlw(mm)
     6473                 surf_usm_v(l)%rad_net_l(m)  = surf_usm_v(l)%rad_net(m)
     6474                 surf_usm_v(l)%rad_lw_dif(m) = surfinlwdif(mm)
     6475                 surf_usm_v(l)%rad_lw_ref(m) = surfinlw(mm) - surfinlwdif(mm)
     6476                 surf_usm_v(l)%rad_lw_res(m) = surfinl(mm)
     6477                 mm = mm + 1
     6478              ENDDO
     6479!
     6480!--           land
     6481              DO  m = surf_lsm_v(l)%start_index(j,i),                       &
     6482                      surf_lsm_v(l)%end_index(j,i)
     6483                 surf_lsm_v(l)%rad_sw_in(m)  = surfinsw(mm)
     6484                 surf_lsm_v(l)%rad_sw_out(m) = surfoutsw(mm)
     6485                 surf_lsm_v(l)%rad_sw_dir(m) = surfinswdir(mm)
     6486                 surf_lsm_v(l)%rad_sw_dif(m) = surfinswdif(mm)
     6487                 surf_lsm_v(l)%rad_sw_ref(m) = surfinsw(mm) - surfinswdir(mm) -     &
     6488                                               surfinswdif(mm)
     6489                 surf_lsm_v(l)%rad_sw_res(m) = surfins(mm)
     6490                 surf_lsm_v(l)%rad_lw_in(m)  = surfinlw(mm)
     6491                 surf_lsm_v(l)%rad_lw_out(m) = surfoutlw(mm)
     6492                 surf_lsm_v(l)%rad_net(m)    = surfinsw(mm) - surfoutsw(mm) +       &
     6493                                               surfinlw(mm) - surfoutlw(mm)
     6494                 surf_lsm_v(l)%rad_lw_dif(m) = surfinlwdif(mm)
     6495                 surf_lsm_v(l)%rad_lw_ref(m) = surfinlw(mm) - surfinlwdif(mm)
     6496                 surf_lsm_v(l)%rad_lw_res(m) = surfinl(mm)
     6497                 mm = mm + 1
     6498              ENDDO
     6499           ENDDO
     6500        ENDDO
     6501     ENDDO
    64026502
    64036503     DO  i = 1, nsurfl
    6404         m  = surfl(im,i)
    64056504        d  = surfl(id, i)
    6406 !
    6407 !--     (1) Urban surfaces
    6408 !--     upward-facing
    6409         IF ( surfl(1,i) == iup_u )  THEN
    6410            surf_usm_h%rad_sw_in(m)  = surfinsw(i)
    6411            surf_usm_h%rad_sw_out(m) = surfoutsw(i)
    6412            surf_usm_h%rad_sw_dir(m) = surfinswdir(i)
    6413            surf_usm_h%rad_sw_dif(m) = surfinswdif(i)
    6414            surf_usm_h%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) -        &
    6415                                       surfinswdif(i)
    6416            surf_usm_h%rad_sw_res(m) = surfins(i)
    6417            surf_usm_h%rad_lw_in(m)  = surfinlw(i)
    6418            surf_usm_h%rad_lw_out(m) = surfoutlw(i)
    6419            surf_usm_h%rad_net(m)    = surfinsw(i) - surfoutsw(i) +          &
    6420                                       surfinlw(i) - surfoutlw(i)
    6421            surf_usm_h%rad_net_l(m)  = surf_usm_h%rad_net(m)
    6422            surf_usm_h%rad_lw_dif(m) = surfinlwdif(i)
    6423            surf_usm_h%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)
    6424            surf_usm_h%rad_lw_res(m) = surfinl(i)
    6425 !
    6426 !--     northward-facding
    6427         ELSEIF ( surfl(1,i) == inorth_u )  THEN
    6428            surf_usm_v(0)%rad_sw_in(m)  = surfinsw(i)
    6429            surf_usm_v(0)%rad_sw_out(m) = surfoutsw(i)
    6430            surf_usm_v(0)%rad_sw_dir(m) = surfinswdir(i)
    6431            surf_usm_v(0)%rad_sw_dif(m) = surfinswdif(i)
    6432            surf_usm_v(0)%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) -     &
    6433                                          surfinswdif(i)
    6434            surf_usm_v(0)%rad_sw_res(m) = surfins(i)
    6435            surf_usm_v(0)%rad_lw_in(m)  = surfinlw(i)
    6436            surf_usm_v(0)%rad_lw_out(m) = surfoutlw(i)
    6437            surf_usm_v(0)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
    6438                                          surfinlw(i) - surfoutlw(i)
    6439            surf_usm_v(0)%rad_net_l(m)  = surf_usm_v(0)%rad_net(m)
    6440            surf_usm_v(0)%rad_lw_dif(m) = surfinlwdif(i)
    6441            surf_usm_v(0)%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)
    6442            surf_usm_v(0)%rad_lw_res(m) = surfinl(i)
    6443 !
    6444 !--     southward-facding
    6445         ELSEIF ( surfl(1,i) == isouth_u )  THEN
    6446            surf_usm_v(1)%rad_sw_in(m)  = surfinsw(i)
    6447            surf_usm_v(1)%rad_sw_out(m) = surfoutsw(i)
    6448            surf_usm_v(1)%rad_sw_dir(m) = surfinswdir(i)
    6449            surf_usm_v(1)%rad_sw_dif(m) = surfinswdif(i)
    6450            surf_usm_v(1)%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) -     &
    6451                                          surfinswdif(i)
    6452            surf_usm_v(1)%rad_sw_res(m) = surfins(i)
    6453            surf_usm_v(1)%rad_lw_in(m)  = surfinlw(i)
    6454            surf_usm_v(1)%rad_lw_out(m) = surfoutlw(i)
    6455            surf_usm_v(1)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
    6456                                          surfinlw(i) - surfoutlw(i)
    6457            surf_usm_v(1)%rad_net_l(m)  = surf_usm_v(1)%rad_net(m)
    6458            surf_usm_v(1)%rad_lw_dif(m) = surfinlwdif(i)
    6459            surf_usm_v(1)%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)
    6460            surf_usm_v(1)%rad_lw_res(m) = surfinl(i)
    6461 !
    6462 !--     eastward-facing
    6463         ELSEIF ( surfl(1,i) == ieast_u )  THEN
    6464            surf_usm_v(2)%rad_sw_in(m)  = surfinsw(i)
    6465            surf_usm_v(2)%rad_sw_out(m) = surfoutsw(i)
    6466            surf_usm_v(2)%rad_sw_dir(m) = surfinswdir(i)
    6467            surf_usm_v(2)%rad_sw_dif(m) = surfinswdif(i)
    6468            surf_usm_v(2)%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) -     &
    6469                                          surfinswdif(i)
    6470            surf_usm_v(2)%rad_sw_res(m) = surfins(i)
    6471            surf_usm_v(2)%rad_lw_in(m)  = surfinlw(i)
    6472            surf_usm_v(2)%rad_lw_out(m) = surfoutlw(i)
    6473            surf_usm_v(2)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
    6474                                          surfinlw(i) - surfoutlw(i)
    6475            surf_usm_v(2)%rad_net_l(m)  = surf_usm_v(2)%rad_net(m)
    6476            surf_usm_v(2)%rad_lw_dif(m) = surfinlwdif(i)
    6477            surf_usm_v(2)%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)
    6478            surf_usm_v(2)%rad_lw_res(m) = surfinl(i)
    6479 !
    6480 !--     westward-facding
    6481         ELSEIF ( surfl(1,i) == iwest_u )  THEN
    6482            surf_usm_v(3)%rad_sw_in(m)  = surfinsw(i)
    6483            surf_usm_v(3)%rad_sw_out(m) = surfoutsw(i)
    6484            surf_usm_v(3)%rad_sw_dir(m) = surfinswdir(i)
    6485            surf_usm_v(3)%rad_sw_dif(m) = surfinswdif(i)
    6486            surf_usm_v(3)%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) -     &
    6487                                          surfinswdif(i)
    6488            surf_usm_v(3)%rad_sw_res(m) = surfins(i)
    6489            surf_usm_v(3)%rad_lw_in(m)  = surfinlw(i)
    6490            surf_usm_v(3)%rad_lw_out(m) = surfoutlw(i)
    6491            surf_usm_v(3)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
    6492                                          surfinlw(i) - surfoutlw(i)
    6493            surf_usm_v(3)%rad_net_l(m)  = surf_usm_v(3)%rad_net(m)
    6494            surf_usm_v(3)%rad_lw_dif(m) = surfinlwdif(i)
    6495            surf_usm_v(3)%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)
    6496            surf_usm_v(3)%rad_lw_res(m) = surfinl(i)
    6497 !
    6498 !--     (2) land surfaces
    6499 !--     upward-facing
    6500         ELSEIF ( surfl(1,i) == iup_l )  THEN
    6501            surf_lsm_h%rad_sw_in(m)  = surfinsw(i)
    6502            surf_lsm_h%rad_sw_out(m) = surfoutsw(i)
    6503            surf_lsm_h%rad_sw_dir(m) = surfinswdir(i)
    6504            surf_lsm_h%rad_sw_dif(m) = surfinswdif(i)
    6505            surf_lsm_h%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) -        &
    6506                                          surfinswdif(i)
    6507            surf_lsm_h%rad_sw_res(m) = surfins(i)
    6508            surf_lsm_h%rad_lw_in(m)  = surfinlw(i)
    6509            surf_lsm_h%rad_lw_out(m) = surfoutlw(i)
    6510            surf_lsm_h%rad_net(m)    = surfinsw(i) - surfoutsw(i) +          &
    6511                                       surfinlw(i) - surfoutlw(i)
    6512            surf_lsm_h%rad_lw_dif(m) = surfinlwdif(i)
    6513            surf_lsm_h%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)
    6514            surf_lsm_h%rad_lw_res(m) = surfinl(i)
    6515 !
    6516 !--     northward-facding
    6517         ELSEIF ( surfl(1,i) == inorth_l )  THEN
    6518            surf_lsm_v(0)%rad_sw_in(m)  = surfinsw(i)
    6519            surf_lsm_v(0)%rad_sw_out(m) = surfoutsw(i)
    6520            surf_lsm_v(0)%rad_sw_dir(m) = surfinswdir(i)
    6521            surf_lsm_v(0)%rad_sw_dif(m) = surfinswdif(i)
    6522            surf_lsm_v(0)%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) -     &
    6523                                          surfinswdif(i)
    6524            surf_lsm_v(0)%rad_sw_res(m) = surfins(i)
    6525            surf_lsm_v(0)%rad_lw_in(m)  = surfinlw(i)
    6526            surf_lsm_v(0)%rad_lw_out(m) = surfoutlw(i)
    6527            surf_lsm_v(0)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
    6528                                          surfinlw(i) - surfoutlw(i)
    6529            surf_lsm_v(0)%rad_lw_dif(m) = surfinlwdif(i)
    6530            surf_lsm_v(0)%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)
    6531            surf_lsm_v(0)%rad_lw_res(m) = surfinl(i)
    6532 !
    6533 !--     southward-facding
    6534         ELSEIF ( surfl(1,i) == isouth_l )  THEN
    6535            surf_lsm_v(1)%rad_sw_in(m)  = surfinsw(i)
    6536            surf_lsm_v(1)%rad_sw_out(m) = surfoutsw(i)
    6537            surf_lsm_v(1)%rad_sw_dir(m) = surfinswdir(i)
    6538            surf_lsm_v(1)%rad_sw_dif(m) = surfinswdif(i)
    6539            surf_lsm_v(1)%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) -     &
    6540                                          surfinswdif(i)
    6541            surf_lsm_v(1)%rad_sw_res(m) = surfins(i)
    6542            surf_lsm_v(1)%rad_lw_in(m)  = surfinlw(i)
    6543            surf_lsm_v(1)%rad_lw_out(m) = surfoutlw(i)
    6544            surf_lsm_v(1)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
    6545                                          surfinlw(i) - surfoutlw(i)
    6546            surf_lsm_v(1)%rad_lw_dif(m) = surfinlwdif(i)
    6547            surf_lsm_v(1)%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)
    6548            surf_lsm_v(1)%rad_lw_res(m) = surfinl(i)
    6549 !
    6550 !--     eastward-facing
    6551         ELSEIF ( surfl(1,i) == ieast_l )  THEN
    6552            surf_lsm_v(2)%rad_sw_in(m)  = surfinsw(i)
    6553            surf_lsm_v(2)%rad_sw_out(m) = surfoutsw(i)
    6554            surf_lsm_v(2)%rad_sw_dir(m) = surfinswdir(i)
    6555            surf_lsm_v(2)%rad_sw_dif(m) = surfinswdif(i)
    6556            surf_lsm_v(2)%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) -     &
    6557                                          surfinswdif(i)
    6558            surf_lsm_v(2)%rad_sw_res(m) = surfins(i)
    6559            surf_lsm_v(2)%rad_lw_in(m)  = surfinlw(i)
    6560            surf_lsm_v(2)%rad_lw_out(m) = surfoutlw(i)
    6561            surf_lsm_v(2)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
    6562                                          surfinlw(i) - surfoutlw(i)
    6563            surf_lsm_v(2)%rad_lw_dif(m) = surfinlwdif(i)
    6564            surf_lsm_v(2)%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)
    6565            surf_lsm_v(2)%rad_lw_res(m) = surfinl(i)
    6566 !
    6567 !--     westward-facing
    6568         ELSEIF ( surfl(1,i) == iwest_l )  THEN
    6569            surf_lsm_v(3)%rad_sw_in(m)  = surfinsw(i)
    6570            surf_lsm_v(3)%rad_sw_out(m) = surfoutsw(i)
    6571            surf_lsm_v(3)%rad_sw_dir(m) = surfinswdir(i)
    6572            surf_lsm_v(3)%rad_sw_dif(m) = surfinswdif(i)
    6573            surf_lsm_v(3)%rad_sw_ref(m) = surfinsw(i) - surfinswdir(i) -     &
    6574                                          surfinswdif(i)
    6575            surf_lsm_v(3)%rad_sw_res(m) = surfins(i)
    6576            surf_lsm_v(3)%rad_lw_in(m)  = surfinlw(i)
    6577            surf_lsm_v(3)%rad_lw_out(m) = surfoutlw(i)
    6578            surf_lsm_v(3)%rad_net(m)    = surfinsw(i) - surfoutsw(i) +       &
    6579                                          surfinlw(i) - surfoutlw(i)
    6580            surf_lsm_v(3)%rad_lw_dif(m) = surfinlwdif(i)
    6581            surf_lsm_v(3)%rad_lw_ref(m) = surfinlw(i) - surfinlwdif(i)
    6582            surf_lsm_v(3)%rad_lw_res(m) = surfinl(i)
    6583         ENDIF
    65846505!
    65856506!--     RTM coupling terms
     
    69426863       INTEGER(iwp)        ::  k_topo               !< vertical index indicating
    69436864                                                    !< topography top for given (j,i)
     6865       INTEGER(iwp)        ::  icol                 !< flat column number (in (y,x) fortran order)
    69446866       INTEGER(iwp)        ::  isurf, ipcgb, imrt
    69456867       INTEGER(iwp)        ::  nzptl, nzubl, nzutl
     
    69656887!--    for any upward-facing wall (see bit 12).
    69666888       nzubl = MINVAL( topo_top_ind(nys:nyn,nxl:nxr,0) )
    6967        nzutl = MAXVAL( topo_top_ind(nys:nyn,nxl:nxr,0) )
     6889       nzutl = MAXVAL( topo_top_ind(nys:nyn,nxl:nxr,5) )
    69686890
    69696891       nzubl = MAX( nzubl, nzb )
     
    70776999!--    Number of horizontal surfaces including land- and roof surfaces in both USM and LSM. Note that
    70787000!--    All horizontal surface elements are already counted in surface_mod.
    7079        startland = 1
    7080        nsurfl    = surf_usm_h%ns + surf_lsm_h%ns
    7081        endland   = nsurfl
    7082        nlands    = endland - startland + 1
    7083 
     7001       nsurfl    = surf_usm_h%ns + surf_lsm_h%ns + surf_def_h(1)%ns
    70847002!
    70857003!--    Number of vertical surfaces in both USM and LSM. Note that all vertical surface elements are
    70867004!--    already counted in surface_mod.
    7087        startwall = nsurfl+1
    70887005       DO  i = 0,3
    70897006          nsurfl = nsurfl + surf_usm_v(i)%ns + surf_lsm_v(i)%ns
    70907007       ENDDO
    7091        endwall = nsurfl
    7092        nwalls  = endwall - startwall + 1
    7093        dirstart = (/ startland, startwall, startwall, startwall, startwall /)
    7094        dirend = (/ endland, endwall, endwall, endwall, endwall /)
    70957008
    70967009!--    fill gridpcbl and pcbl
     
    71967109       ENDIF
    71977110
    7198 !--    add horizontal surface elements (land and urban surfaces)
    7199 !--    TODO: add urban overhanging surfaces (idown_u)
     7111       ALLOCATE( surfl_col_start(0:nnx*nny-1), surfg_col_start(0:(nx+1)*(ny+1)) )
     7112!
     7113!--    Add horizontal and vertical surface elements (land and urban surfaces)
     7114!--    ordered by x,y column (y most varying)
     7115!--    TODO: remove the hard coding of l = 0 to l = idirection
     7116       icol = 0
    72007117       DO i = nxl, nxr
    72017118           DO j = nys, nyn
     7119!
     7120!--           Save column start
     7121              surfl_col_start(icol) = isurf + 1
     7122              icol = icol + 1
     7123!
     7124!--           Horizontal surfaces
    72027125              DO  m = surf_usm_h%start_index(j,i), surf_usm_h%end_index(j,i)
    72037126                 k = surf_usm_h%k(m)
    72047127                 isurf = isurf + 1
    7205                  surfl(:,isurf) = (/iup_u,k,j,i,m/)
    7206                  IF ( rad_angular_discretization ) THEN
    7207                     gridsurf(iup_u,k,j,i) = isurf
    7208                  ENDIF
     7128                 surfl(:,isurf) = (/iup,k,j,i/)
    72097129              ENDDO
    72107130
     
    72127132                 k = surf_lsm_h%k(m)
    72137133                 isurf = isurf + 1
    7214                  surfl(:,isurf) = (/iup_l,k,j,i,m/)
    7215                  IF ( rad_angular_discretization ) THEN
    7216                     gridsurf(iup_u,k,j,i) = isurf
    7217                  ENDIF
     7134                 surfl(:,isurf) = (/iup,k,j,i/)
    72187135              ENDDO
    72197136
    7220            ENDDO
    7221        ENDDO
    7222 
    7223 !--    add vertical surface elements (land and urban surfaces)
    7224 !--    TODO: remove the hard coding of l = 0 to l = idirection
    7225        DO i = nxl, nxr
    7226            DO j = nys, nyn
     7137              DO  m = surf_def_h(1)%start_index(j,i), surf_def_h(1)%end_index(j,i)
     7138                 k = surf_def_h(1)%k(m)
     7139                 isurf = isurf + 1
     7140                 surfl(:,isurf) = (/idown,k,j,i/)
     7141              ENDDO
     7142!
     7143!--           Vertical surfaces
    72277144              l = 0
    72287145              DO  m = surf_usm_v(l)%start_index(j,i), surf_usm_v(l)%end_index(j,i)
    72297146                 k = surf_usm_v(l)%k(m)
    72307147                 isurf = isurf + 1
    7231                  surfl(:,isurf) = (/inorth_u,k,j,i,m/)
    7232                  IF ( rad_angular_discretization ) THEN
    7233                     gridsurf(inorth_u,k,j,i) = isurf
    7234                  ENDIF
     7148                 surfl(:,isurf) = (/inorth,k,j,i/)
    72357149              ENDDO
    72367150              DO  m = surf_lsm_v(l)%start_index(j,i), surf_lsm_v(l)%end_index(j,i)
    72377151                 k = surf_lsm_v(l)%k(m)
    72387152                 isurf = isurf + 1
    7239                  surfl(:,isurf) = (/inorth_l,k,j,i,m/)
    7240                  IF ( rad_angular_discretization ) THEN
    7241                     gridsurf(inorth_u,k,j,i) = isurf
    7242                  ENDIF
     7153                 surfl(:,isurf) = (/inorth,k,j,i/)
    72437154              ENDDO
    72447155
     
    72477158                 k = surf_usm_v(l)%k(m)
    72487159                 isurf = isurf + 1
    7249                  surfl(:,isurf) = (/isouth_u,k,j,i,m/)
    7250                  IF ( rad_angular_discretization ) THEN
    7251                     gridsurf(isouth_u,k,j,i) = isurf
    7252                  ENDIF
     7160                 surfl(:,isurf) = (/isouth,k,j,i/)
    72537161              ENDDO
    72547162              DO  m = surf_lsm_v(l)%start_index(j,i), surf_lsm_v(l)%end_index(j,i)
    72557163                 k = surf_lsm_v(l)%k(m)
    72567164                 isurf = isurf + 1
    7257                  surfl(:,isurf) = (/isouth_l,k,j,i,m/)
    7258                  IF ( rad_angular_discretization ) THEN
    7259                     gridsurf(isouth_u,k,j,i) = isurf
    7260                  ENDIF
     7165                 surfl(:,isurf) = (/isouth,k,j,i/)
    72617166              ENDDO
    72627167
     
    72657170                 k = surf_usm_v(l)%k(m)
    72667171                 isurf = isurf + 1
    7267                  surfl(:,isurf) = (/ieast_u,k,j,i,m/)
    7268                  IF ( rad_angular_discretization ) THEN
    7269                     gridsurf(ieast_u,k,j,i) = isurf
    7270                  ENDIF
     7172                 surfl(:,isurf) = (/ieast,k,j,i/)
    72717173              ENDDO
    72727174              DO  m = surf_lsm_v(l)%start_index(j,i), surf_lsm_v(l)%end_index(j,i)
    72737175                 k = surf_lsm_v(l)%k(m)
    72747176                 isurf = isurf + 1
    7275                  surfl(:,isurf) = (/ieast_l,k,j,i,m/)
    7276                  IF ( rad_angular_discretization ) THEN
    7277                     gridsurf(ieast_u,k,j,i) = isurf
    7278                  ENDIF
     7177                 surfl(:,isurf) = (/ieast,k,j,i/)
    72797178              ENDDO
    72807179
     
    72837182                 k = surf_usm_v(l)%k(m)
    72847183                 isurf = isurf + 1
    7285                  surfl(:,isurf) = (/iwest_u,k,j,i,m/)
    7286                  IF ( rad_angular_discretization ) THEN
    7287                     gridsurf(iwest_u,k,j,i) = isurf
    7288                  ENDIF
     7184                 surfl(:,isurf) = (/iwest,k,j,i/)
    72897185              ENDDO
    72907186              DO  m = surf_lsm_v(l)%start_index(j,i), surf_lsm_v(l)%end_index(j,i)
    72917187                 k = surf_lsm_v(l)%k(m)
    72927188                 isurf = isurf + 1
    7293                  surfl(:,isurf) = (/iwest_l,k,j,i,m/)
    7294                  IF ( rad_angular_discretization ) THEN
    7295                     gridsurf(iwest_u,k,j,i) = isurf
    7296                  ENDIF
     7189                 surfl(:,isurf) = (/iwest,k,j,i/)
    72977190              ENDDO
    72987191           ENDDO
     
    74067299       surf = surfl
    74077300#endif
    7408 
    7409 !--
     7301!
     7302!--    Populate gridsurf with reverse global indices (->surf)
     7303       IF ( rad_angular_discretization )  THEN
     7304          DO  isurf = 1, nsurfl
     7305             gridsurf(surfl(id,isurf),surfl(iz,isurf), &
     7306                      surfl(iy,isurf),surfl(ix,isurf)) = isurf + surfstart(myid)
     7307          ENDDO
     7308       ENDIF
     7309!
     7310!--    Gather column start indices
     7311#if defined( __parallel )
     7312       CALL MPI_Allgather(surfl_col_start, nnx*nny, MPI_INTEGER,             &
     7313                          surfg_col_start, nnx*nny, MPI_INTEGER, comm2d, ierr)
     7314       IF ( ierr /= 0 ) THEN
     7315          WRITE(9,*) 'Error MPI_AllGather1b:', ierr, surfl_col_start, surfg_col_start
     7316          FLUSH(9)
     7317       ENDIF
     7318!
     7319!--    Convert local indices (->surfl) to global (->surf)
     7320       DO  i = 0, numprocs-1
     7321          surfg_col_start(i*nnx*nny:(i+1)*nnx*nny-1) =                      &
     7322                surfg_col_start(i*nnx*nny:(i+1)*nnx*nny-1) + surfstart(i)
     7323       ENDDO
     7324#else
     7325       surfg_col_start(0:(nx+1)*(ny+1)-1) = surfl_col_start(0:(nx+1)*(ny+1)-1)
     7326#endif
     7327       surfg_col_start((nx+1)*(ny+1)) = nsurf+1
     7328!
    74107329!--    allocation of the arrays for direct and diffusion radiation
    74117330       IF ( debug_output )  CALL debug_message( 'allocation of radiation arrays', 'info' )
     
    74837402        REAL(wp), DIMENSION(:), ALLOCATABLE           :: vffrac0       !< dtto (original values)
    74847403        REAL(wp), DIMENSION(:), ALLOCATABLE           :: ztransp       !< array of transparency in z steps
    7485         INTEGER(iwp)                                  :: lowest_free_ray !< index into zdirs
    74867404        INTEGER(iwp), DIMENSION(:), ALLOCATABLE       :: itarget       !< face indices of detected obstacles
    74877405        INTEGER(iwp)                                  :: itarg0, itarg1
     
    75027420        INTEGER(iwp)                                  :: max_track_len !< maximum 2d track length
    75037421#if defined( __parallel )
    7504         INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET:: nzterrl_l
    7505         INTEGER(iwp), DIMENSION(:,:), POINTER         :: nzterrl
     7422        INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET:: nzterrtl_l
     7423        INTEGER(iwp), DIMENSION(:), ALLOCATABLE,TARGET:: nzterrbl_l
     7424        INTEGER(iwp), DIMENSION(:,:), POINTER         :: nzterrtl
     7425        INTEGER(iwp), DIMENSION(:,:), POINTER         :: nzterrbl
    75067426        INTEGER(iwp)                                  :: minfo
    75077427        REAL(wp), DIMENSION(:), POINTER, SAVE         :: lad_s_rma       !< fortran 1D pointer
     
    75407460
    75417461!--     initialize temporary terrain and plant canopy height arrays (global 2D array!)
    7542         ALLOCATE( nzterr(0:(nx+1)*(ny+1)-1) )
     7462        ALLOCATE( nzterrt(0:(nx+1)*(ny+1)-1) )
     7463        ALLOCATE( nzterrb(0:(nx+1)*(ny+1)-1) )
    75437464#if defined( __parallel )
    7544         !ALLOCATE( nzterrl(nys:nyn,nxl:nxr) )
    7545         ALLOCATE( nzterrl_l((nyn-nys+1)*(nxr-nxl+1)) )
    7546         nzterrl(nys:nyn,nxl:nxr) => nzterrl_l(1:(nyn-nys+1)*(nxr-nxl+1))
    7547         nzterrl = topo_top_ind(nys:nyn,nxl:nxr,0)
    7548         CALL MPI_AllGather( nzterrl_l, nnx*nny, MPI_INTEGER, &
    7549                             nzterr, nnx*nny, MPI_INTEGER, comm2d, ierr )
     7465        ALLOCATE( nzterrtl_l((nyn-nys+1)*(nxr-nxl+1)) )
     7466        ALLOCATE( nzterrbl_l((nyn-nys+1)*(nxr-nxl+1)) )
     7467        nzterrtl(nys:nyn,nxl:nxr) => nzterrtl_l(1:(nyn-nys+1)*(nxr-nxl+1))
     7468        nzterrbl(nys:nyn,nxl:nxr) => nzterrbl_l(1:(nyn-nys+1)*(nxr-nxl+1))
     7469        nzterrtl = topo_top_ind(nys:nyn,nxl:nxr,5)
     7470        nzterrbl = topo_top_ind(nys:nyn,nxl:nxr,0)
     7471        CALL MPI_AllGather( nzterrtl_l, nnx*nny, MPI_INTEGER, &
     7472                            nzterrt, nnx*nny, MPI_INTEGER, comm2d, ierr )
    75507473        IF ( ierr /= 0 ) THEN
    7551             WRITE(9,*) 'Error MPI_AllGather1:', ierr, SIZE(nzterrl_l), nnx*nny, &
    7552                        SIZE(nzterr), nnx*nny
     7474            WRITE(9,*) 'Error MPI_AllGather1t:', ierr, SIZE(nzterrtl_l), nnx*nny, &
     7475                       SIZE(nzterrt), nnx*nny
    75537476            FLUSH(9)
    75547477        ENDIF
    7555         DEALLOCATE(nzterrl_l)
     7478        CALL MPI_AllGather( nzterrbl_l, nnx*nny, MPI_INTEGER, &
     7479                            nzterrb, nnx*nny, MPI_INTEGER, comm2d, ierr )
     7480        IF ( ierr /= 0 ) THEN
     7481            WRITE(9,*) 'Error MPI_AllGather1b:', ierr, SIZE(nzterrbl_l), nnx*nny, &
     7482                       SIZE(nzterrb), nnx*nny
     7483            FLUSH(9)
     7484        ENDIF
     7485        DEALLOCATE(nzterrtl_l)
     7486        DEALLOCATE(nzterrbl_l)
    75567487#else
    7557         nzterr = RESHAPE( topo_top_ind(nys:nyn,nxl:nxr,0), (/(nx+1)*(ny+1)/) )
     7488        nzterrt = RESHAPE( topo_top_ind(nys:nyn,nxl:nxr,5), (/(nx+1)*(ny+1)/) )
     7489        nzterrb = RESHAPE( topo_top_ind(nys:nyn,nxl:nxr,0), (/(nx+1)*(ny+1)/) )
    75587490#endif
    75597491        IF ( plant_canopy )  THEN
     
    77017633            !--Select a proper half-sphere for 2D raytracing
    77027634            SELECT CASE ( td )
    7703                CASE ( iup_u, iup_l )
     7635               CASE ( iup )
    77047636                  az0 = 0._wp
    77057637                  naz = raytrace_discrete_azims
     
    77087640                  nzn = raytrace_discrete_elevs / 2
    77097641                  zns = pi / 2._wp / REAL(nzn, wp)
    7710                CASE ( isouth_u, isouth_l )
     7642               CASE ( idown )
     7643                  az0 = 0._wp
     7644                  naz = raytrace_discrete_azims
     7645                  azs = 2._wp * pi / REAL(naz, wp)
     7646                  zn0 = pi / 2._wp
     7647                  nzn = raytrace_discrete_elevs / 2
     7648                  zns = pi / 2._wp / REAL(nzn, wp)
     7649               CASE ( isouth )
    77117650                  az0 = pi / 2._wp
    77127651                  naz = raytrace_discrete_azims / 2
     
    77157654                  nzn = raytrace_discrete_elevs
    77167655                  zns = pi / REAL(nzn, wp)
    7717                CASE ( inorth_u, inorth_l )
     7656               CASE ( inorth )
    77187657                  az0 = - pi / 2._wp
    77197658                  naz = raytrace_discrete_azims / 2
     
    77227661                  nzn = raytrace_discrete_elevs
    77237662                  zns = pi / REAL(nzn, wp)
    7724                CASE ( iwest_u, iwest_l )
     7663               CASE ( iwest )
    77257664                  az0 = pi
    77267665                  naz = raytrace_discrete_azims / 2
     
    77297668                  nzn = raytrace_discrete_elevs
    77307669                  zns = pi / REAL(nzn, wp)
    7731                CASE ( ieast_u, ieast_l )
     7670               CASE ( ieast )
    77327671                  az0 = 0._wp
    77337672                  naz = raytrace_discrete_azims / 2
     
    77447683
    77457684            ALLOCATE ( zdirs(1:nzn), zcent(1:nzn), zbdry(0:nzn), vffrac(1:nzn*naz), &
    7746                        ztransp(1:nzn*naz), itarget(1:nzn*naz) )   !FIXME allocate itarget only
    7747                                                                   !in case of rad_angular_discretization
     7685                       ztransp(1:nzn*naz), itarget(1:nzn*naz) )
    77487686
    77497687            itarg0 = 1
     
    77517689            zcent(:) = (/( zn0+(REAL(izn,wp)-.5_wp)*zns, izn=1, nzn )/)
    77527690            zbdry(:) = (/( zn0+REAL(izn,wp)*zns, izn=0, nzn )/)
    7753             IF ( td == iup_u  .OR.  td == iup_l )  THEN
     7691            IF ( td == iup )  THEN
    77547692               vffrac(1:nzn) = (COS(2 * zbdry(0:nzn-1)) - COS(2 * zbdry(1:nzn))) / 2._wp / REAL(naz, wp)
    77557693!
     
    77597697               ENDDO
    77607698!--            sum of whole vffrac equals 1, verified
     7699            ELSEIF ( td == idown )  THEN
     7700               vffrac(1:nzn) = -(COS(2 * zbdry(0:nzn-1)) - COS(2 * zbdry(1:nzn))) / 2._wp / REAL(naz, wp)
     7701               DO iaz = 1, naz-1
     7702                  vffrac(iaz*nzn+1:(iaz+1)*nzn) = vffrac(1:nzn)
     7703               ENDDO
     7704!--            sum of whole vffrac equals 1, verified
    77617705            ENDIF
    77627706!
     
    77647708            DO iaz = 1, naz
    77657709               azmid = az0 + (REAL(iaz, wp) - .5_wp) * azs
    7766                IF ( td /= iup_u  .AND.  td /= iup_l )  THEN
     7710               IF ( td /= iup  .AND.  td /= idown )  THEN
    77677711                  az2 = REAL(iaz, wp) * azs - pi/2._wp
    77687712                  az1 = az2 - azs
     
    77837727                                    surfstart(myid) + isurflt, facearea(td),  &
    77847728                                    vffrac(itarg0:itarg1), .TRUE., .TRUE.,    &
    7785                                     .FALSE., lowest_free_ray,                 &
    7786                                     ztransp(itarg0:itarg1),                   &
     7729                                    .FALSE., ztransp(itarg0:itarg1),          &
    77877730                                    itarget(itarg0:itarg1))
    77887731
    7789                skyvf(isurflt) = skyvf(isurflt) + &
    7790                                 SUM(vffrac(itarg0:itarg0+lowest_free_ray-1))
    7791                skyvft(isurflt) = skyvft(isurflt) + &
    7792                                  SUM(ztransp(itarg0:itarg0+lowest_free_ray-1) &
    7793                                              * vffrac(itarg0:itarg0+lowest_free_ray-1))
    7794 
     7732               skyvf(isurflt) = skyvf(isurflt) +                    &
     7733                                SUM(vffrac(itarg0:itarg1),          &
     7734                                    MASK=(itarget(itarg0:itarg1)<0))
     7735               skyvft(isurflt) = skyvft(isurflt) +                                   &
     7736                                 SUM(ztransp(itarg0:itarg1) * vffrac(itarg0:itarg1), &
     7737                                     MASK=(itarget(itarg0:itarg1)<0))
     7738!
    77957739!--            Save direct solar transparency
    77967740               j = MODULO(NINT(azmid/                                          &
    77977741                               (2._wp*pi)*raytrace_discrete_azims-.5_wp, iwp), &
    77987742                          raytrace_discrete_azims)
    7799 
    7800                DO k = 1, raytrace_discrete_elevs/2
    7801                   i = dsidir_rev(k-1, j)
    7802                   IF ( i /= -1  .AND.  k <= lowest_free_ray )  &
    7803                      dsitrans(isurflt, i) = ztransp(itarg0+k-1)
    7804                ENDDO
    7805 
     7743!
     7744!--            For down direction there is no direct irradiance, otherwise az0=0
     7745               IF ( td /= idown )  THEN
     7746                  DO  k = 1, raytrace_discrete_elevs/2
     7747                     i = dsidir_rev(k-1, j)
     7748                     IF ( i /= -1  .AND.  itarget(itarg0+k-1) < 0 )  &
     7749                        dsitrans(isurflt, i) = ztransp(itarg0+k-1)
     7750                  ENDDO
     7751               ENDIF
    78067752!
    78077753!--            Advance itarget indices
     
    78207766               itarg0 = 1
    78217767               DO WHILE ( itarg0 <= nzn*naz )
    7822                   IF ( itarget(itarg0) /= -1 )  EXIT
     7768                  IF ( itarget(itarg0) >= 0 )  EXIT
    78237769                  itarg0 = itarg0 + 1
    78247770               ENDDO
     
    79837929              CALL raytrace_2d(ta, yxdir, nzn, zdirs,                                &
    79847930                                   -999, -999._wp, vffrac, .FALSE., .FALSE., .TRUE., &
    7985                                    lowest_free_ray, ztransp, itarget)
     7931                                   ztransp, itarget)
    79867932
    79877933!--           Save direct solar transparency
     
    79917937              DO  k = 1, raytrace_discrete_elevs/2
    79927938                 i = dsidir_rev(k-1, j)
    7993                  IF ( i /= -1  .AND.  k <= lowest_free_ray ) &
     7939                 IF ( i /= -1  .AND.  itarget(k) < 0 ) &
    79947940                    dsitransc(ipcgb, i) = ztransp(k)
    79957941              ENDDO
     
    80107956           zns = pi / REAL(nzn, wp)
    80117957           ALLOCATE ( zdirs(1:nzn), zcent(1:nzn), zbdry(0:nzn), vffrac(1:nzn*naz), vffrac0(1:nzn), &
    8012                       ztransp(1:nzn*naz), itarget(1:nzn*naz) )   !FIXME allocate itarget only
    8013                                                                  !in case of rad_angular_discretization
     7958                      ztransp(1:nzn*naz), itarget(1:nzn*naz) )
    80147959
    80157960           zcent(:) = (/( zn0+(REAL(izn,wp)-.5_wp)*zns, izn=1, nzn )/)
     
    80598004                 CALL raytrace_2d(ta, yxdir, nzn, zdirs,                         &
    80608005                                  -999, -999._wp, vffrac(itarg0:itarg1), .TRUE., &
    8061                                   .FALSE., .TRUE., lowest_free_ray,              &
    8062                                   ztransp(itarg0:itarg1),                        &
     8006                                  .FALSE., .TRUE., ztransp(itarg0:itarg1),       &
    80638007                                  itarget(itarg0:itarg1))
    80648008
    80658009!--              Sky view factors for MRT
    80668010                 mrtsky(imrt) = mrtsky(imrt) + &
    8067                                   SUM(vffrac(itarg0:itarg0+lowest_free_ray-1))
     8011                                  SUM(vffrac(itarg0:itarg1),          &
     8012                                      MASK=(itarget(itarg0:itarg1)<0))
    80688013                 mrtskyt(imrt) = mrtskyt(imrt) + &
    8069                                    SUM(ztransp(itarg0:itarg0+lowest_free_ray-1) &
    8070                                                * vffrac(itarg0:itarg0+lowest_free_ray-1))
     8014                                   SUM(ztransp(itarg0:itarg1) * vffrac(itarg0:itarg1), &
     8015                                       MASK=(itarget(itarg0:itarg1)<0))
    80718016!--              Direct solar transparency for MRT
    80728017                 j = MODULO(NINT(azmid/                                         &
     
    80758020                 DO  k = 1, raytrace_discrete_elevs/2
    80768021                    i = dsidir_rev(k-1, j)
    8077                     IF ( i /= -1  .AND.  k <= lowest_free_ray ) &
     8022                    IF ( i /= -1  .AND.  itarget(itarg0+k-1) < 0 ) &
    80788023                       mrtdsit(imrt, i) = ztransp(itarg0+k-1)
    80798024                 ENDDO
     
    80928037              itarg0 = 1
    80938038              DO WHILE ( itarg0 <= nzn*naz )
    8094                  IF ( itarget(itarg0) /= -1 )  EXIT
     8039                 IF ( itarget(itarg0) >= 0 )  EXIT
    80958040                 itarg0 = itarg0 + 1
    80968041              ENDDO
     
    81838128
    81848129!--     deallocate temporary global arrays
    8185         DEALLOCATE(nzterr)
     8130        DEALLOCATE( nzterrt, nzterrb )
    81868131
    81878132        IF ( plant_canopy )  THEN
     
    86078552!--             in the array nzterr and plantt and id of the coresponding processor
    86088553                CALL radiation_calc_global_offset( box(3), box(2), 0, 1, offs_glob=ig )
    8609                 IF ( box(1) <= nzterr(ig) )  THEN
     8554                IF ( box(1) <= nzterrb(ig) )  THEN
    86108555                    visible = .FALSE.
    86118556                    RETURN
     
    87038648   SUBROUTINE raytrace_2d(origin, yxdir, nrays, zdirs, iorig, aorig, vffrac,  &
    87048649                              calc_svf, create_csf, skip_1st_pcb,             &
    8705                               lowest_free_ray, transparency, itarget)
     8650                              transparency, itarget)
    87068651      IMPLICIT NONE
    87078652
     
    87168661      LOGICAL, INTENT(in)                    ::  create_csf    !< whether to create canopy sink factors
    87178662      LOGICAL, INTENT(in)                    ::  skip_1st_pcb  !< whether to skip first plant canopy box during raytracing
    8718       INTEGER(iwp), INTENT(out)              ::  lowest_free_ray !< index into zdirs
    87198663      REAL(wp), DIMENSION(nrays), INTENT(OUT) ::  transparency !< transparencies of zdirs paths
    8720       INTEGER(iwp), DIMENSION(nrays), INTENT(OUT) ::  itarget  !< global indices of target faces for zdirs
    8721 
    8722       INTEGER(iwp), DIMENSION(nrays)         ::  target_procs
    8723       REAL(wp)                               ::  horizon       !< highest horizon found after raytracing (z/hdist)
     8664      INTEGER(iwp), DIMENSION(nrays), INTENT(OUT) ::  itarget  !< global indices of target faces for zdirs or <0 for sky
     8665
     8666      LOGICAL                                ::  is_mixed_col  !< whether current column contains full-3D geometry
     8667      INTEGER(iwp)                           ::  lowest_free_ray  !< index into zdirs
     8668      INTEGER(iwp)                           ::  lowest_mixed_ray !< index into zdirs
     8669      REAL(wp)                               ::  full_horizon  !< highest full horizon found after raytracing (z/hdist)
    87248670      INTEGER(iwp)                           ::  i, k, l, d
     8671      INTEGER(iwp)                           ::  iray         !< index into zdirs
     8672      INTEGER(iwp)                           ::  isurf        !< index into surf(l)
    87258673      INTEGER(iwp)                           ::  seldim       !< dimension to be incremented
    87268674      REAL(wp), DIMENSION(2)                 ::  yxorigin     !< horizontal copy of origin (y,x)
     
    87298677      REAL(wp)                               ::  nextdist     !< end of current crossing
    87308678      REAL(wp)                               ::  crmid        !< midpoint of crossing
    8731       REAL(wp)                               ::  horz_entry   !< horizon at entry to column
    8732       REAL(wp)                               ::  horz_exit    !< horizon at exit from column
     8679      REAL(wp)                               ::  horz_full_entry  !< full horizon at entry to column
     8680      REAL(wp)                               ::  horz_full_exit   !< full horizon at exit from column
     8681      REAL(wp)                               ::  horz_mixed_entry !< mixed horizon at entry to column
     8682      REAL(wp)                               ::  horz_mixed_exit  !< mixed horizon at exit from column
    87338683      REAL(wp)                               ::  bdydim       !< boundary for current dimension
    87348684      REAL(wp), DIMENSION(2)                 ::  crossdist    !< distances to boundary for dimensions
     
    87388688      INTEGER(iwp), DIMENSION(2)             ::  dimdelta     !< dimension direction = +- 1
    87398689      INTEGER(iwp)                           ::  ip           !< number of processor where gridbox reside
    8740       INTEGER(iwp)                           ::  ig           !< 1D index of gridbox in global 2D array
     8690      INTEGER(iwp)                           ::  ip_last      !< previous number of processor where gridbox reside
     8691      INTEGER(iwp)                           ::  ig           !< 1D index of grid column in global 2D array
     8692      INTEGER(iwp)                           ::  ig_last      !< 1D index of previous column in global 2D array
    87418693      INTEGER(iwp)                           ::  maxboxes     !< max no of CSF created
    87428694      INTEGER(iwp)                           ::  nly          !< maximum  plant canopy height
     
    87468698      INTEGER(iwp)                           ::  zb1
    87478699      INTEGER(iwp)                           ::  nz
    8748       INTEGER(iwp)                           ::  iz
     8700      INTEGER(iwp)                           ::  kz
    87498701      INTEGER(iwp)                           ::  zsgn
    87508702      INTEGER(iwp)                           ::  lastdir      !< wall direction before hitting this column
     
    87558707      INTEGER(iwp)                           ::  wcount       !< RMA window item count
    87568708      INTEGER(MPI_ADDRESS_KIND)              ::  wdisp        !< RMA window displacement
     8709      INTEGER(iwp), DIMENSION(:), ALLOCATABLE::  vert_col     !< vertical faces from previous to current column
     8710      INTEGER(iwp), DIMENSION(:), ALLOCATABLE::  up_col       !< upward oriented surfaces in current column
     8711      INTEGER(iwp), DIMENSION(:), ALLOCATABLE::  down_col     !< downward oriented surfaces in current column
    87578712#endif
    87588713
     
    87698724      yxorigin(:) = origin(2:3)
    87708725      transparency(:) = 1._wp !-- Pre-set the all rays to transparent before reducing
    8771       horizon = -HUGE(1._wp)
     8726      full_horizon = -HUGE(1._wp)
     8727      lowest_mixed_ray = nrays
    87728728      lowest_free_ray = nrays
    8773       IF ( rad_angular_discretization  .AND.  calc_svf )  THEN
     8729      IF ( rad_angular_discretization )  THEN
    87748730         ALLOCATE(target_surfl(nrays))
    87758731         target_surfl(:) = -1
    87768732         lastdir = -999
    87778733         lastcolumn(:) = -999
     8734         ALLOCATE( vert_col(nz_urban_b:nz_urban_t) )
     8735         ALLOCATE( up_col(nz_urban_b:nz_urban_t) )
     8736         ALLOCATE( down_col(nz_urban_b:nz_urban_t) )
    87788737      ENDIF
    87798738
     
    88068765      ENDIF
    88078766
     8767      ip_last = -1
     8768      ig_last = -1
    88088769      lastdist = 0._wp
    88098770
     
    88398800
    88408801!--         calculate index of the grid with global indices (column(1),column(2))
    8841 !--         in the array nzterr and plantt and id of the coresponding processor
    8842             CALL radiation_calc_global_offset( column(2), column(1), 0, 1, offs_glob=ig )
    8843 
    8844             IF ( lastdist == 0._wp )  THEN
    8845                horz_entry = -HUGE(1._wp)
     8802!--         in the array nzterrt/b and plantt and id of the coresponding processor
     8803            CALL radiation_calc_global_offset( column(2), column(1), 0, 1, &
     8804                                               iproc=ip, offs_glob=ig )
     8805
     8806            IF ( ip_last < 0 )  THEN
     8807               horz_full_entry = -HUGE(1._wp)
     8808               horz_mixed_entry = -HUGE(1._wp)
    88468809            ELSE
    8847                horz_entry = (REAL(nzterr(ig), wp) + .5_wp - origin(1)) / lastdist
     8810               horz_full_entry = (REAL(nzterrb(ig), wp) + .5_wp - origin(1)) / lastdist
     8811               horz_mixed_entry = (REAL(nzterrt(ig), wp) + .5_wp - origin(1)) / lastdist
    88488812            ENDIF
    8849             horz_exit = (REAL(nzterr(ig), wp) + .5_wp - origin(1)) / nextdist
    8850 
    8851             IF ( rad_angular_discretization  .AND.  calc_svf )  THEN
    8852 !
    8853 !--            Identify vertical obstacles hit by rays in current column
    8854                DO WHILE ( lowest_free_ray > 0 )
    8855                   IF ( zdirs(lowest_free_ray) > horz_entry )  EXIT
     8813            horz_full_exit = (REAL(nzterrb(ig), wp) + .5_wp - origin(1)) / nextdist
     8814            horz_mixed_exit = (REAL(nzterrt(ig), wp) + .5_wp - origin(1)) / nextdist
     8815            is_mixed_col = ( nzterrt(ig) /= nzterrb(ig) )
     8816
     8817            IF ( rad_angular_discretization )  THEN
     8818!
     8819!--            Identify vertical full obstacles hit by rays in current column,
     8820!--            mixed rays need to be checked if they are already obstructed.
     8821               DO WHILE ( lowest_mixed_ray > lowest_free_ray )
     8822                  IF ( zdirs(lowest_mixed_ray) > horz_full_entry )  EXIT
    88568823!
    88578824!--               This may only happen after 1st column, so lastdir and lastcolumn are valid
    8858                   CALL request_itarget(lastdir,                                         &
    8859                         CEILING(-0.5_wp + origin(1) + zdirs(lowest_free_ray)*lastdist), &
    8860                         lastcolumn(1), lastcolumn(2),                                   &
    8861                         target_surfl(lowest_free_ray), target_procs(lowest_free_ray))
    8862                   lowest_free_ray = lowest_free_ray - 1
     8825                  IF ( target_surfl(lowest_mixed_ray) < 0 )  THEN
     8826                     CALL request_itarget(lastdir,                                                 &
     8827                           CEILING(-0.5_wp + origin(1) + zdirs(lowest_mixed_ray)*lastdist),        &
     8828                           lastcolumn(1), lastcolumn(2),                                           &
     8829                           target_surfl(lowest_mixed_ray))
     8830                  ENDIF
     8831                  lowest_mixed_ray = lowest_mixed_ray - 1
    88638832               ENDDO
    88648833!
    8865 !--            Identify horizontal obstacles hit by rays in current column
    8866                DO WHILE ( lowest_free_ray > 0 )
    8867                   IF ( zdirs(lowest_free_ray) > horz_exit )  EXIT
    8868                   CALL request_itarget(iup_u, nzterr(ig)+1, column(1), column(2), &
    8869                                        target_surfl(lowest_free_ray),           &
    8870                                        target_procs(lowest_free_ray))
    8871                   lowest_free_ray = lowest_free_ray - 1
     8834!--            Identify vertical full obstacles hit by rays in current column,
     8835!--            free rays need no individual checks.
     8836               DO WHILE ( lowest_mixed_ray > 0 )
     8837                  IF ( zdirs(lowest_mixed_ray) > horz_full_entry )  EXIT
     8838!
     8839!--               This may only happen after 1st column, so lastdir and lastcolumn are valid
     8840                  CALL request_itarget(lastdir,                                                    &
     8841                        CEILING(-0.5_wp + origin(1) + zdirs(lowest_mixed_ray)*lastdist),           &
     8842                        lastcolumn(1), lastcolumn(2),                                              &
     8843                        target_surfl(lowest_mixed_ray))
     8844                  lowest_mixed_ray = lowest_mixed_ray - 1
    88728845               ENDDO
    8873             ENDIF
    8874 
    8875             horizon = MAX(horizon, horz_entry, horz_exit)
     8846               IF ( lowest_free_ray > lowest_mixed_ray )  lowest_free_ray = lowest_mixed_ray
     8847!
     8848!--            Identify targets for vertical mixed obstacles.
     8849!--            lowest_mixed_ray now points to bottom of vertical mixed obstacles.
     8850               IF ( is_mixed_col  .AND.  ip_last >= 0 )  THEN
     8851!
     8852!--               Load vertical surfaces belonging to previous column
     8853                  vert_col(:) = -999
     8854                  DO  isurf = surfg_col_start(ig_last), surfg_col_start(ig_last+1)-1
     8855                     IF ( surf(id, isurf) == lastdir )  THEN
     8856                        vert_col(surf(iz, isurf)) = isurf
     8857                     ENDIF
     8858                  ENDDO
     8859!
     8860!--               Previously mixed rays need to be checked whether they are obstructed
     8861                  DO  iray = lowest_mixed_ray, lowest_free_ray-1, -1
     8862                     IF ( zdirs(iray) > horz_mixed_entry )  EXIT
     8863                     IF ( target_surfl(iray) >= 0 )  CYCLE
     8864                     target_surfl(iray) = vert_col(CEILING(-0.5_wp +     &
     8865                                        origin(1) + zdirs(iray)*lastdist)) ! contains -999 if missing surface
     8866                  ENDDO
     8867!
     8868!--               Previously free rays cannot be obstructed yet
     8869                  DO  iray = lowest_free_ray, 1, -1
     8870                     IF ( zdirs(iray) > horz_mixed_entry )  THEN
     8871!
     8872!--                     Extend mixed rays by raising the lowest_free ray
     8873!--                     (remains unchanged if this is hit in 1st iteration)
     8874                        lowest_free_ray = iray
     8875                        EXIT
     8876                     ENDIF
     8877                     target_surfl(iray) = vert_col(CEILING(-0.5_wp +     &
     8878                                        origin(1) + zdirs(iray)*lastdist)) ! contains -999 if missing surface
     8879                  ENDDO
     8880               ENDIF ! end of mixed horizon
     8881!
     8882!--            Identify horizontal full obstacles hit by rays in current column,
     8883!--            mixed rays need to be checked if they are already obstructed.
     8884               DO WHILE ( lowest_mixed_ray > lowest_free_ray )
     8885                  IF ( zdirs(lowest_mixed_ray) > horz_full_exit )  EXIT
     8886                  IF ( target_surfl(lowest_mixed_ray) < 0 )  THEN
     8887                     CALL request_itarget(iup, nzterrb(ig)+1, column(1), column(2), &
     8888                                          target_surfl(lowest_mixed_ray))
     8889                  ENDIF
     8890                  lowest_mixed_ray = lowest_mixed_ray - 1
     8891               ENDDO
     8892!
     8893!--            Identify horizontal full obstacles hit by rays in current column,
     8894!--            free rays need no individual checks.
     8895               DO WHILE ( lowest_mixed_ray > 0 )
     8896                  IF ( zdirs(lowest_mixed_ray) > horz_full_exit )  EXIT
     8897                  CALL request_itarget(iup, nzterrb(ig)+1, column(1), column(2), &
     8898                                       target_surfl(lowest_mixed_ray))
     8899                  lowest_mixed_ray = lowest_mixed_ray - 1
     8900               ENDDO
     8901               IF ( lowest_free_ray > lowest_mixed_ray )  lowest_free_ray = lowest_mixed_ray
     8902!
     8903!--            Identify targets for horizontal mixed obstacles.
     8904!--            lowest_mixed_ray now points _above_ horizontal full obstacles.
     8905               IF ( is_mixed_col )  THEN
     8906!
     8907!--               Load horizontal surfaces corresponding to current column
     8908                  up_col(:) = -999
     8909                  down_col(:) = -999
     8910                  DO  isurf = surfg_col_start(ig), surfg_col_start(ig+1)-1
     8911                     SELECT CASE ( surf(id, isurf) )
     8912                     CASE ( iup )
     8913                        up_col(surf(iz, isurf)) = isurf
     8914                     CASE ( idown )
     8915                        down_col(surf(iz, isurf)) = isurf
     8916                     ENDSELECT
     8917                  ENDDO
     8918!
     8919!--               Previously mixed rays need to be checked whether they are obstructed
     8920                  DO  iray = lowest_mixed_ray, lowest_free_ray-1, -1
     8921                     IF ( zdirs(iray) > MAX(horz_mixed_entry, horz_mixed_exit) )  EXIT
     8922                     IF ( target_surfl(iray) >= 0 )  CYCLE
     8923                     IF ( zdirs(iray) <= 0._wp )  THEN
     8924!
     8925!--                     Downward pointed ray, cycle k down from entry to exit,
     8926!--                     search for upward oriented faces
     8927                        DO  k = FLOOR(0.5_wp + origin(1) + zdirs(iray)*lastdist), &
     8928                               CEILING(-0.5_wp + origin(1) + zdirs(iray)*nextdist)+1, -1
     8929                           target_surfl(iray) = up_col(k) ! contains -999 if missing surface
     8930                           IF ( target_surfl(iray) >= 0 )  EXIT
     8931                        ENDDO
     8932                     ELSE
     8933!
     8934!--                     Upward pointed ray, cycle k up from entry to exit,
     8935!--                     search for downward oriented faces
     8936                        DO  k = CEILING(-0.5_wp + origin(1) + zdirs(iray)*lastdist), &
     8937                               FLOOR(0.5_wp + origin(1) + zdirs(iray)*nextdist)-1
     8938                           target_surfl(iray) = down_col(k) ! contains -999 if missing surface
     8939                           IF ( target_surfl(iray) >= 0 )  EXIT
     8940                        ENDDO
     8941                     ENDIF
     8942                  ENDDO
     8943!
     8944!--               Previously free rays cannot be obstructed yet
     8945                  DO  iray = lowest_free_ray, 1, -1
     8946                     IF ( zdirs(iray) > MAX(horz_mixed_entry, horz_mixed_exit) )  THEN
     8947!
     8948!--                     Extend mixed rays by raising the lowest_free ray
     8949!--                     (remains unchanged if this is hit in 1st iteration)
     8950                        lowest_free_ray = iray
     8951                        EXIT
     8952                     ENDIF
     8953                     IF ( zdirs(iray) <= 0._wp )  THEN
     8954!
     8955!--                     Downward pointed ray, cycle k down from entry to exit,
     8956!--                     search upward oriented faces
     8957                        DO  k = FLOOR(0.5_wp + origin(1) + zdirs(iray)*lastdist), &
     8958                               CEILING(-0.5_wp + origin(1) + zdirs(iray)*nextdist)+1, -1
     8959                           target_surfl(iray) = up_col(k) ! contains -999 if missing surface
     8960                           IF ( target_surfl(iray) >= 0 )  EXIT
     8961                        ENDDO
     8962                     ELSE
     8963!
     8964!--                     Upward pointed ray, cycle k up from entry to exit,
     8965!--                     search downward oriented faces
     8966                        DO  k = CEILING(-0.5_wp + origin(1) + zdirs(iray)*lastdist), &
     8967                               FLOOR(0.5_wp + origin(1) + zdirs(iray)*nextdist)-1
     8968                           target_surfl(iray) = down_col(k) ! contains -999 if missing surface
     8969                           IF ( target_surfl(iray) >= 0 )  EXIT
     8970                        ENDDO
     8971                     ENDIF
     8972                  ENDDO
     8973               ENDIF ! end of mixed horizon
     8974            ENDIF ! rad_angular_discretization
     8975            full_horizon = MAX(full_horizon, horz_full_entry, horz_full_exit)
    88768976
    88778977            IF ( plant_canopy )  THEN
     
    88838983         IF ( nextdist + eps >= distance )  EXIT
    88848984
    8885          IF ( rad_angular_discretization  .AND.  calc_svf )  THEN
     8985         IF ( rad_angular_discretization )  THEN
    88868986!
    88878987!--         Save wall direction of coming building column (= this air column)
    88888988            IF ( seldim == 1 )  THEN
    88898989               IF ( dimdelta(seldim) == 1 )  THEN
    8890                   lastdir = isouth_u
     8990                  lastdir = isouth
    88918991               ELSE
    8892                   lastdir = inorth_u
     8992                  lastdir = inorth
    88938993               ENDIF
    88948994            ELSE
    88958995               IF ( dimdelta(seldim) == 1 )  THEN
    8896                   lastdir = iwest_u
     8996                  lastdir = iwest
    88978997               ELSE
    8898                   lastdir = ieast_u
     8998                  lastdir = ieast
    88998999               ENDIF
    89009000            ENDIF
    89019001            lastcolumn = column
    89029002         ENDIF
     9003         ip_last = ip
     9004         ig_last = ig
    89039005         lastdist = nextdist
    89049006         dimnext(seldim) = dimnext(seldim) + dimdelta(seldim)
     
    89219023!--               For fixed view resolution, we need plant canopy even for rays
    89229024!--               to opposing surfaces
    8923                   lowest_lad = nzterr(ig) + 1
     9025                  lowest_lad = nzterrb(ig) + 1
    89249026               ELSE
    89259027!
    8926 !--               We only need LAD for rays directed above horizon (to sky)
     9028!--               We only need LAD for rays directed above full horizon (to sky)
    89279029                  lowest_lad = CEILING( -0.5_wp + origin(1) +            &
    8928                                     MIN( horizon * rt2_track_dist(i-1),  & ! entry
    8929                                          horizon * rt2_track_dist(i)   ) ) ! exit
     9030                                    MIN( full_horizon * rt2_track_dist(i-1),  & ! entry
     9031                                         full_horizon * rt2_track_dist(i)   ) ) ! exit
    89309032               ENDIF
    89319033!
     
    89709072      ENDIF ! plant_canopy
    89719073
    8972       IF ( rad_angular_discretization  .AND.  calc_svf )  THEN
     9074      IF ( rad_angular_discretization )  THEN
    89739075#if defined( __parallel )
    89749076!--      wait for all gridsurf requests to complete
     
    89799081         ENDIF
    89809082#endif
    8981 !
    8982 !--      recalculate local surf indices into global ones
    8983          DO i = 1, nrays
    8984             IF ( target_surfl(i) == -1 )  THEN
    8985                itarget(i) = -1
    8986             ELSE
    8987                itarget(i) = target_surfl(i) + surfstart(target_procs(i))
    8988             ENDIF
     9083         itarget(:) = target_surfl(:)
     9084         DEALLOCATE( target_surfl )
     9085      ELSE
     9086         iray = nrays
     9087         DO WHILE ( iray > 0  .AND.  zdirs(iray) <= full_horizon )
     9088            itarget(iray) = HUGE(1_iwp)
     9089            iray = iray - 1
    89899090         ENDDO
    8990 
    8991          DEALLOCATE( target_surfl )
    8992 
    8993       ELSE
    8994          itarget(:) = -1
     9091         itarget(1:iray) = -1
    89959092      ENDIF ! rad_angular_discretization
    89969093
     
    90419138!--               cause less moiree.
    90429139                  IF ( .NOT. rad_angular_discretization )  THEN
    9043                      IF ( zdirs(k) <= horizon )  CYCLE
     9140                     ! TODO: remove tradiational discretization (incompatible with full 3D)
     9141                     IF ( zdirs(k) <= full_horizon )  CYCLE
    90449142                  ENDIF
    90459143
     
    90529150                     nz = 2
    90539151                     rt2_dist(nz) = SQRT(dxxyy)
    9054                      iz = CEILING(-.5_wp + zorig, iwp)
     9152                     kz = CEILING(-.5_wp + zorig, iwp)
    90559153                  ELSE
    90569154                     zexit = MIN(MAX(origin(1) + zdirs(k) * rt2_track_dist(i), zbottom), ztop)
     
    90629160                     qdist = rt2_dist(nz) / (zexit-zorig)
    90639161                     rt2_dist(2:nz-1) = (/( ((REAL(l, wp) + .5_wp) * zsgn - zorig) * qdist , l = zb0, zb1 )/)
    9064                      iz = zb0 * zsgn
     9162                     kz = zb0 * zsgn
    90659163                  ENDIF
    90669164
    90679165                  DO  l = 2, nz
    9068                      IF ( rt2_track_lad(iz, i) > 0._wp )  THEN
    9069                         curtrans = exp(-ext_coef * rt2_track_lad(iz, i) * (rt2_dist(l)-rt2_dist(l-1)))
     9166                     IF ( rt2_track_lad(kz, i) > 0._wp )  THEN
     9167                        curtrans = exp(-ext_coef * rt2_track_lad(kz, i) * (rt2_dist(l)-rt2_dist(l-1)))
    90709168
    90719169                        IF ( create_csf )  THEN
     
    90749172                           acsf(ncsfl)%itx = rt2_track(2,i)
    90759173                           acsf(ncsfl)%ity = rt2_track(1,i)
    9076                            acsf(ncsfl)%itz = iz
     9174                           acsf(ncsfl)%itz = kz
    90779175                           acsf(ncsfl)%isurfs = iorig
    90789176                           acsf(ncsfl)%rcvf = (1._wp - curtrans)*transparency(k)*vffrac(k)
     
    90819179                        transparency(k) = transparency(k) * curtrans
    90829180                     ENDIF
    9083                      iz = iz + zsgn
     9181                     kz = kz + zsgn
    90849182                  ENDDO ! l = 1, nz - 1
    90859183               ENDDO ! k = 1, nrays
    90869184            ENDDO ! i = 1, ntrack
    9087 
    9088             transparency(1:lowest_free_ray) = 1._wp !-- Reset rays above horizon to transparent (see NOTE 6778)
     9185!
     9186!--         Reset rays above horizon to transparent (see NOTE 6778)
     9187            WHERE( itarget < 0 )  transparency = 1.0_wp
    90899188         ENDIF
    90909189
     
    90989197!
    90999198!--            See NOTE 6778 above
    9100                IF ( zdirs(k) <= horizon )  CYCLE
     9199               IF ( itarget(k) >= 0 )  CYCLE
    91019200
    91029201               zexit = origin(1) + zdirs(k) * rt2_track_dist(i-1)
     
    91089207                  nz = 2
    91099208                  rt2_dist(nz) = SQRT(dxxyy)
    9110                   iz = NINT(zexit, iwp)
     9209                  kz = NINT(zexit, iwp)
    91119210               ELSE
    91129211                  zorig = MIN(MAX(origin(1) + zdirs(k) * rt2_track_dist(i), zbottom), ztop)
     
    91189217                  qdist = rt2_dist(nz) / (zexit-zorig)
    91199218                  rt2_dist(2:nz-1) = (/( ((REAL(l, wp) + .5_wp) * zsgn - zorig) * qdist , l = zb0, zb1 )/)
    9120                   iz = zb0 * zsgn
     9219                  kz = zb0 * zsgn
    91219220               ENDIF
    91229221
    91239222               DO  l = 2, nz
    9124                   IF ( rt2_track_lad(iz, i) > 0._wp )  THEN
    9125                      curtrans = exp(-ext_coef * rt2_track_lad(iz, i) * (rt2_dist(l)-rt2_dist(l-1)))
     9223                  IF ( rt2_track_lad(kz, i) > 0._wp )  THEN
     9224                     curtrans = exp(-ext_coef * rt2_track_lad(kz, i) * (rt2_dist(l)-rt2_dist(l-1)))
    91269225
    91279226                     IF ( create_csf )  THEN
     
    91309229                        acsf(ncsfl)%itx = rt2_track(2,i)
    91319230                        acsf(ncsfl)%ity = rt2_track(1,i)
    9132                         acsf(ncsfl)%itz = iz
    9133                         IF ( itarget(k) /= -1 ) STOP 1 !FIXME remove after test
     9231                        acsf(ncsfl)%itz = kz
    91349232                        acsf(ncsfl)%isurfs = -1
    91359233                        acsf(ncsfl)%rcvf = (1._wp - curtrans)*transparency(k)*aorig*vffrac(k)
     
    91389236                     transparency(k) = transparency(k) * curtrans
    91399237                  ENDIF
    9140                   iz = iz + zsgn
     9238                  kz = kz + zsgn
    91419239               ENDDO ! l = 1, nz - 1
    91429240            ENDDO ! k = 1, nrays
     
    91449242      ENDIF ! plant_canopy
    91459243
    9146       IF ( .NOT. (rad_angular_discretization  .AND.  calc_svf) )  THEN
    9147 !
    9148 !--      Just update lowest_free_ray according to horizon
    9149          DO WHILE ( lowest_free_ray > 0 )
    9150             IF ( zdirs(lowest_free_ray) > horizon )  EXIT
    9151             lowest_free_ray = lowest_free_ray - 1
    9152          ENDDO
    9153       ENDIF
    9154 
    91559244   CONTAINS
    91569245
    9157       SUBROUTINE request_itarget( d, z, y, x, isurfl, iproc )
     9246      SUBROUTINE request_itarget( d, z, y, x, isurfl )
    91589247
    91599248         INTEGER(iwp), INTENT(in)            ::  d, z, y, x
    91609249         INTEGER(iwp), TARGET, INTENT(out)   ::  isurfl
    9161          INTEGER(iwp), INTENT(out)           ::  iproc
     9250         INTEGER(iwp)                        ::  iproc
    91629251
    91639252#if defined( __parallel )
     
    91819270!--      set index target_surfl(i)
    91829271         isurfl = gridsurf(d,z,y,x)
    9183          iproc  = 0  ! required to avoid compile error about unused variable in serial mode
    91849272#endif
    91859273
     
    92999387
    93009388!-- first check: are the two surfaces directed in the same direction
    9301         IF ( (d==iup_u  .OR.  d==iup_l )                             &
    9302              .AND. (d2==iup_u  .OR. d2==iup_l) ) RETURN
    9303         IF ( (d==isouth_u  .OR.  d==isouth_l ) &
    9304              .AND.  (d2==isouth_u  .OR.  d2==isouth_l) ) RETURN
    9305         IF ( (d==inorth_u  .OR.  d==inorth_l ) &
    9306              .AND.  (d2==inorth_u  .OR.  d2==inorth_l) ) RETURN
    9307         IF ( (d==iwest_u  .OR.  d==iwest_l )     &
    9308              .AND.  (d2==iwest_u  .OR.  d2==iwest_l ) ) RETURN
    9309         IF ( (d==ieast_u  .OR.  d==ieast_l )     &
    9310              .AND.  (d2==ieast_u  .OR.  d2==ieast_l ) ) RETURN
     9389        IF ( d==iup     .AND.  d2==iup    ) RETURN
     9390        IF ( d==isouth  .AND.  d2==isouth ) RETURN
     9391        IF ( d==inorth  .AND.  d2==inorth ) RETURN
     9392        IF ( d==iwest   .AND.  d2==iwest  ) RETURN
     9393        IF ( d==ieast   .AND.  d2==ieast  ) RETURN
    93119394
    93129395!-- second check: are surfaces facing away from each other
    93139396        SELECT CASE (d)
    9314             CASE (iup_u, iup_l)                     !< upward facing surfaces
     9397            CASE (iup)                  !< upward facing surfaces
    93159398                IF ( z2 < z ) RETURN
    9316             CASE (isouth_u, isouth_l)               !< southward facing surfaces
     9399            CASE (isouth)               !< southward facing surfaces
    93179400                IF ( y2 > y ) RETURN
    9318             CASE (inorth_u, inorth_l)               !< northward facing surfaces
     9401            CASE (inorth)               !< northward facing surfaces
    93199402                IF ( y2 < y ) RETURN
    9320             CASE (iwest_u, iwest_l)                 !< westward facing surfaces
     9403            CASE (iwest)                !< westward facing surfaces
    93219404                IF ( x2 > x ) RETURN
    9322             CASE (ieast_u, ieast_l)                 !< eastward facing surfaces
     9405            CASE (ieast)                !< eastward facing surfaces
    93239406                IF ( x2 < x ) RETURN
    93249407        END SELECT
    93259408
    93269409        SELECT CASE (d2)
    9327             CASE (iup_u)                            !< ground, roof
     9410            CASE (iup)                  !< ground, roof
    93289411                IF ( z < z2 ) RETURN
    9329             CASE (isouth_u, isouth_l)               !< south facing
     9412            CASE (isouth)               !< south facing
    93309413                IF ( y > y2 ) RETURN
    9331             CASE (inorth_u, inorth_l)               !< north facing
     9414            CASE (inorth)               !< north facing
    93329415                IF ( y < y2 ) RETURN
    9333             CASE (iwest_u, iwest_l)                 !< west facing
     9416            CASE (iwest)                !< west facing
    93349417                IF ( x > x2 ) RETURN
    9335             CASE (ieast_u, ieast_l)                 !< east facing
     9418            CASE (ieast)                !< east facing
    93369419                IF ( x < x2 ) RETURN
    93379420            CASE (-1)
     
    993810021    INTEGER(iwp) ::  l, m !< index of current surface element
    993910022   
    9940     INTEGER(iwp)                                       :: ids, idsint_u, idsint_l, isurf
     10023    INTEGER(iwp)                                       :: ids, idsint, isurf
    994110024    CHARACTER(LEN=varnamelength)                       :: var
    994210025
     
    995110034           IF ( TRIM(var(k-j+1:k)) == TRIM(dirname(i)) )  THEN
    995210035               ids = i
    9953                idsint_u = dirint_u(ids)
    9954                idsint_l = dirint_l(ids)
     10036               idsint = dirint(ids)
    995510037               var = var(:k-j)
    995610038               EXIT
     
    1039910481          CASE ( 'rtm_rad_net' )
    1040010482!--           array of complete radiation balance
    10401               DO isurf = dirstart(ids), dirend(ids)
    10402                  IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10483              DO isurf = 1, nsurfl
     10484                 IF ( surfl(id,isurf) == idsint )  THEN
    1040310485                    surfradnet_av(isurf) = surfradnet_av(isurf) +               &
    1040410486                                           surfinsw(isurf) - surfoutsw(isurf) + &
     
    1040910491          CASE ( 'rtm_rad_insw' )
    1041010492!--           array of sw radiation falling to surface after i-th reflection
    10411               DO isurf = dirstart(ids), dirend(ids)
    10412                   IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10493              DO isurf = 1, nsurfl
     10494                  IF ( surfl(id,isurf) == idsint )  THEN
    1041310495                      surfinsw_av(isurf) = surfinsw_av(isurf) + surfinsw(isurf)
    1041410496                  ENDIF
     
    1041710499          CASE ( 'rtm_rad_inlw' )
    1041810500!--           array of lw radiation falling to surface after i-th reflection
    10419               DO isurf = dirstart(ids), dirend(ids)
    10420                   IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10501              DO isurf = 1, nsurfl
     10502                  IF ( surfl(id,isurf) == idsint )  THEN
    1042110503                      surfinlw_av(isurf) = surfinlw_av(isurf) + surfinlw(isurf)
    1042210504                  ENDIF
     
    1042510507          CASE ( 'rtm_rad_inswdir' )
    1042610508!--           array of direct sw radiation falling to surface from sun
    10427               DO isurf = dirstart(ids), dirend(ids)
    10428                   IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10509              DO isurf = 1, nsurfl
     10510                  IF ( surfl(id,isurf) == idsint )  THEN
    1042910511                      surfinswdir_av(isurf) = surfinswdir_av(isurf) + surfinswdir(isurf)
    1043010512                  ENDIF
     
    1043310515          CASE ( 'rtm_rad_inswdif' )
    1043410516!--           array of difusion sw radiation falling to surface from sky and borders of the domain
    10435               DO isurf = dirstart(ids), dirend(ids)
    10436                   IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10517              DO isurf = 1, nsurfl
     10518                  IF ( surfl(id,isurf) == idsint )  THEN
    1043710519                      surfinswdif_av(isurf) = surfinswdif_av(isurf) + surfinswdif(isurf)
    1043810520                  ENDIF
     
    1044110523          CASE ( 'rtm_rad_inswref' )
    1044210524!--           array of sw radiation falling to surface from reflections
    10443               DO isurf = dirstart(ids), dirend(ids)
    10444                   IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10525              DO isurf = 1, nsurfl
     10526                  IF ( surfl(id,isurf) == idsint )  THEN
    1044510527                      surfinswref_av(isurf) = surfinswref_av(isurf) + surfinsw(isurf) - &
    1044610528                                          surfinswdir(isurf) - surfinswdif(isurf)
     
    1045110533          CASE ( 'rtm_rad_inlwdif' )
    1045210534!--           array of sw radiation falling to surface after i-th reflection
    10453               DO isurf = dirstart(ids), dirend(ids)
    10454                   IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10535              DO isurf = 1, nsurfl
     10536                  IF ( surfl(id,isurf) == idsint )  THEN
    1045510537                      surfinlwdif_av(isurf) = surfinlwdif_av(isurf) + surfinlwdif(isurf)
    1045610538                  ENDIF
     
    1045910541          CASE ( 'rtm_rad_inlwref' )
    1046010542!--           array of lw radiation falling to surface from reflections
    10461               DO isurf = dirstart(ids), dirend(ids)
    10462                   IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10543              DO isurf = 1, nsurfl
     10544                  IF ( surfl(id,isurf) == idsint )  THEN
    1046310545                      surfinlwref_av(isurf) = surfinlwref_av(isurf) + &
    1046410546                                          surfinlw(isurf) - surfinlwdif(isurf)
     
    1046810550          CASE ( 'rtm_rad_outsw' )
    1046910551!--           array of sw radiation emitted from surface after i-th reflection
    10470               DO isurf = dirstart(ids), dirend(ids)
    10471                   IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10552              DO isurf = 1, nsurfl
     10553                  IF ( surfl(id,isurf) == idsint )  THEN
    1047210554                      surfoutsw_av(isurf) = surfoutsw_av(isurf) + surfoutsw(isurf)
    1047310555                  ENDIF
     
    1047610558          CASE ( 'rtm_rad_outlw' )
    1047710559!--           array of lw radiation emitted from surface after i-th reflection
    10478               DO isurf = dirstart(ids), dirend(ids)
    10479                   IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10560              DO isurf = 1, nsurfl
     10561                  IF ( surfl(id,isurf) == idsint )  THEN
    1048010562                      surfoutlw_av(isurf) = surfoutlw_av(isurf) + surfoutlw(isurf)
    1048110563                  ENDIF
     
    1048410566          CASE ( 'rtm_rad_ressw' )
    1048510567!--           array of residua of sw radiation absorbed in surface after last reflection
    10486               DO isurf = dirstart(ids), dirend(ids)
    10487                   IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10568              DO isurf = 1, nsurfl
     10569                  IF ( surfl(id,isurf) == idsint )  THEN
    1048810570                      surfins_av(isurf) = surfins_av(isurf) + surfins(isurf)
    1048910571                  ENDIF
     
    1049210574          CASE ( 'rtm_rad_reslw' )
    1049310575!--           array of residua of lw radiation absorbed in surface after last reflection
    10494               DO isurf = dirstart(ids), dirend(ids)
    10495                   IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10576              DO isurf = 1, nsurfl
     10577                  IF ( surfl(id,isurf) == idsint )  THEN
    1049610578                      surfinl_av(isurf) = surfinl_av(isurf) + surfinl(isurf)
    1049710579                  ENDIF
     
    1069610778          CASE ( 'rtm_rad_net' )
    1069710779!--           array of complete radiation balance
    10698               DO isurf = dirstart(ids), dirend(ids)
    10699                   IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10780              DO isurf = 1, nsurfl
     10781                  IF ( surfl(id,isurf) == idsint )  THEN
    1070010782                      surfradnet_av(isurf) = surfradnet_av(isurf) / REAL( average_count_3d, kind=wp )
    1070110783                  ENDIF
     
    1070410786          CASE ( 'rtm_rad_insw' )
    1070510787!--           array of sw radiation falling to surface after i-th reflection
    10706               DO isurf = dirstart(ids), dirend(ids)
    10707                   IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10788              DO isurf = 1, nsurfl
     10789                  IF ( surfl(id,isurf) == idsint )  THEN
    1070810790                      surfinsw_av(isurf) = surfinsw_av(isurf) / REAL( average_count_3d, kind=wp )
    1070910791                  ENDIF
     
    1071210794          CASE ( 'rtm_rad_inlw' )
    1071310795!--           array of lw radiation falling to surface after i-th reflection
    10714               DO isurf = dirstart(ids), dirend(ids)
    10715                   IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10796              DO isurf = 1, nsurfl
     10797                  IF ( surfl(id,isurf) == idsint )  THEN
    1071610798                      surfinlw_av(isurf) = surfinlw_av(isurf) / REAL( average_count_3d, kind=wp )
    1071710799                  ENDIF
     
    1072010802          CASE ( 'rtm_rad_inswdir' )
    1072110803!--           array of direct sw radiation falling to surface from sun
    10722               DO isurf = dirstart(ids), dirend(ids)
    10723                   IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10804              DO isurf = 1, nsurfl
     10805                  IF ( surfl(id,isurf) == idsint )  THEN
    1072410806                      surfinswdir_av(isurf) = surfinswdir_av(isurf) / REAL( average_count_3d, kind=wp )
    1072510807                  ENDIF
     
    1072810810          CASE ( 'rtm_rad_inswdif' )
    1072910811!--           array of difusion sw radiation falling to surface from sky and borders of the domain
    10730               DO isurf = dirstart(ids), dirend(ids)
    10731                   IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10812              DO isurf = 1, nsurfl
     10813                  IF ( surfl(id,isurf) == idsint )  THEN
    1073210814                      surfinswdif_av(isurf) = surfinswdif_av(isurf) / REAL( average_count_3d, kind=wp )
    1073310815                  ENDIF
     
    1073610818          CASE ( 'rtm_rad_inswref' )
    1073710819!--           array of sw radiation falling to surface from reflections
    10738               DO isurf = dirstart(ids), dirend(ids)
    10739                   IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10820              DO isurf = 1, nsurfl
     10821                  IF ( surfl(id,isurf) == idsint )  THEN
    1074010822                      surfinswref_av(isurf) = surfinswref_av(isurf) / REAL( average_count_3d, kind=wp )
    1074110823                  ENDIF
     
    1074410826          CASE ( 'rtm_rad_inlwdif' )
    1074510827!--           array of sw radiation falling to surface after i-th reflection
    10746               DO isurf = dirstart(ids), dirend(ids)
    10747                   IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10828              DO isurf = 1, nsurfl
     10829                  IF ( surfl(id,isurf) == idsint )  THEN
    1074810830                      surfinlwdif_av(isurf) = surfinlwdif_av(isurf) / REAL( average_count_3d, kind=wp )
    1074910831                  ENDIF
     
    1075210834          CASE ( 'rtm_rad_inlwref' )
    1075310835!--           array of lw radiation falling to surface from reflections
    10754               DO isurf = dirstart(ids), dirend(ids)
    10755                   IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10836              DO isurf = 1, nsurfl
     10837                  IF ( surfl(id,isurf) == idsint )  THEN
    1075610838                      surfinlwref_av(isurf) = surfinlwref_av(isurf) / REAL( average_count_3d, kind=wp )
    1075710839                  ENDIF
     
    1076010842          CASE ( 'rtm_rad_outsw' )
    1076110843!--           array of sw radiation emitted from surface after i-th reflection
    10762               DO isurf = dirstart(ids), dirend(ids)
    10763                   IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10844              DO isurf = 1, nsurfl
     10845                  IF ( surfl(id,isurf) == idsint )  THEN
    1076410846                      surfoutsw_av(isurf) = surfoutsw_av(isurf) / REAL( average_count_3d, kind=wp )
    1076510847                  ENDIF
     
    1076810850          CASE ( 'rtm_rad_outlw' )
    1076910851!--           array of lw radiation emitted from surface after i-th reflection
    10770               DO isurf = dirstart(ids), dirend(ids)
    10771                   IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10852              DO isurf = 1, nsurfl
     10853                  IF ( surfl(id,isurf) == idsint )  THEN
    1077210854                      surfoutlw_av(isurf) = surfoutlw_av(isurf) / REAL( average_count_3d, kind=wp )
    1077310855                  ENDIF
     
    1077610858          CASE ( 'rtm_rad_ressw' )
    1077710859!--           array of residua of sw radiation absorbed in surface after last reflection
    10778               DO isurf = dirstart(ids), dirend(ids)
    10779                   IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10860              DO isurf = 1, nsurfl
     10861                  IF ( surfl(id,isurf) == idsint )  THEN
    1078010862                      surfins_av(isurf) = surfins_av(isurf) / REAL( average_count_3d, kind=wp )
    1078110863                  ENDIF
     
    1078410866          CASE ( 'rtm_rad_reslw' )
    1078510867!--           array of residua of lw radiation absorbed in surface after last reflection
    10786               DO isurf = dirstart(ids), dirend(ids)
    10787                   IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     10868              DO isurf = 1, nsurfl
     10869                  IF ( surfl(id,isurf) == idsint )  THEN
    1078810870                      surfinl_av(isurf) = surfinl_av(isurf) / REAL( average_count_3d, kind=wp )
    1078910871                  ENDIF
     
    1135411436
    1135511437    CHARACTER (len=varnamelength)                   :: var, surfid
    11356     INTEGER(iwp)                                    :: ids,idsint_u,idsint_l,isurf,isvf,isurfs,isurflt,ipcgb
     11438    INTEGER(iwp)                                    :: ids,idsint,isurf,isvf,isurfs,isurflt,ipcgb
    1135711439    INTEGER(iwp)                                    :: is, js, ks, istat
    1135811440
     
    1137811460           IF ( TRIM(var(k-j+1:k)) == TRIM(dirname(i)) )  THEN
    1137911461              ids = i
    11380               idsint_u = dirint_u(ids)
    11381               idsint_l = dirint_l(ids)
     11462              idsint = dirint(ids)
    1138211463              var = var(:k-j)
    1138311464              EXIT
     
    1159611677      CASE ( 'rtm_rad_net' )
    1159711678!--     array of complete radiation balance
    11598          DO isurf = dirstart(ids), dirend(ids)
    11599             IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     11679         DO isurf = 1, nsurfl
     11680            IF ( surfl(id,isurf) == idsint )  THEN
    1160011681               IF ( av == 0 )  THEN
    1160111682                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = &
     
    1160911690      CASE ( 'rtm_rad_insw' )
    1161011691!--      array of sw radiation falling to surface after i-th reflection
    11611          DO isurf = dirstart(ids), dirend(ids)
    11612             IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     11692         DO isurf = 1, nsurfl
     11693            IF ( surfl(id,isurf) == idsint )  THEN
    1161311694               IF ( av == 0 )  THEN
    1161411695                 local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinsw(isurf)
     
    1162111702      CASE ( 'rtm_rad_inlw' )
    1162211703!--      array of lw radiation falling to surface after i-th reflection
    11623          DO isurf = dirstart(ids), dirend(ids)
    11624             IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     11704         DO isurf = 1, nsurfl
     11705            IF ( surfl(id,isurf) == idsint )  THEN
    1162511706               IF ( av == 0 )  THEN
    1162611707                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinlw(isurf)
     
    1163311714      CASE ( 'rtm_rad_inswdir' )
    1163411715!--      array of direct sw radiation falling to surface from sun
    11635          DO isurf = dirstart(ids), dirend(ids)
    11636             IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     11716         DO isurf = 1, nsurfl
     11717            IF ( surfl(id,isurf) == idsint )  THEN
    1163711718               IF ( av == 0 )  THEN
    1163811719                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinswdir(isurf)
     
    1164511726      CASE ( 'rtm_rad_inswdif' )
    1164611727!--      array of difusion sw radiation falling to surface from sky and borders of the domain
    11647          DO isurf = dirstart(ids), dirend(ids)
    11648             IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     11728         DO isurf = 1, nsurfl
     11729            IF ( surfl(id,isurf) == idsint )  THEN
    1164911730               IF ( av == 0 )  THEN
    1165011731                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinswdif(isurf)
     
    1165711738      CASE ( 'rtm_rad_inswref' )
    1165811739!--      array of sw radiation falling to surface from reflections
    11659          DO isurf = dirstart(ids), dirend(ids)
    11660             IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     11740         DO isurf = 1, nsurfl
     11741            IF ( surfl(id,isurf) == idsint )  THEN
    1166111742               IF ( av == 0 )  THEN
    1166211743                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = &
     
    1167011751      CASE ( 'rtm_rad_inlwdif' )
    1167111752!--      array of difusion lw radiation falling to surface from sky and borders of the domain
    11672          DO isurf = dirstart(ids), dirend(ids)
    11673             IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     11753         DO isurf = 1, nsurfl
     11754            IF ( surfl(id,isurf) == idsint )  THEN
    1167411755               IF ( av == 0 )  THEN
    1167511756                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinlwdif(isurf)
     
    1168211763      CASE ( 'rtm_rad_inlwref' )
    1168311764!--      array of lw radiation falling to surface from reflections
    11684          DO isurf = dirstart(ids), dirend(ids)
    11685             IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     11765         DO isurf = 1, nsurfl
     11766            IF ( surfl(id,isurf) == idsint )  THEN
    1168611767               IF ( av == 0 )  THEN
    1168711768                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinlw(isurf) - surfinlwdif(isurf)
     
    1169411775      CASE ( 'rtm_rad_outsw' )
    1169511776!--      array of sw radiation emitted from surface after i-th reflection
    11696          DO isurf = dirstart(ids), dirend(ids)
    11697             IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     11777         DO isurf = 1, nsurfl
     11778            IF ( surfl(id,isurf) == idsint )  THEN
    1169811779               IF ( av == 0 )  THEN
    1169911780                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfoutsw(isurf)
     
    1170611787      CASE ( 'rtm_rad_outlw' )
    1170711788!--      array of lw radiation emitted from surface after i-th reflection
    11708          DO isurf = dirstart(ids), dirend(ids)
    11709             IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     11789         DO isurf = 1, nsurfl
     11790            IF ( surfl(id,isurf) == idsint )  THEN
    1171011791               IF ( av == 0 )  THEN
    1171111792                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfoutlw(isurf)
     
    1171811799      CASE ( 'rtm_rad_ressw' )
    1171911800!--      average of array of residua of sw radiation absorbed in surface after last reflection
    11720          DO isurf = dirstart(ids), dirend(ids)
    11721             IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     11801         DO isurf = 1, nsurfl
     11802            IF ( surfl(id,isurf) == idsint )  THEN
    1172211803               IF ( av == 0 )  THEN
    1172311804                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfins(isurf)
     
    1173011811      CASE ( 'rtm_rad_reslw' )
    1173111812!--      average of array of residua of lw radiation absorbed in surface after last reflection
    11732          DO isurf = dirstart(ids), dirend(ids)
    11733             IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     11813         DO isurf = 1, nsurfl
     11814            IF ( surfl(id,isurf) == idsint )  THEN
    1173411815               IF ( av == 0 )  THEN
    1173511816                  local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = surfinl(isurf)
     
    1183811919!
    1183911920!--      sky view factor
    11840          DO isurf = dirstart(ids), dirend(ids)
    11841             IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     11921         DO isurf = 1, nsurfl
     11922            IF ( surfl(id,isurf) == idsint )  THEN
    1184211923               local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = skyvf(isurf)
    1184311924            ENDIF
     
    1184711928!
    1184811929!--      sky view factor
    11849          DO isurf = dirstart(ids), dirend(ids)
    11850             IF ( surfl(id,isurf) == ids )  THEN
     11930         DO isurf = 1, nsurfl
     11931            IF ( surfl(id,isurf) == idsint )  THEN
    1185111932               local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = skyvft(isurf)
    1185211933            ENDIF
     
    1186611947
    1186711948            IF ( surf(ix,isurfs) == is  .AND.  surf(iy,isurfs) == js  .AND. surf(iz,isurfs) == ks  .AND. &
    11868                  (surfl(id,isurflt) == idsint_u .OR. surfl(id,isurflt) == idsint_l ) ) THEN
     11949                 surfl(id,isurflt) == idsint ) THEN
    1186911950!
    1187011951!--            correct source surface
     
    1187611957!
    1187711958!--      surface albedo
    11878          DO isurf = dirstart(ids), dirend(ids)
    11879             IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     11959         DO isurf = 1, nsurfl
     11960            IF ( surfl(id,isurf) == idsint )  THEN
    1188011961               local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = albedo_surf(isurf)
    1188111962            ENDIF
     
    1188511966!
    1188611967!--      surface emissivity, weighted average
    11887          DO isurf = dirstart(ids), dirend(ids)
    11888             IF ( surfl(id,isurf) == idsint_u  .OR.  surfl(id,isurf) == idsint_l )  THEN
     11968         DO isurf = 1, nsurfl
     11969            IF ( surfl(id,isurf) == idsint )  THEN
    1188911970               local_pf(surfl(ix,isurf),surfl(iy,isurf),surfl(iz,isurf)) = emiss_surf(isurf)
    1189011971            ENDIF
  • palm/trunk/SOURCE/urban_surface_mod.f90

    r4630 r4653  
    362362               force_radiation_call,                                                               &
    363363               id,                                                                                 &
    364                ieast_l,                                                                            &
    365                ieast_u,                                                                            &
    366                inorth_l,                                                                           &
    367                inorth_u,                                                                           &
    368                isouth_l,                                                                           &
    369                isouth_u,                                                                           &
    370                iup_l,                                                                              &
    371                iup_u,                                                                              &
    372                iwest_l,                                                                            &
    373                iwest_u,                                                                            &
     364               ieast,                                                                              &
     365               inorth,                                                                             &
     366               isouth,                                                                             &
     367               iup,                                                                                &
     368               iwest,                                                                              &
    374369               nz_urban_b,                                                                         &
    375370               nz_urban_t,                                                                         &
     
    13431338    INTEGER(iwp), PARAMETER                            ::  nd = 5                                  !< number of directions
    13441339    CHARACTER(LEN=6), DIMENSION(0:nd-1), PARAMETER     ::  dirname = (/ '_roof ', '_south', '_north', '_west ', '_east ' /)
    1345     INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER         ::  dirint = (/ iup_u, isouth_u, inorth_u, iwest_u, ieast_u /)
     1340    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER         ::  dirint = (/ iup, isouth, inorth, iwest, ieast /)
    13461341
    13471342
     
    24602455
    24612456    INTEGER(iwp), PARAMETER                            ::  nd = 5  !< number of directions
    2462     INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER         ::  dirint =  (/    iup_u, isouth_u, inorth_u,  iwest_u,  ieast_u /)  !<
     2457    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER         ::  dirint =  (/      iup,   isouth,   inorth,    iwest,    ieast /)  !<
    24632458    INTEGER(iwp), DIMENSION(0:nd-1), PARAMETER         ::  diridx =  (/       -1,        1,        0,        3,        2 /)
    24642459                                                           !< index for surf_*_v: 0:3 = (North, South, East, West)
     
    25312526!
    25322527!--       Array of surface height (z)
    2533           IF ( idsint == iup_u )  THEN
     2528          IF ( idsint == iup )  THEN
    25342529             DO  m = 1, surf_usm_h%ns
    25352530                i = surf_usm_h%i(m)
     
    25512546!
    25522547!--       Surface category
    2553           IF ( idsint == iup_u )  THEN
     2548          IF ( idsint == iup )  THEN
    25542549             DO  m = 1, surf_usm_h%ns
    25552550                i = surf_usm_h%i(m)
     
    25712566!
    25722567!--       Transmissivity window tiles
    2573           IF ( idsint == iup_u )  THEN
     2568          IF ( idsint == iup )  THEN
    25742569             DO  m = 1, surf_usm_h%ns
    25752570                i = surf_usm_h%i(m)
     
    25922587!--       Array of sensible heat flux from surfaces
    25932588          IF ( av == 0 )  THEN
    2594              IF ( idsint == iup_u )  THEN
     2589             IF ( idsint == iup )  THEN
    25952590                DO  m = 1, surf_usm_h%ns
    25962591                   i = surf_usm_h%i(m)
     
    26092604             ENDIF
    26102605          ELSE
    2611              IF ( idsint == iup_u )  THEN
     2606             IF ( idsint == iup )  THEN
    26122607                DO  m = 1, surf_usm_h%ns
    26132608                   i = surf_usm_h%i(m)
     
    26322627!--       Array of latent heat flux from surfaces
    26332628          IF ( av == 0 )  THEN
    2634              IF ( idsint == iup_u )  THEN
     2629             IF ( idsint == iup )  THEN
    26352630                DO  m = 1, surf_usm_h%ns
    26362631                   i = surf_usm_h%i(m)
     
    26492644             ENDIF
    26502645          ELSE
    2651              IF ( idsint == iup_u )  THEN
     2646             IF ( idsint == iup )  THEN
    26522647                DO  m = 1, surf_usm_h%ns
    26532648                   i = surf_usm_h%i(m)
     
    26712666!--       Array of latent heat flux from vegetation surfaces
    26722667          IF ( av == 0 )  THEN
    2673              IF ( idsint == iup_u )  THEN
     2668             IF ( idsint == iup )  THEN
    26742669                DO  m = 1, surf_usm_h%ns
    26752670                   i = surf_usm_h%i(m)
     
    26882683             ENDIF
    26892684          ELSE
    2690              IF ( idsint == iup_u )  THEN
     2685             IF ( idsint == iup )  THEN
    26912686                DO  m = 1, surf_usm_h%ns
    26922687                   i = surf_usm_h%i(m)
     
    27102705!--       Array of latent heat flux from surfaces with liquid
    27112706          IF ( av == 0 )  THEN
    2712              IF ( idsint == iup_u )  THEN
     2707             IF ( idsint == iup )  THEN
    27132708                DO  m = 1, surf_usm_h%ns
    27142709                   i = surf_usm_h%i(m)
     
    27272722             ENDIF
    27282723          ELSE
    2729              IF ( idsint == iup_u )  THEN
     2724             IF ( idsint == iup )  THEN
    27302725                DO  m = 1, surf_usm_h%ns
    27312726                   i = surf_usm_h%i(m)
     
    27492744!--       Array of heat flux from ground (land, wall, roof)
    27502745          IF ( av == 0 )  THEN
    2751              IF ( idsint == iup_u )  THEN
     2746             IF ( idsint == iup )  THEN
    27522747                DO  m = 1, surf_usm_h%ns
    27532748                   i = surf_usm_h%i(m)
     
    27662761             ENDIF
    27672762          ELSE
    2768              IF ( idsint == iup_u )  THEN
     2763             IF ( idsint == iup )  THEN
    27692764                DO  m = 1, surf_usm_h%ns
    27702765                   i = surf_usm_h%i(m)
     
    27882783!--       Array of heat flux from window ground (land, wall, roof)
    27892784          IF ( av == 0 )  THEN
    2790              IF ( idsint == iup_u )  THEN
     2785             IF ( idsint == iup )  THEN
    27912786                DO  m = 1, surf_usm_h%ns
    27922787                   i = surf_usm_h%i(m)
     
    28052800             ENDIF
    28062801          ELSE
    2807              IF ( idsint == iup_u )  THEN
     2802             IF ( idsint == iup )  THEN
    28082803                DO  m = 1, surf_usm_h%ns
    28092804                   i = surf_usm_h%i(m)
     
    28272822!--       Array of heat flux from green ground (land, wall, roof)
    28282823          IF ( av == 0 )  THEN
    2829              IF ( idsint == iup_u )  THEN
     2824             IF ( idsint == iup )  THEN
    28302825                DO  m = 1, surf_usm_h%ns
    28312826                   i = surf_usm_h%i(m)
     
    28442839             ENDIF
    28452840          ELSE
    2846              IF ( idsint == iup_u )  THEN
     2841             IF ( idsint == iup )  THEN
    28472842                DO  m = 1, surf_usm_h%ns
    28482843                   i = surf_usm_h%i(m)
     
    28662861!--       Array of heat flux from indoor ground (land, wall, roof)
    28672862          IF ( av == 0 )  THEN
    2868              IF ( idsint == iup_u )  THEN
     2863             IF ( idsint == iup )  THEN
    28692864                DO  m = 1, surf_usm_h%ns
    28702865                   i = surf_usm_h%i(m)
     
    28832878             ENDIF
    28842879          ELSE
    2885              IF ( idsint == iup_u )  THEN
     2880             IF ( idsint == iup )  THEN
    28862881                DO  m = 1, surf_usm_h%ns
    28872882                   i = surf_usm_h%i(m)
     
    29052900!--       Array of heat flux from indoor window ground (land, wall, roof)
    29062901          IF ( av == 0 )  THEN
    2907              IF ( idsint == iup_u )  THEN
     2902             IF ( idsint == iup )  THEN
    29082903                DO  m = 1, surf_usm_h%ns
    29092904                   i = surf_usm_h%i(m)
     
    29222917             ENDIF
    29232918          ELSE
    2924              IF ( idsint == iup_u )  THEN
     2919             IF ( idsint == iup )  THEN
    29252920                DO  m = 1, surf_usm_h%ns
    29262921                   i = surf_usm_h%i(m)
     
    29442939!--       Surface temperature for surfaces
    29452940          IF ( av == 0 )  THEN
    2946              IF ( idsint == iup_u )  THEN
     2941             IF ( idsint == iup )  THEN
    29472942                DO  m = 1, surf_usm_h%ns
    29482943                   i = surf_usm_h%i(m)
     
    29612956             ENDIF
    29622957          ELSE
    2963              IF ( idsint == iup_u )  THEN
     2958             IF ( idsint == iup )  THEN
    29642959                DO  m = 1, surf_usm_h%ns
    29652960                   i = surf_usm_h%i(m)
     
    29832978!--       Surface temperature for window surfaces
    29842979          IF ( av == 0 )  THEN
    2985              IF ( idsint == iup_u )  THEN
     2980             IF ( idsint == iup )  THEN
    29862981                DO  m = 1, surf_usm_h%ns
    29872982                   i = surf_usm_h%i(m)
     
    30012996
    30022997          ELSE
    3003              IF ( idsint == iup_u )  THEN
     2998             IF ( idsint == iup )  THEN
    30042999                DO  m = 1, surf_usm_h%ns
    30053000                   i = surf_usm_h%i(m)
     
    30253020!--       Surface temperature for green surfaces
    30263021          IF ( av == 0 )  THEN
    3027              IF ( idsint == iup_u )  THEN
     3022             IF ( idsint == iup )  THEN
    30283023                DO  m = 1, surf_usm_h%ns
    30293024                   i = surf_usm_h%i(m)
     
    30433038
    30443039          ELSE
    3045              IF ( idsint == iup_u )  THEN
     3040             IF ( idsint == iup )  THEN
    30463041                DO  m = 1, surf_usm_h%ns
    30473042                   i = surf_usm_h%i(m)
     
    30673062!--       Near surface temperature for whole surfaces
    30683063          IF ( av == 0 )  THEN
    3069              IF ( idsint == iup_u )  THEN
     3064             IF ( idsint == iup )  THEN
    30703065                DO  m = 1, surf_usm_h%ns
    30713066                   i = surf_usm_h%i(m)
     
    30863081
    30873082          ELSE
    3088              IF ( idsint == iup_u )  THEN
     3083             IF ( idsint == iup )  THEN
    30893084                DO  m = 1, surf_usm_h%ns
    30903085                   i = surf_usm_h%i(m)
     
    31093104!--       Wall temperature for  iwl layer of walls and land
    31103105          IF ( av == 0 )  THEN
    3111              IF ( idsint == iup_u )  THEN
     3106             IF ( idsint == iup )  THEN
    31123107                DO  m = 1, surf_usm_h%ns
    31133108                   i = surf_usm_h%i(m)
     
    31263121             ENDIF
    31273122          ELSE
    3128              IF ( idsint == iup_u )  THEN
     3123             IF ( idsint == iup )  THEN
    31293124                DO  m = 1, surf_usm_h%ns
    31303125                   i = surf_usm_h%i(m)
     
    31483143!--       Window temperature for iwl layer of walls and land
    31493144          IF ( av == 0 )  THEN
    3150              IF ( idsint == iup_u )  THEN
     3145             IF ( idsint == iup )  THEN
    31513146                DO  m = 1, surf_usm_h%ns
    31523147                   i = surf_usm_h%i(m)
     
    31653160             ENDIF
    31663161          ELSE
    3167              IF ( idsint == iup_u )  THEN
     3162             IF ( idsint == iup )  THEN
    31683163                DO  m = 1, surf_usm_h%ns
    31693164                   i = surf_usm_h%i(m)
     
    31873182!--       Green temperature for  iwl layer of walls and land
    31883183          IF ( av == 0 )  THEN
    3189              IF ( idsint == iup_u )  THEN
     3184             IF ( idsint == iup )  THEN
    31903185                DO  m = 1, surf_usm_h%ns
    31913186                   i = surf_usm_h%i(m)
     
    32043199             ENDIF
    32053200          ELSE
    3206              IF ( idsint == iup_u )  THEN
     3201             IF ( idsint == iup )  THEN
    32073202                DO  m = 1, surf_usm_h%ns
    32083203                   i = surf_usm_h%i(m)
     
    32263221!--       Soil water content for  iwl layer of walls and land
    32273222          IF ( av == 0 )  THEN
    3228              IF ( idsint == iup_u )  THEN
     3223             IF ( idsint == iup )  THEN
    32293224                DO  m = 1, surf_usm_h%ns
    32303225                   i = surf_usm_h%i(m)
     
    32373232             ENDIF
    32383233          ELSE
    3239              IF ( idsint == iup_u )  THEN
     3234             IF ( idsint == iup )  THEN
    32403235                DO  m = 1, surf_usm_h%ns
    32413236                   i = surf_usm_h%i(m)
Note: See TracChangeset for help on using the changeset viewer.